Magnetic exchange interactions in yttrium iron garnet: A fully relativistic ﬁrst-principles investigation

Magnetic isotropic and Dzyaloshinskii-Moriya interactions in yttrium iron garnet have been obtained by ab initio fully relativistic calculations. The calculated coupling constants are in agreement with available experimental data. Using linear spin-wave theory, we are able to reproduce the experimental magnon spectrum including the spin-wave gap and stiffness. The way to calculate the exchange coupling constants using the Korringa-Kohn-Rostoker formalism for large magnetic systems such as complex oxides is discussed in detail. DOI: 10.1103/PhysRevB.


I. INTRODUCTION
Yttrium iron garnet (YIG) is well known as a ferrimagnetic insulator with the largest spin-wave lifetime and relatively high Curie temperature [1].Thus it is arguably the most important material used in microwave and recent spintronic devices [2].While this has been known for some time, a resurgence in interest was caused after the discovery of electrically injected spin waves into YIG and detection of their transmission over macroscopic distances [3].It was also shown that magnons, excited electrically, can diffuse over long distances even at room temperatures [4].The magnon polarization was also directly measured for the first time in YIG, giving a direct proof of its ferrimagnetic nature [5].
Models of the magnon spectrum of YIG have mostly been based on neutron scattering, covering only three of the lowest energy modes out of 20 distinct magnon branches [6].The first model covering all 20 magnon branches was obtained by timeof-flight inelastic neutron scattering [7].It was shown that the magnetic interactions were long ranged and more complex [7] than previously understood [1,6,8,9].The full magnon dispersions, investigated using neutron scattering, were reproduced by a model with three nearest-neighbor exchange integrals [10].With that being mentioned, the exchange constants reported in literature vary considerably, since they depend on Published by the American Physical Society under the terms of the Creative Commons Attribution 4.0 International license.Further distribution of this work must maintain attribution to the author(s) and the published article's title, journal citation, and DOI.Funded by Bibsam.the experimental technique as well as on the fitting procedure [1,[6][7][8][9][10].
The complexity of YIG's crystal structure makes accurate ab initio calculations of exchange coupling constants a formidable task.As a consequence, ab initio calculated magnetic exchange couplings are scarce.In previous investigations [11], the exchange constants were determined by fitting total energies of different collinear spin configurations to the Heisenberg model.Including nearest-neighbor (NN) and next-nearest-neighbor (NNN) interactions, it was found that the spin-wave spectrum of YIG is very sensitive to small changes in the exchange interactions and small changes appear to give dramatically different spectra [11].The most recent investigation was based on the noninteracting transverse spin susceptibility from quasiparticle self-consistent GW (QSGW) band structure calculations [12].This state-of-the-art calculation demonstrates a parameter-free multiscale modeling of YIG from first principles.
No fully relativistic ab initio results have been reported, and thus there has been no theoretical investigation of the antisymmetric Dzyaloshinskii-Moriya interaction (DMI) [13,14].It should be noted that there has been one attempt to include DMI in modeling of the magnon spectra in bulk YIG [7].Models including DMI were found to destabilize the magnetic structure for arbitrarily small perturbations when fitting to the experimental magnon spectrum.Nevertheless, DMI is of vital importance especially in systems with broken inversion symmetry, where it could be responsible for the formation of chiral magnetic structures.This is the case in deposited ultrathin YIG films, where DMI-induced chiral propagation of spin waves was recently observed by different experimental groups [15,16].In this paper, the isotropic and Dzyaloshinskii-Moriya interactions in YIG have been calculated within the framework of the fully relativistic Korringa-Kohn-Rostoker (KKR) Green's function formalism.The critical temperature has been estimated using the Metropolis Monte Carlo technique and the magnon spectrum modeled by linear spin-wave theory.The paper is organized as follows: Sec.II describes the methodology and computational details of the calculations.Section III presents the calculated results and discusses the physics of magnetic interactions in YIG.The main conclusions are stated in Sec.IV.

