Spin-coupling-induced Improper Polarizations and Latent Magnetization in Multiferroic BiFeO3

Multiferroic BiFeO3 (BFO) that exhibits a gigantic off-centering polarization (OCP) is the most extensively studied material among all multiferroics. In addition to this gigantic OCP, the BFO having R3c structural symmetry is expected to exhibit a couple of parasitic improper polarizations owing to coexisting spin-polarization coupling mechanisms. However, these improper polarizations are not yet theoretically quantified. Herein, we show that there exist two distinct spin-coupling-induced improper polarizations in the R3c BFO on the basis of the Landau-Lifshitz-Ginzburg theory: ΔP LF arising from the Lifshitz gradient coupling in a cycloidal spin-density wave, and ΔP ms originating from the biquadratic magnetostrictive interaction. With the help of ab initio calculations, we have numerically evaluated magnitudes of these improper polarizations, in addition to the estimate of all three relevant coupling constants. We further predict that the magnetic susceptibility increases substantially upon the transition from the bulk R3c BFO to the homogeneous canted spin state in a constrained epitaxial film, which satisfactorily accounts for the experimental observation. The present study will help us understand the magnetoelectric coupling and shed light on design of BFO-based materials with improved multiferroic properties.

room-temperature exchange coupling between the BFO AFM order and the Co overlayer with ~90° in-plane Co-moment rotation upon single-step ferroelectric switching of the monodomain BFO. This has important consequences for practical, low power non-volatile ME devices utilizing BFO 20 .
In spite of extensive studies on the R3c BFO, however, possible ME coupling mechanisms and associated improper polarizations are not yet quantitatively resolved or lucidly explained. Herein, we show unequivocally that there exist two distinct spin-coupling-induced improper polarizations in the R3c BFO on the basis of the Landau-Lifshitz-Ginzburg phenomenological theory. These are: (i) a small parasitic improper polarization (ΔP LF ) originating from the Lifshitz exchange coupling, which can be equated with the S i × S j -type vector coupling induced polarization (ΔP DM ) caused by the reverse Dzyaloshinskii-Moriya (DM) interaction, and (ii) a S i · S j -type scalar coupling induced polarization (ΔP ms ) caused by the biquadratic magnetostrictive exchange interaction. With the help of ab initio density-functional theory (DFT) calculations, we have evaluated magnitudes of these improper polarizations, in addition to the numerical estimate of all three relevant ME coupling constants. We have further predicted that the magnetic susceptibility increases substantially upon the transition from the bulk R3c BFO to the homogeneous canted spin state in a constrained epitaxial film. These studies will help us comprehensively understand the ME coupling mechanisms in the R3c BFO and shed light on design of BFO-based materials with improved multiferroic properties.

