First-principles approach to the dynamic magnetoelectric couplings for the non-reciprocal directional dichroism in BiFeO₃^*

The FE polarization ${{\bf{P}}}^{{\rm{FE}}}\parallel {\bf{z}}^{\prime}$ emerging below ${T}_{{\rm{C}}}$ can take eight different orientations along the four cubic diagonals ${\bf{z}}^{\prime} =\langle 1,1,1\rangle$. For a given ${\bf{z}}^{\prime} =[1,1,1]$, the three possible orientations for the ${{\bf{x}}}_{m}^{\prime }$ cycloidal modulation wavevectors are ${{\bf{x}}}_{1}^{\prime }=[1,-1,0]$, ${{\bf{x}}}_{2}^{\prime }=[1,0,-1]$, and ${{\bf{x}}}_{3}^{\prime }=[0,1,-1]$ with corresponding ${\bf{y}}{^{\prime} }_{m}={\bf{z}}^{\prime} \times {{\bf{x}}}_{m}^{\prime }$. In magnetic domain m, the cycloidal ordering wavevectors are

$$\begin{eqnarray}&&{{\bf{Q}}}_{m}={{\bf{Q}}}_{0}+\displaystyle \frac{2\sqrt{2}\pi \delta }{a}{{\bf{x}}}_{m}^{\prime },\end{eqnarray} \tag{ 1 }$$

where ${{\bf{Q}}}_{0}=(\pi /a)(1,1,1)$ is the wavevector for simple G-type antiferromagnet and $a=3.96\;\mathring{\rm A}$ is the pseudo-cubic lattice constant. Hence, the ordering wavevectors are ${{\bf{Q}}}_{1}=(2\pi /a)(0.5+\delta ,0.5-\delta ,0.5)$, ${{\bf{Q}}}_{2}=(2\pi /a)(0.5+\delta ,0.5,0.5-\delta )$, and ${{\bf{Q}}}_{3}=(2\pi /a)(0.5,0.5+\delta ,0.5-\delta )$. In terms of $\delta \ll 1$, the cycloidal period is $\lambda =a/(\sqrt{2}\delta )$.

FE and AF distortions create the DM interactions ${{\boldsymbol{D}}}_{{\rm{FE}}}$ and ${{\boldsymbol{D}}}_{{\rm{AF}}}$. Including all magnetic anisotropies produced by these distortions, the spin Hamiltonian can be written

