Magnetic ground state of SrRuO3 thin film and applicability of standard first-principles approximations to metallic magnetism

A systematic first-principles study has been performed to understand the magnetism of thin film SrRuO3 which lots of research efforts have been devoted to but no clear consensus has been reached about its ground state properties. The relative t2g level difference, lattice distortion as well as the layer thickness play together in determining the spin order. In particular, it is important to understand the difference between two standard approximations, namely LDA and GGA, in describing this metallic magnetism. Landau free energy analysis and the magnetization-energy-ratio plot clearly show the different tendency of favoring the magnetic moment formation, and it is magnified when applied to the thin film limit where the experimental information is severely limited. As a result, LDA gives a qualitatively different prediction from GGA in the experimentally relevant region of strain whereas both approximations give reasonable results for the bulk phase. We discuss the origin of this difference and the applicability of standard methods to the correlated oxide and the metallic magnetic systems.

Here we perform a systematic and comparative investigation of SRO113 based on the standard first-principles computation methods. We elucidate the controversies in the previous calculations 16,[20][21][22][23][24] and partly resolve the unsettled issues including the thickness and strain dependent transition of magnetic and electronic properties [15][16][17][18][20][21][22][23][24] . The role of strain, correlation and lattice distortions are clarified. We further analyze the applicability of primary approximations, namely, LDA (local density approximation) and GGA (generalized gradient approximation), to metallic magnetism. Landau free energy analysis and the magnetization-energy-ratio plot clearly show the different tendency of favoring the magnetic moment formation. In case of thin film, as a result, LDA gives a qualitatively different prediction from GGA in the experimentally relevant region of strain whereas both approximations give reasonable results for the bulk phase. We discuss the origin of this difference and the applicability of standard methods to the correlated oxide and the metallic magnetic systems. Our work carries important information for future study especially when the experimental information is severely limited just as in the thin film or heterostructure.

