First-principles approach to the dynamic magnetoelectric couplings for the non-reciprocal directional dichroism in BiFeO3

Due to the complicated magnetic and crystallographic structures of BiFeO3, its magnetoelectric (ME) couplings and microscopic model Hamiltonian remain poorly understood. By employing a first-principles approach, we uncover all possible ME couplings associated with the spin-current (SC) and exchange-striction (ES) polarizations, and construct an appropriate Hamiltonian for the long-range spin-cycloid in BiFeO3. First-principles calculations are used to understand the microscopic origins of the ME couplings. We find that inversion symmetries broken by ferroelectric and antiferroelectric distortions induce the SC and the ES polarizations, which cooperatively produce the dynamic ME effects in BiFeO3. A model motivated by first principles reproduces the absorption difference of counter-propagating light beams called non-reciprocal directional dichroism. The current paper focuses on the spin-driven (SD) polarizations produced by a dynamic electric field, i.e. the dynamic ME couplings. Due to the inertial properties of Fe, the dynamic SD polarizations differ significantly from the static SD polarizations. Our systematic approach can be generally applied to any multiferroic material, laying the foundation for revealing hidden ME couplings on the atomic scale and for exploiting optical ME effects in the next generation of technological devices such as optical diodes.

static magneto-capacitance measurements are not feasible and to type-I multiferroics like BiFeO 3 where static magneto-capacitance measurements are often hindered by the large preexisting FE polarization. BiFeO 3 has two distinctive structural distortions that remove inversion centers and couple to the electric component of light. One is the FE distortion G - 4 [111] that breaks global inversion symmetry (IS). The other is the antiferroelectric (AF) octahedral rotation + R 4 [111] that breaks the local IS between nearest neighbor spins. Using a first-principles approach tied to a microscopic Hamiltonian, we demonstrate that all ME couplings are microscopically driven by distinct combinations of these two inherent structural distortions.
The first-principles approach described in this paper has already laid the foundation for two previous studies of BiFeO 3 . This approach was used [18,19] to predict the dynamic NDD observed in BiFeO 3 even at room temperature. As discussed in section 3, four spin-current (SC) polarizations~Ś S i j associated with FE and AF distortions cooperatively induce the strong NDD in BiFeO 3 . This approach was also used [20] to predict that the static SD polarization −3 μC cm −2 in BiFeO 3 points opposite to the preexisting FE polarization. A record high among all known multiferroics, this SD polarization is produced by the ES contribution~· S S i j discussed in section 4.

Microscopic spin-cycloid model for BiFeO 3
The FE polarization ¢ P z FE  emerging below T C can take eight different orientations along the four cubic diagonals ¢ = á 2 . FE and AF distortions create the DM interactions D FE and D AF . Including all magnetic anisotropies produced by these distortions, the spin Hamiltonian can be written where á ñ i j , and á ñ¢ i j , represent nearest and next-nearest neighbor spins, respectively. This is the most general Hamiltonian that includes the allowed distortions in R3c BiFeO 3 but neglecting exchange anisotropy terms of the form which are usually small for transition metal ions with half-filled d-shell such as Fe 3+ . Moreover, due to the long wavelength of the cycloid, exchange anisotropy can be effectively absorbed into the single-ion anisotropy (SIA) > K 0, which favors spin alignment along ¢ z . All terms in this Hamitlonian are also essential to explain the spin modes of BiFeO 3 observed using THz spectroscopy [21,22].
Since the FE distortion is uniform, the D FE sum is translation invariant. Due to the translation-odd + R 4 [111] AF octahedral rotation, the D AF sum contains the coefficient -( ) 1 n i , which alternates from one hexagonal layer = ¢ · n a z R 3 i i to the next. Simplified forms for the DM terms  FE SC and  AF SC are given in appendix A. By ignoring the higher cycloidal harmonics but including the tilt [23] τ produced by D AF , the classicsl spin state can be approximated [24] as so that the spins on each hexagonal layer depend only on the integer = ¢ · r a x R 2 . The weak FM moment m = ¢ S M y 2 0 B 0 of the canted antiferromagnetic phase above H c is related to the tilt by [21] t = S S sin 0 . For [5,25]