$$\begin{eqnarray}&&{ \mathcal H }=\ {{ \mathcal H }}_{{\rm{FE}}}^{{\rm{SC}}}\ +\ {{ \mathcal H }}_{{\rm{AF}}}^{{\rm{SC}}}\ +{{ \mathcal H }}^{{\rm{EX}}}+{{ \mathcal H }}^{{\rm{SIA}}},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}_{{\rm{FE}}}^{{\rm{SC}}}=\displaystyle \sum _{\langle i,j\rangle }{{\boldsymbol{D}}}_{{\rm{FE}}}\cdot ({{\bf{S}}}_{i}\times {{\bf{S}}}_{j}),\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}_{{\rm{AF}}}^{{\rm{SC}}}=\displaystyle \sum _{\langle i,j\rangle }\;{(-1)}^{{n}_{i}}{{\boldsymbol{D}}}_{{\rm{AF}}}\cdot ({{\bf{S}}}_{i}\times {{\bf{S}}}_{j}),\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}^{{\rm{EX}}}=-{J}_{1}\displaystyle \sum _{\langle i,j\rangle }{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}-{J}_{2}\displaystyle \sum _{\langle i,j\rangle ^{\prime} }{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{{ \mathcal H }}^{{\rm{SIA}}}=-K\displaystyle \sum _{i}{({{\bf{S}}}_{i}\cdot {\bf{z}}^{\prime} )}^{2},\end{eqnarray} \tag{ 6 }$$

where $\langle i,j\rangle$ and $\langle i,j\rangle ^{\prime}$ represent nearest and next-nearest neighbor spins, respectively. This is the most general Hamiltonian that includes the allowed distortions in R3c BiFeO₃ but neglecting exchange anisotropy terms of the form

$$\begin{eqnarray*}&&{J}_{x^{\prime} }({S}_{{ix}^{\prime} }{S}_{{jx}^{\prime} }+{S}_{{iy}^{\prime} }{S}_{{jy}^{\prime} })+{J}_{z^{\prime} }{S}_{{iz}^{\prime} }{S}_{{jz}^{\prime} },\end{eqnarray*}$$

which are usually small for transition metal ions with half-filled d-shell such as Fe³⁺. Moreover, due to the long wavelength of the cycloid, exchange anisotropy can be effectively absorbed into the single-ion anisotropy (SIA) $K\gt 0$, which favors spin alignment along ${\bf{z}}^{\prime}$. All terms in this Hamitlonian are also essential to explain the spin modes of BiFeO₃ observed using THz spectroscopy [21, 22].

Since the FE distortion is uniform, the ${{\boldsymbol{D}}}_{{\rm{FE}}}$ sum is translation invariant. Due to the translation-odd ${{\rm{R}}}_{4}^{+}$[111] AF octahedral rotation, the ${{\boldsymbol{D}}}_{{\rm{AF}}}$ sum contains the coefficient ${(-1)}^{{n}_{i}}$, which alternates from one hexagonal layer ${n}_{i}=\sqrt{3}{\bf{z}}^{\prime} \cdot {{\bf{R}}}_{i}/a$ to the next. Simplified forms for the DM terms ${{ \mathcal H }}_{{\rm{FE}}}^{{\rm{SC}}}$ and ${{ \mathcal H }}_{{\rm{AF}}}^{{\rm{SC}}}$ are given in appendix A.

By ignoring the higher cycloidal harmonics but including the tilt [23] τ produced by ${{\boldsymbol{D}}}_{{\rm{AF}}}$, the classicsl spin state can be approximated [24] as

$$\begin{eqnarray}&&{S}_{x^{\prime} }({\bf{R}})=S{(-1)}^{n+1}\mathrm{cos}\tau \mathrm{sin}(2\pi \delta r),\end{eqnarray} \tag{ 7 }$$

$$\begin{eqnarray}&&{S}_{y^{\prime} }({\bf{R}})=S\mathrm{sin}\tau \mathrm{sin}(2\pi \delta r),\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{S}_{z^{\prime} }({\bf{R}})=S{(-1)}^{n+1}\mathrm{cos}(2\pi \delta r),\end{eqnarray} \tag{ 9 }$$

so that the spins on each hexagonal layer depend only on the integer $r=\sqrt{2}{{\bf{x}}}^{\prime }\cdot {\bf{R}}/a$. The weak FM moment ${{\bf{M}}}_{0}=2{\mu }_{{\rm{B}}}{S}_{0}{\bf{y}}^{\prime}$ of the canted antiferromagnetic phase above H_c is related to the tilt by [21] $\mathrm{sin}\tau ={S}_{0}/S$. For [5, 25] ${M}_{0}=0.03\;{\mu }_{{\rm{B}}}$, $\tau =0.006$ or ${0.34}^{\circ }$. By comparison, the local spin-density approximation (LSDA)+U (U = 5 eV) calculations described in the next section yield ${M}_{0}=0.029\;{\mu }_{{\rm{B}}}$. Because higher harmonics are neglected, averages taken with the tilted cycloid in zero magnetic field introduce a very small error of order 10⁻⁵. Quantum fluctuations about the classical spin state are expected to be small for $S=5/2$ Fe³⁺ ions.

First-principles calculations were performed using density functional theory (DFT) from the VASP code within LSDA+U. The Hubbard U = 5 eV and the exchange ${J}_{{\rm{H}}}=0\;{\rm{eV}}$ parameters were optimized for Fe³⁺ in BiFeO₃ [26, 27]. We employed projector augmented wave potentials [28, 29]. To integrate over the Brillouin zone, we constructed a supercell made of 2 × 2 × 2 perovskite units (40 atoms, 8 f.u.) and a 3 × 3 × 3 Monkhorst–Pack k-points mesh. The DM interactions ${{\boldsymbol{D}}}_{{\rm{FE}}}$ and ${{\boldsymbol{D}}}_{{\rm{AF}}}$ were evaluated with 4 × 2 × 2 units (80 atoms, 16 f.u.) and a 1 × 3 × 3 Monkhorst–Pack mesh. The wave functions were expanded with plane waves up to an energy cutoff of 500 eV. To calculate exchange interactions (J_n), we applied four different magnetic configurations (G-AFM, C-AFM, A-AFM and FM). We estimated ${{\boldsymbol{D}}}_{{\rm{FE}}}$ and ${{\boldsymbol{D}}}_{{\rm{AF}}}$ by replacing all except four of the Fe³⁺ cations with Al³⁺ [26] in the 80 atom unit cell. As shown in table 1, the LSDA+U results are in excellent agreement with recent neutron-scattering measurements [22].

Table 1.  Calculated magnetic interaction parameters (meV) compared to spin model results fitted to neutron-scattering measurements [22]. ${D}_{{\rm{AF}}}$ splits into two components parallel (A = 0.042) and perpendicular (B = 0.075) to spin bond direction. The components A and B are explained in appendix A.2.

meV	J1	${D}_{{\rm{FE}}}$	${D}_{{\rm{AF}}}$	K
LSDA+U	−6.1	0.089	0.042, 0.075	3.5 × 10−3
Neutron	−5.3	0.103	0.064	4.1 × 10−3

After obtaining the exchange, DM, and SIA interactions, we calculate their derivatives with respect to an applied electric field parallel to a cartesian direction. A dielectric constant $\epsilon =90$ is used to estimate the SD polarizations when the electric field is perpendicular to the rhombohedral axis [30]. To simulate atomic displacements driven by the applied field E_α, we evaluate the lowest-frequency polar eigenvector from the dynamical matrix by forcibly moving the atoms incrementally from the ground state R3c structure. The resulting energy difference between the two structures is divided by the induced electric polarization ${P}_{\alpha }^{{\rm{ind}}}$. The major difference in the polar eigenvectors obtained from the dynamic and the force-constant matrices arises from the Fe–O–Fe bond angle. The eigenvectors of the dynamic matrix reduce the bond-angle while the eigenvectors of the force-constant matrix raises that angle (see appendix B). These opposing tendencies produce distinct ME behaviors for dynamic and static electric fields. Although this paper evaluates the SD polarizations for dynamic electric fields, the general formalism is applicable in both the static and dynamic limits.

The cross products ${{\bf{S}}}_{i}\times {{\bf{S}}}_{j}$ modulate the Fe–O–Fe bond angle and produce the SC polarizations [31]. These SC polarizations are simply obtained from electric-field derivatives of Dzyaloshinskii–Moriya interaction Hamiltonian. In BiFeO₃, FE and AF distortions generate SC polarizations ${{\bf{P}}}_{{\rm{FE}}}^{{\rm{SC}}}$ and ${{\bf{P}}}_{{\rm{AF}}}^{{\rm{SC}}}$ associated with the electric-field derivatives of the DM interactions ${{\boldsymbol{D}}}_{{\rm{FE}}}$ and ${{\boldsymbol{D}}}_{{\rm{AF}}}$. These are calculated using the procedure described in [20].

The first SC polarization is induced by the response of the FE distortion to an external electric field:

$$\begin{eqnarray}&&{P}_{\mathrm{FE},\gamma }^{{\rm{SC}}}=-\displaystyle \frac{1}{N}\displaystyle \frac{\partial {{ \mathcal H }}_{{\rm{FE}}}^{{\rm{SC}}}}{\partial {E}_{\gamma }}=-\displaystyle \frac{1}{N}\displaystyle \sum _{u,\langle i,j{\rangle }^{{\bf{u}}}}\displaystyle \frac{\partial {{\bf{D}}}_{{\rm{FE}}}^{u}}{\partial {E}_{\gamma }}\cdot ({{\bf{S}}}_{i}\times {{\bf{S}}}_{j}),\end{eqnarray} \tag{ 10 }$$

where $\langle i,j{\rangle }^{{\bf{u}}}$ is a sum over nearest neighbors with ${{\bf{R}}}_{j}-{{\bf{R}}}_{i}=a{\bf{u}}$ and ${\bf{u}}={\bf{x}}$, ${\bf{y}}$, or ${\bf{z}}$ cubic axis. The electric-field derivatives of the DM interactions ${{\bf{f}}}^{u\gamma }=\partial {{\boldsymbol{D}}}_{{\rm{FE}}}^{u}/\partial {E}_{\gamma }$ are given in appendix C and table 2. While the derivative ${{\bf{f}}}^{\alpha \alpha }$ of ${{\boldsymbol{D}}}_{{\rm{FE}}}^{\alpha }$ between spins ${{\bf{S}}}_{j}$ and ${{\bf{S}}}_{i}$ with ${{\bf{R}}}_{j}-{{\bf{R}}}_{i}$ parallel to the electric field is parallel to ${{\boldsymbol{D}}}_{{\rm{FE}}}^{\alpha }$, the derivative ${{\bf{f}}}^{\alpha \beta }$ ($\beta \ne \alpha$) of ${{\boldsymbol{D}}}_{{\rm{FE}}}^{\alpha }$ between spins with ${{\bf{R}}}_{j}-{{\bf{R}}}_{j}$ perpendicular to the electric field is perpendicular to ${{\boldsymbol{D}}}_{{\rm{FE}}}^{\alpha }$, as shown in figure 1.

Zoom In Zoom Out Reset image size

Table 2.  SD polarizations from ES, SC and SIA. Shown are the calculated (LSDA+U) electric-field derivatives of J₁, ${{\boldsymbol{D}}}_{{\rm{FE}}},{{\boldsymbol{D}}}_{{\rm{AF}}}$, and K. The upper left and right scripts denote the directions of the spin bond and electric field, respectively. ${f}_{\beta }^{\alpha \alpha }=-{f}_{\gamma }^{\alpha \alpha }$, ${f}_{\gamma }^{\alpha \beta }=-{f}_{\gamma }^{\beta \alpha }$, and ${a}_{\beta }^{\alpha \alpha }={a}_{\gamma }^{\alpha \alpha }$ by R3c symmetry as in appendix C. α, β, and γ are in ascending order so that ${\epsilon }_{\alpha \beta \gamma }=1$. All units are nC cm⁻².

	SC polarization from ${{\boldsymbol{D}}}_{{\rm{FE}}}$			SC polarization from ${{\boldsymbol{D}}}_{{\rm{AF}}}$		ES polarization from J1
	${f}_{\beta }^{\alpha \alpha }$	${f}_{\gamma }^{\alpha \beta }$	${f}_{\beta }^{\alpha \beta }$	${a}_{\alpha }^{\alpha \alpha }$ + 2${a}_{\beta }^{\alpha \beta }$	${a}_{\beta }^{\alpha \alpha }+{a}_{\alpha }^{\alpha \beta }+{a}_{\gamma }^{\alpha \beta }$	${C}_{{\rm{AF}}}$	${C}_{{\rm{FE}}}$
LSDA+U	9	17	14	17	−19	−250	−350
NDD	36	29	29	28	−7.2	—	—

In the lab reference frame $\{x,y,z\}$, regrouping terms for domain 2 with ${{\bf{x}}}^{\prime }=[1,0,-1]$ using equations(7)–(9) yields

$$\begin{eqnarray}&&{{\bf{P}}}_{{\rm{FE}}}^{{\rm{SC}}}=-\displaystyle \frac{1}{N}\;{\underline{{\rm{\Lambda }}}}^{{\rm{FE}}}\cdot \displaystyle \sum _{\langle i,j{\rangle }^{{\bf{x}}}}({{\bf{S}}}_{i}\times {{\bf{S}}}_{j}),\end{eqnarray} \tag{ 11 }$$

where

$$\begin{eqnarray}&&{\underline{{\rm{\Lambda }}}}^{{\rm{FE}}}=\left(\begin{array}{l}{{\bf{f}}}^{{xx}}\\ {{\bf{f}}}^{{xy}}\\ {{\bf{f}}}^{{xz}}\end{array}\right)-\left(\begin{array}{l}{{\bf{f}}}^{{zx}}\\ {{\bf{f}}}^{{zy}}\\ {{\bf{f}}}^{{zz}}\end{array}\right)=\left(\begin{array}{ccc}-h & f-g & -f\\ g & 2h & g\\ -f & f-g & -h\end{array}\right)\end{eqnarray} \tag{ 12 }$$

with $f={f}_{\beta }^{\alpha \alpha }$, $g={f}_{\beta }^{\alpha \beta }$, and $h={f}_{\gamma }^{\alpha \beta }$ given in table 2. Fits to the NDD [19] described in section 6 imply that g = h.

The second SC polarization alternates in sign due to the alternating AF rotations along [111]:

$$\begin{eqnarray}&&{P}_{{\rm{AF}},\gamma }^{{\rm{SC}}}=-\displaystyle \frac{1}{N}\displaystyle \frac{\partial {{ \mathcal H }}_{{\rm{AF}}}^{{\rm{SC}}}}{\partial {E}_{\gamma }}=-\displaystyle \frac{1}{N}\displaystyle \sum _{u,\langle i,j{\rangle }^{{\bf{u}}}}{(-1)}^{{n}_{i}}\displaystyle \frac{{{\boldsymbol{D}}}_{{\rm{AF}}}^{u}}{\partial {E}_{\gamma }}\cdot ({{\bf{S}}}_{i}\times {{\bf{S}}}_{j}).\end{eqnarray} \tag{ 13 }$$

The SC polarization components ${{\bf{a}}}^{u\gamma }=\partial {{\boldsymbol{D}}}_{{\rm{AF}}}^{u}/\partial {E}_{\gamma }$ are evaluated in table 2. While the derivative ${{\bf{a}}}^{\alpha \alpha }$ of ${{\boldsymbol{D}}}_{{\rm{AF}}}^{\alpha }$ between spins ${{\bf{S}}}_{i}$ and ${{\bf{S}}}_{j}$ with ${{\bf{R}}}_{j}-{{\bf{R}}}_{i}$ parallel to the electric field is nearly anti-parallel to ${{\boldsymbol{D}}}_{{\rm{AF}}}^{\alpha }$, the derivative ${{\bf{a}}}^{\alpha \beta }$ ($\beta \ne \alpha$) of ${{\boldsymbol{D}}}_{{\rm{AF}}}^{\alpha }$ between spins with ${{\bf{R}}}_{j}-{{\bf{R}}}_{i}$ perpendicular to the electric field is perpendicular to ${{\boldsymbol{D}}}_{{\rm{AF}}}^{\alpha }$, as shown in figure 1.

Appendix D shows that the SC polarization can be simplified as

$$\begin{eqnarray}&&{{\bf{P}}}_{{\rm{AF}}}^{{\rm{SC}}}=\displaystyle \frac{1}{\sqrt{3}N}\;{\underline{{\rm{\Lambda }}}}^{{\rm{AF}}}\cdot \displaystyle \sum _{u,\langle i,j{\rangle }^{{\bf{u}}}}{(-1)}^{{n}_{i}}({{\bf{S}}}_{i}\times {{\bf{S}}}_{j}),\end{eqnarray} \tag{ 14 }$$

where

$$\begin{eqnarray}{\underline{{\rm{\Lambda }}}}^{{\rm{AF}}}={{\bf{a}}}^{{xx}}+{{\bf{a}}}^{{yy}}+{{\bf{a}}}^{{zz}}=\left(\begin{array}{ccc}s & t & t\\ t & s & t\\ t & t & s\end{array}\right)\end{eqnarray} \tag{ 15 }$$

with $s={a}_{\alpha }^{\alpha \alpha }+2{a}_{\beta }^{\alpha \beta }$ and $t={a}_{\beta }^{\alpha \alpha }+{a}_{\alpha }^{\alpha \beta }+{a}_{\gamma }^{\alpha \beta }$ given in table 2.

The absence of an inversion center between neighboring spin sites induces the ES bond polarizations. Since the scalar product ${{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}$ is altered by external perturbations such as temperature, electric field, or magnetic field, FE and AF distortions each generates its own ES polarization.

For symmetric exchange couplings, ES is dominated by the response of the nearest-neighbor interaction J₁:

$$\begin{eqnarray}&&{{ \mathcal H }}_{{\rm{ex}}}=-\displaystyle \sum _{\langle i,j\rangle }{J}_{1}\ {{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}=-\displaystyle \sum _{u,\langle i,j{\rangle }^{{\bf{u}}}}{J}_{1}^{u}\ {{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}.\end{eqnarray} \tag{ 16 }$$

The two ES polarizations ${{\bf{P}}}_{{\rm{FE}}}^{{\rm{ES}}}$ and ${{\bf{P}}}_{{\rm{AF}}}^{{\rm{ES}}}$ are closely related to one another. The electric-field derivatives $\underline{{\rm{\Gamma }}}$ are given in the cubic coordinate system by

$$\begin{eqnarray}&&{P}_{{\rm{FE}},\alpha }^{{\rm{ES}}}=-\displaystyle \frac{1}{N}\displaystyle \frac{\partial {{ \mathcal H }}_{{\rm{ex}}}}{\partial {E}_{\alpha }}=\displaystyle \sum _{\beta }{{\rm{\Gamma }}}_{\alpha \beta }^{{\rm{FE}}}\;{W}_{1\beta },\end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray}{\underline{{\rm{\Gamma }}}}^{{\rm{FE}}}=\left(\begin{array}{ccc}{C}_{\parallel } & {C}_{\perp } & {C}_{\perp }\\ {C}_{\perp } & {C}_{\parallel } & {C}_{\perp }\\ {C}_{\perp } & {C}_{\perp } & {C}_{\parallel }\end{array}\right),\end{eqnarray} \tag{ 18 }$$

$$\begin{eqnarray}&&{W}_{1u}=\displaystyle \frac{1}{N}\displaystyle \sum _{\langle i,j{\rangle }^{{\bf{u}}}}{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j},\end{eqnarray} \tag{ 19 }$$

where ${C}_{\perp }=\partial {J}_{1}^{\beta }/\partial {E}_{\alpha }$ ($\beta \ne \alpha$) and ${C}_{\parallel }=\partial {J}_{1}^{\alpha }/\partial {E}_{\alpha }$ for spin bonds perpendicular and parallel to the electric field, respectively.

Because the AF octahedral rotation is perpendicular to ${\bf{z}}^{\prime}$, the ES polarization associated with AF rotations is also perpendicular to ${\bf{z}}^{\prime}$ with

$$\begin{eqnarray}&&{{\boldsymbol{P}}}_{{\rm{AF}}}^{{\rm{ES}}}={C}_{{\rm{AF}}}\;{\bf{z}}^{\prime} \times {{\bf{W}}}_{2},\end{eqnarray} \tag{ 20 }$$

$$\begin{eqnarray}&&{P}_{\mathrm{AF},\alpha }^{{\rm{ES}}}=\displaystyle \sum _{\beta }{{\rm{\Gamma }}}_{\alpha \beta }^{{\rm{AF}}}\;{W}_{2\beta },\end{eqnarray} \tag{ 21 }$$

$$\begin{eqnarray}&&{W}_{2u}=\displaystyle \frac{1}{N}\sum _{\langle i,j{\rangle }^{{\bf{u}}}}{(-1)}^{{n}_{i}}\;{{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j},\end{eqnarray} \tag{ 22 }$$

$$\begin{eqnarray}{\underline{{\rm{\Gamma }}}}^{{\rm{AF}}}=\left(\begin{array}{ccc}0 & -({C}_{\parallel }-{C}_{\perp }) & {C}_{\parallel }-{C}_{\perp }\\ {C}_{\parallel }-{C}_{\perp } & 0 & -({C}_{\parallel }-{C}_{\perp })\\ -({C}_{\parallel }-{C}_{\perp }) & {C}_{\parallel }-{C}_{\perp } & 0\end{array}\right).\end{eqnarray} \tag{ 23 }$$

Unlike W_1u, W_2u alternates in sign due to the opposite AF rotations on adjacent hexagonal layers.

The first ES polarization parallel to ${\bf{z}}^{\prime}$ with coefficient ${C}_{{\rm{FE}}}=2{C}_{\perp }+{C}_{\parallel }$ modulates the FE polarization that already breaks IS above ${T}_{{\rm{N}}}$. The second ES polarization perpendicular to ${\bf{z}}^{\prime}$ has coefficient ${C}_{{\rm{AF}}}={C}_{\perp }-{C}_{\parallel }$. The AF rotations affect the bonds between nearest-neighbor spins in the plane normal to ${\bf{z}}^{\prime}$ because each oxygen moves along the directions $[0,-1,1]$, $[1,0,-1]$, and $[-1,1,0]$, perpendicular to ${\bf{z}}^{\prime}$. Thus, the second ES polarization is associated with atomic displacements perpendicular to ${\bf{z}}^{\prime}$ and parallel to the AF rotation.

Figure 2 demonstrates the strong anisotropy in the response of magnetic exchange to an electric field. While ${C}_{\perp }$ arises from the change in Fe–O–Fe bond angle due to a polar distortion, ${C}_{\parallel }$ arises from bond contraction. As shown, ${C}_{\parallel }$ is much more sensitive to an electric field than ${C}_{\perp }$. Since the ME anisotropy ${C}_{{\rm{AF}}}={C}_{\parallel }-{C}_{\perp }$ produces an ES polarization, the AF rotation angle is changed by the spin ordering. In particular, the negative sign of ${C}_{{\rm{AF}}}=-250$ nC cm⁻² indicates that the rotation angle increases with the dot product ${{\bf{S}}}_{i}\cdot {{\bf{S}}}_{j}$ because oxygen atoms moving in the AF plane have a negative effective charge ${Z}_{{\rm{O}}}^{*}({\rm{DFT}})=-3.3e$.

Zoom In Zoom Out Reset image size

The anisotropic ES polarization components ${C}_{\perp }$ and ${C}_{\parallel }$ cooperatively induce the ES polarization along ${\bf{z}}^{\prime}$ under the IS broken by the FE polarization. In contrast to our previous study [20] on the response to a static electric field (${C}_{{\rm{FE}}}=215$ nC cm⁻²), we obtain a negative ${C}_{{\rm{FE}}}=-350$ nC cm⁻² in a dynamic electric field. Appendix B describes the different eigenvectors of the dynamic and force-constant matrices. While Fe moves upward with respect to oxygen in the static regime, Fe moves downward in the dynamic regime because its mass is much larger than that of oxygen. Therefore, a static ${\bf{E}}$ increases the bond angle of Fe–O–Fe (positive ${C}_{{\rm{FE}}}$) but a dynamic ${\bf{E}}$ decreases the bond angle (negative ${C}_{{\rm{FE}}}$) due to the Goodenough–Kanamori rules [32].

We recently predicted [20] that the static SD polarization of BiFeO₃ is about −3 μC cm⁻² along a cubic diagonal opposite to the FE polarization emerging below ${T}_{{\rm{C}}}$. The electronic plus atomic contribution to the SD polarization is −1.3 μC cm⁻² and the lattice-deformation contribution is −1.7 μC cm⁻², which were slightly underestimated (−1.0 and −1.3, respectively) in previous literature [33, 34]. The total SD polarization (−3 μC cm⁻²) is higher than observed in any other known multiferroic material [20].

The most stringent test yet for the microscopic model proposed above is its ability to predict the NDD which is the asymmetry ${\rm{\Delta }}\;\alpha (\omega )$ in the absorption $\alpha (\omega )$ of light when the direction of light propagation is reversed. The absorption of THz light is given by $\alpha (\omega )=(2\omega /c){\rm{Im}}N{(\omega )}_{{ij}}$, where [35, 36]

$$\begin{eqnarray}&&N{(\omega )}_{{ij}}\approx \sqrt{({\epsilon }_{{ii}}+{\chi }_{{ii}}^{{ee}}(\omega ))(1+{\chi }_{{jj}}^{{mm}}(\omega ))}\pm {\chi }_{{ji}}^{{me}}(\omega )\end{eqnarray} \tag{ 24 }$$

is the complex refractive index for a linearly polarized beam, ${\underline{\chi }}^{{ee}}$, ${\underline{\chi }}^{{mm}}$ and ${\underline{\chi }}^{{me}}$ are the dielectric, magnetic, and magnetoelectric susceptibility tensors describing the dynamical response of the spin system [13, 15, 17, 35] and $\underline{\epsilon }$ is the dielectric constant related to the charge response. Subscripts i and j are fixed by the electric ${\bf{e}}$ and magnetic ${\bf{h}}$ polarization directions, respectively. The second term, which depends on the light propagation direction and produces NDD, is separated from the mean absorption by writing $N{(\omega )}_{{ij}}=\bar{N}(\omega )\pm {\chi }_{{ji}}^{{me}}(\omega )$.

Summing over the spin-wave modes n at the cycloidal ordering wavevector ${\bf{Q}}$, ${\rm{\Delta }}\;\alpha (\omega )=(4\omega /c){\rm{Im}}{\chi }_{{ji}}^{{me}}(\omega )$ is given by

$$\begin{eqnarray}&&{\rm{\Delta }}\;\alpha (\omega )=\displaystyle \sum _{n}{A}_{n}\;\delta (\omega -{\omega }_{n}),\end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray}&&{A}_{n}={NX}{\omega }_{n}\;{\rm{Re}}\;{\rm{}}\{{\rho }_{n0}{\mu }_{0n}\},\end{eqnarray} \tag{ 26 }$$

$$\begin{eqnarray}&&{\rho }_{0n}=\langle 0| {{\bf{P}}}^{{\rm{SD}}}\cdot {\bf{e}}/{ \mathcal V }| {\text{}}n\rangle ,\ \ \ \ \ \ \ \ \ \ \ \ {\mu }_{0n}=\langle 0| {\bf{M}}\cdot {\bf{h}}/{\mu }_{{\rm{B}}}| {\text{}}n\rangle ,\end{eqnarray} \tag{ 27 }$$

where ${\bf{M}}=(2{\mu }_{{\rm{B}}}/N){\sum }_{i}{{\bf{S}}}_{i}$ is the magnetization, ${ \mathcal V }={a}^{3}$ is the volume per Fe site, ${{\bf{P}}}^{{\rm{SD}}}/{ \mathcal V }=({{\bf{P}}}^{{\rm{ES}}}+{{\bf{P}}}^{{\rm{SC}}})/{ \mathcal V }$ is the net SD polarization given in units of, nC cm⁻², and $X=(4\pi {\mu }_{{\rm{B}}}/{\hslash })\;\mathrm{nC}\;{{\rm{cm}}}^{-2}=0.1388\;{{\rm{cm}}}^{-1}$.

For each field orientation, the integrated weight of every spectroscopic peak at ${\omega }_{n}$ is compared with the measured values, thereby eliminating estimates of the individual peak widths. Experimental results for the NDD with field along ${\bf{m}}=[1,-1,0]$ are plotted in figure 3(a) for ${\bf{e}}=[1,-1,0]$. Fits to the NDD are based on the plotted 2, 4, 6, 8, 10, and 12 T data sets. For each data set, we evaluate the integrated weights for the eight modes [22] ${{\rm{\Psi }}}_{0}$, ${{\rm{\Phi }}}_{1}^{(1)}$, ${{\rm{\Psi }}}_{1}^{(1,2)}$, ${{\rm{\Phi }}}_{2}^{(1,2)}$, and ${{\rm{\Psi }}}_{2}^{(1,2)}$ between roughly 12 and 35 cm⁻¹.

Zoom In Zoom Out Reset image size

Comparing figures 3(a) and (b), the NDD for ${\bf{m}}=[1,-1,0]$ is dominated by the two sets of SC polarizations ${{\bf{P}}}_{{\rm{FE}}}^{{\rm{SC}}}$ and ${{\bf{P}}}_{{\rm{AF}}}^{{\rm{SC}}}$. Table 2 indicates that the fitting results are not significantly changed by including the ES polarizations. Figure 3(c) attempts to fit the experimental data using the ES polarizations alone. Clearly, the ES polarizations by themselves cannot produce the observed NDD.

Figures 3(a) and (b) indicate that the various components of the SC polarizations in BiFeO₃ are captured by first-principles calculations and that the NDD is not strongly affected by the ES terms This selectivity originates from the spin dynamics of this nearly collinear antiferromagnet. Due to the very small SIA on the S = 5/2 Fe³⁺ spins, each magnon mode can be described as the pure precession of Fe³⁺ spins: the oscillating component $\delta {{\bf{S}}}_{i}^{\omega }$ of the spin on site i is perpendicular to its equilibrium direction ${{\bf{S}}}_{i}^{0}$. Since neighboring spins in the long-wavelength spin cycloid of BiFeO₃ are close to collinear, a dynamic polarization is effectively induced by SC terms such as ${{\bf{S}}}_{i}^{0}\times \delta {{\bf{S}}}_{i+1}^{\omega }$. However, the dynamic polarization generated by ES terms ${{\bf{S}}}_{i}^{0}\cdot \delta {{\bf{S}}}_{i+1}^{\omega }$ is almost zero. The spin stretching modes observed in strongly anisotropic magnets [35, 38] do not appear in BiFeO₃.

Zoom In Zoom Out Reset image size

Recent work [37] explains the observed static polarization perpendicular to ${\bf{z}}^{\prime}$ by the ${\bf{z}}^{\prime} \times {{\bf{P}}}_{{\rm{FE}}}^{{\rm{SC}}}$ term proportional to h − f. Although the fitting and LSDA+U values for $h-f\equiv {f}_{\gamma }^{\alpha \beta }-{f}_{\beta }^{\alpha \alpha }\sim \pm 10$ nC cm⁻² in table 2 are an order of magnitude smaller than required by that work, keep in mind that the SC parameters given in table 2 were evaluated or fitted for a dynamic electric field.

Although DFT calculations underestimate the ME coefficients compared to the NDD fitting results in table 2, they nicely demonstrate which of the symmetry-allowed ME couplings are relevant and which are negligibly small. Combining the two methodologies allows a more unambiguous determination of these coupling parameters. There are several possible explanations for the difference between the results obtained from DFT calculations and the NDD fitting. First, a larger dielectric constant $\epsilon$ could produce better agreement between DFT and NDD since the SD polarizations are proportional to $\epsilon$ through

$$\begin{eqnarray}&&{P}_{\alpha }^{{\rm{SD}}}=\displaystyle \frac{\partial { \mathcal H }}{\partial {E}_{\alpha }}=\displaystyle \frac{\partial {P}_{\alpha }^{{\rm{ind}}}}{\partial {E}_{\alpha }}\displaystyle \frac{\partial { \mathcal H }}{\partial {P}_{\alpha }^{{\rm{ind}}}}\approx \epsilon \displaystyle \frac{\partial { \mathcal H }}{\partial {P}_{\alpha }^{{\rm{ind}}}}.\end{eqnarray} \tag{ 28 }$$

Second, higher-frequency polar modes not considered here also can affect NDD. Third, a smaller Hubbard U will increase the SD polarizations and improve the agreement with the experimental fits. Fourth, magnon modes were observed between $\nu =15$ and 40 cm⁻¹ while we calculated the SC coupling constants in the dynamical limit. The crossover frequency ${\omega }_{{\rm{c}}}$ between static and dynamical behavior lies between 0 and the polar phonon at $\omega =78$ cm⁻¹. If ${\omega }_{{\rm{c}}}$ lies in the middle of the measured frequencies, then the SC fitting parameters may differ from the dynamical couplings obtained from LSDA +U.

In order to study the ME couplings in complex multiferroic systems, first-principles calculations must be anchored to the right microscopic Hamiltonian. With two sets of SC polarizations derived from the two distinct structural distortions, BiFeO₃ is a good example of how our atomistic approach works for complex materials. This paper calculated only the ionic displacement contribution to the ME coupling which is typically larger than the purely electronic contribution [34, 39]. The lattice deformation contribution to the SD polarization was discussed in our previous work [20].

The higher-frequency polar modes contribute to the electric-field induced displacement. Their contributions are proportional to the mode strength ${Z}^{2}/{\omega }^{2}$, where Z is the mode effective charge and ω its frequency. From our dynamical matrix calculations, the mode strengths of the higher frequency modes are less than 30% smaller than the strength of the lowest mode. Therefore, the lowest mode makes dominates the electric-field induced polar displacement.

The advantages (large FE polarization, high ${T}_{{\rm{C}}}$, and high ${T}_{{\rm{N}}}$) of BiFeO₃ are also major obstacles to understanding the ME couplings that produce the SD polarizations below ${T}_{{\rm{N}}}$. Leakage currents and the preexisting large FE polarization at high temperatures have hampered magneto-capacitance measurements and hidden the SD polarizations. Although recent neutron-scattering measurements [10] imply a large ES polarization, most other ME polarizations are unknown. We show that NDD measurements combined with first-principles calculations based on a microscopic model reveal the hidden SC polarizations. In particular, this approach allows us to disentangle the delicate SC polarizations and the hidden ES polarizations associated with AF rotations. We envision that intrinsic methods such as NDD will reveal hidden ME couplings in many materials and rekindle the investigation of type-I multiferroics.

We acknowledge discussions with H Kim, E Bousquet, Nobuo Furukawa, S Miyahara, J Musfeldt, U Nagel, S Okamoto, S Bordács and T Rõõm. Research sponsored by the US Department of Energy, Office of Basic Energy Sciences, Materials Sciences and Engineering Division. IK was supported by the Hungarian Research Fund OTKA K 108918. JHL was in part supported by the year of 2015 Research Fund (1.150132.01) of the UNIST (Ulsan National Institute of Science and Technology). We also thank Hee Taek Yi and Sang-Wook Cheong for preparation of the BiFeO₃ sample.

A.1. FE-induced Dzyaloshinskii–Moriya interaction

Since the FE vectors ${{\boldsymbol{D}}}_{{\rm{FE}}}^{u}$ are given by (0, ${D}_{{\rm{FE}}}$, $-{D}_{{\rm{FE}}}$), ($-{D}_{{\rm{FE}}}$, ${D}_{{\rm{FE}}}$, 0), and (${D}_{{\rm{FE}}}$, $-{D}_{{\rm{FE}}}$, 0) between nearest spins along ${\bf{x}}$, ${\bf{y}}$, and ${\bf{z}}$, respectively, the FE-induced DM interaction can be transformed as:

$$\begin{eqnarray}{{ \mathcal H }}_{{\rm{FE}}}^{{\rm{SC}}} & = & \displaystyle \sum _{{{\bf{R}}}_{i},{{\bf{R}}}_{j}={{\bf{R}}}_{i}+a{\bf{u}}}{{\boldsymbol{D}}}_{{\rm{FE}}}^{u}\cdot ({{\bf{S}}}_{i}\times {{\bf{S}}}_{j})\\ & = & {D}_{1}\;\displaystyle \sum _{{{\bf{R}}}_{i},{{\bf{R}}}_{j}={{\bf{R}}}_{i}+a{\bf{u}}}({\bf{z}}^{\prime} \times {\bf{u}})\cdot ({{\bf{S}}}_{i}\times {{\bf{S}}}_{j}),\end{eqnarray} \tag{ A1 }$$

where ${D}_{1}={D}_{{\rm{FE}}}\approx 0.154$ meV is now larger by $\sqrt{2}$ than in previous work [21].

A.2. AF-induced Dzyaloshinskii–Moriya interaction

The AF interactions ${{\boldsymbol{D}}}_{{\rm{AF}}}^{u}$ along ${\bf{x}}$, ${\bf{y}}$, and ${\bf{z}}$ can be written

$$\begin{eqnarray}&&{{\boldsymbol{D}}}_{{\rm{AF}}}^{x}=B({\bf{y}}+{\bf{z}})+A{\bf{x}},\end{eqnarray} \tag{ A2 }$$

$$\begin{eqnarray}&&{{\boldsymbol{D}}}_{{\rm{AF}}}^{y}=B({\bf{z}}+{\bf{x}})+A{\bf{y}},\end{eqnarray} \tag{ A3 }$$

$$\begin{eqnarray}&&{{\boldsymbol{D}}}_{{\rm{AF}}}^{z}=B({\bf{x}}+{\bf{y}})+A{\bf{z}}.\end{eqnarray} \tag{ A4 }$$

For domain 2 with ${{\bf{x}}}^{\prime }=[1,0,-1]$,

$$\begin{eqnarray}{{ \mathcal H }}_{{\rm{AF}}}^{{\rm{SC}}} & = & \displaystyle \sum _{{{\bf{R}}}_{i},{{\bf{R}}}_{j}={{\bf{R}}}_{i}+a{\bf{u}}}\;{(-1)}^{{n}_{i}}{{\boldsymbol{D}}}_{{\rm{AF}}}^{u}\cdot ({{\bf{S}}}_{i}\times {{\bf{S}}}_{j})\\ & = & \sqrt{3}\displaystyle \sum _{{{\bf{R}}}_{i}}^{\prime }{\bf{z}}^{\prime} \cdot \{B\ {{\bf{S}}}_{r}\times ({{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{x}}}+2{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{y}}}+{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{z}}})\\ & & +A\ {{\bf{S}}}_{{{\bf{R}}}_{i}}\times ({{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{x}}}+{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{z}}})\}\\ & & +\displaystyle \sum _{{{\bf{R}}}_{i}}^{\prime }{\bf{y}}^{\prime} \cdot \{(B-A)\ {{\bf{S}}}_{{{\bf{R}}}_{i}}\times ({{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{x}}}-2{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{y}}}+{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{z}}})\}\\ & \approx & \sqrt{3}\displaystyle \sum _{{{\bf{R}}}_{i}}^{\prime }{\bf{z}}^{\prime} \cdot \{B\ {{\bf{S}}}_{{{\bf{R}}}_{i}}\times ({{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{x}}}+2{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{y}}}+{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{z}}})\\ & & +A\ {{\bf{S}}}_{{{\bf{R}}}_{i}}\times ({{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{x}}}+{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{z}}})\},\\ & \approx & \sqrt{3}(4B+2A)\displaystyle \sum _{{{\bf{R}}}_{i}}^{\prime }{\bf{z}}^{\prime} \cdot ({{\bf{S}}}_{{{\bf{R}}}_{i}}\times {{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{y}}}),\end{eqnarray} \tag{ A5 }$$

