week ending 7 FEBRUARY 2014
PHYSICAL REVIEW LETTERS
PRL 112, 057202 (2014)
Prediction of a Novel Magnetoelectric Switching Mechanism in Multiferroics 1
Yurong Yang,1,2 Jorge Íñiguez,3 Ai-Jie Mao,4 and L. Bellaiche1
Physics Department and Institute for Nanoscience and Engineering, University of Arkansas, Fayetteville, Arkansas 72701, USA 2 Physics Department, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China 3 Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Spain 4 Institute of Atomic and Molecular Physics, Sichuan University, Chengdu 610065, China (Received 30 May 2013; published 5 February 2014) We report a first-principles study of the recently predicted Pmc21 phase of the multiferroic BiFeO3 material, revealing a novel magnetoelectric effect that makes it possible to control magnetism with an electric field. The effect can be viewed as a two-step process: Switching the polarization first results in the change of the sense of the rotation of the oxygen octahedra, which in turn induces the switching of the secondary magnetic order parameter. The first step is governed by an original trilinear-coupling energy between polarization, octahedral tilting, and an antiferroelectric distortion. The second step is controlled by another trilinear coupling, this one involving the predominant and secondary magnetic orders as well as the oxygen octahedral tilting. In contrast with other trilinear-coupling effects in the literature, the present ones occur in a simple ABO3 perovskite and involve a large polarization. DOI: 10.1103/PhysRevLett.112.057202
PACS numbers: 75.85.+t, 77.55.Nv, 77.80.Fm
Multiferroics have experienced a resurgence of interest because of their inherent cross coupling between electric and magnetic degrees of freedom [1–7]. In particular, a systematic control of the magnitude and crystallographic direction of magnetic order parameters by an electric field is attractive for the design of original devices, and fundamentally interesting. One way to attain such control is via the existence of an energy term ΔE that (1) directly couples magnetic moments and electric polarization P and (2) is linear in P. Familiar examples of such an energy term are specific Dzyaloshinskii-Moriya interactions [8–10] of the form ΔE ∼ P · ðL × MÞ [10,11], the Lifshitz invariants ΔE ∼ P · ½Lð∇ · LÞ þ L × ð∇ × LÞ and ΔE ∼ P · ½Mð∇ · MÞ − ðM · ∇ÞMÞ [12,13], where L and M are the antiferromagnetic vector and magnetization, respectively. Another example is the spin-current model [14–16] for which ΔE ∼ ðP × eij Þ · ðmi × mj Þ, where eij is the unit vector joining site i to site j and where mi and mj are the magnetic moments located at these sites i and j. The pioneering work of Ref. [17] also shows that there is another way to control magnetic quantities by electric polarization. This can be thought of as indirect since it consists of two steps. A trilinear “structural-only” coupling between polarization and two different oxygen octahedral tilting modes (also termed antiferrodistortive or AFD motion) first leads to the reversal of one octahedral tilt when switching the polarization. Then, this reversal can result in the change of direction of a magnetic order parameter (via a second energy that was not provided in Ref. [17] and that should couple magnetic moments and octahedral tilt). It is important to know if other two-step processes, allowing a systematic control of magnetic properties by electrical polarization, can exist. In particular, can the first energy also involve another kind of structural quantity, in addition to 0031-9007=14=112(5)=057202(5)
polarization and AFD motions, namely, antiferroelectricity? Another question to resolve is to determine the expression of the second energy that couples magnetic degrees and the polarization-switchable structural property. We provide answers to these issues via first-principles calculations. A novel two-step magnetoelectric effect is discovered. Its first energy is an original term that trilinearly couples polarization, one AFD distortion, and one antiferroelectric mode. This structural-only energy results in the change of the sense of the rotation of the oxygen octahedra when polarization is reversed. The second energy is another coupling involving predominant and secondary magnetic order parameters, as well as the AFD mode. It forces the secondary magnetic order parameter to switch its direction, once oxygen octahedral tilting is reversed by the change of sign of the polarization. The designation of the phenomenon as having two steps is made to indicate that two different energies are involved, even if the switchings of the two physical properties resulting from the existence of these two energies can be simultaneous. Density-functional calculations within the local spin density approximation are performed for BiFeO3 (BFO) systems, using the Vienna ab initio simulation package (VASP) [18]. As in Refs. [19–21], we also include a Hubbard correction for Fe ions, adopting the self-consistent value of U ¼ 3.8 eV [22,23] (note that this choice of U provided a good agreement between theory and experiment even for films under a large magnitude of strain; see, e.g., Ref. [24] and references therein). Spin-orbit and noncollinear magnetism are taken into account, and the projected augmented wave method and an energy cutoff of 500 eV are employed. We focus on the orthorhombic Pmc21 phase of (001) BFO films that has recently been found for large enough tensile strain [20] (note that this state has also been predicted in BaMnO3 , EuTiO3 , or CaTiO3 , and that increasing the temperature should lead to a
057202-1
© 2014 American Physical Society
PRL 112, 057202 (2014)
PHYSICAL REVIEW LETTERS
decrease of the critical strain at which this phase becomes the ground state). A 40-atom supercell is used to model it, along with a 3×3×3 k-point mesh. Its lattice vectors are a1 ¼2aIP x, a2 ¼ 2aIP y, and a3 ¼ aIP ½δ1 x þ δ2 y þ ð2 þ δ3 Þz. aIP is the in-plane lattice constant. It is chosen here to be 4.135 Å, which leads to a tensile strain of about 6% in BFO (note that BFO films experiencing strains having a magnitude even larger than 6% have been experimentally grown [25]). x, y, and z are unit vectors along the [100], [010], and [001] pseudocubic directions, respectively. The δ1 , δ2 , and δ3 variables and the atomic positions are relaxed to minimize the total energy. As sketched in Fig. 1, this Pmc21 state is characterized by (1) a polarization (P) lying along an in-plane h110i direction, (2) a Bi-driven antiferroelectric (AFE) vector (A) that is associated with the M point of the cubic, five-atom first ¯ Brillouin zone and that is oriented along an in-plane h110i direction that is perpendicular to the polarization, and (3) by an in-phase oxygen octahedral tilting (that is also associated with the M point of the first Brillouin zone) about the out-of-plane [001] direction. Such a tilting is represented by a pseudovector ω for which the magnitude provides the value of the tilting angle while its direction is either parallel or antiparallel to the out-of-plane direction—indicating if this tilting is counterclockwise versus clockwise, respectively. The superposition of P and A results in cations moving in a zigzag pattern [20], as shown in the top panels of Fig. 1. The polarization is calculated from the Bloch representation of the modern theory of polarization [26]. The strength of the antiferroelectricity vector is estimated by A¼ðdBi;A =dBi;P ÞP, where dBi;A repre¯ sents the atomic motion of Bi ions (along h110i) associated with antiferroelectricity, while dBi;P is the Bi displacement (along h110i) associated with the polarization P. Let us denote by þP and −P the cases in which P is parallel or antiparallel to the [110] direction, respectively. Similarly, once we choose a specific Bi site to define the origin of the AFE vector [27], þA and −A represent the AFE ¯ and [110], ¯ vectors along [110] respectively. Finally, once we select an Fe site to define the origin of the AFD motion [27], þω and −ω indicate that the oxygen octahedral tiltings (about a specific line joining Fe atoms along [001]) are counterclockwise and clockwise, respectively. When considering the eight cases resulting from the combination of the positive and negative signs of P, A, and ω, we found that the ground state of the Pmc21 state is fourfold degenerate; namely, the (þP, þA, þω), (þP, −A, −ω), (−P, −A, þω), and (−P, þA, −ω) all have the same energy, as well as identical magnitude for their P, A, and ω vectors. The other four states (−P, þA, þω), (−P, −A, −ω), (þP, −A, þω), and (þP, þA, −ω) are unstable and transform to one of the four stable equilibrium states during relaxation. These results suggest the existence of a very peculiar physical energy that couples P, A, and ω in a very precise manner. The Supplemental Material [28] provides the expression of this energy term: ΔE1 ¼ −C1 PAω; where C1 is a positive coupling parameter.
(1)
week ending 7 FEBRUARY 2014
FIG. 1 (color online). Schematization of properties in the equilibrium (þP, þA, þω) state (left panels) and (−P, þA, −ω) phase (right panels). The upper panels display the zigzag motions of Bi atoms in an (x, y) plane. The middle panels show (i) oxygen octahedral tiltings (represented by curled arrows), (ii) magnetic dipoles associated with the weak A-type antiferromagnetic vector (purple arrows) and with the predominant G-type AFM vector (represented by green arrows), and (iii) Bi displacements that are associated with AFE (light blue arrows). The lettering “Fe” (respectively, “Bi”) indicates the iron (respectively, bismuth) ion whose location serves as the origin for the definition of the ω, G, and AAFM (respectively, A) vectors [27]. The lower panels schematize the directions of P, the AFE vector, the AFD vector, the G-type AFM vector, and the weak A-type AFM vector.
As characteristics of ferroelectrics, the application of an electric field that is opposed to an initial polarization direction should switch the polarization. According to Eq. (1), this switch should make the oxygen octahedra rotate from counterclockwise to clockwise, if the AFE vector maintains its direction. Equation (1) further suggests that reverting this electric field should then allow the oxygen octahedra to tilt again in a counterclockwise manner and thus hints towards a systematic control of the oxygen octahedral tilting direction by an electric field. To test such a possibility, Fig. 2(a) reports the total energy resulting from the variation of the polarization from its positive equilibrium value of þ1.12 C=m2 to its opposite, negative equilibrium value (along ½1¯ 1¯ 0), and vice versa. Figures 2(b) and 2(c) show the AFE vector and out-of-plane, in-phase octahedral tilting resulting from the atomic relaxation occurring for each considered value of the polarization selected during the pathways (þP → −P) and (−P → þP). They reveal that one can indeed go from the equilibrium (þP, þA, þω) state to the equilibrium (−P, þA,
057202-2
PRL 112, 057202 (2014)
PHYSICAL REVIEW LETTERS
FIG. 2 (color online). Behavior of various quantities as a function of P, during the ðþP; þA; þωÞ → ð−P; þA; −ωÞ and ð−P; þA; −ωÞ → ðþP; þA; þωÞ pathways. (a),(b),(c) The total energy, AFE vector, and in-phase octahedral tilting, respectively. The insets of (b) schematize the evolution of the resulting Bi displacements along these paths, with the red, blue, and black arrows representing the polarization-induced, AFE-driven, and total Bi displacements, respectively. In (c), symbols represent data from first-principles calculations while the blue (respectively, red) solid line shows the fitting of these data by fourth order polynomials of P for positive (respectively, negative) ω.
−ω) phase (and vice versa) when continuously changing the magnitude and sign of the polarization. Such a transformation occurs via a pathway having an energy barrier of 0.36 eV per formula unit [29], along which the AFE vector does not vary too much. This implies, as shown in the insets of Fig. 2(b), that the overall cation displacement gradually evolves from a zigzag pattern (when A and P are close to each other in magnitude, that is, for a large magnitude of P) to a pure AFE pattern (which occurs when P vanishes). The overall symmetry of the structure always remains Pmc21 during the whole pathway, i.e., even when the polarization vanishes. The considered pathway is only one possibility out of the many that could occur in an actual material, but is very sound. Indeed, due to the persistence of the AFE distortion throughout, this path involves a quasicontinuous rotation of the electric dipoles associated with individual Bi cations, which are always in an off-centered, low-energy position. We also considered the transformation from (þP, þA, þω) to the equivalent state (−P, −A, þω) and found that this possibility, which involves the switching of the AFE vector, has a higherenergy barrier of 0.84 eV=f.u. associated with it. Figure 2 also shows that two solutions (P, þA, þω) and (P, þA, −ω) exist when P ranges from −0.60 to 0.60 C=m2 ,
week ending 7 FEBRUARY 2014
with (P, þA, þω) [respectively, (P, þA, −ω)] having a lower energy for positive (respectively, negative) polarization—as is consistent with Eq. (1). On the other hand, only one of these two solutions is found for larger magnitudes of P. Such a result can be understood by Eq. (1) that indicates that the difference in energy between the (P, þA, þω) and (P, þA, −ω) states is −2C1 PAω and thus increases with the magnitude of P. For large values of P, the state with higher energy thus becomes unstable, and only the state with lower energy [that is, (þP, þA, þω) for large positive P or (−P, þA, −ω) for negative P having a large magnitude] can persist. Figure 2(c) further shows that the behavior of the z component of the AFD vector ω as a function of P is quite complex, both when ω is positive or negative. As also shown in Fig. 2(c), such behaviors can be fitted by a polynomial offourth order in P (that is derived in the Supplemental Material [28]). Trilinear couplings involving P have also been proposed for the layered Ca3 Mn2 O7 [17], PbTiO3 =SrTiO3 superlattices [31], and strained CaTiO3 [32]. However, the other two quantities appearing in the trilinear couplings of these systems are either two different AFD motions [17,31] or two different antipolarmodes[32].TheproposedEq.(1)isthereforeoriginal since three different physical quantities (polarization, AFE, and AFD motions) are intertwined. Moreover, the polarizations involved in the trilinear couplings of Refs. [17,31,32] are small (about 10−6 to 0.2 C=m2 ), while the polarization of our equilibrium Pmc21 phase is much larger (namely, 1.12 C=m2 ). Moreover, our footnote or Ref. [33] provides a discussion about the nature of this latter polarization (that is, proper versus improper). Let us now investigate if the polarization-induced switching between (þP, þA, þω) and (−P, þA, −ω) can have some consequence on magnetism. For that, let us first choose a G-type antiferromagnetic (AFM) vector G being aligned along [110] to demonstrate magnetoelectric switching [35]. Our calculations predict that, in that case, a weak A-type AFM ¯ vector (to be denoted by AAFM ) is created along ½110 if the system is in the (þP, þA, þω) equilibrium state, while it is ¯ direction if the system is in the (−P, along the opposite ½110 þA, −ω) equilibrium state (see Fig. 1) [27]. One should thus be able to reverse the direction of the weak AFM vector AAFM when applying an external electric field [since such a field should induce the switching between (þP, þA, þω) and (−P, þA, −ω) states]. To shed more light on such an effect, Fig. 3 shows the evolution of the projection of AAFM onto ¯ for the paths considered in Fig. 2, that is, as a function of ½110 P varying between þ1.12 and −1.12 C=m2 . This projection is denoted as AAFM . Continuously decreasing P from its largest studied positive value to its largest-in-magnitude negative value leads to the following behaviors: (i) the negative AAFM first decreases in strength until P is close to þ0.75 C=m2 ; (ii) the negative AAFM then increases in magnitude until exhibiting a maximum for small polarizations; and (iii) the AAFM remains negative and decreases again in magnitude to suddenly jump for a negative P that is in between −0.6 and −0.8 C=m2 and then adopts a positive
057202-3
PRL 112, 057202 (2014)
PHYSICAL REVIEW LETTERS
week ending 7 FEBRUARY 2014
confirm such an exciting effect. Moreover, the discovered two-step control of the secondary magnetic order parameter is not restricted to BFO, since Eqs. (1) and (2) emerge from symmetry considerations that are not material specific.
FIG. 3 (color online). Evolution of the A-type AFM vector AAFM as a function of P. The square symbols show the values predicted from first principles. The solid lines represent the fitting of these data by AAFM ¼ −0.00051Gω þ 0.059GAPþ 0.00011GωP2 − 0.064GAP3 , where G is the G-type AFM vector (that has a magnitude of around 4μB ), A is the polarization-dependent AFE vector displayed in Fig. 2(b), and ω is the in-phase oxygen octahedral tilting [that also depends on P, as shown in Fig. 2(c)].
value for P less than −0.8 C=m2 . The change in sign of AAFM occurs precisely when ω also reverts its sign [see Fig. 2(c)] along the equilibrium (þP, þA, þω) → equilibrium (−P, þA, −ω) pathway. Figure 3 also shows that AAFM behaves in a mirroring fashion when studying the reverse ð−P; þ A; −ωÞ → ðþP; þ A; þ ωÞ pathway. Throughout these two pathways, we numerically found that the G-type antiferromagnetic vector not only remains large in magnitude (close to 4 bohr magneton) but also stays oriented along [110]. To understand these effects, one has to recall that Ref. [36] suggested an energy ΔE2 ¼ −C2 ω · ðG × AAFM Þ;
(2)
where C2 is a positive coefficient. In other words, ω, G, and AAFM want to form a direct triad, as schematized in Fig. 1. Such an equation explains why the combined existence of a G-type AFM vector and in-phase oxygen octahedral tilting results in a weak A-type AFM vector. It also explains why the polarization-induced switching of ω [such a switch is driven by Eq. (1)] leads to a reversal of AAFM (see Fig. 1). The control of the AAFM vector by an electric field is thus made possible via a two-step process [described by Eqs. (1) and (2)], both of which involve AFD motions. Such a two-step process represents a novel magnetoelectric effect. As is also detailed in the Supplemental Material [28], the nontrivial behavior of AAFM as a function of polarization (Fig. 3) can be reproduced from the energy terms discussed above. The proposed control of magnetism by polarization is not restricted to the switching of the A-type AFM vector. For instance, let us imagine a Pmc21 state for which the predominant magnetic ordering is of C-type. Reference [36] tells us that Eq. (2) would have to be replaced by ΔE2 ¼ −C2 ω · ðCAFM × FÞ, where CAFM is the C-type AFM vector and F is the (weak) ferromagnetic vector. Combining this latter equation with Eq. (1) implies that a net magnetic moment should also be controllable by polarization, in that case. Further simulations we conducted indeed
This work is supported by DOE grant ER-46612. We also acknowledge NNSF of China (Grants No. 51032002, No. 11374162, and No. 11104190) and 863 Program of China (Grant No. 2011AA050526), and MINECO-Spain (Grants No. MAT2010-18113 and No. CSD2007-00041). Some computations were also made possible thanks to the MRI grant 0722625, MRI-R2 grant 0959124, and CITRAIN grant 0918970 from NSF, the ONR grant N0001407-1-0825 (DURIP) and a Challenge grant from the Department of Defense.
[1] M. Fiebig, J. Phys. D 38, R123 (2005). [2] Z. J. Huang, Y. Cao, Y. Y. Sun, Y. Y. Xue, and C. W. Chu, Phys. Rev. B 56, 2623 (1997). [3] T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, and Y. Tokura, Nature (London) 426, 55 (2003). [4] T. Lottermoser, T. Lonkai, U. Amann, D. Hohlwein, J. Ihringer, and M. Fiebig, Nature (London) 430, 541 (2004). [5] C. Ederer and N. A. Spaldin, Phys. Rev. B 74, 020401 (2006). [6] C. Ederer and C. J. Fennie, J. Phys. Condens. Matter 20, 434219 (2008). [7] C. J. Fennie, Phys. Rev. Lett. 100, 167203 (2008). [8] I. Dzyaloshinsky, J. Phys. Chem. Solids 4, 241 (1958). [9] T. Moriya, Phys. Rev. Lett. 4, 228 (1960). [10] J. F. Scott and R. Blinc, J. Phys. Condens. Matter 23, 113202 (2011). [11] S. Prosandeev, I. A. Kornev, and L. Bellaiche, Phys. Rev. B 83, 020102(R) (2011). [12] I. Sosnowska and A. Zvezdin, J. Magn. Magn. Mater. 140– 144, 167 (1995). [13] M. Mostovoy, Phys. Rev. Lett. 96, 067601 (2006). [14] H. Katsura, N. Nagaosa, and A. V. Balatsky, Phys. Rev. Lett. 95, 057205 (2005). [15] D. Rahmedov, D. Wang, J. Íñiguez, and L. Bellaiche, Phys. Rev. Lett. 109, 037207 (2012). [16] A. Raeliarijaona, S. Singh, H. Fu, and L. Bellaiche, Phys. Rev. Lett. 110, 137205 (2013). [17] N. A. Benedek and C. J. Fennie, Phys. Rev. Lett. 106, 107204 (2011). [18] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). [19] S. Prosandeev, I. A. Kornev, and L. Bellaiche, Phys. Rev. Lett. 107, 117602 (2011). [20] Y. Yang, W. Ren, M. Stengel, X. H. Yan, and L. Bellaiche, Phys. Rev. Lett. 109, 057602 (2012). [21] Y. Yang, M. Stengel, W. Ren, X. H. Yan, and L. Bellaiche, Phys. Rev. B 86, 144114 (2012). [22] V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein, J. Phys. Condens. Matter 9, 767 (1997). [23] I. A. Kornev, S. Lisenkov, R. Haumont, B. Dkhil, and L. Bellaiche, Phys. Rev. Lett. 99, 227602 (2007). [24] Z. Chen, S. Prosandeev, Z. L. Luo, W. Ren, Y. Qi, C. W. Huang, L. You, C. Gao, I. A. Kornew, T. Wu, et al., Phys. Rev. B 84, 094116 (2011).
057202-4
PRL 112, 057202 (2014)
PHYSICAL REVIEW LETTERS
[25] R. J. Zeches, M. D. Rossell, J. X. Zhang, A. J. Hatt, Q. He, C.-H. Yang, A. Kumar, C. H. Wang, A. Melville, C. Adamo, et al., Science 326, 977 (2009). [26] R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993). P [27] The AFE vector A is defined as A¼1=N j dBi;A;j ð−1Þlj þmj , where the sum runs over all the N Bi sites j, and where dBi;A;j is the displacement of the Bi ion at site j that is induced by antiferroelectricity. The site j is located at alat ðlj x þ mj y þ nj zÞ from a Bi site chosen as the origin, where alat is the fiveatom cubic lattice constant and where lj , mj , and nj are all integers. Similarly, the in-phase AFD vector ω is defined as P ω ¼ 1=N k ωk ð−1Þlk þmk , where the sum runs over all the N Fe sites k, and where ωk is the vector characterizing the tilting centered on the Fe site k. The latter site is located at alat ðlk x þ mk y þ nk zÞ from an Fe site selected as the origin, with lk , mk , and nk being integers. The specific Bi and Fe atoms that have been chosen as the origins to define A and ω, respectively, are indicated in the middle panels of Fig. 1. Note that thePAFM vectors AAFM andP G are defined by AAFM ¼1=N k μk ð−1Þnk and G ¼ 1=N k μk ð−1Þlk þmk þnk , where μk is the magnetic moment at Fe site k. [28] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.112.057202 for additional details and derivations related to the energy of Eq. (1), as well as, to the dependence of the antiferrodistortive vector and (weak) A-type antiferromagnetic vector on polarization. [29] Note that dividing this barrier of 0.36 eV per formula unit by the spontaneous electric dipole moment (that is, the product between the polarization and the unit volume) of the investigated system in its equilibrium phase yields an electric field of 0.8 × 109 V=m. Such a field has the same order of magnitude as the one recently experimentally applied in BFO films (that is, 109 V=m) [30]. This implies that the barrier of 0.36 eV per formula unit can be overcome by a realistic electric field. Moreover, this barrier should decrease when increasing the temperature, rendering the electric field needed to overcome the barrier even smaller and easier to experimentally achieve at finite temperature. [30] H. Yamada et al., ACS Nano 7, 5385 (2013). [31] E. Bousquet, M. Dawber, N. Stucki, C. Lichtensteiger, P. Hermet, S. Gariglio, J.-M. Triscone, and P. Ghosez, Nature (London) 452, 732 (2008). [32] Q. Zhou and K. M. Rabe, arXiv:1306.1839v1. [33] Computing the total energy as a function of polarization, once starting from the P4=mmm phase (which is obtained from the ideal cubic structure by imposing the epitaxial constraint) and setting the AFE and AFD vectors to both vanish, constitutes what we will call “case 1.” It leads to a symmetric double well with an energy gain of 0.78 eV per formula unit, with the energy-versus-polarization curve being nicely fitted by a polynomial solely involving the second and fourth powers of the polarization. This is the strongest instability of the reference structure, since letting only the AFE (respectively, AFD) motions occur results in a energy gain of 0.67 eV (respectively, 0.44 eV) per formula unit. As a result, our system can be classified as a “proper ferroelectric” if one uses the P4=mmm phase as a reference. Let us now consider another reference structure, that is, the Pmc21 state, for which the polarization is zero but for which the AFE and AFD
week ending 7 FEBRUARY 2014
vectors are optimized to minimize the total energy. Such a new reference structure has four equivalent minima—that are ðþA; þωÞ, ð−A; þωÞ, ðþA; −ωÞ, and ð−A; −ωÞ—with an energy gain of 0.92 eV per formula unit with respect to the P4=mmm phase that is larger than the instability due to polarization only. Moreover, starting from the ðþA; þωÞ minimum and allowing the polarization to vary from negative to positive values (while freezing the A and ω values) constitutes what we call “case 2” in the following. It leads to an energy-versus-polarization single well having an energy gain of 1.05 eV per formula unit with respect to the P4=mmm phase, with the energy being nearly linear around P ¼ 0 [as is consistent with the trilinear coupling of Eq. (1)]. Fitting this energy-versus-polarization curve by a polynomial up to fourth order in polarization yields a negative linear coefficient [as is also consistent with the trilinear coupling of Eq. (1)], while the quadratic parameter is still negative but has a magnitude around 19 times smaller than that in case 1. Let us finally consider another Pmc21 structure as a reference that is the one for which the polarization still vanishes but the AFE and AFD vectors are those of the equilibrium Pmc21 state (for which P, A, and ω are all optimized to minimize the total energy). This second Pmc21 reference structure also has four equivalent minima (due to the possible positive and negative signs of A and ω) with an energy gain of 0.76 eV per formula unit with respect to the P4=mmm phase that is now slightly smaller than the instability due to polarization only. Starting from the ðþA; þωÞ minima and letting the polarization vary from negative to positive values (while freezing the A and ω values) constitutes what we we denote as “case 3.” It results in an energy-versus-polarization curve that is still nearly linear around P ¼ 0, but such a curve now has a double well that is asymmetric (the lowest well has an energy gain of 1.29 eV per formula unit with respect to the P4=mmm phase). Such a change in shape of the energy-versus-polarization curve with respect to case 2 arises from the fact that its fitting by a polynomial up to fourth order in polarization still yields a negative linear coefficient as in case 2, but the magnitude of the negative quadratic parameter is now about 6 times larger than that in case 2. This dramatic change of the magnitude of the quadratic term between case 1, case 2, and case 3 originates from the sum between the (negative) b0 coefficient and the (positive) eω2 þ fA2 terms introduced in Eq. (8) of the Supplemental Material [28], since these latter terms depend on the magnitude of the AFD and AFE vectors (note that case 3 has a much lower ω than case 2 due to the biquadratic repulsion between AFD motions and polarization). All the results indicated here may imply that the traditional concept of proper-versus-improper ferroelectricity can be rather subtle in our case and may have to be revisited or generalized [34]. They also suggest that, in systems like this one, it may be mandatory to study temperature-dependent properties in order to determine the character of the ferroelectricity. [34] M. Stengel, C. J. Fennie, and P. Ghosez, Phys. Rev. B 86, 094112 (2012). [35] Additional first-principles calculations indicate that the total energy is quasidegenerate (namely, within 0.7 meV=5 atoms) when placing G along several high-symmetry directions, with a minute preference for the G-type AFM vector to lie along [001]. [36] L. Bellaiche, Z. Gui, and I. A. Kornev, J. Phys. Condens. Matter 24, 312201 (2012).
057202-5
PHYSICAL REVIEW LETTERS
PRL 112, 057202 (2014)
Prediction of a Novel Magnetoelectric Switching Mechanism in Multiferroics 1
Yurong Yang,1,2 Jorge Íñiguez,3 Ai-Jie Mao,4 and L. Bellaiche1
Physics Department and Institute for Nanoscience and Engineering, University of Arkansas, Fayetteville, Arkansas 72701, USA 2 Physics Department, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China 3 Institut de Ciència de Materials de Barcelona (ICMAB-CSIC), Campus UAB, 08193 Bellaterra, Spain 4 Institute of Atomic and Molecular Physics, Sichuan University, Chengdu 610065, China (Received 30 May 2013; published 5 February 2014) We report a first-principles study of the recently predicted Pmc21 phase of the multiferroic BiFeO3 material, revealing a novel magnetoelectric effect that makes it possible to control magnetism with an electric field. The effect can be viewed as a two-step process: Switching the polarization first results in the change of the sense of the rotation of the oxygen octahedra, which in turn induces the switching of the secondary magnetic order parameter. The first step is governed by an original trilinear-coupling energy between polarization, octahedral tilting, and an antiferroelectric distortion. The second step is controlled by another trilinear coupling, this one involving the predominant and secondary magnetic orders as well as the oxygen octahedral tilting. In contrast with other trilinear-coupling effects in the literature, the present ones occur in a simple ABO3 perovskite and involve a large polarization. DOI: 10.1103/PhysRevLett.112.057202
PACS numbers: 75.85.+t, 77.55.Nv, 77.80.Fm
Multiferroics have experienced a resurgence of interest because of their inherent cross coupling between electric and magnetic degrees of freedom [1–7]. In particular, a systematic control of the magnitude and crystallographic direction of magnetic order parameters by an electric field is attractive for the design of original devices, and fundamentally interesting. One way to attain such control is via the existence of an energy term ΔE that (1) directly couples magnetic moments and electric polarization P and (2) is linear in P. Familiar examples of such an energy term are specific Dzyaloshinskii-Moriya interactions [8–10] of the form ΔE ∼ P · ðL × MÞ [10,11], the Lifshitz invariants ΔE ∼ P · ½Lð∇ · LÞ þ L × ð∇ × LÞ and ΔE ∼ P · ½Mð∇ · MÞ − ðM · ∇ÞMÞ [12,13], where L and M are the antiferromagnetic vector and magnetization, respectively. Another example is the spin-current model [14–16] for which ΔE ∼ ðP × eij Þ · ðmi × mj Þ, where eij is the unit vector joining site i to site j and where mi and mj are the magnetic moments located at these sites i and j. The pioneering work of Ref. [17] also shows that there is another way to control magnetic quantities by electric polarization. This can be thought of as indirect since it consists of two steps. A trilinear “structural-only” coupling between polarization and two different oxygen octahedral tilting modes (also termed antiferrodistortive or AFD motion) first leads to the reversal of one octahedral tilt when switching the polarization. Then, this reversal can result in the change of direction of a magnetic order parameter (via a second energy that was not provided in Ref. [17] and that should couple magnetic moments and octahedral tilt). It is important to know if other two-step processes, allowing a systematic control of magnetic properties by electrical polarization, can exist. In particular, can the first energy also involve another kind of structural quantity, in addition to 0031-9007=14=112(5)=057202(5)
polarization and AFD motions, namely, antiferroelectricity? Another question to resolve is to determine the expression of the second energy that couples magnetic degrees and the polarization-switchable structural property. We provide answers to these issues via first-principles calculations. A novel two-step magnetoelectric effect is discovered. Its first energy is an original term that trilinearly couples polarization, one AFD distortion, and one antiferroelectric mode. This structural-only energy results in the change of the sense of the rotation of the oxygen octahedra when polarization is reversed. The second energy is another coupling involving predominant and secondary magnetic order parameters, as well as the AFD mode. It forces the secondary magnetic order parameter to switch its direction, once oxygen octahedral tilting is reversed by the change of sign of the polarization. The designation of the phenomenon as having two steps is made to indicate that two different energies are involved, even if the switchings of the two physical properties resulting from the existence of these two energies can be simultaneous. Density-functional calculations within the local spin density approximation are performed for BiFeO3 (BFO) systems, using the Vienna ab initio simulation package (VASP) [18]. As in Refs. [19–21], we also include a Hubbard correction for Fe ions, adopting the self-consistent value of U ¼ 3.8 eV [22,23] (note that this choice of U provided a good agreement between theory and experiment even for films under a large magnitude of strain; see, e.g., Ref. [24] and references therein). Spin-orbit and noncollinear magnetism are taken into account, and the projected augmented wave method and an energy cutoff of 500 eV are employed. We focus on the orthorhombic Pmc21 phase of (001) BFO films that has recently been found for large enough tensile strain [20] (note that this state has also been predicted in BaMnO3 , EuTiO3 , or CaTiO3 , and that increasing the temperature should lead to a
057202-1
© 2014 American Physical Society
PRL 112, 057202 (2014)
PHYSICAL REVIEW LETTERS
decrease of the critical strain at which this phase becomes the ground state). A 40-atom supercell is used to model it, along with a 3×3×3 k-point mesh. Its lattice vectors are a1 ¼2aIP x, a2 ¼ 2aIP y, and a3 ¼ aIP ½δ1 x þ δ2 y þ ð2 þ δ3 Þz. aIP is the in-plane lattice constant. It is chosen here to be 4.135 Å, which leads to a tensile strain of about 6% in BFO (note that BFO films experiencing strains having a magnitude even larger than 6% have been experimentally grown [25]). x, y, and z are unit vectors along the [100], [010], and [001] pseudocubic directions, respectively. The δ1 , δ2 , and δ3 variables and the atomic positions are relaxed to minimize the total energy. As sketched in Fig. 1, this Pmc21 state is characterized by (1) a polarization (P) lying along an in-plane h110i direction, (2) a Bi-driven antiferroelectric (AFE) vector (A) that is associated with the M point of the cubic, five-atom first ¯ Brillouin zone and that is oriented along an in-plane h110i direction that is perpendicular to the polarization, and (3) by an in-phase oxygen octahedral tilting (that is also associated with the M point of the first Brillouin zone) about the out-of-plane [001] direction. Such a tilting is represented by a pseudovector ω for which the magnitude provides the value of the tilting angle while its direction is either parallel or antiparallel to the out-of-plane direction—indicating if this tilting is counterclockwise versus clockwise, respectively. The superposition of P and A results in cations moving in a zigzag pattern [20], as shown in the top panels of Fig. 1. The polarization is calculated from the Bloch representation of the modern theory of polarization [26]. The strength of the antiferroelectricity vector is estimated by A¼ðdBi;A =dBi;P ÞP, where dBi;A repre¯ sents the atomic motion of Bi ions (along h110i) associated with antiferroelectricity, while dBi;P is the Bi displacement (along h110i) associated with the polarization P. Let us denote by þP and −P the cases in which P is parallel or antiparallel to the [110] direction, respectively. Similarly, once we choose a specific Bi site to define the origin of the AFE vector [27], þA and −A represent the AFE ¯ and [110], ¯ vectors along [110] respectively. Finally, once we select an Fe site to define the origin of the AFD motion [27], þω and −ω indicate that the oxygen octahedral tiltings (about a specific line joining Fe atoms along [001]) are counterclockwise and clockwise, respectively. When considering the eight cases resulting from the combination of the positive and negative signs of P, A, and ω, we found that the ground state of the Pmc21 state is fourfold degenerate; namely, the (þP, þA, þω), (þP, −A, −ω), (−P, −A, þω), and (−P, þA, −ω) all have the same energy, as well as identical magnitude for their P, A, and ω vectors. The other four states (−P, þA, þω), (−P, −A, −ω), (þP, −A, þω), and (þP, þA, −ω) are unstable and transform to one of the four stable equilibrium states during relaxation. These results suggest the existence of a very peculiar physical energy that couples P, A, and ω in a very precise manner. The Supplemental Material [28] provides the expression of this energy term: ΔE1 ¼ −C1 PAω; where C1 is a positive coupling parameter.
(1)
week ending 7 FEBRUARY 2014
FIG. 1 (color online). Schematization of properties in the equilibrium (þP, þA, þω) state (left panels) and (−P, þA, −ω) phase (right panels). The upper panels display the zigzag motions of Bi atoms in an (x, y) plane. The middle panels show (i) oxygen octahedral tiltings (represented by curled arrows), (ii) magnetic dipoles associated with the weak A-type antiferromagnetic vector (purple arrows) and with the predominant G-type AFM vector (represented by green arrows), and (iii) Bi displacements that are associated with AFE (light blue arrows). The lettering “Fe” (respectively, “Bi”) indicates the iron (respectively, bismuth) ion whose location serves as the origin for the definition of the ω, G, and AAFM (respectively, A) vectors [27]. The lower panels schematize the directions of P, the AFE vector, the AFD vector, the G-type AFM vector, and the weak A-type AFM vector.
As characteristics of ferroelectrics, the application of an electric field that is opposed to an initial polarization direction should switch the polarization. According to Eq. (1), this switch should make the oxygen octahedra rotate from counterclockwise to clockwise, if the AFE vector maintains its direction. Equation (1) further suggests that reverting this electric field should then allow the oxygen octahedra to tilt again in a counterclockwise manner and thus hints towards a systematic control of the oxygen octahedral tilting direction by an electric field. To test such a possibility, Fig. 2(a) reports the total energy resulting from the variation of the polarization from its positive equilibrium value of þ1.12 C=m2 to its opposite, negative equilibrium value (along ½1¯ 1¯ 0), and vice versa. Figures 2(b) and 2(c) show the AFE vector and out-of-plane, in-phase octahedral tilting resulting from the atomic relaxation occurring for each considered value of the polarization selected during the pathways (þP → −P) and (−P → þP). They reveal that one can indeed go from the equilibrium (þP, þA, þω) state to the equilibrium (−P, þA,
057202-2
PRL 112, 057202 (2014)
PHYSICAL REVIEW LETTERS
FIG. 2 (color online). Behavior of various quantities as a function of P, during the ðþP; þA; þωÞ → ð−P; þA; −ωÞ and ð−P; þA; −ωÞ → ðþP; þA; þωÞ pathways. (a),(b),(c) The total energy, AFE vector, and in-phase octahedral tilting, respectively. The insets of (b) schematize the evolution of the resulting Bi displacements along these paths, with the red, blue, and black arrows representing the polarization-induced, AFE-driven, and total Bi displacements, respectively. In (c), symbols represent data from first-principles calculations while the blue (respectively, red) solid line shows the fitting of these data by fourth order polynomials of P for positive (respectively, negative) ω.
−ω) phase (and vice versa) when continuously changing the magnitude and sign of the polarization. Such a transformation occurs via a pathway having an energy barrier of 0.36 eV per formula unit [29], along which the AFE vector does not vary too much. This implies, as shown in the insets of Fig. 2(b), that the overall cation displacement gradually evolves from a zigzag pattern (when A and P are close to each other in magnitude, that is, for a large magnitude of P) to a pure AFE pattern (which occurs when P vanishes). The overall symmetry of the structure always remains Pmc21 during the whole pathway, i.e., even when the polarization vanishes. The considered pathway is only one possibility out of the many that could occur in an actual material, but is very sound. Indeed, due to the persistence of the AFE distortion throughout, this path involves a quasicontinuous rotation of the electric dipoles associated with individual Bi cations, which are always in an off-centered, low-energy position. We also considered the transformation from (þP, þA, þω) to the equivalent state (−P, −A, þω) and found that this possibility, which involves the switching of the AFE vector, has a higherenergy barrier of 0.84 eV=f.u. associated with it. Figure 2 also shows that two solutions (P, þA, þω) and (P, þA, −ω) exist when P ranges from −0.60 to 0.60 C=m2 ,
week ending 7 FEBRUARY 2014
with (P, þA, þω) [respectively, (P, þA, −ω)] having a lower energy for positive (respectively, negative) polarization—as is consistent with Eq. (1). On the other hand, only one of these two solutions is found for larger magnitudes of P. Such a result can be understood by Eq. (1) that indicates that the difference in energy between the (P, þA, þω) and (P, þA, −ω) states is −2C1 PAω and thus increases with the magnitude of P. For large values of P, the state with higher energy thus becomes unstable, and only the state with lower energy [that is, (þP, þA, þω) for large positive P or (−P, þA, −ω) for negative P having a large magnitude] can persist. Figure 2(c) further shows that the behavior of the z component of the AFD vector ω as a function of P is quite complex, both when ω is positive or negative. As also shown in Fig. 2(c), such behaviors can be fitted by a polynomial offourth order in P (that is derived in the Supplemental Material [28]). Trilinear couplings involving P have also been proposed for the layered Ca3 Mn2 O7 [17], PbTiO3 =SrTiO3 superlattices [31], and strained CaTiO3 [32]. However, the other two quantities appearing in the trilinear couplings of these systems are either two different AFD motions [17,31] or two different antipolarmodes[32].TheproposedEq.(1)isthereforeoriginal since three different physical quantities (polarization, AFE, and AFD motions) are intertwined. Moreover, the polarizations involved in the trilinear couplings of Refs. [17,31,32] are small (about 10−6 to 0.2 C=m2 ), while the polarization of our equilibrium Pmc21 phase is much larger (namely, 1.12 C=m2 ). Moreover, our footnote or Ref. [33] provides a discussion about the nature of this latter polarization (that is, proper versus improper). Let us now investigate if the polarization-induced switching between (þP, þA, þω) and (−P, þA, −ω) can have some consequence on magnetism. For that, let us first choose a G-type antiferromagnetic (AFM) vector G being aligned along [110] to demonstrate magnetoelectric switching [35]. Our calculations predict that, in that case, a weak A-type AFM ¯ vector (to be denoted by AAFM ) is created along ½110 if the system is in the (þP, þA, þω) equilibrium state, while it is ¯ direction if the system is in the (−P, along the opposite ½110 þA, −ω) equilibrium state (see Fig. 1) [27]. One should thus be able to reverse the direction of the weak AFM vector AAFM when applying an external electric field [since such a field should induce the switching between (þP, þA, þω) and (−P, þA, −ω) states]. To shed more light on such an effect, Fig. 3 shows the evolution of the projection of AAFM onto ¯ for the paths considered in Fig. 2, that is, as a function of ½110 P varying between þ1.12 and −1.12 C=m2 . This projection is denoted as AAFM . Continuously decreasing P from its largest studied positive value to its largest-in-magnitude negative value leads to the following behaviors: (i) the negative AAFM first decreases in strength until P is close to þ0.75 C=m2 ; (ii) the negative AAFM then increases in magnitude until exhibiting a maximum for small polarizations; and (iii) the AAFM remains negative and decreases again in magnitude to suddenly jump for a negative P that is in between −0.6 and −0.8 C=m2 and then adopts a positive
057202-3
PRL 112, 057202 (2014)
PHYSICAL REVIEW LETTERS
week ending 7 FEBRUARY 2014
confirm such an exciting effect. Moreover, the discovered two-step control of the secondary magnetic order parameter is not restricted to BFO, since Eqs. (1) and (2) emerge from symmetry considerations that are not material specific.
FIG. 3 (color online). Evolution of the A-type AFM vector AAFM as a function of P. The square symbols show the values predicted from first principles. The solid lines represent the fitting of these data by AAFM ¼ −0.00051Gω þ 0.059GAPþ 0.00011GωP2 − 0.064GAP3 , where G is the G-type AFM vector (that has a magnitude of around 4μB ), A is the polarization-dependent AFE vector displayed in Fig. 2(b), and ω is the in-phase oxygen octahedral tilting [that also depends on P, as shown in Fig. 2(c)].
value for P less than −0.8 C=m2 . The change in sign of AAFM occurs precisely when ω also reverts its sign [see Fig. 2(c)] along the equilibrium (þP, þA, þω) → equilibrium (−P, þA, −ω) pathway. Figure 3 also shows that AAFM behaves in a mirroring fashion when studying the reverse ð−P; þ A; −ωÞ → ðþP; þ A; þ ωÞ pathway. Throughout these two pathways, we numerically found that the G-type antiferromagnetic vector not only remains large in magnitude (close to 4 bohr magneton) but also stays oriented along [110]. To understand these effects, one has to recall that Ref. [36] suggested an energy ΔE2 ¼ −C2 ω · ðG × AAFM Þ;
(2)
where C2 is a positive coefficient. In other words, ω, G, and AAFM want to form a direct triad, as schematized in Fig. 1. Such an equation explains why the combined existence of a G-type AFM vector and in-phase oxygen octahedral tilting results in a weak A-type AFM vector. It also explains why the polarization-induced switching of ω [such a switch is driven by Eq. (1)] leads to a reversal of AAFM (see Fig. 1). The control of the AAFM vector by an electric field is thus made possible via a two-step process [described by Eqs. (1) and (2)], both of which involve AFD motions. Such a two-step process represents a novel magnetoelectric effect. As is also detailed in the Supplemental Material [28], the nontrivial behavior of AAFM as a function of polarization (Fig. 3) can be reproduced from the energy terms discussed above. The proposed control of magnetism by polarization is not restricted to the switching of the A-type AFM vector. For instance, let us imagine a Pmc21 state for which the predominant magnetic ordering is of C-type. Reference [36] tells us that Eq. (2) would have to be replaced by ΔE2 ¼ −C2 ω · ðCAFM × FÞ, where CAFM is the C-type AFM vector and F is the (weak) ferromagnetic vector. Combining this latter equation with Eq. (1) implies that a net magnetic moment should also be controllable by polarization, in that case. Further simulations we conducted indeed
This work is supported by DOE grant ER-46612. We also acknowledge NNSF of China (Grants No. 51032002, No. 11374162, and No. 11104190) and 863 Program of China (Grant No. 2011AA050526), and MINECO-Spain (Grants No. MAT2010-18113 and No. CSD2007-00041). Some computations were also made possible thanks to the MRI grant 0722625, MRI-R2 grant 0959124, and CITRAIN grant 0918970 from NSF, the ONR grant N0001407-1-0825 (DURIP) and a Challenge grant from the Department of Defense.
[1] M. Fiebig, J. Phys. D 38, R123 (2005). [2] Z. J. Huang, Y. Cao, Y. Y. Sun, Y. Y. Xue, and C. W. Chu, Phys. Rev. B 56, 2623 (1997). [3] T. Kimura, T. Goto, H. Shintani, K. Ishizaka, T. Arima, and Y. Tokura, Nature (London) 426, 55 (2003). [4] T. Lottermoser, T. Lonkai, U. Amann, D. Hohlwein, J. Ihringer, and M. Fiebig, Nature (London) 430, 541 (2004). [5] C. Ederer and N. A. Spaldin, Phys. Rev. B 74, 020401 (2006). [6] C. Ederer and C. J. Fennie, J. Phys. Condens. Matter 20, 434219 (2008). [7] C. J. Fennie, Phys. Rev. Lett. 100, 167203 (2008). [8] I. Dzyaloshinsky, J. Phys. Chem. Solids 4, 241 (1958). [9] T. Moriya, Phys. Rev. Lett. 4, 228 (1960). [10] J. F. Scott and R. Blinc, J. Phys. Condens. Matter 23, 113202 (2011). [11] S. Prosandeev, I. A. Kornev, and L. Bellaiche, Phys. Rev. B 83, 020102(R) (2011). [12] I. Sosnowska and A. Zvezdin, J. Magn. Magn. Mater. 140– 144, 167 (1995). [13] M. Mostovoy, Phys. Rev. Lett. 96, 067601 (2006). [14] H. Katsura, N. Nagaosa, and A. V. Balatsky, Phys. Rev. Lett. 95, 057205 (2005). [15] D. Rahmedov, D. Wang, J. Íñiguez, and L. Bellaiche, Phys. Rev. Lett. 109, 037207 (2012). [16] A. Raeliarijaona, S. Singh, H. Fu, and L. Bellaiche, Phys. Rev. Lett. 110, 137205 (2013). [17] N. A. Benedek and C. J. Fennie, Phys. Rev. Lett. 106, 107204 (2011). [18] G. Kresse and D. Joubert, Phys. Rev. B 59, 1758 (1999). [19] S. Prosandeev, I. A. Kornev, and L. Bellaiche, Phys. Rev. Lett. 107, 117602 (2011). [20] Y. Yang, W. Ren, M. Stengel, X. H. Yan, and L. Bellaiche, Phys. Rev. Lett. 109, 057602 (2012). [21] Y. Yang, M. Stengel, W. Ren, X. H. Yan, and L. Bellaiche, Phys. Rev. B 86, 144114 (2012). [22] V. I. Anisimov, F. Aryasetiawan, and A. I. Lichtenstein, J. Phys. Condens. Matter 9, 767 (1997). [23] I. A. Kornev, S. Lisenkov, R. Haumont, B. Dkhil, and L. Bellaiche, Phys. Rev. Lett. 99, 227602 (2007). [24] Z. Chen, S. Prosandeev, Z. L. Luo, W. Ren, Y. Qi, C. W. Huang, L. You, C. Gao, I. A. Kornew, T. Wu, et al., Phys. Rev. B 84, 094116 (2011).
057202-4
PRL 112, 057202 (2014)
PHYSICAL REVIEW LETTERS
[25] R. J. Zeches, M. D. Rossell, J. X. Zhang, A. J. Hatt, Q. He, C.-H. Yang, A. Kumar, C. H. Wang, A. Melville, C. Adamo, et al., Science 326, 977 (2009). [26] R. D. King-Smith and D. Vanderbilt, Phys. Rev. B 47, 1651 (1993). P [27] The AFE vector A is defined as A¼1=N j dBi;A;j ð−1Þlj þmj , where the sum runs over all the N Bi sites j, and where dBi;A;j is the displacement of the Bi ion at site j that is induced by antiferroelectricity. The site j is located at alat ðlj x þ mj y þ nj zÞ from a Bi site chosen as the origin, where alat is the fiveatom cubic lattice constant and where lj , mj , and nj are all integers. Similarly, the in-phase AFD vector ω is defined as P ω ¼ 1=N k ωk ð−1Þlk þmk , where the sum runs over all the N Fe sites k, and where ωk is the vector characterizing the tilting centered on the Fe site k. The latter site is located at alat ðlk x þ mk y þ nk zÞ from an Fe site selected as the origin, with lk , mk , and nk being integers. The specific Bi and Fe atoms that have been chosen as the origins to define A and ω, respectively, are indicated in the middle panels of Fig. 1. Note that thePAFM vectors AAFM andP G are defined by AAFM ¼1=N k μk ð−1Þnk and G ¼ 1=N k μk ð−1Þlk þmk þnk , where μk is the magnetic moment at Fe site k. [28] See Supplemental Material at http://link.aps.org/ supplemental/10.1103/PhysRevLett.112.057202 for additional details and derivations related to the energy of Eq. (1), as well as, to the dependence of the antiferrodistortive vector and (weak) A-type antiferromagnetic vector on polarization. [29] Note that dividing this barrier of 0.36 eV per formula unit by the spontaneous electric dipole moment (that is, the product between the polarization and the unit volume) of the investigated system in its equilibrium phase yields an electric field of 0.8 × 109 V=m. Such a field has the same order of magnitude as the one recently experimentally applied in BFO films (that is, 109 V=m) [30]. This implies that the barrier of 0.36 eV per formula unit can be overcome by a realistic electric field. Moreover, this barrier should decrease when increasing the temperature, rendering the electric field needed to overcome the barrier even smaller and easier to experimentally achieve at finite temperature. [30] H. Yamada et al., ACS Nano 7, 5385 (2013). [31] E. Bousquet, M. Dawber, N. Stucki, C. Lichtensteiger, P. Hermet, S. Gariglio, J.-M. Triscone, and P. Ghosez, Nature (London) 452, 732 (2008). [32] Q. Zhou and K. M. Rabe, arXiv:1306.1839v1. [33] Computing the total energy as a function of polarization, once starting from the P4=mmm phase (which is obtained from the ideal cubic structure by imposing the epitaxial constraint) and setting the AFE and AFD vectors to both vanish, constitutes what we will call “case 1.” It leads to a symmetric double well with an energy gain of 0.78 eV per formula unit, with the energy-versus-polarization curve being nicely fitted by a polynomial solely involving the second and fourth powers of the polarization. This is the strongest instability of the reference structure, since letting only the AFE (respectively, AFD) motions occur results in a energy gain of 0.67 eV (respectively, 0.44 eV) per formula unit. As a result, our system can be classified as a “proper ferroelectric” if one uses the P4=mmm phase as a reference. Let us now consider another reference structure, that is, the Pmc21 state, for which the polarization is zero but for which the AFE and AFD
week ending 7 FEBRUARY 2014
vectors are optimized to minimize the total energy. Such a new reference structure has four equivalent minima—that are ðþA; þωÞ, ð−A; þωÞ, ðþA; −ωÞ, and ð−A; −ωÞ—with an energy gain of 0.92 eV per formula unit with respect to the P4=mmm phase that is larger than the instability due to polarization only. Moreover, starting from the ðþA; þωÞ minimum and allowing the polarization to vary from negative to positive values (while freezing the A and ω values) constitutes what we call “case 2” in the following. It leads to an energy-versus-polarization single well having an energy gain of 1.05 eV per formula unit with respect to the P4=mmm phase, with the energy being nearly linear around P ¼ 0 [as is consistent with the trilinear coupling of Eq. (1)]. Fitting this energy-versus-polarization curve by a polynomial up to fourth order in polarization yields a negative linear coefficient [as is also consistent with the trilinear coupling of Eq. (1)], while the quadratic parameter is still negative but has a magnitude around 19 times smaller than that in case 1. Let us finally consider another Pmc21 structure as a reference that is the one for which the polarization still vanishes but the AFE and AFD vectors are those of the equilibrium Pmc21 state (for which P, A, and ω are all optimized to minimize the total energy). This second Pmc21 reference structure also has four equivalent minima (due to the possible positive and negative signs of A and ω) with an energy gain of 0.76 eV per formula unit with respect to the P4=mmm phase that is now slightly smaller than the instability due to polarization only. Starting from the ðþA; þωÞ minima and letting the polarization vary from negative to positive values (while freezing the A and ω values) constitutes what we we denote as “case 3.” It results in an energy-versus-polarization curve that is still nearly linear around P ¼ 0, but such a curve now has a double well that is asymmetric (the lowest well has an energy gain of 1.29 eV per formula unit with respect to the P4=mmm phase). Such a change in shape of the energy-versus-polarization curve with respect to case 2 arises from the fact that its fitting by a polynomial up to fourth order in polarization still yields a negative linear coefficient as in case 2, but the magnitude of the negative quadratic parameter is now about 6 times larger than that in case 2. This dramatic change of the magnitude of the quadratic term between case 1, case 2, and case 3 originates from the sum between the (negative) b0 coefficient and the (positive) eω2 þ fA2 terms introduced in Eq. (8) of the Supplemental Material [28], since these latter terms depend on the magnitude of the AFD and AFE vectors (note that case 3 has a much lower ω than case 2 due to the biquadratic repulsion between AFD motions and polarization). All the results indicated here may imply that the traditional concept of proper-versus-improper ferroelectricity can be rather subtle in our case and may have to be revisited or generalized [34]. They also suggest that, in systems like this one, it may be mandatory to study temperature-dependent properties in order to determine the character of the ferroelectricity. [34] M. Stengel, C. J. Fennie, and P. Ghosez, Phys. Rev. B 86, 094112 (2012). [35] Additional first-principles calculations indicate that the total energy is quasidegenerate (namely, within 0.7 meV=5 atoms) when placing G along several high-symmetry directions, with a minute preference for the G-type AFM vector to lie along [001]. [36] L. Bellaiche, Z. Gui, and I. A. Kornev, J. Phys. Condens. Matter 24, 312201 (2012).
057202-5