Results
Bulk. SRO113 is a ferromagnetic metal (FM-M) with Curie temperature of 160 K and magnetic moment (M) of 1.1-1.7 μ B /f.u. 13 . Orthorhombic crystal structure is stabilized below 850 K with GdFeO 3 -type distortion 13 . To investigate SRO113, we first optimized the lattice parameters and the result is in good agreement with the previous studies (Table 1) 20,21,30,31 . A well-known trend of the underestimated/overestimated lattice parameters by LDA/GGA is clearly seen. The deviation from experimental lattice values is within ±1.2%. Rotation and tilting angles calculated by LDA and GGA are slightly overestimated (within ~1.6-2.2°) in comparison with the experimental value of 8.61° and 8.60°, respectively.
Ferromagnetic ground state is well reproduced by both LDA and GGA with the optimized and experimental lattice parameters. The calculated moments are summarized in Table 1 being consistent with the previous studies 20, 30,31 . About 70% of total moment resides at Ru-site and the other amount comes from O-site due to the hybridization 20,25 . GGA moment is larger than LDA by about factor of two. Here we note this unusually large difference of the calculated moments by two different exchange-correlation (XC) functionals. It may indicate either that GGA overestimates the ferromagnetism or that LDA underestimates it. This difference becomes critical when one tries to predict the ferromagnetic instability in the thin films for which the experimental detections are largely limited. This point will be highlighted and further discussed later in this paper.
To have further understanding, we performed the fixed spin moment (FSM) calculations which have not been reported before in these materials. Our result is summarized in Fig. 1. The aforementioned magnetic behavior, namely the enhanced (suppressed) ferromagnetism by GGA (LDA), is again clearly manifested in this plot (see blue and red lines in Fig. 1). Note that this is not just the effect from the different lattice constants given by GGA and LDA optimization. Even when the same crystal structures were used, GGA total energy more favors the ferromagnetic solution than LDA; compare the blue-solid line (filled triangles) with the yellow-dashed line (open circles), and compare the green-dotted line (open triangles) with the red-dashed-dotted lines (filled circles). It is also noted that the GGA-optimized crystal structure is more favorable to ferromagnetism than the LDA structure; compare the blue solid line (filled triangles) with the green dotted lines (open triangles), and compare the yellow-dashed line (open circles) with the red-dashed-dotted lines (filled circles).
The ferromagnetism in ruthenates has been understood traditionally based on the Stoner model 13,16,25 . Therefore, it would be fruitful to estimate the Stoner parameter (I) from our FSM calculation results. Landau free energy can be written as, where E(M) is the total energy as a function of magnetic moment (Fig. 1). The uniform spin susceptibility is given by   Table 1. The LDA result is in good agreement with the previous study 25 . Both LDA and GGA results satisfy the Stoner's criterion of ferromagnetism and the calculated magnetic moments are in reasonable range consistent with experiments. The calculated I and I N(E F ) values by GGA are larger than those of LDA by ~17% and ~40%, respectively. Namely, the stronger preference for FM solution by GGA (than LDA) can also be seen in these parameters. We once again emphasize that this cannot be attributed to the lattice effect. It is rather directly related to the parameterization of XC functional itself and thus requires further understanding at the more fundamental level. Our result demonstrates that the special care needs to be paid in ruthenates and related systems when one uses first-principles methods even within the standard XC approximations.
Undistorted thin films: 1a × 1a lateral unitcell. Now let us turn our attention to the thin film. The previous calculations and experimental results are not quite consistent with one another regarding, for example, the thickness-dependent magnetic states on SrTiO 3 substrate 13, 15-18, 20-24 . Since SRO113 thin film and its magnetic property have been studied often within 1a × 1a lateral unitcell (uc) 16,20 , we focus on 1a × 1a films before considering the effect of tilting/rotational distortions of RuO 6 octahedra. Our results are summarized in Fig. 2. The calculated magnetic moments by GGA (red circles) and LDA (blue triangles) are presented as a function of strain. Figure 2(a-c) shows the result of the 3-uc, 2-uc, and 1-uc thick SRO113 films, respectively. In the 1a × 1a lateral uc, antiferromagnetic spin order can not be simulated and the finite moment indicates the ferromagnetic order as in the previous studies 16,20 . It is noted that GGA more prefers the ferromagnetic solution than LDA, which is the same tendency as observed in the calculations of bulk SRO113.
One general feature found in Fig. 2 is that the ferromagnetism is suppressed by compressive strain. Both LDA and GGA predict the zero moment or paramagnetic phase under the large enough compressive strain. Importantly, however, the difference between GGA and LDA leads to a different prediction for the existence of ferromagnetic order and the critical value of strain below (above) which the magnetic moment or ferromagnetic order vanishes (sets in). In GGA, the ferromagnetic order is robust for 3-uc, 2-uc, and 1-uc films, which is significantly different from the LDA predictions of ≤−0.5 (3-uc)-0.0% (2-uc).
Note that, in the experimentally accessible strain range, GGA gives a qualitatively different prediction from LDA especially for the case of 1-uc. Vertical dashed lines indicate the strain values corresponding to the available substrates. Within the −2-+1.5% range, LDA gives zero (or very small) moment while GGA predicts magnetic phase and the moment is quite large. It should be noted that the zero strain in our calculations is set to be the LDA/GGA-optimized lattice parameters obtained by each XC functional and therefore two values at a given strain value represent the different in-plane lattice parameters; larger in GGA case and smaller in LDA.
As an example, let us first consider the mono-layer SRO113 grown on NdGaO 3 substrate (dark-blue vertical lines) 32 . Due to the lattice mismatch, this corresponds to the −1.81% of compressive strain situation. For this case, GGA predicts the ferromagnetic spin order while LDA solution is nonmagnetic. The same is true for the case of Sr 2 RuO 4 substrate 33 which corresponds to −1.55% compressive strain. Our result therefore raises a question about the predictive power of standard first-principles methodology. On the one hand, it shows that the special care needs to be paid in the first-principles study of correlated oxide thin films and/or heterostructures even with the standard XC energy functionals. On the other, our result calls the further method development which can better describe the electronic correlations and the structural properties. We note that in this kind of systems Another notable example is SRO113 film on SrTiO 3 substrate (−0.38% compressive strain). Contrary to the previous experimental reports [15][16][17][18] , there is no clear magnetic-to-nonmagnetic and metal-to-insulator (MI) transition in both of our GGA and LDA results. In 1-uc case, GGA result is ferromagnetic with ~1.4 μ B /Ru while a significantly smaller magnetic moment of ~0.2 μ B /Ru is found in LDA (Fig. 2(c)). It may seem inconsistent not only with experiments but also with the previous LDA 16, 20 calculations since they report the magnetic to non-magnetic transition as the SRO layer thickness is reduced although the critical thickness is still under debate [15][16][17][18]20 . It should be noted however that in the previous LDA 16 studies, the experimental SrTiO 3 lattice parameter (3.905 Å) is taken for the calculations. According to our optimized lattice parameters, this corresponds to +0.69% tensile strain for the case of LDA and −2.08% compressive strain for GGA, both of which are significantly different from the experimental situation of −0.38% compressive strain. In fact, our results of +0.69% tensile and −0.38% compressive strain are consistent with the previous calculations 16,20 . At the correct strain of SRO113 on SrTiO 3 (−0.38%), the LDA magnetic moment gradually decreases from 3-uc to 1-uc which may be understood as a magnetic to non-magnetic transition observed in experiments. In GGA, however, ferromagnetism survives down to 1-uc thickness. This behavior is double checked by other codes, namely, OpenMX and ecalj.
It is clear that understanding the SRO113 thin film based on 1a × 1a is quite limited. The tilting and rotation distortions can play important roles in determining the electronic and magnetic property.