where the primed sum over ${{\bf{R}}}_{i}$ is restricted to either n_i odd or even hexagonal layers. Based on the tilted cycloid of equations (7)–(9) the ${\bf{z}}^{\prime}$ term dominates because ${{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{x}}}-2{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{y}}}+{{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{z}}}$ is of order ${\delta }^{2}\sim 2\times {10}^{-5}$.

Previously, the second DM term was written

$$\begin{eqnarray}{{ \mathcal H }}_{{\rm{AF}}}^{{\rm{SC}}} & = & {D}_{2}\displaystyle \sum _{{{\bf{R}}}_{i},{{\bf{R}}}_{j}={{\bf{R}}}_{i}+a{\bf{u}}}\;{(-1)}^{{n}_{i}}{\bf{z}}^{\prime} \cdot ({{\bf{S}}}_{i}\times {{\bf{S}}}_{j})\\ & = & 2\sqrt{3}{D}_{2}\displaystyle \sum _{{{\bf{R}}}_{i}}^{\prime }{\bf{z}}^{\prime} \cdot ({{\bf{S}}}_{{{\bf{R}}}_{i}}\times {{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{x}}}+{{\bf{S}}}_{{{\bf{R}}}_{i}}\times {{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{y}}}+{{\bf{S}}}_{{{\bf{R}}}_{i}}\times {{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{z}}})\\ & \approx & 6\sqrt{3}{D}_{2}\displaystyle \sum _{{{\bf{R}}}_{i}}^{\prime }{\bf{z}}^{\prime} \cdot ({{\bf{S}}}_{{{\bf{R}}}_{i}}\times {{\bf{S}}}_{{{\bf{R}}}_{i}+a{\bf{y}}}),\end{eqnarray} \tag{ A6 }$$