Theoretical Analysis
Improper Polarization and Invariant caused by the Reverse DM Interaction. As mentioned previously, the R3c BFO is represented by a pseudo-cubic unit cell with its proper polarization along the cubic [111] c direction (i.e., [001] h in the hexagonal setting; Fig. 1). Let us define [111] c = [001] h as the z-direction. On the other hand, the incommensurate cycloidal spin structure with a periodicity of 620 Å suggests the appearance of a small improper polarization in the R3c BFO via the reverse Dzyloshinskii-Moriya (DM) coupling [21][22][23] . The spin-density wave (SDW) associated with this incommensurate spin cycloid is characterized by the propagation vector Q along the [110] h direction in the hexagonal setting 18,19,24,25 . Let us define [110] h as the x-direction (Fig. 1). We will first examine the magnitude of this reverse DM coupling-induced parasitic polarization which is directly linked to the spin cycloid before presenting the relevant thermodynamic potential of the bulk R3c BFO. The induced polarization (ΔP DM ) by the reverse DM interaction is expressed by the following form: where ê ij denotes a unit vector connecting the two neighboring magnetic (spin) moments, M i and M j , at the sites i and j, respectively. Thus, we have to evaluate the position-dependent . Thus, M i (≡M(0)) is given by On the other hand, the two neighboring canted magnetizations, m 1 and m 2 , at the site j, that are δ away from x = 0 are g ive n by , w he re ϕ ϕ ϕ ′ ≡ + Δ and ϕ ϕ ϕ ″ ≡ − Δ with the variation of the spin-canting angle (Δϕ) associated with the transla- y . This expression is readily transformed into the following form using elementary trigonometric algebra: Plugging Eqs (2) and (3) into Eq. (1) yields the following expression for the improper polarization induced by the reverse DM coupling: The last expression of P DM is valid if ϕ Δ ≈ cos( ) 1. As defined previously, ϕ Δ denotes the variation of the spin-canting angle associated with the translation from one Fe-site to the nearest-neighbor Fe-site along the x-axis. Thus, Δϕ = 360° × (5.58 Å/620 Å) = 3.24°, where 5.58 Å does correspond to the distance between the two nearest-neighbor Fe ions along the x-axis 17 . This predicts that ϕ Δ = .
≈ cos( ) 0 9984 1 and supports the validity of the last expression in Eq. (4). We will quantitatively examine the validity of this proposition of ignoring the y-component of the reverse DM coupling induced polarization in 'Discussion' section. It is interesting to notice that P DM is location-independent but depends on the DM coupling strength (d DM ), canting angle ϕ ( ), and the periodicity of SDW through ϕ Δ . Thus, the improper polarization induced by the reverse DM coupling does uniformly polarizes along the z-axis, i.e., along [001] h [ Fig. 2].
Within a continuum approximation for magnetic properties, the DM interaction responsible for this cycloidal modulation of spin moments in the R3c BFO can be expressed by inhomogeneous invariants, so-called 'Lifshitz  invariants' in the free-energy density 26 . Since the leading terms in the DM interaction are linear with respect to first spatial derivatives of magnetization in an antisymmetric mathematical form, the Lifshitz invariant for the R3c BFO can be written as where γ s and γ′ s denote the Lifshitz relativistic P-L exchange-coupling constants. Since the R3c BFO is antiferromagnetic (AFM), we used the magnitude of the AFM Néel vector, L, instead of the net magnetization order parameter, M, [≡ + m m ] 1 2 , in the description of Δf LF . Herein, the DM vector is replaced by γ P s since the polarization couples to gradients of M or L, thereby inducing an inhomogeneous cycloidal spin configuration. In the next section, we will use this form of Δf LF in constructing the thermodynamic potential for the R3c BFO.
Landau-Lifshitz-Ginzburg Thermodynamic Potential. Before presenting the free-energy density of the multiferroic R3c BiFeO 3 (BFO), we have first examined the free-energy density of the ferroelectric subsystem, Δf P ( ), on the basis of Landau-Ginzburg-Devonshire phenomenological theory [27][28][29] . As mentioned previously, the R3c BFO is represented by a pseudo-cubic unit cell with its proper polarization along the cubic [111] c direction (i.e., [001] h in the hexagonal setting; Fig. 1). Then, the free-energy density (thermodynamic potential) of the ferroelectric subsystem can be expanded on the basis of a paraelectric prototypic cell having cubic P m m 3 symmetry: where χ p denotes the dielectric stiffness, and ξ ξ ′ ″ , p p are high-order stiffness coefficients. In Eq. (6), the three polarization components along the three orthogonal cubic directions are denoted by P x ,P y , and P z . Then, the proper ferroelectric polarization (P) along the pseudo-cubic [111] c is given by = + + P P P P x y z 2 2 2 2 . Substituting this relation into Eq. (6) and suitably rearranging the resulting relation, one can obtain the following relation: where ξ‴ p is defined by ξ ξ ξ ≡ ″ − ′ . . Substituting this relation into Eq. (7), one can immediately obtain the following expression: Considering the Lifshitz invariant for the cycloidal modulation of spin moments [Eq. (5)] and the free-energy density for the ferroelectric subsystem [Eq. (8)], one can write down the Landau-Lifshitz-Ginzburg thermodynamic potential of the single crystalline R3c BFO in terms of two independent order parameters, P and L, where L is an AFM Néel vector describing the staggered sublattice magnetization. The model free-energy density Δf ( ) with respect to the paraphrase where where P denotes the magnitude of the total ferroelectric polarization (proper + improper) developed along the hexagonal c-axis, i.e., [001] h , or, equivalently, along [111] c of the pseudo-cubic unit cell ( Fig. 1) 13,14 . According to our Berry-phase calculations, P ≈ P z and P z is as high as ~90 μC/cm 2 for the undoped BFO having the R3c space-group symmetry 17  not independent of each other but are interrelated by M 2 − L 2 = 4(m 1 · m 2 ). Thus, the magnetostrictive coupling invariant in Eq. (9) can be rewritten using the following biquadratic form: Lifshitz Invariant associated with Cycloidal Spin-density Wave. As described previously, the single crystalline R3c BFO is characterized by the incommensurate SDW with the propagation vector Q along the [110] h spiral direction in hexagonal setting 18,19,24 . As schematically shown in Fig. 3, the x-location-dependent Néel vector [L(x)] forms a continuously varying cycloidal vector on the x-z plane with its propagation direction Q along x (=[110] h ). Thus, the Néel vector is given by In obtaining the last expression, we omitted terms representing the square of Néel-vector gradients. According to Mostovoy 33 , these nonlinear square terms do not practically contribute to the uniform improper polarization caused by the Lifshitz P-L exchange coupling.

