# Spin-Induced Polarizations and Non-Reciprocal Directional Dichroism of the Room-Temperature Multiferroic BiFeO^{1}^{1}1Copyright notice: This manuscript has been authored by UT-Battelle, LLC under Contract No. DE-AC05-00OR22725 with the U.S. Department of Energy. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a non-exclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this manuscript, or allow others to do so, for United States Government purposes. The Department of Energy will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan (http://energy.gov/downloads/doe-public-access-plan).

###### Abstract

A microscopic model for the room-temperature multiferroic BiFeO that includes two Dzyaloshinskii-Moriya interactions and single-ion anisotropy along the ferroelectric polarization predicts both the zero-field spectroscopic modes as well as their splitting and evolution in a magnetic field. Due to simultaneously broken time-reversal and spatial-inversion symmetries, the absorption of light changes as the magnetic field or the direction of light propagation is reversed. We discuss three physical mechanisms that may contribute to this absorption asymmetry known as non-reciprocal directional dichroism: the spin current, magnetostriction, and single-ion anisotropy. We conclude that the non-reciprocal directional dichroism in BiFeO is dominated by the spin-current polarization and is insensitive to the magnetostriction and easy-axis anisotropy. With three independent spin-current parameters, our model accurately describes the non-reciprocal directional dichroism observed for magnetic field along . Since some modes are almost transparent to light traveling in one direction but opaque for light traveling in the opposite direction, BiFeO can be used as a room-temperature optical diode at certain frequencies in the GHz to THz range. Our work demonstrates that an analysis of the non-reciprocal directional dichroism spectra based on an effective spin model supplemented by first-principles calculations can produce a quantitative microscopic theory of the magnetoelectric couplings in multiferroic materials.

###### pacs:

75.25.-j, 75.30.Ds, 75.50.Ee, 78.30.-j## I Introduction

BiFeO is the only material known to exhibit multiferroic behavior at room temperature. Because its ferroelectric (FE) transition temperature teague70 () K is significantly higher than its Néel transition temperature sosnowska82 () K, BiFeO is a type I multiferroic. Although the non-magnetic FE polarization lebeugle07 () C/cm is much larger than the magnetic contribution kadomtseva04 (); tokunaga10 (); park11 (); lee13 () induced by the distorted spin cycloid sosnowska82 (); lebeugle08 (); rama11a (); herrero10 (); sosnowska11 (), the magnetic domain distribution of BiFeO can be manipulated by an applied electric field lebeugle08 (); slee08 ().

A great deal has been learned about BiFeO since the first single crystals became available for inelastic neutron scattering jeong12 (); matsuda12 (); xu12 (), Raman scattering cazayous08 (); rov10 (), and THz spectroscopy talbayev11 (); nagel13 () measurements. It is now understood that two sets of interactions control the cycloid of BiFeO: two Dzyaloshinskii-Moriya (DM) interactions produced by broken inversion symmetry and a single-ion anisotropy sosnowska95 () (ANI) along the direction of the FE polarization . Whereas the DM interaction sosnowska82 () perpendicular to is responsible for the long 62 nm cycloidal period, the DM interaction kadomtseva04 (); ed05 (); pyatakov09 (); ohoyama11 () along is responsible for a small cycloidal tilt pyatakov09 (). Above the critical magnetic field , the cycloidal tilt develops into the weak ferromagnetic (FM) moment tokunaga10 (); park11 (); rama11b () of a G-type antiferromagnet (AF) that is isosymmetrically canted by an antiferrodistortive rotation (R[1,1,1]) of the structure ed05 ().

Inelastic neutron scattering measurements jeong12 (); matsuda12 (); xu12 () were used to extract the AF nearest- and next-nearest neighbor exchange interactions Hes () meV and meV between the Fe spins on the pseudo-cubic unit cell sketched in Fig.1(a) with lattice constant . However, those measurements lacked the sensitivity to resolve the ordering wavevectors on either side of the G-type AF wavevector at , where is inversely proportional to the cycloidal period . Recent neutron scattering measurements jeong14 () with higher precision were able to distinguish the two cycloidal ordering wavevectors and found that Hes () meV and meV. But even those measurements lacked the precision to obtain , which was set to zero.