which also uses equations (7)–(9). Therefore, ${D}_{2}=(A+2B)/3=0.064$ meV, which is in excellent agreement with previous work [21].

As noted in section 4, ${{\rm{C}}}_{{\rm{FE}}}$ is negative from the eigenmode of the dynamic matrix while it is positive from the eigenmode of the force-constant matrix [20]. This difference originates from the opposite changes to the Fe–O–Fe bond angle: the bond angle increases in the static limit while it decreases in the dynamic limit.

Defining ${{\bf{f}}}^{u\gamma }=\partial {{\boldsymbol{D}}}_{{\rm{FE}}}^{u}/\partial {E}_{\gamma }$,

$$\begin{eqnarray}&&{{\boldsymbol{D}}}_{{\rm{FE}}}^{x}=(0,D,-D),\ \ \ \ \ \ {{\boldsymbol{D}}}_{{\rm{FE}}}^{y}=(-D,0,D),\ \ \ \ \ \ {{\boldsymbol{D}}}_{{\rm{FE}}}^{z}=(D,-D,0)\end{eqnarray} \tag{ C1 }$$

$$\begin{eqnarray}&&{{\bf{f}}}^{{xx}}=(0,f,-f),\ \ \ \ \ \ {{\bf{f}}}^{{yx}}=(-g,0,-h),\ \ \ \ \ \ {{\bf{f}}}^{{zx}}=(g,h,0),\end{eqnarray} \tag{ C2 }$$