it can be shown that the Lifshitz invariant is given by
In obtaining the last expression, we used the equality that since both L x and L z are independent of the z-coordinate [ Fig. 3] and that P x = 0 as the SDW-propagation direction (Q) is parallel to x. Combining this result with the two relations given in Eq. (10), one can obtain the following expression of the Lifshitz invariant: where L o (≡ L = |m 1 − m 2 |) denotes the magnitude of Néel vector and is given by R. de Sousa and J. E. Moore 24 used another form of the Lifshitz exchange coupling for the R3c BFO, which is apparently different from Δf LF presented in Eq. (11). It is given by After carrying out several steps for mathematical rearrangements, one can show that In obtaining the above relation, we used the equality that . As mentioned previously, P x = 0 since the SDW-propagation direction (Q) is parallel to x. Then, combining Eq. (15) with Eq. (10) immediately yields This indicates that the two apparently different forms of the Lifshitz invariant [i.e., Eqs (11) and (14)] are equal to each other under the condition of cycloidal spin ordering confined in an x-z plane with the translational symmetry along ŷ and the propagation direction along x which is parallel to [110] h .

Free-energy Minimization for Deducing Two Distinct Improper Polarizations.
Incorporating the results of Eq. (13) for the Lifshitz invariant into Eq. (9), we obtain a simplified form of the thermodynamic potential.
where we used the notation P for the proper off-centering polarization, P z , since P z ≈ P. Then, consider the free-energy functional for a finite volume, ΔF = ∫dvΔf(P,L). We impose the following equality for equilibrium: . Since the P-L cross-coupling is sufficiently weak, one can establish: Let us first consider the equilibrium off-centering (proper) polarization. One can immediately establish the following equality by imposing Eq. (18) to Eq. (17): If the Lifshitz-coupling term were absent, one would obtain the following relation from Eq. (19): Then, one immediately obtains the following expression for the equilibrium polarization (P o ) in the absence of the Lifshitz exchange coupling: where P eq(0) denotes the equilibrium (proper) off-centering polarization in the absence of any intrinsic ME coupling and ΔP ms (≡ ΔP ms ) represents a small improper polarization caused by the magnetostrictive exchange coupling. Thus, P o is comprised of two distinct terms, namely, = +Δ we establish the following relations from Eq. (20): , where T c denotes the ferroelectric Curie temperature (≈1100 K for the R3c BFO) 11,34 . Thus, χ p < 0 below T c . Let the total equilibrium polarization that satisfies Eq. (19) be P eq . We then establish In the above equation, ΔP LF (≡ ΔP LF ) appears due to the last term in the parenthesis of Eq. (19) and thus represents a small improper polarization arising from the Lifshitz coupling.
Owing to the Lifshitz invariant, one cannot obtain a correct analytic solution of P eq from Eq. (19). We thus treat ΔP LF as a small perturbation to P o and obtain a reasonably accurate solution for ΔP LF . Substituting Eq. As discussed previously, the terms inside the first parenthesis are equal to 0. Neglecting the term containing ΔP ( ) LF 2 , we eventually derive the following expressions of ΔP LF : This equality was previously discussed in conjunction with Eq. (19).
On the other hand, the improper polarization caused by the magnetostrictive exchange coupling (ΔP ms ) is given in Eq. The above equation demonstrates that the improper polarization caused by the magnetostrictive interaction belongs to S i · S j -type scalar coupling. According to Eq. (25), ΔP ms is inversely proportional to P o . In contrast, ΔP LF is inversely proportional to the square of P o . Since |cos θ| ≈ 1, the ratio of these two antiparallel improper polarizations is obtained from Eqs (24)  where P is nearly equal to P z or P o . The subscript 'q' appeared in Δf q emphasizes that the magnetostrictive interaction is represented by a biquadratic coupling of the form, P 2 M 2 or P 2 L 2 . According to φ 4 -expansion adopted in Eq. (8) [or Eq. (17)], the Landau coefficient ξ p should be positive. Thus, Eq. (26) tells us that the biquadratic P-M cross-coupling thermodynamically stabilizes the R3c BFO system if γ q < 0. On the contrary, the biquadratic magnetostrictive coupling destabilizes the system with a concomitant decrease in P o (i.e., ΔP ms < 0) if γ q > 0. According to the experimental result reported by S. Lee et al. 25 , the latter case (γ q > 0) is applicable to the R3c BFO. Our theoretical estimate also supports this conclusion ('Discussion' section). Since ΔP LF appeared in Eq. (24) is equal to ΔP DM that is given in Eq. (4), one can derive the following analytical expression for the reverse DM interaction coefficient (d DM ) which is a measure of the strength of the reverse DM coupling: where κ G denotes the Ginzburg gradient-energy coefficient. The last term inside the parenthesis of Eq. (28) can be rewritten as Then, we obtain the following type Euler-Lagrange equation from Eq. (28): The Ginzburg gradient term in the above equation is comprised of three distinct terms, namely, Substituting Eq. (32) into Eq. (28) yields the following expression for the equilibrium magnitude of the AFM Néel vector: Thus, L eq(0) denotes the equilibrium magnitude of the Néel vector in the absence of any coupling (i.e., γ q = γ s = κ G = 0). Then, the equilibrium magnetic remanence (M eq ) is related to L eq via the following relation: 2 (0) and M eq(0) denotes the equilibrium magnetic remanence in the absence of any coupling. Thus, the Lifshitz exchange coupling enhances M eq if γ < 0 s which corresponds to thermodynamically favorable Lifshitz coupling [Eq. (13)]. According to our theoretical estimate, γ s is indeed negative as described in 'Discussion' section. In contrast, the Ginzburg gradient term always suppresses M eq as κ > . 0 G The Ginzburg gradient energy can be computed by considering the space average of the gradient term in Eq.
It can be shown readily that On the contrary, (∇L y ) 2 = 0 since L y = 0. Thus, the space average of the Ginzburg gradient term is Since κ G > 0, the Ginzburg gradient term always increases the free-energy density.

