Spin-current driven Dzyaloshinskii-Moriya interaction in the multiferroic BiFeO3 from first-principles

The electrical control of magnons opens up new ways to transport and process information for logic devices. In magnetoelectrical multiferroics, the Dzyaloshinskii-Moriya (DM) interaction directly allow for such a control and, hence, is of major importance. We determine the origin and the strength of the (converse) spin current DM interaction in the R3c bulk phase of the multiferroic BiFeO3 based on density functional theory. Our data supports only the existence of one DM interaction contribution originating from the spin current model. By exploring then magnon dispersion in the full Brillouin Zone, we show that the exchange is isotropic, but the DM interaction and anisotropy prefer any propagation and any magnetization direction within the full (111) plane. Our work emphasizes the significance of the asymmetric potential induced by the spin current over the structural asymmetry induced by the anionic octahedron in multiferroics such as BiFeO3.

Introduction BiFeO 3 (BFO) is one of the few singlephase multiferroics which exhibit a large spontaneous polarization and a long range magnetic order at roomtemperature. BFO has an antiferromagnetic (AFM) texture that can be approximated locally by a G-type order in its R3c ground state with a Néel temperature of 643 K [1]. R3c BFO exhibits a polarization of about 90 µC/cm 2 along the pseudo-cubic < 111 > symmetry equivalent directions with a Curie temperature of 1123 K [2][3][4][5].
In bulk, the AFM order is modified by the Dzyaloshinskii-Moryia (DM) interaction [6,7] creating a spin spiral that propagates along the [110] direction (known as type-I cycloid), with magnetic moments lying in the plane formed by the polarization and the propagation direction. The periodicity of this spin spiral is 62 nm [8,9]. Additionally, another propagation direction along [112] has been reported in BFO, referred to as type-II cycloid [10][11][12][13][14].
The AFM spin spiral couples to the polarization via the magnetoelectric (ME) effect. The ME effect can have multiple origins [15]. The DM interaction directly couples atomic displacements to the AFM spin spiral and opens up magnons electrical control [16]. This effect is at the center of a new logic device which intents to electrically control out-of-equilibrium spin spirals -also called magnons -in BFO to transport and process information [17,18]. The exploration of the stability mechanisms of different types of spin spirals is of paramount importance and has been subjected to a lot of research.
A phenomenological model based on Lifschitz invariant established that the DM interaction stabilizing the spin spiral had the form α(P ) · L i × L j where P is the polarization and L is the AFM vector [19]. This phenomenological model was completed by a microscopic model based on tight binding approximation which attributed the presence of a non-zero polarization to the presence of a spin spiral. The hybridization between the d -orbitals of the metal ions and the p-orbitals of the Oxygen results in the polarization P ∝ e ij ×(S i ×S j ) where S i and S j are spins on site i and j, respectively and e ij is the unit vector between site i and j [20]. The link between this polarization and the presence of the DM interaction was explicitly written by Rahmedov et al. where the spin spiral in BFO was explained based on the presence of two chiral interactions of different symmetries [21]. The first term D wFM ∝ (ω i − ω j ) · (S i × S j ) couples the oxygen octahedra tilts ω to the magnetic texture and stabilizes the magnetic moment in the (111)-plane perpendicular to the tilts rotation vectors. This DM contribution induces the weak ferromagnetic moment in BFO [22]. The second term has the form D SC ∝ (u i × e ij ) · (S i × S j ) and couples the polarization direction u i to the magnetic moments and favors the stabilization of the magnetic moments perpendicular to the (111) plane. By varying these energy terms, several types of spin cycloids have been predicted in BFO [23].
To compare the energies of these different spin spirals, both DM contributions must be computed from DFT. D wFM was first predicted to create the weak magnetic moment in BFO. Therefore it should lie along the [111] direction [22]. This DM contribution was obtained from DFT calculations (146 µeV [24] up to 304 µeV [25]) in relative good agreement with the experimental value (163 µeV [26]). Note that a full parametrization of BFO in the R3c phase obtained from experimental measurements is given in Refs. 27 and 28. They found D SC = 110 µeV and D wFM = 50 µeV which are significantly lower than previous measurements and calculations. To our surprise and to the best of our knowledge, D SC has not yet been obtained from DFT calculations which does not allow to discuss the stability of the different spin spirals.
Here, we determine the origin and the strength of the spin current DM interaction in BiFeO 3 from DFT and  show that in this type of multiferroics, the DM interaction originates from the asymmetric potential within the Fe cations rather than the structural distortions induced by the O anionic octahedra. We calculate the energies of spin spirals E(q) for different propagation directions q in the full pseudo-cubic BZ which shows an isotropic exchange interaction. The spin-orbit coupling (SOC) contribution is fully quenched for spin spirals propagating along the polarization direction which only agrees with the spin current model [20,21] and excludes the Levy and Fert model [29]. Finally, all magnetic interactions e.g. the magnetic exchange, the spin current model and the anisotropy lead to a degeneracy of spin cycloids within the (111) plane which suggest that both type-I and type-II cycloid could coexist in bulk R3c BFO.
Methods We have calculated the energy dispersions E(q) of flat spin spiral states. These are the general solution of the Heisenberg model on a periodic lattice and can be characterized by the spin spiral vector q. This vector determines the propagation direction of the spin spiral. A magnetic moment S i at an atom position r i is given by S i = m cos(q · r i ), sin(q · r i ), 0 where m is the magnitude of the magnetic moment (see supplemental material [30]). We choose q along several high symmetry directions of the rhombohedral BZ. For simplification, we present all data within the pseudo-cubic Brillouin zone, shown in Figure 1 (a). The different directions are drawn as colored arrows, starting from the R point. Every high symmetry point in the BZ is connected to a certain collinear state -sketched in Fig. 1 (b).
Results The calculated energy dispersion (points) obtained without SOC is presented in Fig. 2 (a) where the paths are R → X → R → M → R → Γ → R (cf. lower x axis). On the upper x axis, the respective directions according to Fig. 1 can be seen. Without SOC, the dispersion shows an energy minimum at the R point (Gtype AFM state) [31] and the Γ point (the FM state) the [110] [111] [111] [112] [112] [210] [210] [100] [100] [111] [111] highest energy. The energy dispersion is mapped onto an extended Heisenberg model H ex = − ij J ij (S i · S j ) to determine the strength of the Heisenberg exchange interaction parameters J ij [lines in panel (a) and (b)] beyond nearest neighbors (see supplemental material [30] for detailed information). The exchange between nearest neighbors is dominant ∼ −27 meV/Fe atom capturing the large energy differences between the high symmetry points, but seven neighbors are necessary to describe the curvature around the R points (values can be found in Table I). A closeup is shown in Fig. 2 (b). The exchange interaction results in a good description of the DFT calculations. Note that the energy differences are very small and a numerical error on this energy scale is expected. Nevertheless, all the energies (data and fit) around the collinear state at the R points have the same curvature, e.g. the same effective exchange interaction [32]. Hence, the exchange in BFO is fully isotropic.
To determine the strength of the Dzyaloshinskii-Moriya (DM) interaction, we apply spin-orbit coupling (SOC) in first order perturbation [33] to every calculated point of Figure 2 [30]. The resulting energy contribution due to SOC, ∆E SOC is illustrated in Figure 3. To quantify the strength of the DM interaction, we apply the model of Levy and Fert (LF model) [29] which is sketched in panel (a) of Fig. 3. The LF model is typically explained as a superexchange mechanism between two magnetic atoms S i , S j via a third non-magnetic atom that holds a large spin-orbit coupling. Due to the direction dependent scattering of the electrons, non-collinear magnetic structures with a specific sense of rotation (clockwise or counterclockwise) are preferred. Using the symmetry of the DM vector D ij ∝ Ri×Rj |Ri×Rj | where R i,j are the position vectors of magnetic moments S i,j with respect to the atom of large SOC, it is possible to fit every direction in panel (c) separately and to determine the DM interaction (cf. [30]). This is the reason, why in other publications, the magnitude for the DM vector can be obtained from first principles. However, since our calculations include different high symmetry directions including the [111] direction, a combined mapping of the DM vector to our data is not reproducing the DFT calculations [see red line in Fig. 3 (c)]. Hence, the model of Levy and Fert cannot be used to describe the DM interaction in BFO. Note that the DM interaction D wFM which couples the oxygen octahedra tilts to the magnetic structure and explains the weak FM canting coincides with the LF model. It would give rise to an energy contribution of SOC in the [111] direction. However, our tests reveal completely vanishing energies in the direction of polarization.
In the literature, another term is shown to describe the non-collinear canting in BiFeO 3 , the so called (converse) spin-current model [13,20,21,34]. This magnetoelectric effect is based on the spin current acting on non-collinear magnets and leads to the polarization P ∝ e ij × (S i × S j ). Following the convention of Ref. 21, the SC DM interaction is described via where u is the unit vector in the direction of polarization, e ij the unit vector connecting magnetic moments S i,j at sites i, j. The parameter C ij describes the strength of the DM interaction. The difference between the DM vector D ij and the SC vector C ij (u × e ij ) lies within the symmetry. Depending on the symmetry of the system however, both vectors can coincide, e.g. in two-dimensional interfaces. In BiFeO 3 both vectors differ -most noticeably in the [111] direction -and hence the two models lead to different results modeling the DM interaction. Mapping Eq. (1) to the energy contribution ∆E SOC gives rise to the blue lines in Fig. 3 (c) and (d). The SC model describes the DFT results almost perfectly. Here, the fit contains three neighbors, where the contributions of the second and the third neighbors are a minor correction to the first neighbor (cf. Tab. I). Note that even an effective nearest neighbor approximation for C ij , cf.
C eff shows a qualitatively very good agreement with the DFT calculations (cf. supplemental material [30]).
In panel (d) of Figure 3, we see the energy differences between the different pseudo-cubic propagation directions of spin spirals around the ground state. Both data and SC model exhibit the steepest slope in the [110] and [112] directions meaning the DM interaction for these two directions is the largest. Furthermore, for both directions, points and lines show complete energy degeneracy.
[110] and [112] referring typically to the type-I and type-II cycloids, respectively, lie within the (111) plane. Our tests reveal that any direction within the (111) plane is equally preferred by the spin-current DM interaction. That means, including exchange and spin-current DM interaction restricts the possibility of spin cycloid propagation directions in BiFeO 3 bulk to any direction of the two-dimensional (111) plane. This is in accordance with Ref. 28, where it is also stated that q can point along any direction in the hexagonal (111) plane.
For a detailed insight into the origin of the DM interaction in BFO, we calculated the element resolved energy contribution due to SOC (Fig. 4). Here, the gray points show the total energy contribution from Fig. 3 (c) and the colored points represent the elements of BiFeO 3 (lines serve as guide to the eye). In panels (b), (c), (d) the energy contributions due to Bi, Fe and O, respectively are presented. Since in the rhombohedral unit cell, two Bi, two Fe and six O atoms are used, a sketch of BiFeO 3 in cubic representation is shown in (a) with the same color code for each element as in the graphs. Even though Bi is the element of large SOC and the O atoms break the inversion symmetry, their contribution to the total DM interaction in BFO is negligible [panels (b) and (d)]. This emphasizes that the model of Levy and Fert relying on an heavy metallic element with large SOC and structural asymmetry does not provide a good description for BFO. The total DM interaction -as a sum of all contributions -is represented by the whole contribution of both Fe atoms. This is similar to the Rashba DM interaction observed in on 3d unsupported monolayers (UMLs) under an external electric field [36] and in Graphene/Co(0001) due to potential gradient between graphene and Co [37]. It creates an internal asymmetry of the potential leading to a non-vanishing DM interaction even though the structure itself does not hold an asymmetry. In BFO, the internal electric field creates the same asymmetry within the potential of the Fe atoms and hence the DM interaction is driven by the spin current of the non-collinear magnetic structure.
As mentioned above, the exchange is isotropic within the full 3D BZ, whereas the spin current driven DM interaction narrows the spin cycloid directions down to the [110] [111] [111] [112] [112] [210] [210] [100] [100] [111] [111]  [20,21]. In systems with a polarization P, the spin-current vector u × ei,j is perpendicular to the direction of polarization. (c) Energy contribution due to spinorbit coupling ∆ESOC to the energy dispersion of spin spirals along different directions within the pseudo-cubic Brillouin zone [ Fig. 1  shown that the small easy (111) plane anisotropy in R3c BFO is a result of both strong out-of-plane anisotropy (in [111] direction) driven by the Bi-Fe ferroelectric displacement in competition with a strong easy plane anisotropy (perpendicular to the [111] direction) stemming from the oxygen octahedra tilts. Here, we calculate the anisotropy energy using a 360 • rotation within the (111) plane. The magnetocrystalline anisotropy energy (MAE) is not changing for any direction in this plane [30]. In bulk BFO, the MAE does not prefer a certain magnetization direction within the (111) plane.
[111]   Conclusion Our work demonstrates based on density functional theory that the Dzyaloshinskii-Moriya interaction in R3c bulk BiFeO 3 is governed by the spin-current model stemming from the non-collinear antiferromagnetic structure. The off-centered displacement of Fe and Bi induces an asymmetric shape of the internal potential in the Fe atoms and consequently, the Fe atoms carry almost the whole contribution to the DM interaction. This effect is not restricted to BFO in particular and should occur in other multiferroic materials. In the case of BFO, by including exchange, DM interaction and anisotropy, we show that any spin cycloid propagation direction in the (111) plane are energetically degenerate. This could explain an effect of continuously rotating cycloids in BFO that has recently been observed [39] and explained as a surface effect. Based on our results, this effect should also be measurable in bulk.