$$\begin{eqnarray}&&{{\bf{f}}}^{{xy}}=(0,g,h),\ \ \ \ \ \ {{\bf{f}}}^{{yy}}=(-f,0,f),\ \ \ \ \ \ {{\bf{f}}}^{{zy}}=(-h,-g,0),\end{eqnarray} \tag{ C3 }$$

$$\begin{eqnarray}&&{{\bf{f}}}^{{xz}}=(0,-h,-g),\ \ \ \ \ \ {{\bf{f}}}^{{yz}}=(h,0,g),\ \ \ \ \ \ {{\bf{f}}}^{{zz}}=(f,-f,0),\end{eqnarray} \tag{ C4 }$$

where $f\equiv {f}_{\beta }^{\alpha \alpha },g\equiv {f}_{\beta }^{\alpha \beta }$, and $h\equiv {f}_{\gamma }^{\alpha \beta }$.

Defining ${{\bf{a}}}^{u\gamma }=\partial {{\boldsymbol{D}}}_{{\rm{AF}}}^{u}/\partial {E}_{\gamma }$,

$$\begin{eqnarray}&&{{\boldsymbol{D}}}_{{\rm{FE}}}^{x}=(A,B,B),\ \ \ \ \ \ {{\boldsymbol{D}}}_{{\rm{FE}}}^{y}=(B,A,B),\ \ \ \ \ \ {{\boldsymbol{D}}}_{{\rm{FE}}}^{z}=(B,B,A),\end{eqnarray} \tag{ C5 }$$