Discussion
Estimate of the Two distinct Improper Polarizations. Having theoretically identified the two spin-coupling-induced improper polarizations in the R3c BFO (i.e., Δ Δ P P , ms LF ), we now focus on the numerical estimate of these values with the help of ab initio density-functional theory (DFT) calculations and experimental measurements. For this, let us first consider the improper polarization caused by the magnetostrictive exchange coupling, ΔP ms . As given in Eq. (20), this induced polarization is defined by Δ ≡ − P P P Let us estimate the improper polarization (ΔP LF ) arising from the Lifshitz gradient coupling. For this, we have to first consider the slowly varying spin reorientation with the periodicity of 620 Å along [110] h or x-axis. The distance between the two neighboring Fe-spin sites along the [110] h SDW propagation axis is 5.580 Å. Thus, the spin-rotation angle (Δφ) between the two neighboring Fe sites in the x-z plane is given by Δφ = 360° × (5.58/620) = 3.24°. We have imposed the spin-orientation structure and performed ab initio calculations by adopting a 2 × 2 × 1 supercell. The estimated ab initio value of ΔP LF is ~15 nC/cm 2 . However, this computed value is substantially smaller than the experimental value of 36 nC/cm 2 along [001] h 36 . Considering reliability of our ab initio value, we will adopt this experimental value in the evaluation of the Lifshitz coupling constant, γ s , in the next section. We schematically depict these two distinct improper polarizations with their directions in Fig. 4 and these can be summarized by the following ratio: −ΔP ms :+ΔP LF = 20:36 = 5:9.  ms LF and ΔP n in the R3c BFO. Among these three improper polarizations, ΔP LF is parallel to P eq(0) but ΔP ms is anti-parallel to P eq(0) . In contrast, ΔP n is perpendicular to P eq(0) . coefficient ξ p can be derived in terms of P eq and Δf ( ) , p eq where Δ ( ) f p eq denotes the equilibrium free-energy density of the ferroelectric subsystem with respect to that of the paraelectric reference system. This can be written explicitly as To estimate the Lifshitz exchange-coupling constant (γ s ), we have reconsidered Eq. (24). As mentioned previously, we have adopted the experimental value of 3.6 × 10 −4 (C/m 2 ) as ΔP LF . In addition to this, we use the following values for the estimate of γ s : P o = 0.863 (C/m 2 ), ξ = + . × ⋅ + J m C 1 175 10 ( / ),  . Polarization-electric field (P−E) hysteresis loop of the preferential [111] c -oriented BFO film (400-nm-thick) obtained at a measuring frequency of 1 kHz at 300 K. As shown, the net switching polarization (2P r ) of the [111] c -preferentially grown BFO film is ~180 μC/cm 2 , which indicates that the remanent polarization (P r ) is ~90 μC/cm 2 . According to the XRD pattern (Fig. 5), the in-plane film strain in the [111] coriented BFO layer is fully relaxed at the thickness of 400 nm. Thus, the measured P r value (~90 μC/cm 2 ) can be treated as a bulk P r value of the R3c BFO. We used the following values to estimate the above ratio: Δϕ = 3.24°, ϕ = 0.7°, d DM = +4.98 × 10 +44 (T 2 /J · V · m 2 ), and , where χ denotes the y-component spin-canting angle along the x-axis. In contrast, ϕ designates the x-component spin-canting angle along the z-axis [ Fig. 2]. According to our previous ab initio calculations 17 , χ = 0.203° along the x-axis, [110] h . Plugging all these values into Eq. (39), one obtains that ΔP DM (y)/ΔP DM (z) is as small as 1.65 × 10 −47 . This clearly justifies our previous proposition that the y-component value of the reverse DM coupling-induced polarization is completely negligible as compared with the corresponding z-component value.