By contrast, the frequencies of the spin-wave (SW) modes at the ordering wavevector can be precisely measured with Raman scattering cazayous08 () and THz spectroscopy talbayev11 (); nagel13 (). The parameters , , and were estimated by fitting the frequencies fishman13a () of the four observed zero-field THz modes. With no remaining adjustable parameters, that same model predicted fishman13b () the evolution and activation of the THz modes nagel13 () in a magnetic field along .

We now use this microscopic model to predict the asymmetry in the absorption of light when the direction of the magnetic field or, equivalently, the direction of light propagation is reversed. Called non-reciprocal directional dichroism (NDD), absorption asymmetry was first observed by Hopfield and Thomas hop60 () over 50 years ago in CdS. Much more recently, the precise symmetry requirements for NDD in magnetic materials were systematically investigated by Szaller et al. [szaller13, ]. Strong NDD is expected for the spin excitations of multiferroic materials when both time reversal and spatial inversion symmetries are broken by the spin state. Both the magnetic and electric components of THz radiation can excite SWs in multiferroic materials. The NDD exhibited by simultaneously electric- and magnetic-dipole active excitations has been extensively studied in BaCoGeO [kezsmarki11, ; bordacs12, ; miyahara11, ; kezsmarki14, ], SrCoSiO [kezsmarki14, ], CaCoSiO [kezsmarki14, ], GdTbMnO [takahashi13, ], and EuYMnO [takahashi12, ].

Because the cycloidal spin state is produced by the competition between DM, exchange, and ANI interactions, three distinct physical mechanisms can produce NDD in BiFeO: the spin current (SC) driven by the DM interactions, magnetostriction (MS) or the electric-field induced changes in the exchange interactions, and the electric-field induced changes in the ANI. Remarkably, the dynamical magnetoelectric coupling governing the NDD in BiFeO is dominated by the two sets of SC polarizations associated with and . Qualitatively, the SC dominates the magnetoelectric coupling in BiFeO because spin fluctuations are transverse to the almost collinear, cycloidal spin state . Since but (for nearby sites and ) and the ANI is extremely weak, spin fluctuations more strongly affect the SC-induced polarization than the MS- and ANI-induced polarizations.

As a fraction of the total light absorption at a given wavelength, NDD is most pronounced for a mode with fluctuations out of the cycloidal plane at 15.5 cm. At this wavenumber, BiFeO is almost transparent for light traveling in one direction but opaque for light traveling in the opposite direction. Therefore, BiFeO can be used as an optical diode that operates up to room temperature.

Despite the success of our model describing the NDD for magnetic field along , several questions remain open. Although our model predicts NDD to be absent for light propagating along , a static magnetic field along , and THz electric-field orientation or , weak NDD has been observed for a magnetic field along under those conditions. An optical misalignment of the THz electric- and magnetic-field vectors and may be responsible for this effect. In addition, the mean absorption (the absorption averaged over positive and negative magnetic fields) is not as accurately predicted by our model as the NDD.

This paper complements a recent work istun () that presents detailed experimental results for both the individual absorptions and the NDD. We have divided this paper into six sections. Section II presents our microscopic model and Section III presents the predicted mode frequencies. Section IV describes the three polarization mechanisms and presents results for the magnetization and polarization matrix elements, with symmetry relations provided by Local Spin-Density Approximation (LSDA)+ calculations. Results for the NDD are presented in Section V. Section VI contains a discussion and conclusion. While Appendix A summarizes the experimental details, Appendices B, C, and D treat the SC-, MS-, and ANI-induced polarizations, respectively. For convenience, the unit vectors used in this paper are given in Table I.

