The multiferroic phase of DyFeO$_{3}$:an ab--initio study

By performing accurate ab-initio density functional theory calculations, we study the role of $4f$ electrons in stabilizing the magnetic-field-induced ferroelectric state of DyFeO$_{3}$. We confirm that the ferroelectric polarization is driven by an exchange-strictive mechanism, working between adjacent spin-polarized Fe and Dy layers, as suggested by Y. Tokunaga [Phys. Rev. Lett, \textbf{101}, 097205 (2008)]. A careful electronic structure analysis suggests that coupling between Dy and Fe spin sublattices is mediated by Dy-$d$ and O-$2p$ hybridization. Our results are robust with respect to the different computational schemes used for $d$ and $f$ localized states, such as the DFT+$U$ method, the Heyd-Scuseria-Ernzerhof (HSE) hybrid functional and the GW approach. Our findings indicate that the interaction between the $f$ and $d$ sublattice might be used to tailor ferroelectric and magnetic properties of multiferroic compounds.


I. INTRODUCTION
Multiferroics (MFs) are compounds where long-range magnetic and dipolar orders coexist [1]. There is plenty of fascinating physics in these materials, due to the strong entanglement of spin-charge-orbital degrees of freedom [2,3] and a great potential for technological applications has already been recognized [4][5][6]. The coupling between magnetism and ferroelectricity can be used for sensing applications, but also for memory devices where data is typically stored as magnetic information and read out electrically. Recently, several manganese and iron oxides have been shown to possess strong coupling; however, ferroelectricity in these materials is rather weak and only the electrical polarization can be switched by a magnetic field (but not viceversa, a limitation for many applications). Within this framework, DyFeO 3 is a very interesting compound [7,8] because it shows large ferroelectric polarization combined with a strong magnetoelectric coupling. Furthermore, the ferroelectric polarization in DyFeO 3 is induced by the peculiar magnetic structure, obtained upon applying an external magnetic field.
DyFeO 3 belongs to the the class of perovkite oxides, such as RMO 3 (where R is a rare-earth ion and M is a transition metal atom). From the computational point of view, their description poses serious problems. These compounds are often (strongly) correlated materials, involving d− and f − electronic charge with significant spatial localization. First-principles density functional theory (DFT) calculations in the most commonly applied local density approximation (LDA) or generalized gradient approximation (GGA) have to face well known deficiencies: the nonlocality of the screened exchange interaction is not well taken into account and the electrostatic self-interaction is not entirely compensated [13]. Since semilocal functionals tend to delocalize f -states, the f -electrons are often kept "frozen" in the core, and the origin of multiferroicity is generally attributed to the spin-charge-orbital degrees of freedom of the M sublattice. Although this "standard" approach helped in clarifying many mechanisms leading to ferroelectricity in MFs [14], the influence of f electrons on multiferroicity has not been extensively investigated yet by ab-initio calculations, despite several experiments clearly point out f -electrons to play an important role in MF properties of RMO 3 compounds [15,16]. This is especially so for DyFeO 3 , where the electric polarization appears at a very low temperature, corresponding to ordering of Dy spins [9][10][11]. Recall that Dy is expected to have a 3+ oxidation state, i.e. a f 9 configuration [12].
In this study, stimulated by the recent experimental demonstration of a magnetic-field-induced ferroelectric (FE) phase of DyFeO 3 (DFO) [9][10][11], we performed ab-initio simulations in order to show the key role played by f electrons in stabilizing spin-driven ferroelectricity. Our calculations support the experimental study of Ref. 10,11 and, at the same time, quantify the polarization and shed light on the microscopic mechanism (i.e. based on electronic structure analysis) of the origin of ferroelectricity. To the best of our knowledge, this is the first ab-initio study of the above-mentioned effect. Our simulations were mainly carried out within a DFT+U approach for localized electrons. In addition, we used the Heyd-Scuseria-Ernzerhof (HSE) screened hybrid functional [17,18], which has been shown to improve the description of d-and f -electron systems [19][20][21] over LDA or GGA. Finally, our calculations were benchmarked by single shot GW using HSE wave functions (G 0 W 0 @HSE); this treatment is expected to give a very accurate description of the electronic ground state [22].