Release of the Latent Magnetization in a Constrained Thin Film.
According to the study of Bai and co-workers 37 , the epitaxial film-constraint induces the destruction of a spatially modulated cycloidal spin structure in the bulk R3c BFO, releasing a latent AFM component locked within the cycloid. This corresponds to a transition from the incommensurately modulated cycloidal spin state to the homogeneously canted spin state in a constrained film with the onset thickness of ~150 nm 38 . Ryu and co-workers 38 further showed that the release of a latent AFM magnetization associated with the transition to the homogeneous spin state accompanies with a pronounced increase in the magnetic susceptibility (χ m ) in epitaxially constrained BFO thin films. We will quantitatively account for these observations by using Eq. (34). For an epitaxially constrained thin film having a homogeneous spin structure, we impose that γ = 0 s and κ G = 0 by considering the disappearance of a spatially modulated cycloidal spin structure 39 . For a constrained epitaxial thin film, Eq. (34) thus reads: where M eq f ( ) denotes the equilibrium magnetic remanence of the constrained epitaxial film. On the other hand, M eq(0) denotes the equilibrium magnetic remanence in the absence of any ME coupling. Thus, Using this relation and Eq. (34) for M eq b ( ) , one can derive the following relation for the latent magnetization released by the transition to the homogeneously canted spin state in a constrained film: where M eq(b) denotes the equilibrium magnetic remanence of the bulk R3c BFO and γ s < 0. In obtaining Eq. (40) we used the following equality: . Thus, we have estimated all four coupling constants needed for the Landau-Lifshitz-Ginzburg treatment of the R3c BFO.
The saturation magnetization (M s ) or magnetic susceptibility (χ M ), in general, can be readily estimated from the M-H hysteresis curve. On the contrary, ΔM eq is too small 30,38 to quantitatively discuss this effect in terms of γ s and κ G . Thus, we have examined a variation in χ M . In doing this, we first consider the thermodynamic potential of the AFM subsystem. One can write the following relation by exploiting Eq. (17) for Δf(L), Eq. (13)  The term inside a small bracket, 2γ s PQ + κ G Q 2 , is non-zero for a bulk crystal but is zero for a constrained epitaxial film. If this term is positive, the inverse susceptibility decreases or equivalently, M s increases upon the transition to the homogeneous spin state in a constrained epitaxial film. On the contrary, the reverse is true if this term is have further predicted that the magnetic susceptibility (χ m ) increases substantially upon the transition to the homogeneous spin state in a constrained epitaxial BFO film, which accounts for the experimental observation well.

Computational Methods
To obtain material parameters needed to quantitatively estimate three distinct polarizations and coupling constants, we have performed first-principles DFT calculations on the basis of the generalized gradient approximation (GGA) 40 and the GGA+U method 41 implemented with the projector augmented-wave (PAW) method 42 using the Vienna ab initio simulation package (VASP) 43,44 . The Hubbard U eff of 4.5 eV was chosen on the basis of empirical corrections. We explicitly treated five valence electrons for Bi (6s 2 6p 3 ), eight for Fe (3d 6 4s 2 ), and six for oxygen (2s 2 2p 4 ). Actual DFT calculations were performed using (i) a 4 × 4 × 3 Monkhorst-Pack k-point mesh 45 centered at Γ for the R3c structure, (ii) a 500 eV plane-wave cutoff, and (iii) the tetrahedron method with the Blöchl corrections for the Brillouin-zone integrations 46 . Structural optimizations were basically performed for the 30-atoms cell which corresponds to a hexagonal unit cell. In contrast, we adopted a 2 × 2 × 1 hexagonal supercell (containing 4 unit cells) to evaluate ΔP LF . The ions were relaxed until the Hellmann-Feynman forces on them were less than 0.01 eV/Å.