## Ii Microscopic model

In a magnetic field , the spin state and SW excitations of BiFeO are evaluated from the microscopic Hamiltonian

(1) |

where , , or connects with its nearest neighbor . Since the unit vector points along a cubic diagonal parallel to the FE polarization , the sum has the form proposed by Katsura et al. katsura05 (). The hexagonal layers normal to are separated by and are labeled by the integer . Consequently, the sum alternates sign from one hexagonal layer to the next. Notice that the local DM interactions and are, respectively, perpendicular and parallel to .

, , | Pseudo-cubic laboratory reference frame |
---|---|

Orientation of the electric polarization | |

along one of the cubic diagonals | |

, , | Cycloidal reference frame |

, , | Cycloidal reference frame for domain |

, , or | |

Orientation of the static magnetic field | |

Local single-ion ANI axis | |

Direction of light propagation | |

Orientation of the THz electric field | |

Orientation of the THz magnetic field |

There are eight possible orientations for along the four cubic diagonals. For every possible , the three magnetic domains have different and . When (all unit vectors in Table I are assumed normalized to 1), the possible orientations for the axis are , , and with corresponding . These three magnetic domains have cycloidal ordering wavevectors

(2) |

Hence, the ordering wavevectors for different domains are , , and . In terms of , the period of the cycloid in zero field is nm.

As mentioned above, the DM interactions and only couple nearest-neighbor sites. In a previous formulation fishman13a (); fishman13b () of this microscopic model, coupled next-neighbor sites within the same hexagonal layer. Due to the very long cycloidal period of BiFeO, the equilibrium and dynamical properties of these two Hamiltonians are the same up to errors of order . Specifically, earlier predictions for the SW mode frequencies fishman13a (); fishman13b () and critical magnetic field fishman13c () are unchanged. However, the earlier DM interaction is now multiplied by . Because the nearest-neighbor DM interactions are much larger than those between next-neighbor spins, the Hamiltonian above provides a close connection with recent first-principles calculations ed05 (); junun ().

Since the and terms in depend only on , is independent of the magnetic domain. For a specific domain , the first SC term can be written , where

(3) |

(4) |

where is a sum over nearest neighbors with . These relations assume that the spins on each hexagonal layer depend only on the integer . So for domain 2, . The cross products in Eqs.(3) and (4) couple spins with indices and on neighboring layers.

The second SC term proportional to can be written , where

(5) |

Like , also couples neighboring spins on neighboring layers.

The nearest- and next-nearest neighbor exchange interactions Hes () meV and meV were obtained from recent inelastic neutron scattering measurements jeong12 (); matsuda12 (); xu12 () between 5.5 meV and 72 meV. On the other hand, the small interactions , , and that control the cycloid can be obtained from THz spectroscopy measurements talbayev11 (); fishman13a () below 5.5 meV (44.3 cm) in zero magnetic field.

We have neglected the broken spatial symmetry between the exchange interactions due to the rhombohedral distortion. While all interactions must remain the same due to the rotational symmetry about , may reflect the rhombohedral distortion. For example, next-nearest neighbors separated by and may experience slightly different exchange interactions, denoted by and in Fig.1(a), because while . However, based on the excellent agreement between theory and experiment for the mode frequencies reported in Section III and because is already so small, we expect this exchange anisotropy to have a very minor effect on the NDD.

For a given set of interaction parameters, the spin state of BiFeO is obtained by minimizing the energy over a set of variational parameters fishman13b (). Fixing , where is an integer, the energy is minimized over the variational parameters on a unit cell with sites along and two hexagonal layers. The spin state on layer is assumed to be identical to the spin state on layer . The wavevector parameter is determined as a function of field by an additional minimization loop over . In zero field, and . We verify that the corresponding spin state provides at least a metastable minimum of the energy by checking that the classical forces on each spin vanish.

Ignoring the cycloidal harmonics produced by and but including the tilt pyatakov09 () produced by , the spin state in zero field can be approximated by