First-principles method
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 0 eV H parameters were optimized for Fe 3+ in BiFeO 3 [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 D FE and D 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 D FE and D AF by replacing all except four of the Fe 3+ cations with Al 3+ [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].
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 =  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 a P 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.

SC polarizations
The cross productsŚ S i 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 3 , FE and AF distortions generate SC polarizations P FE SC and P AF SC associated with the electric-field derivatives of the DM interactions D FE and D 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: j j perpendicular to the electric field is perpendicular to a D FE , as shown in figure 1. Table 1. Calculated magnetic interaction parameters (meV) compared to spin model results fitted to neutron-scattering measurements [22]. D 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.
In the lab reference frame { } x y z , , , regrouping terms for domain 2 with ¢ = - 1, 0, 1 using equations are evaluated in table 2. While the derivative aa a of a D AF between spins S i and S j with -R R j i parallel to the electric field is nearly anti-parallel to a D AF , the derivative ab a (b a ¹ ) of a D AF between spins with -R R j i perpendicular to the electric field is perpendicular to a D AF , as shown in figure 1.
Appendix D shows that the SC polarization can be simplified as  where

ES polarizations
The absence of an inversion center between neighboring spin sites induces the ES bond polarizations. Since the scalar product · S S i 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 1 : The two ES polarizations P FE ES and P AF ES are closely related to one another. The electric-field derivatives G are given in the cubic coordinate system by for spin bonds perpendicular and parallel to the electric field, respectively.
Because the AF octahedral rotation is perpendicular to ¢ z , the ES polarization associated with AF rotations is also perpendicular to ¢ z with Unlike W 1u , W 2u alternates in sign due to the opposite AF rotations on adjacent hexagonal layers. The first ES polarization parallel to ¢ z with coefficient  modulates the FE polarization that already breaks IS above T N . The second ES polarization perpendicular to ¢ z has coefficient The AF rotations affect the bonds between nearest-neighbor spins in the plane normal to ¢ z because each oxygen moves along the directions -[ ] 0, 1, 1 , -[ ] 1, 0, 1 , and -[ ] 1, 1, 0 , perpendicular to ¢ z . Thus, the second ES polarization is associated with atomic displacements perpendicular to ¢ z and parallel to the AF rotation. Figure 2 demonstrates the strong anisotropy in the response of magnetic exchange to an electric field. Whilê C arises from the change in Fe-O-Fe bond angle due to a polar distortion, C  arises from bond contraction. As shown, C  is much more sensitive to an electric field thanĈ . Since the ME anisotropy = -Ĉ C C AF  produces an ES polarization, the AF rotation angle is changed by the spin ordering. In particular, the negative sign of = -C 250 AF nC cm −2 indicates that the rotation angle increases with the dot product · S S i j because oxygen atoms moving in the AF plane have a negative effective charge = - The anisotropic ES polarization componentsĈ and C  cooperatively induce the ES polarization along ¢ z under the IS broken by the FE polarization. In contrast to our previous study [20] on the response to a static electric field ( = C 215 FE nC cm −2 ), we obtain a negative = -C 350 FE nC cm −2 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 E increases the bond angle of Fe-O-Fe (positive C FE ) but a dynamic E decreases the bond angle (negative C FE ) due to the Goodenough-Kanamori rules [32].
We recently predicted [20] that the static SD polarization of BiFeO 3 is about −3 μC cm −2 along a cubic diagonal opposite to the FE polarization emerging below T C . The electronic plus atomic contribution to the SD polarization is −1.3 μC cm −2 and the lattice-deformation contribution is −1.7 μC cm −2 , which were slightly underestimated (−1.0 and −1.3, respectively) in previous literature [33,34]. The total SD polarization (−3 μC cm −2 ) is higher than observed in any other known multiferroic material [20].

Origin of NDD
The most stringent test yet for the microscopic model proposed above is its ability to predict the NDD which is the asymmetry a w D ( )in the absorption a w ( ) of light when the direction of light propagation is reversed. The absorption of THz light is given by a w is the complex refractive index for a linearly polarized beam, c ee , c mm and c me are the dielectric, magnetic, and magnetoelectric susceptibility tensors describing the dynamical response of the spin system [13,15,17,35] and  is the dielectric constant related to the charge response. Subscripts i and j are fixed by the electric e and magnetic 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 Summing over the spin-wave modes n at the cycloidal ordering wavevector Q, a w w c w Im ji me is given by 1, 1, 0 is dominated by the two sets of SC polarizations P FE SC and P AF 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 3 are captured by first-principles calculations and that the NDD is not strongly affected by the ES terms This selectivity originates    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  could produce better agreement between DFT and NDD since the SD polarizations are proportional to  through 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 n = 15 and 40 cm −1 while we calculated the SC coupling constants in the dynamical limit. The crossover frequency w c between static and dynamical behavior lies between 0 and the polar phonon at w = 78 cm −1 . If w c lies in the middle of the measured frequencies, then the SC fitting parameters may differ from the dynamical couplings obtained from LSDA +U.

Discussion
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 3 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 w Z 2 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 C , and high T N ) of BiFeO 3 are also major obstacles to understanding the ME couplings that produce the SD polarizations below T 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. meV is now larger by 2 than in previous work [21].

A.2. AF-induced Dzyaloshinskii-Moriya interaction
The AF interactions D u AF along x, y, and z can be written where the primed sum over R i is restricted to either n i odd or even hexagonal layers. Based on the tilted cycloid of equations (7) which also uses equations (7)- (9). Therefore,