$$\begin{eqnarray}&&{{\bf{a}}}^{{xx}}=(a,b,b),\ \ \ \ \ \ {{\bf{a}}}^{{yx}}=(d,c,e),\ \ \ \ \ \ {{\bf{a}}}^{{zx}}=(d,e,c),\end{eqnarray} \tag{ C6 }$$

$$\begin{eqnarray}&&{{\bf{a}}}^{{xy}}=(c,d,e),\ \ \ \ \ \ {{\bf{a}}}^{{yy}}=(b,a,b),\ \ \ \ \ \ {{\bf{a}}}^{{zy}}=(e,d,c),\end{eqnarray} \tag{ C7 }$$

$$\begin{eqnarray}&&{{\bf{a}}}^{{xz}}=(c,e,d),\ \ \ \ \ \ {{\bf{a}}}^{{yz}}=(e,c,d),\ \ \ \ \ \ {{\bf{a}}}^{{zz}}=(b,b,a),\end{eqnarray} \tag{ C8 }$$

where $a\equiv {a}_{\alpha }^{\alpha \alpha },b\equiv {a}_{\beta }^{\alpha \alpha },c\equiv {a}_{\alpha }^{\alpha \beta },d\equiv {a}_{\beta }^{\alpha \beta }$, and $e\equiv {a}_{\gamma }^{\alpha \beta }$.