(6) | |||||

(7) | |||||

(8) |

This tilted cycloid is plotted in Fig.1(b). Averages over this state are readily performed using , , and . In zero field, averages over the tilted cycloid are fairly accurate because fishman13a () even harmonics like vanish and . Corrections to the averages are then of order .

For comparison, the spin state of the canted AF at zero field can be simply written in terms of the canting angle within the coordinate system as

(9) |

on hexagonal layer . Recall that fishman13a () where is the weak FM moment of the AF phase along above . Whereas susceptibility measurements tokunaga10 () indicate that , a recent neutron-scattering study rama11b () suggests that equivalent to . By contrast, LSDA+ ( eV) junun () gives , in agreement with the former experimental result. Note that is a linear function of and of .

We now adopt a different approach to estimate . The three parameters , , and are fixed by two conditions: the period of the cycloid must match the measured period and and the frequencies of the four predicted SW modes in zero field must match the measured frequencies fishman13a (). A third condition is provided by the dependence of the predicted critical field on . As shown in Fig.2, the measured critical field of 18.8 T for requires that Hes () , corresponding to or 0.45. While meV is virtually independent of , linearly increases with . Figure 2 indicates that increases almost quadratically with from a value of meV when . Corresponding to , we obtain meV and meV. A somewhat smaller value meV was given in Ref.[jeong14, ], which took .

With other parameters fixed and , a value of smaller than about 0.079 meV would stabilize a different canted AF phase above with spins tilted above and below the plane due to the dominant single-ion ANI. Hence, the coplanar AF phase of Eq.(9) is barely stabilized by the second DM interaction.

## Iii Spectroscopic mode frequencies

Using the parameters given above, predicts the evolution of the modes with magnetic field nagel13 (); fishman13b (); istun () for all orientations . SW modes at the ordering wavevector can be labeled deSousa08 () as in-cycloidal-plane modes and out-of-cycloidal-plane modes. In an extended zone scheme, those mode frequencies are plotted versus for wavevector in Fig.3(a). For simplicity, and denote both the modes and their frequencies. Neglecting higher spin harmonics, and . It follows that .

Higher harmonics generated by the tilt and ANI split each mode with into two labeled or . For the predicted parameters of BiFeO, those modes are plotted versus wavevector in Fig.3(b). While the modes are strongly affected by the spin harmonics, the former mode scheme remains fairly accurate for . Because the splitting of the low-frequency modes was not considered, recent Raman studies cazayous08 (); rov10 () misidentified the observed modes with some out-of-plane modes mistaken for in-plane modes and vice versa.

Despite the substantial splitting of and , is only slightly larger than . The nearly degenerate and modes cannot be separated by THz measurements talbayev11 (); nagel13 () in zero field.

In Fig.4, the predicted and measured missing () mode frequencies are plotted versus field for orientations , , and . Experimental data was not available for the THz modes above for the last two field orientations. The experimental results for and are presented here for the first time with experimental details summarized in Appendix A. Because its frequency was too low, was not detected when and . The predicted mode frequencies of the stable domain(s) are presented in the solid curves: domain 1 for and and domains 2 and 3 for . For , the mode that dips below arises from metastable domains 2 and 3, as seen by the agreement with the dashed curve. Hence, metastable domains may survive up to about 10 T.

With , the agreement between experiment and theory is even better than previously reported nagel13 () for with . Nevertheless, that agreement deteriorates somewhat above 12 T, particularly for , when avoided mode crossings strongly affect the mode frequencies. It is possible that the trial spin state is not sophisticated enough at high magnetic fields. For example, the spin state in high magnetic fields may have a periodicity greater than two hexagonal layers.

Above , the canted AF state of Eq.(9) supports only two modes that are labeled and in Fig.4. Because the transition at is first order, the spectroscopic modes change discontinuously at the critical field.