II. METHODOLOGY
The crystal structure of YIG has a bcc cubic cell (space group Ia 3d), which contains 160 atoms (80 atoms in the primitive cell; see Table I and Fig. 1) [17,18].The Y metal ions occupy dodecahedral 24c sites, and the Fe atoms occupy octahedral 16a sites and tetrahedral 24d sites.The Fe ions at the 16a sites (denoted as Fe O ) and the 24d sites (denoted as Fe T ) form two magnetic sublattices that are ferrimagnetically coupled.The structure can also be considered as an interconnection of distorted octahedra, tetrahedra, and dodecahedra with shared O atoms at the corners of the polyhedra.
The structure optimization has been performed using the projector-augmented plane-wave method (PAW) [20] as implemented in the Vienna ab initio simulation package (VASP) [21].The exchange-correlation functional was treated in the local density approximation (LDA) as well as in the generalized gradient approximation (GGA).The cutoff energy for the plane-wave basis was set to 500 eV.The convergence tolerance for the total energy was set to 10 −6 eV/at., while forces were minimized up to 10 −3 eV/Å or less.With these criteria, convergence of the k-point mesh was obtained for a 6 × 6 × 6 grid, in the Monkhorst and Pack scheme [22].
Since 3d electrons already show a noticeable degree of localization, treatment in LDA or GGA may not be fully adequate.To overcome this limitation, an on-site Coulomb correction term was included for the 3d states, as in the LDA+U or GGA+U approach [23].We employ the formulation of Dudarev et al. [24], in which the isotropic screened on-site Coulomb interaction is given by an effective interac-FIG.1. YIG unit cell.The dodecahedrally coordinated Y ions (green) occupy the 24c Wyckoff sites, the octahedrally coordinated Fe O ions (gold) occupy the 16a sites, and the tetrahedrally coordinated Fe T ions (blue) occupy the 24d sites.The oxygen 96h sites (red) are shown.This figure was created using VESTA [19].
tion parameter U eff = U − J, where U and J are the screened Coulomb and exchange parameters, respectively [25].We also explored the usage of hybrid functionals, which can improve the description of the electronic structure of insulators and reproduce experimentally observed band gaps in a wide variety of materials [26].Therefore we also performed PAW-VASP calculations with the revised Heyd-Scuseria-Ernzerhof functional HSE06 [27].
After having performed the structural relaxation of YIG, we investigated the electronic and magnetic properties in more detail.To this aim, we used the fully relativistic spin-polarized KKR method, as implemented in the SPRKKR package [28,29], which is known to give an accurate description of magnetic properties.A characteristic feature of the KKR method is that it allows for the separation of the potential aspects of a material, embodied in the scattering t matrices, from the structural aspects, represented by the structure constants of the underlying lattice, which describe electron propagation between the atoms.
YIG, having 80 atoms in the primitive cell, is an open system.To describe this system in the framework of the KKR formalism within the atomic-sphere approximation (ASA), 32 empty spheres were included in the unit cell.This setup makes it possible to minimize the overlap between all the atomic spheres, filling the space without a large geometry violation.All atomic radii including empty spheres have been rescaled in order to keep the volume fixed.It should be noted that our test calculations without empty spheres could not reproduce the electronic structure of YIG (the density of states).Moreover, there was an artificially induced magnetic moment on the O atoms, m O = 0.23 μ B /at.
For large-scale systems such as YIG, the computational complexity of the simulations can be challenging.However, this issue can be circumvented by replacing the long-range structure constants with more localized ones.This is achieved in the screened KKR (SKKR) or the tight-binding KKR (TB-KKR) method [28][29][30][31].In this approach, the singlesite scattering t matrix of the real system is replaced by its tight-binding (TB) counterpart, t TB , which is derived from a repulsive potential.One can then solve in a first step the corresponding multiple-scattering problem for this new reference system.
In contrast to the free-electron Green's function, this new auxiliary Green's function G TB (r, r , E ) decays very rapidly in space, i.e., with increasing distance |r − r |.As a consequence, the corresponding scattering path operator τ TB nn is essentially zero if the distance between sites n and n is greater than the next-nearest-neighbor distance.By applying the Dyson equation, the multiple-scattering problem for the real system can now be solved by using the Green's function of the new reference system [31].Thus this efficient scheme dealing with the auxiliary TB-reference system allows us to treat large complex magnetic systems.Now we can focus on the calculations of the magnetic exchange interactions for YIG, using the KKR Green's function method in the fully relativistic formulation.As was discussed in the literature (e.g., Refs.[32,33]), the exchange coupling parameters accounting for relativistic effects can be written in a tensorial form, J i j , giving access to the extended form of the Heisenberg Hamiltonian where the single-site anisotropy K is also included.The elements of the exchange coupling tensor in the multiplescattering formalism in Ref. [32] are given by the expression with τ i j (E ) being the matrix of the scattering path operator [34].The matrix elements T i( j),η(ξ ) (E ) in Eq. ( 2) [with η(ξ ) = x, y, z] characterize the perturbation due to infinitesimal rotation of the exchange field on site i( j), B i( j) xc , and have the following form: where Z n ( r, E ) and J n ( r, E ) are functions representing the regular and irregular solutions of the single-site Dirac equation, given in the (κμ) representation.Here, stands for (κμ), with κ and μ being relativistic quantum numbers.β is one of the standard 4 × 4 Dirac matrices, and ση(ξ ) is the corresponding spin Pauli matrix [35].The present version of Eqs.(2) and (3) restricts their use at the moment to a pure local spin-density approximation (LSDA)-based Hamiltonian for which B xc represents the difference in the spin-up and spin-down LSDA exchange-correlation potentials.However, it is straightforward to extend Eqs. ( 2) and (3) to the LDA+U or dynamical mean-field theory (DMFT) case [36].
In an alternative scheme reported by Udvardi et al. [33], the elements of the exchange tensor represented in spherical coordinates are given by the expression with where t i is the single-site scattering matrix at site i, calculated for collinear magnetic structure; R i (θ i , φ i ) is a rotation matrix; and α i = (θ i , φ i ).As one can see, the latter scheme is associated with the perturbation due to the infinitesimal rotation of the single-site scattering t matrix that accounts for all the exchange and correlation effects.This makes this scheme suitable for calculations based on LDA+U or GGA+U for strongly correlated materials (e.g., Refs.[37,38]).Therefore all discussion below concerns the exchange parameters calculated using the approach as formulated in Ref. [33].
Representing the exchange tensor J in the Cartesian frame of reference, the Hamiltonian in Eq. ( 1) can be reduced to a more convenient form given by the expression with the isotropic exchange coupling parameters J i j = 1 3 (J xx i j + J yy i j + J zz i j ) and the DMI vector D i j with the components , where ηξ ζ is the Levi-Civita tensor [39].
The KKR calculations were performed using an angular momentum cutoff of l max = 3 for the expansion of the Green's function.The Brillouin zone was sampled with a 10 × 10 × 10 k-point mesh.The energy integration of the Green's function was done with 32 energy points along a semicircle contour.The convergence of the results was checked by increasing the k-point grid and the energy points.The exchange-correlation potential was modeled in LSDA, using the parametrization of Vosko, Wilk, and Nusair (VWN) [40].