Distorted thin films:
a a 2 × 2 lateral unitcell. The octahedral distortion (tilting and rotation of RuO 6 cage) and the use of enlarged lateral uc can lead to a qualitatively different ground state solution. The 2-uc-thick SRO113 is found to have FM-M ground state in the strain range of about −3 to +2.5%. At ~3%, G-type-like (i.e., all the nearest-neighbor couplings are antiferromagnetic) AF-I is stabilized in both LDA and GGA with Ru-site magnetic moment of 0.74 μ B /Ru (LDA) and 1.51 μ B /Ru (GGA). It is markedly different from the strained bulk SRO113 for which ferromagnetic ground state is fairly robust over a wide range of strain and G-type AF-I state is not stable 30,31 .
To understand the detailed electronic structure and its relation to the magnetism, we focus on the thinnest case (i.e., 1-uc thick) in the remaining part of this subsection. The ground state phase diagram as a function of strain is presented in Fig. 3; (a) LDA and (b) GGA. The first thing to be noted is that both XC functionals predict FM-M and AF-I ground state in the large compressive and large tensile strain region, respectively. The size of magnetic moment is generally larger in GGA than LDA, which is the same feature found in the previous cases of bulk and undistorted thin film.
In the experimentally more accessible range, however, LDA and GGA give qualitatively different solutions. In the GGA phase diagram (Fig. 3(b)), AF-I is the ground state around zero strain region (≤±1% strain) while LDA predicts either FM-M or antiferromagnetic metal (AF-M). The AF-M phase of LDA (which is absent in GGA phase diagram) is attributed to the smaller energy splitting between d xy and d yz,zx states near E F . Thus, for the case of SrTiO 3 32 , DyScO 3 32 or GdScO 3 14 substrate, we have two different predictions regarding the ground state property from two standard XC functionals. It therefore raises a serious question about the predictive power of the current first-principles method to describe the correlated oxide thin film and/or heterostructure for which the experimental information is often limited. Further investigation and development are urgently necessary.
Noticeable is our GGA result of AF-I ground state for SrTiO 3 substrate (see Fig. 3(b)). We emphasize that this AF-I phase has never been achieved before by LDA and GGA 16,[20][21][22] . Since AF-I has previously been obtained only by DFT + U 21,22 or DFT + DMFT (DFT + dynamical mean-field theory) 24 and it can be consistent with the recent experiments reporting an insulating state carrying no net moment 17,18 , it has been discussed that the insulating gap of 1-uc SRO113 is opened by on-site Coulomb correlation 24 . However, Fig. 4(a) clearly shows the band-gap  Table 1). The vertical dashed lines represent the strains realizable by NdGaO 3 (blue), Sr 2 RuO 4 (light blue), SrTiO 3 (orange), DyScO 3 (magenta), and GdScO 3 (red) as substrates.
without the explicit inclusion of U while the gap size is smaller than that of DMFT (E g DMFT = 1.0 eV) 24 . Our result therefore implies that the insulating ground state of mono-layer SRO113 is not just attributed to the local Coulomb physics within Ru-d electrons, but the homogeneous electron approximation can describe the gap opening. It also shows that the careful consideration of lattice effect is important to simulate the experimental situation of complex oxides. A common trend found in both LDA and GGA phase diagram is that the tensile strain makes system antiferromagnetic while the compressive strain favors ferromagnetic. It can be simply understood by defining an energy level difference between Ru-d yz,zx and d xy state:  Fig. 4(b). In the large Δ limit (large tensile strain regime), the separation between ε yz,zx and ε xy is large and it approximately corresponds to the (two-band) half-filled case within the ionic picture of Ru 4+ (see inset of Fig. 4(a)) as in the case of Ca 2 RuO 4 36 . Thus the antiferromagnetic spin order is naturally stabilized via superexchange. On the other hand, the small Δ limit (compressive strain) basically corresponds to the bulk SRO113 regime in which ferromagnetic order is stable.  The inset of Fig. 4(b) shows the calculated Δ GGA as a function of SRO thickness. A decreasing trend is clearly noticed as the number of layers increases. Since this is a qualitatively consistent trend with the magnetic transition observed in the several experiments as a function of layer thickness, our result implies that the Δ plays an important role in the thickness-dependent transition. Therefore it is probably not understood solely from N(E F ) and Stoner's criterion 16 , but rather comparable with the case of SrVO 3 /SrTiO 3 in which the lifted orbital degeneracy plays an important role to induce the MI transition 37,38 . We do not find any systematic trend in N(E F ) over the range of strain. In the sense that the FM-M to AF-I transition occurs in between 2-and 1-uc thickness, our GGA result is not in quite good agreement with experiments which report the transitions in between 5-and 4-uc 15 , or 4-and 3-uc 17, 18 , or 3-to 2-uc 16 . Further inclusion of electron correlations might be needed to correctly describe the Δ and other electronic properties 19,[38][39][40][41][42][43] .
It might be informative to see DFT + U 44 results. In spite of its static Hartree-Fock nature, it has been used to study SRO113 21,22 . With a fixed U value to a recent cRPA (constrained random phase approximation) result (U = 3.5 eV 24 ) and at a −0.5% compressive strain (close to SrTiO 3 value), AF-I phase is found to be the ground state.