II. COMPUTATIONAL DETAILS
Calculations were performed using the projector augmented-wave (PAW) method [23,24] with the Perdew-Burke-Ernzerhof (PBE) GGA functional [25]. We used DFT+U within Dudarev's approach [26] using U eff = U − J=3 eV and 4 eV for Fe-d and Dy-f states, respectively. The energy cutoff was set to 400 eV and a 4×2×4 Monkhorst-Pack grid of k-points was used. We treated the Dy f electrons both as valence and as core states [27]. The Berry phase approach [28,29] was used to calculate the macroscopic polarization P . Non collinear magnetism was treated in accordance with Ref. [30]. Spin-orbit coupling (SOC) was included for the end-point states of the adiabatic path, i.e. paraelectric (PE) and ferroelectric (FE) DyFeO 3 , see below. For hybrid functionals, we used the HSE functional, recently implemented in VASP [31]. GW calculations [32,33] were performed on top of the HSE ground state [22]. The experimental lattice constants for orthorhombic DFO were used (space group P nma, with a=5.596Å, b=7.629Å, c=5.301Å) [7]. Starting from experimental atomic positions, we performed atomic relaxations until residual Hellman-Feynman forces were <0.01 eV/Å.
In figure 1, we show the two FE states of DFO, FE 1 and FE 2 , with opposite polarization. These represent a simplified version of the complex experimental magnetic-field-induced polar states (cfr figure 1 (c) and (d) of Ref. 10), obtained by neglecting the x and z components of Dy and Fe spins. [41]. In both, the FE 1 and the FE 2 cases, the Fe and Dy sublattices show an A-type magnetic structure, i.e. ferromagnetic (FM) and antiferromagnetic (AFM) intralayer and interlayer coupling, respectively, with the spins pointing along the in-plane c axis. In FE 1 , the stacking of the spins along the out-of-plane b axis is ↓Fe-↓Dy-↑Fe-↑Dy; the FE 2 state is obtained from the FE 1 state by rotating Dy spins by 180 • , so that the spin stacking is ↓Fe-↑Dy-↑Fe-↓Dy. In the PE state, the atoms are arranged according to the P nma (D 2h ) space group and the Dy spins are still intra-layer FM coupled, but rotated with respect to the Fe spins by 90 • . In FE 1 , the spins of a Fe layer become parallel to the moments on one of the nearest-neighbor Dy layers and antiparallel to the other: Dy layers should then displace cooperatively towards Fe layers with parallel spins via exchange striction, giving rise to alternating short-long-short-long interlayer Dy-Fe distances. Accordingly, a polarization P along the b axis should be generated. In the FE 2 state, the flip of Dy spins would cause a reversal of P . Finally, in the PE state, each Fe sheet is sandwiched by layers with ⊥ Dy spins, and no interlayer dimerization is expected. This microscopic mechanism, which involves frustrated interactions between rare earth and transition metal ions and its optimization by exchange striction, was proposed in Refs. 10, 11.

III. SWITCHABLE FERROELECTRIC STATES
Our calculations confirm the above interpretation. First of all, the total energies of FE 1 and FE 2 are degenerate.
[42] Furthermore, from the symmetry point-of-view, the rotation of Dy spins from PE to FE 1 or to FE 2 causes a symmetry lowering from the non-polar space group 62 (D 2h ) to a polar space group 33 (C 2v ), giving rise to a polarization along the b axis. Indeed, the distortions lead to an alternate short-long-short-long interlayer distance, with d  electrons play a key role in stabilizing the FE state. In order to prove this, we froze the f electrons in the core and we calculated the electronic ground state using the previously relaxed FE 1 structure. First, the electronic contribution to the polarization, P ele vanishes; then, by allowing the ions to relax, the ionic contribution, P ionic , becomes negligible as well, and the final crystal structure is non-polar. In summary, when treating the f electrons as valence states, the PE state becomes unstable, the D 2h point group symmetry is spontaneously broken to C 2v and the system evolves towards a stable and polar state. If the f electrons are removed from the valence and frozen in the core, the PE state remains stable. This unambiguously confirms that f states are a necessary ingredient for ferroelectricity in DyFeO 3 .