A. Structure and ground-state properties
The optimized atomic positions were obtained by PAW-VASP at the experimental atomic volume, using the GGA in the parametrization by Perdew, Burke, and Ernzerhof (PBE) [41].The results are reported in Table I, in comparison to the experimental crystallographic data [17,18].The calculated Wyckoff positions of Fe O , Fe T , and Y agree very well with the experimental data, and the small discrepancies observed for the Wyckoff positions of O are within the expected computational accuracy.
The relaxed atomic positions were used to calculate the ground-state properties by the PAW and KKR methods.The experimental lattice parameter was used for the KKR calculations including the calculations of the exchange coupling.In the PAW calculations, however, the lattice parameter was optimized.An overview of the calculated properties compared with experimental data is given in Table II.A more complete overview of available experimental [10,[42][43][44][45][46][47] and theoretical [11,[48][49][50][51][52][53][54][55][56] data for YIG can be found in Appendix A.
The theoretical lattice parameters obtained in GGA and HSE06 are in good agreement with the experimental values, considering also the fact that ab initio calculations were performed for 0 K.As one can often observe, LDA-based calculations tend to underestimate the lattice parameter, while GGA leads to a slight overestimation.Overall, the hybrid functional HSE06 gives the best agreement with the available experimental values.This is consistent with the fact that YIG is a ferrimagnetic insulator, and therefore its physical properties are strongly affected by an improved description of the exchange mechanism.It should be noted that for thin garnet films the lattice constants can depend on the thicknesses of the films as well as on the type of substrate [57].Therefore one should not expect a perfect quantitative agreement between bulk crystals and thin films.
LDA calculations predict YIG to be a metal, in contrast to the experimental observation that YIG is an insulator [56].This is because the LDA theory fails to account for the strong intra-atomic correlation of the localized Fe 3d electrons.The inclusion of the intra-atomic correlation by LDA+U or GGA+U allows for the opening of a band gap in YIG.The dependence of the band gap on U eff is illustrated in Fig. 2, as obtained from LDA+U and GGA+U calculations.The values for U eff = 0 eV correspond to the standard LDA and GGA values.The result from HSE06 calculations is also reported.The experimental band gap reported in the literature is around 2.9 eV, as can be inferred from the values of 2.81 eV (magneto-optical spectra at 300 K [46]), 2.89 eV (optical absorption and Faraday rotation at 6 K [46]), and 2.85 eV (electrical measurements, in the band model [45]).Figure 2 illustrates that GGA+U calculations predict a band gap within the experimental range only for very high U eff values.The need for high values of U eff is supported by Ref. [58], in which the optimal values (9.8 eV for Fe O and 9.1 eV for Fe T ) were obtained by constrained density functional theory and linear response theory.The hybrid functional, HSE06, provides a good description of the experimentally observed band gap without external parameters.It should be noted that, because LDA calculations predict YIG to be a metal in contrast to the experimental observation, one could expect magnetic moments in LDA to be smaller than for other schemes.
In order to discuss the magnetic properties of YIG, we present in Table II the calculated spin magnetic moments of Fe, obtained within the VASP and SPRKKR calculations, in comparison with the experimental results.As one can see, the magnetic moments at the Fe O and Fe T sites, obtained within the GGA+U approach with U eff = 7.1 eV, lead to a good agreement with the data obtained by neutron diffraction experiments [55] and based on the measurements of electronic magnetic circular dichroism (EMCD) [47].On the other hand, the calculations which do not include a specific correction for the local Coulomb repulsion substantially underestimate the magnetic moments, indicating a crucial impact of local correlations on the magnetic properties of YIG.The dependence of Fe spin magnetic moment on U eff has therefore been systematically investigated.The corresponding results are shown in Fig. 3  moments are larger in GGA+U than in LDA+U for the same experimental lattice constant.
The orbital contribution to the magnetic moment is almost completely quenched due to crystal field effects.Even so, the small orbital moments show a U eff dependence.For U eff smaller than 4 eV, Fe O has a larger orbital moment which decreases with increasing U eff , while Fe T shows the opposite trend.The orbital moments on both Fe sublattices converge towards the same value for U eff larger than 4 eV.Our calculated values for the orbital magnetic moments are close to previous ab initio calculations [51].The only available experimental data on the orbital magnetism are obtained via EMCD [47] (see Appendix A).

B. Magnetic exchange interactions
First of all, considering the interatomic exchange interactions in YIG, we focus on the interactions between the magnetic atoms, i.e., Fe, characterized by robust local spin magnetic moments.This exchange can be associated with the dominating superexchange mechanism [59], as it often happens in transition-metal oxides, with oxygen p orbitals responsible for the mediation of the Fe-Fe exchange interactions.This type of interaction leads usually to the antiferromagnetic coupling, although it can also result in a ferromagnetic coupling, according to the Goodenough-Kanamori rules [60,61].
We shall now focus on the Fe-Fe interatomic exchange interactions, obtained within the KKR formalism.Describing the geometry for the exchange interaction, Fig. 4 presents the first six nearest neighbors in 1/4 of the ferrimagnetic YIG unit cell.J ad , J dd , J aa , J ad , and J dd correspond to the first five coordination shells, respectively, while J aa and J dd correspond to the sixth shell.
Additional contributions to the Fe-Fe exchange interaction can appear due to the induced spin magnetic moment on the oxygen atoms.In this case, an interaction between the Fe and O magnetic moments occurs, leading to an effective Fe-O-Fe exchange interaction as discussed in Refs.[62,63].Following Refs.[62,63], the effective magnetic exchange interactions,

J eff
i j (which contain, besides the Fe-Fe exchange, the Fe-O interactions as well), can be written as The i and j label Fe sites, k and l label O sites, and J ik represents the Fe-O exchange interactions between Fe site i and O site k.Before we discuss in detail the magnetic exchange interactions in a YIG crystal, one has to stress once again some technical aspects, which are crucial for accurate calculations of the interatomic exchange interactions at long distances.As discussed above, the large size of the unit cell for this material leads to a daunting challenge for accurate computations due to the long-range spatial dependence of the KKR structure constants in real space.Therefore the TB-KKR method was used in the present calculations, which allows an efficient treatment of large systems due to the rapid convergence of the KKR structure constants.Moreover, it should be noted that the inclusion of Fe-O exchange interactions in the so-called effective magnetic exchange interactions [62,63] improved the agreement with the experimental data for the oxygenmediated coupling in YIG and decreased the scatter of values of the exchange interactions.
Having set up a suitable approach to calculate accurate interatomic exchange interactions, we have to decide upon suitable values for the Coulomb interaction parameters.The choice of U eff is a fundamental problem per se but becomes particularly tedious for large systems such as YIG [64,65].In this paper, these parameters have been chosen from the estimation of the critical temperature T C (Fig. 5) and the spin-wave stiffness constant D (Fig. 6) and their comparison with available experimental data.
The critical temperature was obtained by solving for the magnetic spin configuration governed by the isotropic part of Eq. ( 6), using the Metropolis Monte Carlo technique.This approach allows for a proper treatment of thermal fluctuations, without the approximation of a mean-field approach.The simulation box considered for the calculations contained 10 240 atoms (an 8 × 8 × 8 supercell based on the 20-Fe-atom unit cell).Test calculations for larger boxes show that this size is large enough to avoid finite-size effects, in agreement with Ref. [66].It should be mentioned that the effect of the DMI on the Curie temperature and spin-wave stiffness is negligible, while it increases the computational complexity significantly.This is the reason why the results presented here include only the isotropic part of the exchange.[7], Shamoto et al. [10], and Man et al. [71].Theoretical data are from Cherepanov et al. [1] and Xie et al. [11].
The calculated critical temperature as a function of U eff is shown in Fig. 5.The trend observed with respect to U eff is well defined.The larger U eff is, the lower is T C .This is due to the fact that a large U leads to more localized orbitals, which results in larger moments (see Fig. 3) but reduced exchange interactions.With U eff = 4.1 eV, T C is 590 K, in good agreement with the experimental value of 560 K [67].
The spin-wave stiffness constant D was obtained from low-temperature magnon spectra that were modeled by linear spin-wave theory (LST) as implemented in SPINW [68,69] (see also Appendix B).In this approach, the low-temperature spectrum is calculated by diagonalizing the Hamiltonian, which should be identical to the classical solution up to a finite zero-point energy shift at low temperatures.We add a small external field B z = 0.01 T to specify a magnetization (spin quantization) axis.Using a least-squares fit of the quadratic energy dispersion relation, hω = Dq 2 , the spin-wave stiffness could be estimated.Since this approximation is only valid for small q around the point, the interval used in the fitting is decreased until convergence is obtained.The effect of the DMI on the spin-wave stiffness is negligible.Figure 6 shows the spin-wave stiffness constant D as a function of U eff .The stronger the Coulomb interaction is, the softer the spin waves are (see Fig. 9 in Appendix B).This will be discussed in more detail in the next section.For U eff values in the interval 4.1-6.1 eV, the calculated D values are in the range of the experimental values [7,10,70,71] and previously reported theoretical values [1,11,12].
In conclusion, the previous analysis allows us to estimate that U eff = 4.1 eV are the most appropriate Coulomb interaction parameters for YIG.They provide both critical temperature and spin-wave stiffness in agreement with experimental and theoretical values reported in the literature.Moreover, they lead to a good estimate of magnetic moments (Fig. 3) and band gap (Fig. 2) as well.The latter is, however, slightly underestimated.
In the following, we focus on the analysis of the interatomic exchange interactions for the chosen Coulomb interaction parameters.Values obtained by means of the TB-KKR approach, with and without the renormalization due to the Fe-O terms [Eq.( 7)], are reported in Table III and Fig. 7, together with values extracted from experimental data [1,6,7], as well as values acquired from previous calculations [11,12].One should keep in mind that the experimental values have a rather large uncertainty, since they are affected by the type of experiment and by the choice of the empirical model used for extracting the exchange interactions.In addition, one should consider that exchange parameters from different models are compared.Nevertheless, the experimental and calculated results, presented in Table III, are in reasonable agreement for the first six coordination shells [6][7][8]10].This is a further proof that our procedure for setting an appropriate value for U eff provides a good description of the magnetic properties in YIG.
Figure 7 shows that the dominant exchange Fe-Fe interaction is between nearest-neighbor Fe sites.This is an antiferromagnetic coupling and is responsible for results in an antiparallel alignment of the magnetic moments of Fe atoms, placed in the octahedral (a) and tetrahedral (d) sites.The TABLE III.Magnetic exchange interaction parameters (meV).See Fig. 4 for the geometry and notation.Experimental and theoretical data include Princep et al. [7], Shamoto et al. [10], Plant [6], Cherepanov et al. [1], Xie et al. [11], and Barker et al. [12].Positive values of exchange interactions correspond to ferromagnetic ordering, while negative values refer to antiferromagnetic ordering.n is the number of coordination shells.k n is the number of Fe atoms in the nth coordination shell.

This paper
Ref. [11] TB-KKR (eff), NN NNN U -J = 4.1 eV model model ferrimagnetism in this system originates from the unequal numbers of a (minority) and d (majority) Fe sites.For a schematic view of the geometry of the magnetic interactions in YIG, we refer the reader to Fig. 4. Experimental (exp) data are from Plant [6] (1977), Plant [8] (1983), Princep et al. [7], and Shamoto et al. [10].Theoretical data are from Cherepanov et al. [1], Xie et al. [11], and Barker et al. [12].Vertical dashed lines show the coordination shells for Fe atoms.More information can be found in Table III.NS fit refers to the neutron scattering (NS) data fit.
The number of Fe neighbors for each of the first eight coordination shells is given in Table III.Each Fe atom in an octahedral position has six nearest-neighbor Fe atoms in tetrahedral positions, while each Fe atom in a tetrahedral position has four nearest-neighbor Fe atoms in octahedral positions.The second-nearest neighbors are four Fe atoms in tetrahedral sites.Note that the distances for first and second coordination shells are close to each other (see also Fig. 7).
According to Ref. [7], there are two distinct third-nearestneighbor bonds (J aa ), which have identical length but differ in symmetry.Our results show that the exchange coupling constants for the third-nearest-neighbor shell have the same values, i.e., there is no second solution.However, the sixthnearest-neighbor bonds correspond to two types of exchange constants, J aa and J dd , having six and two neighbors, respectively (see Fig. 4).As can be seen from Table III and Fig. 7, these exchange couplings have different values.Note that the magnetic exchange couplings for the sixth coordination shell are small; however, it was shown that even small features of the exchange couplings can be important for describing the experimentally observed magnon spectrum [7].The reason is that exchange coupling interactions are determined by the local environment of the Fe sites.Our results support the conclusion made in Ref. [7].
Using the TB-KKR setup described earlier, we also calculate the Dzyaloshinskii-Moriya interactions (DMIs).The obtained values for the DMI vectors D i j in the z direction are D z 01 = 0.14 meV, D z 02 = 0.10 meV, and D z 03 = 0.03 meV, respectively.The DMIs for the atoms which are farther away than the third coordination shell are negligible.The calculated DMIs for the first neighboring Fe atoms are much smaller than the isotropic exchange interactions (about a few percent), which is consistent with their relativistic origin.While DMI is generally not expected for ordered systems with a center of inversion symmetry, it can occur for complex structures [72].The reason why it appears in YIG is because of the large unit cell with magnetic ions in different octahedral and tetrahedral oxygen surroundings.Hence DMI is the result of local symmetry breaking in YIG due to the different environments.While the antisymmetric exchange interactions are generally expected to be proportional to the symmetric ones [14,72,73], the magnitude of the calculated Fe-Fe DMI for the second Fe coordination shell is rather significant.The slow decrease of the DMI with increasing distance, when compared with the isotropic exchange interactions, has been discussed in the literature [32].
A short comment on the magnetocrystalline anisotropy energy is appropriate.In bulk YIG, due to its cubic structure [74], the magnetocrystalline anisotropy energy is known to be small, and therefore it is not addressed in this paper.However, in thin YIG films, the magnetocrystalline anisotropy can be tuned by strain.In the presence of strain-induced tetragonal distortion, magnetocrystalline anisotropy has been induced in thin YIG films grown on different substrates [75].Moreover, perpendicular magnetic anisotropy was obtained in epitaxial YIG thin films [76].

C. Magnon spectrum
The calculated exchange interactions for the established Coulomb interaction parameters (see Sec. III B) can be used in the SPINW code [68,69] to obtain the spin-wave spectrum.The effective spin Hamiltonian included a bilinear Heisenberg exchange term and a Dzyaloshinskii-Moriya term, while the cubic anisotropy energy was neglected, due to its almost vanishing value.
The calculated spin-wave spectra are shown in Fig. 8, together with experimental data from Refs.[6,7,10,71], taken at various temperatures (10, 83, 100, and 295 K).The top panel contains the spectrum calculated from the imaginary part of the spin-spin correlation function S xy (Q, ω) − S yx (Q, ω).The spin-wave dispersions are plotted as black lines, while the intensity convoluted in energy with a 0.75-THz full width at half maximum (FWHM) Gaussian is shown in red and blue colors.Red and blue indicate the positive and negative chirality (polarization) of the modes, respectively.The intensity (thickness) of the lines gives an estimate of the magnon density.For the sake of comparison, the adiabatic spectrum obtained using exchange couplings extracted from recent experimental studies [7,10] is also presented in the bottom panel.Note that only the low-temperature data should be compared with the current theoretical model.
Let us discuss the dispersion relation of the acoustic mode.At low energies, the acoustic mode has a parabolic behavior.The slope of the calculated acoustic mode (Fig. 8) agrees with the neutron-scattering data for low temperatures from Refs.[10,71]; however, it differs slightly from the results of Ref. [6].The calculated spin-wave stiffness constant D, which is governed by the second derivative of the energy dispersion at the point, is in the range of experimental values for the chosen value of U eff (Fig. 6).It should be FIG.8. Top: calculated spin-wave spectrum compared with experimental data from Plant [6] (295 K, cyan filled squares; 83 K, magenta filled circles), Shamoto et al. [10] (295 K, white filled circles), and Man et al. [71] (10 K, green filled circles).Bottom: adiabatic spin-wave spectrum, obtained using exchange couplings from this paper (black thick lines) and from recent experimental studies by Shamoto et al. [10] (red lines) and Princep et al. [7] (blue lines).Here, r.l.u., reciprocal lattice units.
mentioned that the values for D reported in the literature, obtained by different experimental methods and under different conditions, have a wide scatter.At higher energies, instead, the magnon dispersion starts deviating from the quadratic function, in good agreement with the reported experimental results [10].
As shown in Fig. 8, the highest-frequency mode of the calculated spectrum is close to results obtained in Ref. [10].However, the other high-frequency modes, between 22 and 24 THz, differ from those reported in Ref. [10].For the middle frequency range of the spectrum, 12-19 THz, the calculated spin-wave modes are placed between branches generated from magnetic interactions extracted from experiments [7,10].
The spin-wave gap, which is defined as the energy difference between the two lowest modes (acoustic and first optical mode) at the point, can be extracted from Fig. 8 and is about 9 THz.This value agrees well with the experimental data from Refs.[7,10].On the other hand, the experimental value from Ref. [6] is about 8 THz.This lower value is most likely due to the fact that the measurements were performed at a higher temperature, namely, 83 K.It should also be mentioned that the estimated value is an improvement compared with the latest theoretical investigations [11,12,77].
The spin-wave dispersion spectrum consists of many modes, in contrast to only two modes as might naively be expected based on the cubic crystal structure and the ferrimagnetic ordering of the compound.This magnon spectrum reflects the complexity of the magnetic structure of YIG.
A final comment regards magnon-phonon coupling.Magnetoacoustic phenomena appear much brighter than in usual ferromagnets because the magnetoelastic coupling is enhanced by the exchange interlattice interactions.In YIG, the magnon-phonon coupling is responsible for many phenomena and effects that might be useful for future technologies [4,54,71,[78][79][80][81].This interesting and important issue is beyond the scope of this paper and is left for future investigations.

IV. CONCLUSIONS
Exchange coupling constants and Dzyaloshinskii-Moriya interactions in yttrium iron garnet have been obtained by ab initio fully relativistic calculations.The calculated isotropic coupling constants are in agreement with available experimental data.The sixth-nearest-neighbor bonds correspond to two types of exchange constants between Fe atoms in octahedral (a) and tetrahedral (d) positions.The reason is that there are two distinct sixth-nearest-neighbor bonds, which have identical length but differ in the symmetry of the local environment of the Fe sites.The calculated Dzyaloshinskii-Moriya contributions to the exchange are short ranged.The screened KKR formalism seems to be a very promising tool for the calculations of the exchange coupling constants in complex magnetic oxides.The calculated effective magnetic exchange interactions due to the oxygen-mediated coupling in YIG can improve the agreement with experimental data.It should be noted that the hybrid functional, HSE06, improves the description of the electronic structure of YIG insulators and reproduces the experimentally observed band gap in YIG.
The spin-wave dispersion has been modeled using linear spin-wave theory with the calculated exchange parameters.The inclusion of Dzyaloshinskii-Moriya interaction showed a negligible effect.The calculated spin-wave stiffness is in the range of experimental values.The magnon spectrum agrees with experimental results, and the spin-wave gap is well reproduced.
The present investigation has been performed for the magnetically ordered state, but the description of the magnetic state for YIG at finite temperatures has not been included.The questions of whether the effect of the magnetic disorder on the magnetic exchange interactions is large or not, and how to correctly take into account the correlation effects (and/or define U eff ) call for further investigations.Nevertheless, we believe that the obtained results correctly describe qualitative features of the magnetic exchange interactions in YIG. from the JRG Program at APCTP through the Science and Technology Promotion Fund and Lottery Fund of the Korean Government, as well as from the Korean Local Governments-Gyeongsangbuk-do Province and Pohang City.J.M. would like to acknowledge the project CEDAMNF (CZ.02.1.01/0.0/0.0/15_003/0000358) of the Ministry of Education, Youth and Sports (Czech Republic).S.M. and H.E. acknowledge financial support by the DFG via SFB 1277 (Emergent Relativistic Effects in Condensed Matter -From Fundamental Aspects to Electronic Functionality).

APPENDIX A: LITERATURE DATA FOR YIG
Table IV displays a comprehensive list of available experimental and theoretical data for YIG.The abbreviations used in the table are defined in Sec.III of the main text.The orbital magnetic moment m L values, reported in the EMCD experiment [47], are 0.28 ± 0.03 and 0.31 ± 0.03 μ B /at., for the Fe O and Fe T sites, respectively.This large discrepancy with our results might be due to the different orientation of the magnetic moments in the experiment and ab initio calculations.In the calculations, the Fe orientation of the magnetic moments was taken along the [001] direction and not relaxed.In contrast, the measurements show the magnetic moments' alignment within the 111 plane.

APPENDIX B: MAGNON SPECTRA AS A FUNCTION OF U
We simulated the spin-wave spectra of YIG at low temperatures by means of SPINW [69]. Figure 9 shows the acoustic branch at low energy as a function of U (J = 0.9 eV) for the N--H direction.The larger the U , the softer the lowenergy ferromagnetic acoustic mode is.The calculated spectra were used to estimate the spin-wave stiffness constant D (see Fig. 6).
FIG. 2. Band gap (eV) of YIG as a function of U eff = U − J.The theoretical values were obtained using different exchange-correlation functionals as implemented in VASP and SPRKKR.Experimental values of the band gap are represented by the gray bar (see Table IV in Appendix A).

FIG. 3 .
FIG. 3. Total (m tot = m S + m L , top panel) and orbital (m L , bottom panel) magnetic moments of Fe atoms as a function of U eff = U − J.Note that for PAW calculations, m tot = m S .

FIG. 4 .
FIG. 4. One-fourth of the ferrimagnetic YIG unit cell.The solid lines represent the first six next-nearest-neighbor exchange interactions.The subscripts aa, dd, and ad are Fe O -Fe O , Fe T -Fe T , and Fe O -Fe T interactions, respectively.In the upper left corner, magnetic octahedral and tetrahedral sites are shown, and they compose the majority and minority spin sublattice of the ferrimagnet, respectively.

FIG. 5 .
FIG. 5. Monte Carlo (MC) results for the critical temperature as a function of U eff = U − J.The experimental critical temperature is shown as a black solid line [67].

TABLE I .
Crystal structure of YIG.The experimental structure (upper half of table) and relaxed atomic positions (lower half of table) are given.

TABLE II .
Overview of the ground-state properties of YIG, including lattice constant a(Å) and the total magnetic moment at relevant sites, Fe O and Fe T .Calculated values are reported alongside available experimental data.
a Optimized lattice constant obtained from ab initio calculations.

TABLE IV .
Calculated properties of YIG compared with available data from theoretical and experimental investigations.The calculated lattice constant a (Å), band gap E g (eV), and total magnetic moment (μ B ) are presented.For relativistic SPRKKR calculations, the orbital magnetic moments m L are shown in parentheses.FLAPW, full-potential linearized augmented plane-wave method; PBEsol, revised PBE exchange-correlation density functional.
a Optimized lattice constant from ab initio calculations.