The estimates given above for , , and were based on fits to the four THz modes observed talbayev11 (); nagel13 () in zero field exc (). Experimental data points in Fig.4 indicate that those four modes correspond to / (nearly degenerate), , , and with frequencies 16.2, 20.7, 22.4, and 27.6 cm, respectively.

## Iv Polarization matrix elements

At zero magnetic field, only a few of the SW modes are optically active with finite magnetic-dipole resonance matrix elements , where

(10) |

is the magnetization operator, is the ground state with no SWs, and is the th excited state with a single SW mode at the cycloidal wavevector . At finite magnetic fields, all SW modes also have non-zero matrix elements of the induced electric polarization . The coexistence of the magnetic-dipole and polarization matrix elements is responsible for the NDD observed in the THz absorption spectra for field along . The physical mechanisms that contribute to below can be divided into three classes: SC, MS, and ANI.

For the SC- and MS-induced polarizations, we use LSDA+ calculations junun () to simplify the matrices connecting the induced polarizations with the spin operators. This greatly reduces the number of polarization parameters. In some instances, those matrices were simplified even further, either because some matrix elements were roughly equal or because additional matrix elements had a negligible effect on the NDD. Those additional simplifications are described in Ref.[junun, ]. This section expresses the induced polarizations in the cycloidal reference frame . In the laboratory reference frame , the induced polarizations are given in Appendices B, C, and D.

### iv.1 SC-induced polarizations

The SC-induced polarization is produced by shifts in the O locations due to the hopping of electrons between Fe 3 and O 2 orbitals katsura05 (); jia07 (); kaplan11 (). The first SC-induced polarization is created by the well-known inverse DM interaction katsura05 (); mostovoy06 (); sergienko06 () corresponding to the term in the Hamiltonian. This polarization can be generally written as

(11) |

where was defined by Eqs.(3) and (4). According to Eq.(B), the four nonzero matrix elements of are , , and , where the plus sign is for domain 2 and the minus sign is for domains 1 or 3.

In a simplified version of the first SC-induced polarization with , the diagonal terms and would vanish. Then and would reduce to the usual form katsura05 () for the inverse DM interaction:

(12) |

with , , and so that .

The second SC-induced polarization is associated with the DM interaction :

(13) |

where was defined by Eq.(5). As shown in Appendix B, the coefficient may differ from the and coefficients .

### iv.2 MS-induced polarizations

The first MS-induced polarization is produced by the uniform displacement of Fe with respect to O:

(16) |

(17) |

(18) |

It is easy to show that . For a simple twisted cycloid,

(19) |

The energy uniformly shifts all the nearest-neighbor interactions by .

The second MS-induced polarization can be written deSousa13 ()

(20) |

(21) |

Unlike , alternates sign from one hexagonal layer to the next. The cross product with in Eq.(20) ensures that remains a polar vector polar (). For a simple tilted cycloid in zero field, . The energy shifts the nearest-neighbor exchange interaction by an amount proportional to . For example, the nearest-neighbor exchange between spins at and is shifted by . Appendix C shows that .

The MS-induced polarization associated with next-nearest neighbor sites can be similarly constructed starting with

(22) |

(23) |

where all next-nearest neighbor pairs are double counted with , , and . So for , , , , and . For next nearest neighbors, both and lie on either even or odd layers.

Since , the polarizations associated with and are

(24) |

(25) |

For a simple twisted cycloid,

(26) |

while . The energy uniformly shifts all the next nearest-neighbor interactions by .

Another possible MS-induced polarization is associated with the spin exchange ANI or different exchange couplings for different spin components . Because it is of order , this polarization can be neglected.

### iv.3 ANI-induced polarizations

The ANI-induced polarization , which arises from the spin-dependent hybridization between the Fe ions and their ligands, contains components perpendicular or parallel to . As shown in Appendix D, the perpendicular polarization has two sets of terms associated with the electric-field dependence of the local single-ion ANI axis defined by Eq.(64). The first set is produced by the dependence of the polar angle on the electric field :