A. Adiabatic path
To shed light on the onset of ferroelectricity, we construct an adiabatic path by progressively rotating the Dy spins from 90 to 0 degree, i.e. from the PE to FE 1 state. The results are summarized in figure 2: panel (I) shows the energy difference between the FE 1 and the PE phases, evaluated at the ideal centrosymmetric (CS) ionic structure (in blue) and relaxed configuration (in red); in (II), (III), (IV) we show the electronic, ionic and total FE polarization, respectively, evaluated at the ideal (blue) and relaxed (red) ionic structure; in (IV) the Point-Charge-Model (PCM) estimate of the polarization is also reported (in black); in (V) the change in the interlayer distances along the path (see also Fig. 1) is shown. Panel (I) clearly shows that ferroelectricity in DFO is magnetically induced. In fact, already for the CS structure, the spin rotations give rise to an energy gain (blue bars), further enhanced by ionic relaxations (red bars). Note that the relaxation of the electronic degrees of freedoms accounts for most of the total stabilization energy (cfr blue and red bars). The energy gain increases from left to right in panel (I), i.e. towards a collinear configuration. In parallel, P ele also increases from left to right, even in the ideal CS structure (blue), as expected for magnetically induced MFs [34]; ionic relaxations further increase P ele [see panel (II)]. From panel (III), we see that P ion is opposite to P ele and of the same order of magnitude. However, P ion does not fully compensate P ele [see panel (IV)], giving rise to a total polarization P tot of ∼ 0.20 µCcm −2 for collinear spins. The inclusion of SOC confirms P tot , i.e. 0.18 µCcm −2 . Notably, this value is in good agreement with the estimate given in Ref. [10]. HSE also confirms the value of FE polarization. Furthermore, note that P pcm is not only opposite to P tot , but also smaller in absolute value than P tot . This confirms that electronic degrees of freedom trigger the FE transition. Finally, in panel (V), the evolution of the interlayer distances along the polar b axis clearly shows that the Dy and Fe layers are coupled, giving rise to spin-driven interlayer dimerization.

B. Ionic relaxations and origin of the ferroelectricity
In the PE phase Fe, Dy, O ap (apical oxygens), O eq (equatorial oxygens) occupy the 4a, 4c, 4c, and 8d Wyckoff positions (WPs), respectively. When the symmetry is lowered to C 2v , Dy and O ap change their site symmetry to 4a and the O eq become inequivalent (8d → 4a + 4a). This is readily explained by considering the local spin configuration around O eq s: when the latter are sandwiched by FM coupled Fe and Dy layers, they have two ↑Fe and two ↑Dy atoms as nearest neighbors (these will be labelled O ↑,↑ eq ); when sandwiched by Fe and Dy layers AFM coupled, they have two ↑Fe and two ↓Dy atoms as nearest neighbors (labelled as O ↑,↓ eq ). First of all, all O eq s carry a spin-induced moment parallel to the neighboring Fe atom. Furthermore, by imposing the FE 1 spin configuration on top of the CS ionic structure, O ↑,↑ eq and O ↑,↓ eq become inequivalent : the O ↑,↑ eq has ±0.194 µ B and the O ↑,↓ eq has ±0.207 µ B as induced spin moment. To rule out any numerical artifact on this small difference, we impose the PE spin configuration on top of the CS ionic structure. In this case, all O eq s carry induced spin moments of exactly the same magnitude. In passing, we note that all oxygens remain equivalent when freezing the Dy-f electrons in the core. Furthermore, we performed an analysis of the symmetry breaking distortions [35,36]. The mode decomposition confirms that a FE mode is involved, called GM 4−. In Fig. 3 we show the pattern of atomic displacements with respect to the CS structure. The inequivalency of O eq is subtle: O ↑,↑ eq (O ↑,↓ eq ) move in such a way to decrease (increase) the distance to its neighbor Dy atom. For O ↑,↑ eq , d Dy−O is 2.478Å; for O ↑,↓ eq , d Dy−O is 2.496Å (the same distance in the PE phase is 2.487Å), suggesting that a weak bonding interaction is active between the FM layers, leading to a polarization along b. A useful tool for studying tiny differences in bonding interaction in solid state systems is the electron localization function (ELF) [37,38]. The ELF values lie by definition between zero and one. Values are close to 1, if in the vicinity of one electron no other electron with the same spin may be found, for instance as occurs in bonding pairs. Here, we consider the difference in ELF (DELF) between the situation when f electrons are in the valence and when they are frozen in the core, for the same ionic configuration, i.e. DELF( r)=ELF f val ( r)-ELF fcore ( r). The physical interpretation is as follows: positive values of DELF show up in regions where the electron localization is higher, i.e. the bonding between FM layers is strengthened. In Fig. 3 (right part) we show a positive isosurface of DELF projected into the ab plane. Clearly, it is mainly localized between FM layers and, more specifically, in the region between Dy and O ↑,↑ eq . This points to a bonding interaction between FM layers mediated by O ↑,↑ eq .