For domain 2 with ${{\bf{x}}}^{\prime }=[1,0,-1]$

$$\begin{eqnarray}&&{{ \mathcal T }}_{x}=\displaystyle \frac{1}{\sqrt{3}}{{ \mathcal T }}_{z^{\prime} }-\displaystyle \frac{1}{\sqrt{6}}{{ \mathcal T }}_{y^{\prime} }\ \ \ \ {{ \mathcal T }}_{y}=\displaystyle \frac{\sqrt{6}}{3}{{ \mathcal T }}_{y^{\prime} }+\displaystyle \frac{1}{\sqrt{3}}{{ \mathcal T }}_{z^{\prime} }\ \ \ \ {{ \mathcal T }}_{z}={{ \mathcal T }}_{x},\end{eqnarray} \tag{ D1 }$$

which use

$$\begin{eqnarray}&&{{ \mathcal T }}_{u}\equiv \displaystyle \frac{3}{N}\displaystyle \sum _{i}{(-1)}^{{n}_{i}}({{\bf{S}}}_{i}\times {{\bf{S}}}_{i+u}).\end{eqnarray} \tag{ D2 }$$

The SD polarization associated with ${{\boldsymbol{D}}}_{{\rm{AF}}}$ is

$$\begin{eqnarray}&&{P}_{x}^{{\rm{SC}}}={{\bf{a}}}^{{xx}}\cdot {{ \mathcal T }}_{x}+{{\bf{a}}}^{{yx}}\cdot {{ \mathcal T }}_{y}+{{\bf{a}}}^{{zx}}\cdot {{ \mathcal T }}_{z}\end{eqnarray} \tag{ D3 }$$