Discussion
Our results show that the prediction of the magnetic ground state of complex oxide thin film or heterostructure is a challenge to the current first-principles methodology. Without further information such as structural property and/or magnetic order (usually known from experiment in the case of bulk), it is difficult to make a prediction for thin film SRO113. Note that LDA and GGA are not a bad choice to describe the bulk SRO113 in the sense that both predict the correct FM-M ground state. Although some detailed features in the electronic structure are not well captured by LDA or GGA (e.g., Hubbard-like state ~−1.3 eV below the Fermi energy; see, ref. 45), these standard approximations provide reasonable descriptions of this metallic bulk material. Arguably, they are better than other improved but static methods to take orbital-dependent correlations into serious account such as DFT + U 11 . It is therefore surprising to see such a large difference between LDA and GGA in this system.
We first note that this difficulty is related to the complex nature of correlated oxides in which the spin, orbital, and lattice degree of freedom are tightly coupled to each other. On the one hand, it is the reason why the parameter-free first-principles techniques are strongly requested. On the other hand, this complexity makes the first-principles description be a greater challenge. Note that an improved description of electronic correlations cannot immediately be the right answer unless it can simultaneously give the better description of the other degrees of freedom such as lattice. For example, if we resort to DMFT approximation to get a better electronic structure of ruthenates, we have to describe the force and the lattice relaxation simultaneously at the same level. This development is technically non-trivial [46][47][48] . Similar issues and related intrinsic ambiguities can also be found in other techniques such as SIC-DFT (self-interaction corrected DFT) and hybrid functional.
Here it is instructive to recall that the metallic magnetism has been a challenge to the first-principles calculation [49][50][51][52] . A classical example is bcc Fe for which minimum energy solution of LDA is paramagnetic hcp phase (with the optimized lattice parameter) 49 unlike Co and Ni 53 . In other words, without a priori knowledge of ferromagnetism, the situation of bcc Fe is similar with SRO113 thin film; paramagnetic ground state by LDA and ferromagnetic in GGA. Since the use of experimental lattice parameter gives the correct ferromagnetic solution in both LDA and GGA, bcc Fe is an easier case in the practical sense. In SRO113 thin film, on the other hand, the limitation of computation method is magnified due to no a priori information of structure and magnetic order.
A more detailed nature of enhanced ferromagnetism by GGA functional can be seen in Fig. 5 where the PBE (Perdew-Burke-Zunger) correction of spin-polarized energy density with respect to LDA is presented as a Figure 5. The calculated α XC for c-SRO113 (black crosses) and bcc Fe (red crosses). About 85000 real-space grid points are obtained from the self-consistently converged electronic structure within LDA. Experimental lattice parameters are used. ζ is fixed to 0.02 which is the average value of c-SRO113. The α XC change is negligible in the range of 0 < ζ < 0.5. function of inverse density (r s ) and density gradient (s) for the relative spin-polarization (ζ refer to the spin-polarized part of XC energy density within PBE and LDA, respectively. For the definitions of ε, r s , ζ, s, F XC (r s , ζ, s), and f(ζ), we follow the original PBE 54 and LDA-PZ (Perdew-Zunger) 55 papers. By definition, α XC represents the energy gain by GGA-PBE over LDA. In the yellow-red colored region of Fig. 5, GGA-PBE prefers the magnetic solution more than LDA while in the dark blue region LDA-like solution is favored. About 85000 real-space grid points (using OpenMX code) are plotted in this map for bcc Fe (red-colored points) and c-SRO113 (black-colored points), respectively. It is clear that both bcc Fe and c-SRO113 have large portion of grid points in the region of α XC larger than unity. It is responsible for the larger spin splitting (~ε ζ ∂∆ ∂ / XC ) and thus the moment formation enhanced by PBE correction. Since many interesting aspects of complex oxide films and interfaces are directly or indirectly related to the metallic magnetism [56][57][58][59][60] , it is strongly requested to develop more reliable computation methods. One promising aspect is the same overall trend given by LDA and GGA. For example, the moment is suppressed by compressive strain in the case of 1a × 1a lateral uc ( Fig. 2(a,b)). FM-M and AF-I phase in the large compressive and large tensile strain, respectively, are also consistently reproduced by LDA and GGA (Fig. 3). Further investigation of parameterization hopefully gives a way to improve approximations.

Summary
From a systematic first-principles investigation, we elucidate the controversies in the previous calculations and experiments. A part of issues has been clarified including the thickness-and strain-dependent phase transition as well as the ground state property of monolayer SRO113. The role of strain, correlation and lattice distortions are explored in detail. Applicability of LDA and GGA has been discussed in the realistic material context by using Landau free energy analysis and magnetization-energy-ratio plot. While the overall feature is commonly reproduced by both approximations, they give qualitatively different predictions for the case of thin film in the experimentally relevant region of strain. Our work provides useful information not only for ruthenates and other correlated oxides, but also for further first-principles studies in the related systems such as metallic magnetism.

Methods
We choose LDA parameterized by Perdew and Zunger 55 and GGA by Perdew, Burke, and Ernzerhof 54 as two representative XC functionals. The data obtained by these two functionals is presented as our main results and discussed in a comparative way. We used projector-augmented-wave (PAW) method 61 as implemented in the VASP software package 62 . While VASP-PAW is used to obtain our main results, we also double checked some cases with other methods, confirming that our conclusions are valid. The other codes we used include OpenMX 63 , which is based on the norm-conserving pseudopotential and pseudoatomic orbitals, and ecalj 64 based on the full-potential 'LMTO' or 'PMT' basis 65 .
Tetrahedron method with Blöchl correction has been used for Brillouin zone integration 66 , which we found is important to achieve the enough accuracy. In some cases the use of Gaussian broadening is found to give qualitatively different predictions depending on the broadening value. For bulk structure optimizations, 12 × 12 × 12 and 8 × 8 × 6 k-points were used for 5-atom uc of c-SRO113 and 20-atom uc of o-SRO113, respectively. The force criterion of 1 meV/Å was adopted for structural optimization. This choice of k meshes and the 500 eV energy cutoff were found to be satisfactory from our convergence test. For the electronic structure calculations, however, we took a denser k-grid of 16 × 16 × 16 for c-SRO113 and 12 × 12 × 10 for o-SRO113 to calculate DOS for the given optimized structures. To simulate SRO113 thin films, we considered the 3-uc, 2-uc, and 1-uc slab geometries terminated with SrO layers which is known to be the case in experimental situations 13 . The k-grids used for slab calculations were 12 × 12 × 2 for c-SRO113 and 8 × 8 × 2 for o-SRO113. The vacuum thickness is ≥15 Å which is large enough to simulate the experimental situation. The Ru radius is set to 1.402 Å and the Ru moment dependence on this setting is negligible (~0.01 μ B ).
Epitaxial strain is one of the key control parameters in the thin film growth. It is therefore of crucial importance to determine how to simulate the in-plane lattice parameter corresponding to the experimental situation. Note that this is a non-trivial issue because LDA-and GGA-optimized lattice parameters are both in general different from the experimental values. While the former tends to underestimate the lattice parameter, the latter overestimates (see Table 1). Thus, the naive choice of a XC functional or taking experimental lattice value can easily be misleading. Previous studies are not quite consistent in this regard. In the current study, we take the optimized lattice constant with a given XC functional as the reference value (i.e., zero strain) and define the compressive/tensile strain value with respect to it. This can be the most reasonable way to simulate the experiments and can serve as a solid reference for the future study. This point should be kept in mind when our results are compared with the previous calculation results 16,22 .