C. Electronic structure fingerprint
A careful inspection of the orbital-decomposed magnetic moments reveals that the Dy-5d states are polarized only if the 4f states are in the valence: the 4f states couple locally to the 5d spin moments by intra-atomic 4f -5d exchange interaction [39]. The 5d states are much more extended than 4f electrons, suggesting that the glue that finally couples FM Fe and Dy layers is the interatomic interaction between the Fe 3d states and the Dy 5d states mediated by the oxygen p states. In Fig. 4 we show the local Density of States (DOS) at O ↑,↑ eq atoms and the Dy-d and Dy-f states. In panel (I), a clear interaction between d and f states of Dy and oxygen states is visible (cfr dotted ellipse). Note that the same interaction involves the Fe 3d states as well (not shown in Fig. 4) disappears when freezing the f states in the core (see panel (II)), for which the intra-atomic Dy-f -d and interatomic O-p hybridization disappears. Panel (III) shows the Dy-f DOS calculated using DFT+U , HSE and G 0 W 0 @HSE calculations. In the (relevant) occupied manifold, the chosen U nicely fits the HSE DOS, which, in turn, is rather close to the G 0 W 0 calculations. The corrections beyond a mere DFT+U approach show up in the unoccupied states by opening the band gap; however, this does not change our conclusions, as far as the polarization is concerned. The fact that the Dy and Fe interaction is mediated by O ↑,↑ eq -sp states suggests possible routes to tailor the ferroelectric polarization. For instance, compressive or tensile strain along the b axis should increase or decrease the tilting of the octahedra, favoring or disfavoring the interaction via the intermediate O ↑,↑ eq states, i.e. enhancing or reducing the ferroelectric polarization.

V. CONCLUSIONS
Several results emerge from our study: i) the FE state of DFO is driven by exchange striction, confirming the qualitative explanation given in Refs.10, 11; ii) two degenerate and switchable polar states exist characterized by a sign-reversal of the FE polarization (±P ) and linked by an adiabatic path (connecting the two ferroelectric (FE) states FE 1 and FE 2 to the same paraelectric reference structure) can be obtained through a relative rotation of the direction of Dy spins (with respect to Fe spins); iii) the coupling between Dy and Fe spin sublattices is mediated by Dy-d and O-2p states; iv) the estimated FE polarization is in agreement with experiments; (v) by freezing the f states in the core instead of relaxing them in the valence, we confirm the crucial role played by f electrons in establishing the spin-driven ferroelectricity. More generally, our study suggests that f electrons might play an important role in the ferroelectric properties of other RMO 3 compounds (where 4f electrons are often neglected in theoretical calculations); (v) last, but not least, our electronic structure analysis suggests possible routes to tailor the ferroelectric polarization, owing to the strong Dy-Fe coupling via intermediate equatorial oxygens. Further study is in progress to confirm this expectation.