(27) |

which agrees with the first ANI-induced polarization proposed by deSousa et al. [deSousa13, ].

An additional perpendicular polarization

(28) |

with was proposed in Ref.[deSousa13, ]. However, the cross terms () in Eq.(68) cancel this contribution.

The second set of perpendicular ANI-induced terms is produced by the dependence of the azimuthal angle on :

(29) |

which was not previously proposed.

We also construct the ANI-induced polarization parallel to produced by the electric-field dependence of the constant :

(30) |

which shifts the single-ion ANI by . For a simple tilted cycloid in zero field, includes only a contribution from and is parallel to .

### iv.4 Total induced polarization

With all proposed terms, the net induced polarization in the cycloidal phase is . For the simple tilted cycloid,

(31) |

Of course, the components of perpendicular to do not change the magnitude of the total polarization significantly. The change in polarization from the paramagnetic phase above to the cycloidal phase below is given by Eq.(IV.4). Recently, Lee et al. [lee13, ] observed that has a magnitude of about 400 nC/cm and opposes due to the suppressed displacement of the Fe ions compared to the Bi ions.

By comparison, the induced polarization of the canted AF evaluated using Eq.(9) is given by

(32) |

which has no ANI contribution because the spins are in the plane. So the change in polarization from the AF phase to the cycloidal phase at zero field is given by

(33) |

Despite an early measurement of 1 nC/cm [kadomtseva04, ], the magnitude of the polarization change below extrapolated to zero field has recently been estimated as 40 nC/cm [tokunaga10, ; park11, ].

The Hamiltonian in zero electric field can be simply written in terms of the induced polarizations as

(34) |

Introducing the field dependence of the DM interactions, we find and . Similarly, , , and .

All components of the induced polarization appear in above. Because appears above in the paramagnetic phase, each static magnetically-induced polarization along corresponds to a term in the Hamiltonian. Due to the symmetry lowering associated with , each bilinear spin term that appears in also contributes to an induced polarization parallel to .

Taking in Eq.(IV.1),

(35) |

has no components perpendicular to . Components of the operator perpendicular to would then contribute only to the transition matrix elements . In other words, includes all induced polarizations with static contributions but not induced polarizations with only dynamical contributions . For example, does not appear in because .

We used Eq.(IV.4) to check our numerical results for the matrix elements . Since , the appropriate sum of polarization matrix elements with the field-dependent term must vanish when . We verified that this condition is indeed satisfied for all excited states and magnetic fields.

## V THz Absorption

The absorption of THz light is given by where miyahara11 (); miyahara14 ()

(36) |

is the complex refractive index for a linearly polarized beam, , and are the dielectric, magnetic, and magnetoelectric susceptibility tensors describing the dynamical response of the spin system kezsmarki11 (); bordacs12 (); miyahara11 (); kezsmarki14 () and is the background dielectric constant tensor associated with charge excitations at higher energies. Subscripts and refer to the electric and magnetic 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 SW modes at the cycloidal ordering wavevector , is given by

(37) |

(38) |

(39) |

(40) |

where is the volume per Fe site, is given in units of nC/cm and

(41) |

The THz electric and magnetic fields are polarized in the and directions, respectively.

After expanding for small susceptibilities, we find that is given by

(42) |

(43) |

where

(44) |

(45) |

Notice that .

The dielectric constant depends on the polarization of light. Based on a fit to the interference fringes, and 51.55 for and , respectively.

For each orientation of the static magnetic field and light polarization, the integrated weight of every spectroscopic peak at is compared with the measured values. This eliminates estimates of the individual peak widths. Because the polarization and magnetization matrix elements are generally complex with an arbitrary overall phase that differs for each mode , we can choose to be real. Other magnetization and polarization matrix elements for mode are then either purely real or imaginary. Under reversal of the field orientation, our numerical results indicate that