$$\begin{eqnarray}&&=\displaystyle \frac{1}{\sqrt{6}}(-{{\bf{a}}}^{{xx}}+2{{\bf{a}}}^{{yx}}-{{\bf{a}}}^{{zx}})\cdot {{ \mathcal T }}_{y^{\prime} }+\displaystyle \frac{1}{\sqrt{3}}({{\bf{a}}}^{{xx}}+{{\bf{a}}}^{{yy}}+{{\bf{a}}}^{{zz}})\cdot {{ \mathcal T }}_{z^{\prime} }\end{eqnarray} \tag{ D4 }$$

$$\begin{eqnarray}&&\approx \displaystyle \frac{1}{\sqrt{3}}({{\bf{a}}}^{{xx}}+{{\bf{a}}}^{{yy}}+{{\bf{a}}}^{{zz}})\cdot {{ \mathcal T }}_{z^{\prime} }.\end{eqnarray} \tag{ D5 }$$

Similarly

$$\begin{eqnarray}&&{P}_{y}^{{\rm{SC}}}\approx \displaystyle \frac{1}{\sqrt{3}}({{\bf{a}}}^{{xx}}+{{\bf{a}}}^{{yy}}+{{\bf{a}}}^{{zz}})\cdot {{ \mathcal T }}_{z^{\prime} },\end{eqnarray} \tag{ D6 }$$

$$\begin{eqnarray}&&{P}_{z}^{{\rm{SC}}}\approx \displaystyle \frac{1}{\sqrt{3}}({{\bf{a}}}^{{xx}}+{{\bf{a}}}^{{yy}}+{{\bf{a}}}^{{zz}})\cdot {{ \mathcal T }}_{z^{\prime} }.\end{eqnarray} \tag{ D7 }$$

So in the local frame

$$\begin{eqnarray}&&{P}_{{{\bf{x}}}^{\prime }}^{{\rm{SC}}}={P}_{{\bf{y}}^{\prime} }^{{\rm{SC}}}=0,\end{eqnarray} \tag{ D8 }$$

$$\begin{eqnarray}&&{P}_{{\bf{z}}^{\prime} }^{{\rm{SC}}}=\displaystyle \frac{1}{\sqrt{3}}({P}_{x}^{{\rm{SC}}}+{P}_{y}^{{\rm{SC}}}+{P}_{z}^{{\rm{SC}}})=\displaystyle \frac{1}{3}({{\bf{a}}}^{{xx}}+{{\bf{a}}}^{{yy}}+{{\bf{a}}}^{{zz}})\cdot {{ \mathcal T }}_{z^{\prime} }\end{eqnarray} \tag{ D9 }$$

plus a correction of order ${\delta }^{2}\sim 2\times {10}^{-5}$.

The polarization matrix used to evaluate the NDD is given by

$$\begin{eqnarray}{{\bf{a}}}^{{xx}}+{{\bf{a}}}^{{yy}}+{{\bf{a}}}^{{zz}}=\left(\begin{array}{ccc}a+2\;d & b+c+e & b+c+e\\ b+c+e & a+2\;d & b+c+e\\ b+c+e & b+c+e & a+2\;d\end{array}\right),\end{eqnarray} \tag{ D10 }$$

where $a+2d=17$ nC cm⁻² and $b+c+e=-19$ nC cm⁻² are obtained from first principles as given in table 2. ($a\equiv {a}_{\alpha }^{\alpha \alpha }=4.1$ nC cm⁻², $b\equiv {a}_{\beta }^{\alpha \alpha }=-21$ nC cm⁻², $c\equiv {a}_{\alpha }^{\alpha \beta }=-6.7$ nC cm⁻², $d\equiv {a}_{\beta }^{\alpha \beta }=6.4$ nC cm⁻², and $e\equiv {a}_{\beta }^{\alpha \beta }=8.9$ nC cm⁻²).