Ferromagnetism of LaCoO$_3$ films

We study ferromagnetic ordering and microscopic inhomogeneity in tensile strained LaCoO$_3$ using numerical simulations. We argue that both phenomena originate from effective superexchange interactions between atoms in the high-spin (HS) state mediated by the intermediate-spin excitations. We derive a model of the HS excitation as a bare atomic state dressed by electron and electron-hole fluctuations on the neighbor atoms. We construct a series of approximations to account for electron correlation effects responsible for HS fluctuations and magnetic exchange. The obtained amplitudes and directional dependence of magnetic couplings between the"dressed"HS states show a qualitative agreement with experimental observations and provide a new physical picture of LaCoO$_3$ films.


Introduction
LaCoO 3 is an non-magnetic insulator, which develops a Curie-Weiss magnetic response and eventually transforms to metal with increasing temperature [1][2][3][4]. This peculiar behavior arising from proximity of several atomic multiplets has attracted attention for decades and the debate is not over yet. Strained films of LaCoO 3 were shown to develop a ferromagnetic (FM) order [5][6][7][8][9] while remaining insulating [10]. Small value of saturated magnetization suggests that the FM state is not uniform on the atomic scale. Fujioka et al. [11] observed a commensurate structure with the propagation vector (1/4, -1/4, 1/4). A superstructure with magnetic atoms separated by three non-magnetic ones was observed with transmission electron microscopy [12].
The discussion of bulk LaCoO 3 revolves around the nature of the lowest atomic excited state with the alternatives being the high-spin (HS) or the intermediate-spin (IS) states [13][14][15][16][17]. While spectroscopic data [18,19] as well as crystal-field considerations [20] clearly favor the HS state, recent resonant inelastic x-ray scattering (RIXS) experiments [21,22] show that IS states cannot be neglected despite their single-ion energy of several 100 meV above the low-spin (LS) ground state. A peculiar feature of LaCoO 3 is the short-range ferromagnetic correlation in the middle-temperature range [23,24], which becomes more pronounced in thin films leading to a ferromagnetic phase below 94 K [11]. Should the magnetism of LaCoO 3 originate from the HS state of Co 3+ ions only, an anti-ferromagnetic (AFM) coupling between these is expected [25].
There were several attempts to explain emerging FM long-range order in insulating LaCoO 3 under tensile strain [24,[26][27][28][29]. The LDA+U density functional calculations using Hartree-Fock-like treatment of intra-atomic d-d interaction [27,28] argue in favor of the superexchangebased mechanism for FM ordering. However, a more detailed analysis within the same framework allowing for AFM arrangement of HS shows that the LDA+U treatment can not be considered as conclusive due to diverse results, which strongly depend on the employed interaction parameters (see Ref. [28] and Appendix A for more details).
In this paper, we build a model of LaCoO 3 in the strained tetragonal structure with the aim to explain its ferromagnetism and atomic-scale inhomogeneity. Despite being an insulator with quenched charge fluctuations, the number of local low-energy degrees of freedom involving LS, IS, and HS states is prohibitively large for a direct numerical treatment. Therefore, we construct a series of approximations. Starting from electronic description, we integrate out the charge fluctuations to arrive at an effective model in the Hilbert space spanned by the LS, IS, and HS states [21,30]. We take advantage of the tetragonal crystal field, which lifts the orbital degeneracy of the HS and IS states and thus simplifies the description of the strained LaCoO 3 films in comparison to the bulk material. This model describes a gas of IS excitonic particles with strongly-directional mobility depending on their orbital character. Strong local attraction between IS excitons of different orbital character renders their immobile bound pairs (HS states) to be the lowest excitations of the system. In the next step, we integrate out the IS excitons and formulate a model in terms of HS particles. Dressed in that way HS bi-excitons possess strongly-directional inter-site interaction varying between repulsion along the crystallographic axes and in-plane next-nearest-neighbor attraction in the diagonal direction. As a result, spatially separated chains of HS sites provide an energetically stable building block. We show that dynamics of such chain can be described by the Ising-like model with FM coupling, which provides the lowest level of simplification we arrive at. We further study the arrangement of such chains in a three-dimensional (3D) system.

Theoretical approach
Our approach consists of three successive steps (see also Fig. 1): (A) Density-functional theory (DFT) calculation is performed for the strained LaCoO 3 structure [31]. A tight-binding (d-only) model in Wannier basis spanning the Co d-like bands is constructed [32,33] and augmented with the local Coulomb interaction within the Co 3d shell, which amounts to projecting out the O-2p orbitals.
(B) The 5-orbital Hubbard model is reduced to an effective low-energy model spanning the LS, HS and IS states by employing the Schrieffer-Wolff transformation [34]. This amounts to integrating out the local charge fluctuations from d 6 to d 5 and d 7 configurations. Furthermore, we use the tetragonal symmetry to restrict our model to the lowest HS-orbital singlet and two orbital flavors of the IS multiplet that are coupled. The model at this level describes a gas of hard-core IS excitons, which propagate on the lattice. The IS excitons with different orbital flavors interact via strong local attractive interaction, which gives rise to local bound pairs, bi-excitons.
(C) Analysing the behavior of the model (B), we come to the conclusion that for realistic strength of the attractive IS-IS interaction, the local bi-excitons are stable objects consisting of atomic HS state dressed with IS-IS fluctuations to the nearest-neighbor (nn) sites. The last layer of approximation consists in integrating out the IS states. This way we arrive at a model formulated in terms of HS particles only. These turn out to be immobile in the xy plane and interacting via strongly directional interactions.
In the following, we provide technical details about execution of the above program and point out the special features of the models (A)-(C), which justify our treatment. The numerical results and their discussion is left for the next section.
Step (A) follows a standard approach of constructing the multi-orbital Hubbard model used for correlated materials [35]. DFT computational parameters are summarized in Appendix A. To incorporate the correlation effects originating from electron interactions in the Co 3+ d 6 shell of LaCoO 3 , we introduce the 5-orbital Hubbard Hamiltonian The local part H (i) at consists of the tetragonal crystal field (diagonal in cubic harmonics basis), spin-orbit coupling (SOC), and the Coulomb interaction, while H t describes the inter-atomic hopping. The coupling constants are given in Appendix B, c † iκ (c iκ ) are the fermionic creation (annihilation) operators, i labels the cobalt ion sites, while κ, λ, µ, ν are the combined orbital and spin state indices.
In following, we proceed without spin-orbit interaction H SOC , which we treat perturbatively as described later. This way we can use higher symmetry of the Hamilatonian and apply corrections at the end instead of neglecting small terms at various points of steps (B) and (C). We diagonalize the local Hamiltonian H (i) at and define the low-energy subspace L to be the space spanned by 25 lowest d 6 states: the spin-singlet S = 0 LS state, the spin-triplet S = 1 IS states with the t 1g orbital symmetry (ISx ≡ (y 2 − z 2 ) ⊗ yz, ISŷ ≡ (z 2 − x 2 ) ⊗ zx, ISẑ ≡ (x 2 − y 2 ) ⊗ xy), and the spin S = 2 HS states with the t 2g orbital symmetry (HSx ≡ yz, HSŷ ≡ zx, HSẑ ≡ xy). The strong-coupling expansion consists in the Schrieffer-Wolff transformation to the second order in the hopping H t [36] (we consider the nn processes only). The resulting effective Hamiltonian can be written in the form where the operator T αβ i describes a transition between local states |β and |α on the site i, while e labels the nn bonds.
The number of local degrees of freedom in the above model is still too high to allow for a numerical treatment of even a small cluster. We notice several symmetries that are weakly broken by SOC: (i) Viewing the HS state as a pair of IS excitons, the number of IS excitons of a given orbital character is conserved. This is a simple consequence of the fact that a t 2ghole cannot change its character (nn hopping is only possible between the same t 2g orbitals). This will allow us to work in subspaces indexed by the total number of IS excitons. (ii) The hopping amplitude of ISα along the α-axis is very small. We will assume that this hopping vanishes completely and thus ISα are confined to planes perpendicular to the α-axis. (iii) HS state can be viewed as a bi-exciton, e.g, HSẑ ≡ ISx ⊗ ISŷ. Non-orthogonality of the e g orbitals forming the two IS excitons is accounted for by a correlated hopping terms in the effective Hamiltonian. (iv) The tetragonal crystal field splits HSẑ from the set of HS states to be the lowest state. In the following, we will consider the subspace containing only HSẑ and the states coupled to it, i.e., ISx and ISŷ (see also Fig. 2). This choice is justified a posteriori by noting that for realistic parameters, the lowest excited states are HS-like. (v) The Hamiltoninan (2) without SOC has the SU(2) spin symmetry, which will allow us to further restrict the Hilbert space in numerical calculations.
The coupling constants of the model (2) consist of nn coupling and on-site energies. The nn coupling is determined by the products of hopping elements in H t divided by denominator of the order of F 0 . Changing the parameters of local interaction leads to more or less uniform scaling of all nn couplings, which does not change the studied physics significantly. Therefore, we treat the nn couplings as fixed, see Appendix B for details. The on-site energies of the IS and HS excitations, which are obtained by the strong-coupling expansion and measured relative to the LS background, are less reliable, since they are calculated from differences of large numbers. In particular, the lowest excitation energy is experimentally known to be rather small (< 20 meV). For this reason we fix the on-site energy of the ISx ,ŷ exciton to E IS = 330 meV, which matches well the RIXS data for the bulk LaCoO 3 , and treat the onsite IS-IS attraction U , E HS = 2E IS − U , as adjustable parameter. Since all our calculations are performed for a fixed number of IS excitons, the results do not depend on E IS . This way we can analyze relative stability of various states, while the absolute excitation energy above the singlet ground state is not studied.
The step between models (B) and (C) is executed by performing series of exact diagonl-ization calculations for small clusters. The computational effort is heavily reduced by the confinement of IS excitons into their respective planes (or even lines in calculations for a single xy plane). We start by investigating the structure of a single HS excitation. To this end, we perform calculations in the subspace with one of each ISx and ISŷ exciton. We find that for realistic values of U , the ground state in the 2D case is a HS state located on the intersection of the "life-lines" of the two IS excitons, see Fig. 2(b) and (c), with some admixture of IS-IS configurations. In the 3D case, the "life-lines" turn into "life-planes" and the HS state with the cloud of IS-IS fluctuations can propagate along a vertical line parallel to the z-axis in Fig. 2(b). Next, we study the effect of SOC in the single HS set-up. We find that it can be approximated by replacing the spin S = 2 and the isotropic exchange S i · S j with a pseudospin J = 1 anisotropic J-J exchange with stronger in-plane component.
Finally, we study clusters with pairs of HS bi-excitons, i.e., the subspaces containing two of each ISx and ISŷ excitons. We determine that the lowest states are linear combinations of the dressed HS bi-excitons and use the respective eigenenergies to determine parameters of the model (C) formulated in terms of HS states only. In addition to the bare nearest-neighbor HS-HS repulsion present in the model (B), we obtain longer-range interactions, which arise from the interaction of the IS-IS clouds surrounding the HS cores.

HS-IS model
Upon the Schrieffer-Wolff transformation of the model (1) and restricting to the lowest-energy sector containing HSẑ, the Hamiltonian (2) acquires the form where L † i , b † i,Λs and h † i,s are creation operators of Schwinger bosons representing LS, IS Λ=x,ŷ and HSẑ states, respectively, with the spin projection s and orbital flavor Λ on the lattice site i. The physical states fulfil the hard-core constraint  Table 1: Coupling constants of the effective model (3). Following symmetries are obeyed: J The coefficients are normalization factors times the Clebsch-Gordan coefficients, which ensure the SO(3) spin symmetry. The summation bounds of the spin index are s = −1, 0, 1 and s = −2, . . . , 2 for IS and HS states, respectively. It is also assumed that any expression with s out of the given bounds has zero prefactor. The index λ = x, y, z labels directions on the cubic lattice. In the following, we use the number operators n i,a = s a † i,s a i,s and n a = i n i,a . The coupling constants are given in the Table 1.
The model (3) has several specific features mentioned above: i) The total numbers n bx +n h and n bŷ + n h are conserved. ii) The HS h-particles do not hop directly. iii) The IS bx-and bŷ-particles can hop only within the yz and xz planes, respectively.
While the Hamiltonian (3) looks complicated, its basis structure is actually simple. It describes two types of excitonic particles b † i,Λs (Λ =x,ŷ), which can move on the lattice in mutually orthogonal planes, interact via strong attractive ferromagnetic on-site interaction U and spin-dependent nn interaction. Since the HS state is not exactly a product of two IS states, the hopping and nn interaction amplitudes to some extent vary depending on the initial state. Therefore, we have t for the |LS, IS ↔ |IS, LS process, t c for |HS, LS ↔ |IS, IS process andt for |HS, IS ↔ |IS, HS process. Similarly, the nn terms with J 1 , J 2 , and J 3 correspond to IS-IS, HS-IS, and HS-HS interactions, respectively. The superscript index of J (k) refers to spin-multipole order of the interaction with k = 0 corresponding to the n a n b density-density interaction, k = 1 to S · S dipole-dipole interaction. Symmetry allows also quadrupole k = 2 interactions, which cannot appear in the second-order processes and thus are absent in our treatment.
Finally, we briefly discuss the physics underlying the coupling constants in Table 1. The planar nature of the exciton hopping refers to the planar nature of both the hole and electron part of the exciton. For example, the ISx exciton consists of y 2 − z 2 and yz electron and hole part, respectively, both of which have small hopping amplitude along the x axis (see also Fig. 2). Positive sign of the hopping amplitude corresponds to the excitonic-band maximum at Γ point, in agreement with the experiment [21]. To understand the origin of nn interactions, it is necessary to note that on the LS background the IS and HS (bi-)excitons reduce their energy by virtual hopping of electrons and holes to the nn sites. For the HS bi-excitons these processes happen in all directions, while for IS excitons in the respective "life-planes" only. Placing (bi-)excitons of the same orbital character on the nn sites blocks some of these processes, thus causes an effective repulsion, e.g., J 3x . Should the spins of these nn (bi-)excitons be anti-parallel, the Pauli blocking is less effective, which adds an antiferromagnetic dipole-dipole component J 2xŷ , and J (1) 3x . The interaction on the bonds, along which one or both of the excitons have vanishing hopping amplitudes, is weak, J (k) 1xxx being example of the latter.

Structure of a single HS state
Quantum fluctuations. We begin with numerical analysis of the structure of a single HSẑ excitation. To this end, we investigate the subspace of the total spin S = 2 containing one of each ISx and ISŷ excitons. Since the IS excitons are constrained to mutually perpendicular planes, the 3D problem is reduced to two intersecting 2D planes, each containing one IS excitons. Similarly, the 2D problem of a single xy plane, depicted in Fig. 2(b), reduces to two intersecting 1D lines. We use exact diagonalization for finite systems with linear dimension L = 8 and periodic boundary conditions. The physics is governed by the ratio of the attractive interaction to hopping U/t. In Figure 3, we show the energies as functions of U/t (equivalently, E HS ), as well as the density distributions in the ground state of the 2D and 3D model. For weak U , the ground state of the 3D model consists of delocalized (plane-wave) IS excitons living in their respective life-planes. Between U/t of 4 and 6, a bound HSẑ bi-exciton is formed that can propagate along the intersection of the two planes forming thus a narrow band.
What is a realistic magnitude of the on-site IS-IS interaction U in LaCoO 3 and how accurately should it be determined? While the excitation gap of LaCoO 3 , because of its smallness, is very sensitive to the size of U , the wave functions of the excitations and their mutual interactions are much less sensitive. Setting U to approximate the experimental energỹ E HS of HS excitation thus provides a solid starting point for analysis of interactions between excitons. The corresponding value ofẼ HS in bulk LaCoO 3 is about 10-20 meV [18,19]. The gap for elementary excitations in the strained LaCoO 3 films is smaller or even negative 1 . In both cases, |Ẽ HS | E IS , with the experimental value of E IS between 300-350 meV [21]. We can therefore safely assume thatẼ HS ≈ 0 in our considerations. With some provision for the uncertainty of E IS ,Ẽ HS calculated at E IS = 330 meV puts the realistic U/t to the 8-10 interval, see Fig. 3.
In the determined U/t range, the ground state of the 2D and 3D problems is a HS biexciton with |HS, LS ↔ |IS, IS quantum fluctuations confined to the nn bonds. The ground state energyẼ HS = 2E IS − U − ∆ Q is dominated by the bare HS contribution 2E IS − U with a correction ∆ Q of about 40 meV and 80 meV due to quantum fluctuations in the 2D and 3D cases, respectively. The difference comes from effectively doubling the number of accessible neighbors by going from 2D to 3D. While in the 2D case the HS bi-exciton is confined to the intersection of 1D life-lines of the IS excitons, in the 3D case it can hop along the line of intersection of the corresponding life-planes with a moderate amplitude of +10 meV. Inclusion of SOC lowersẼ HS by approximately 20 meV, an amount that does not affect the presented estimate of U .    Figure 4: Splitting of atomic energies by SOC, where the maximal strength corresponds to ζ = 56 meV. In this limit, the structure of 5 lowest states is shown by plotting weights of the states (a)-(e). The labels and numbers on the top of color bars correspond to the orbital character and spin projection s z in the cubic-harmonics basis (ζ = 0), respectively.
The main result so far is the observation that in an insulator the HS bi-exciton is stable and can be viewed as an atomic HS state dressed with IS-IS fluctuations confined to the nn sites. With this picture in mind we will discuss the effect of SOC in the following.
Spin-orbit coupling. The SOC introduces an on-site mixing between LS and IS, IS and HS states as well as mixing within the HS and IS multiplets. The former causes violation of the conservation law 2n HS + n IS = const. Thanks to the large separation between E IS on one hand and E HS ≈ E LS on the other, the violation is weak and will be neglected in the following. The latter makes the hopping of IS excitons and thus HS bi-excitons weakly three-dimensional, however, the effect is much weaker compared to the inter-site interaction, derived in the next section, and is also neglected.
Thanks to the sizeable tetragonal crystal field, the triplet states retain a dominant HSẑ character. The leading effect of SOC is splitting of the HSẑ quintuplet into a ground state quasi-triplet and an excited doublet with energy separation of about 20 meV, see Fig. 4(a). Projection on the quasi-triplet state, while neglecting the admixture of HSx and HSŷ, approximately amounts to elimination of the S z = ±2 states of the HSẑ quintuplet. We introduce a pseudospin J = 1 to describe the ground state triplet, which leads to a substitution S z → J z , S x,y → √ 3J x,y . The rotationally invariant terms S i · S j obtained without SOC thus undergo a simple transformation As a result, we can proceed with analysis of the model without SOC and include the leading correction due to SOC by introducing the above substitution to the final results.  Table 2: Values for the HS-HS interaction amplitudes in meV units in the effective model (5) at two different amplitudes of U/t.
Evaluating L z and L x,y matrices in the basis of the eigenstates of atomic Hamiltonian with SOC (see Appendix C for details), we estimate the the effective g factors for the lowest J = 1 triplet from the the operator relation M = µ B (2S + L) = µ B gJ. The obtained g values are almost isotropic, g zz ≈ 2.73 and g xx,yy ≈ 2.64.

HS-only model
Having established the stability of HS bi-exciton with respect to HS ↔ IS-IS dissociation, we eliminate IS excitons and formulate a model H C in terms of renormalized HS bi-excitons. These states can be viewed as a local HS core surrounded by IS fluctuations on nn sites. While multi-particle interactions are possible due to non-local structure of the renormalized HS states, we restrict ourselves to pairwise interactions and assume interactions involving simultaneously three or more particles to be less important. To estimate the parameters of H C , we diagonalalize H B in the subspace with two HSẑ bi-excitons in various geometries and analyze the spectrum and characteristics of the lowest eigenstates.
Interaction in the xy plane. To analyze the in-plane HS-HS interactions, we choose 2D geometry, in which the two HSẑ are confined to the xy plane. Neglected IS fluctuations in the z-direction, which do not contribute to xy interactions, amount to about 5% of the HSẑ wavefunction. Hence, the interactions may be overestimated by 5% in 2D geometry with respect to the 3D one, a correction that is well beyond the expected accuracy of our treatment.
The confinement of ISx and ISŷ excitations to lines parallel to y and x axis, respectively, simplifies the calculations significantly: it allows to work in geometries corresponding to the pairs of lines parallel to the axes (with a somewhat different treatment of the situation, where the two HS states appear on the same line), see Fig. 2(e)-(g). The interaction strength for a given HS separation is obtained simply as the difference between the ground-state energy for a given two-HS geometry and twice the energy of a single HS state on a lattice of the same size.
First, we analyze the interactions in the fully FM polarized configuration, S z = 4, in geometry with periodic boundary conditions and the linear size L = 8. The results are summarized in Table 2. The nn interaction is strongly repulsive. It is dominated by the explicit nn HS-HS interaction in H B , see Eq. (3) and Table 1. The AFM spin-dependent part does not overcome the spin-independent repulsion, thus the nn interaction remains strongly repulsive irrespective of spin configuration.
The second-neighbor (nnn) interaction along each of crystallographic axes is also repulsive. The origin of this repulsion is the hard-core constraint on the IS states of the same orbital character in H B , which restricts the virtual IS fluctuations. Since the constraint applies irrespective of spin projection, the interaction is spin-independent. The interactions at longer distances are repulsive for the same reason, but their amplitude is negligible at U/t ≥ 8 due to localization of the IS cloud to nn sites. The interaction between HS bi-excitons that are not on the same line is negligible for the same reason, except for the nnn interaction, which requires a special treatment.
The ground state of the configuration with the two corners of a square plaquette, see Fig. 2(e), is a quasi-doublet, which originates from the HS bi-excitons exchanging the corners of the plaquette. Comparing the energies of the doublet with the energy of single HS bi-exciton as shown in Fig. 5(a), we obtain an attractive diagonal interaction as well as the amplitude of the pair-hopping t p . Both the attractive nnn interaction and the pair-hopping have their origin in the on-site attraction between ISx and ISŷ, which can meet in the "empty" corner of the plaquette. Since the on-site attraction is strongly FM, we expect the diagonal interaction as well as the pair-hopping to be strongly spin-dependent.
To investigate the spin structure of the diagonal nnn interactions, we perform calculations in the subspaces with total spin projections S z =4, 3 and 2. The calculations are performed on a 2 × 2 plaquette with open boundary conditions. The results summarized in Fig. 5 reveal that: i) The spin-independent part of the nnn interaction (and pair-hopping) is small. ii) The spin exchange interaction is FM and well approximated by S i · S j . iii) The pair hopping term has more complicated dependence on S i · S j and is large only for parallel spins. We point out that the mechanism of the attractive FM nnn interaction is similar to the classical 90 • Goodenough-Kanamori superexchange, if one replaces the virtual excursions of electrons to the common anion site with virtual excursion of IS excitons to the common LS site.
According to the listed observations, we introduce an extended Heisenberg model in terms of the dressed HS states only (on the top of the LS vacuum), where τ r (τ ± ) are the Pauli matrices (ladder operators) representing pseudospin-1/2 operators in the (HS,LS) basis, S i are the conventional spin S = 2 operators of the HS state, ij indicates a summation over the nnn sites along two diagonals. While the indices i and j correspond to the sites on one diagonal in a plaquette, the sites i and j are their mirror images placed on the perpendicular one. The first term corresponds a single HS bi-exciton with the excitation energy µ, the second one describes repulsion of HS bi-excitons along the crystallographic axes with characteristic values given in Table 2. The third term is the diagonal nnn interaction and the fourth is the pair-hopping on the square plaquette.
Interaction out of the xy plane. Extraction of the out-of-plane interaction parameters is cumbersome due to the mobility of HS bi-exciton in the z-direction. The repulsion along z-axis has similar characteristics and origin as the repulsion along the x/y-axis, see Table 2.
There is a strong nn repulsion, originating from the bare HS-HS interaction in H B , and a sizable nnn repulsion, while the interaction between more distant neighbors is negligible. The mechanism of attractive diagonal interaction J ex present in xy plane is not active on xz and yz plaquettes. The reason is that the HS-IS repulsion relevant in this case J (k) 2xŷ is much stronger than the repulsion J (k) 2xx relevant for the xy plaquette. Simplified calculations on the 2 × 2 xz (yz) plaquette suggest that the diagonal interaction is an order of magnitude weaker (∼ 0.1 meV). To obtain a reliable spin structure of this weak interaction, explicit inclusion of SOC may be necessary and is not attempted here.

HS order on the lattice
Next, we discuss ordering of the HS bi-excitons on the lattice. So far we have used the fact that the number of excitons per orbital flavor is approximately conserved, thus performed calculations for specific exciton-charge sectors without need to specify the energy of HS biexciton with respect to the LS vacuum, µ in Eq. (5). However, the energy µ is crucial for determination of the global ground state. Experimental estimates of µ in bulk LaCoO 3 yield 10-20 meV [18,37]. To calculate µ with the desired meV accuracy from first principles is unrealistic (unlike the other coupling constants that are effectively ratios of small and large parameters).
Therefore, we base our discussion on the experimentally observed evolution of LaCoO 3 with increasing lateral strain from non-magnetic bulk to ferromagnetic LaCoO 3 /SrTiO 3 films. The strain-induced tetragonal crystal field causes splitting of the HS multiplet into the lowenergy HSẑ singlet and the high-energy HSx ,ŷ doublet, while the energy of the HSẑ excitation above the LS state decreases. As a result, the excitation energy µ of the dressed HSẑ bi-exciton also decreases from its bulk value. Considering the ground state of the model (5) as a function of decreasing µ, we start from an empty (LS) lattice. Upon reaching a critical µ c , the nnn interactions will lead to formation of diagonal ferromagnetic chains separated by at least two empty sites, see Fig. 6, due to the strong nn and nnn repulsive interactions. Neglecting the pair-hopping terms and replacing the Heisenberg exchange with the Ising one, the model (5) reduces to a generalization of the classical Blume-Emery-Griffiths model [38]. While the pair-hopping plays a role in determining µ c and chain separation in a single xy-layer problem, we conjecture that the pair hopping is strongly inhibited in the 3D problem due to strong nn repulsion along the z-axis. Neglecting the pair-hopping, we obtain µ c = 4J ex for a single layer problem. By reducing µ below µ c , the FM chains start to form with a separation dictated by the repulsive interactions.
To proceed with the 3D problem, we assume identical order in all layers, i.e., we only have to determine the type of stacking. Considering nn and nnn repulsion along x, y and z axes yields the 1/3 HS-filling of the lattice with the stacking shown in Fig. 6(b). Note that this order breaks the tetragonal symmetry. Another possible stacking obtained at filling 1/4 is shown in Fig. 6(c). In this case, the chains in the adjacent planes are mutually perpendicular and the tetragonal symmetry is preserved. The (1/4, 1/4, 1/4) periodicity is consistent with the finding of Ref. [11]. We point out that the orders at 1/3 and 1/4 filling are not distinguished by the out-of-plane nnn interactions. The larger spacing between the chains along the z-axis at 1/4 filling may be a consequence of the quantum-mechanical effects that we did not consider in this subsection, i.e., the in-plane pair-hopping or hopping along the z-axis estimated above.
It is clear that lowering of µ favors higher lattice fillings. This can be achieved either by applying more strain or applying a strong magnetic field. This way one can proceed from 1/4 via 1/3 to 1/2 filling of the lattice. To achieve n HS = 1/2, the magnetic energy must overcome the nnn repulsion V per atom (see Table 2). With the present estimates of the g factor and the nnn amplitudes, the necessary field falls to the 90-130 T range. A further increase to the regime n HS > 1/2 requires to overcome the nn repulsion in the 200-300 meV range, thus cannot be achieved with the available magnetic fields.
Finally, we comment on the HS-HS FM exchange reported in some DFT studies, where the local-density approximation (LDA+U) was employed. The corresponding calculations by the same approach (see Appendix A) show that there is no universal 90 • FM exchange mechanism at work and that the sign of the effective coupling depends strongly on the local interaction parameters. Indeed, the setting on Co sublattice in the spin-state ordered phase is not analogous to two cations hybridizing with two orthogonal orbitals on the common anion, for which the 90 • Goodenough-Kanamori rule applies [27]. In the spin-state ordered phase, the two HS ions hybridize with a common x 2 − y 2 orbital on the LS site and, therefore, the AFM exchange should be expected, which is the case for effective interaction below 7 eV (neglecting the structural relaxation).

Relationship to bulk LaCoO 3
An effective Hamiltonian for bulk LaCoO 3 can be constructed along the same lines as the present model for the strained film. Given the dominant HS-HS repulsion on nn bonds, the question arises why does the bulk system remains uniform [39] despite a high concentration of HS excitations [18,37,40] at temperatures that are order of magnitude lower than the nn interaction. We point out that a spin-state ordered state (checkerboard arrangement of HS and LS sites) is expected at these HS concentrations [41][42][43] even if the T = 0 ground state is a pure LS state [21,44,45]. To understand why the bulk LaCoO 3 at elevated temperature does not exhibit the checkerboard spin-state order or a more complex order as in the strained films, while it stays insulating, is one of the central questions concerning this materials. We suggest that the difference between bulk LaCoO 3 and the strained films is the number of available HS states. While the repulsion between the HS excitations with the same orbital character is strong in all directions, the repulsion between HS excitations with different orbital character is weaker 2 and can be further reduced by anti-ferromagnetic orientation of their spin moments [25]. We hypothesize that in bulk LaCoO 3 the orbital degeneracy allows the HS excitations come closer to each other and prevents formation of the spin-state ordered phase. The tetragonal crystal field of the strained films selects the HSẑ state as the only relevant low-energy excitation, which results in the spin-state ordered configuration supporting FM order.

Summary
We have studied the insulating ferromagnetism of LaCoO 3 films associated with the state disproportionation on Co sites. We have constructed a low-energy effective model in two steps of strong-coupling expansion. In the first step, we convert the multi-orbital Hubbard model into a model of strongly interacting mobile IS excitons and localized HS bi-excitons. In the second step, we eliminate the IS excitons to obtain a model of quasi-immobile HS biexcitons. The approach allows us to estimate the relative strength and magnitude of various couplings in the model.
The key observations are the following. The HS state can be viewed as a bound pair of two IS excitons, which undergo virtual hopping to the neighbor sites. These fluctuations mediate HS-HS interactions beyond the nn bonds. We find that HS-HS interactions along x, y, and z axes are repulsive, while the nnn (±1, ±1, 0) interaction is attractive for ferromagnetic orientation of the HS moments. We identify this interaction as the origin of the FM state of strained LaCoO 3 and thus diagonal FM HS chains can be viewed as the basic building blocks of the FM state. The strong nn and nnn repulsion along crystallographic axes is responsible for the separation of HS chains by at least two nominal LS sites.
We point out that the HS-HS interaction mediated by IS hopping is not captured by LDA+U [13,26,27,42] or LDA plus the dynamical mean-field theory [25,43,46,47] calculations, while both of these approaches take into account the superexchange processes mediated by electron hopping. We have also argued that the main effect explaining the qualitatively different behavior of bulk LaCoO 3 and strained films is the orbital degeneracy of the HS multiplet, which is removed in the films by the tetragonal crystal field.

A DFT results
Computational details. We performed DFT calculations based on the exchange-correlation functional of the local density approximation (LDA) with the all-electron full-potential code wien2k [31]. For simulating the strained structure, a = 3.89Å and c = 3.79Å are chosen for the lattice parameters. The basis size was determined by R MT × K max = 7.0, with the muffin-tin radii of 2.5 for La, 1.91 for Co, and 1.65 for O in units of Bohr radius. To simulate the magnetic phases in the LaCoO 3 system, we constructed a 4 × 4 × 4 supercell. Fig. 7 shows three configurations: FM, AFM, and paramagnetic (PM). The AFM order is set in all a, b, and c directions, resulting a structure of the space group C2/m (no. 12). The Brillouin zone was sampled with a regular mesh containing up to 78 irreducible k points, which is sufficient to obtain a converged spin magnetic moment.
The LDA+U approach within the fully-localized limit (FLL) [48] is performed to investigate the magnetic phases of Co sites due to strong correlation effects. To this end, we employed the effective Hubbard parameter defined by Effects of U and SOC. The differences of the total energies of three configurations in a wide range of U ef f are shown in Fig. 8. We have confirmed that the solutions obtained for all U ef f values (ranging from 4 to 8 eV) consist of HS and LS states only, while the IS contributions are negligibly small. As U ef f increases, the magnetically-ordered configurations become more favorable than the PM one. At U ef f < 7 eV, the AFM instability becomes dominant. U ef f values in this range are reasonable: U ef f = 3.8 eV for Ref. [28] and U ef f = 5.4 − 7.0 eV  for Ref. [27] to study the HS-LS mixture. At U ef f > 7 eV, the FM configuration has the lowest energy, but the energy gain with respect to the AFM setting remains small (about 2 meV/f.u.). Fig. 9 shows the energies at different magnetization directions in the presence of SOC. Our results indicate that the easy axis lies on the ab plane and can not be aligned to the c axis. For the AFM phase, the prefferred orientation corresponds to the [110] direction. Compared to Fig. 8(a), SOC results in an increase of the energy difference E P M − E F M/AF M by about 20 meV that agrees well with the calculated energy lowering of HS z states in Fig. 4.

B Coupling constants
The matrix elements determining the free-particle local and hopping terms in Eq. (1) are obtained from the Wannier projection on the Co d-only model. In the absence of SOC, the local matrix is diagonal in the cubic harmonics basis. Denoting the states d xy , d yz , d zx , d x 2 −y 2 , and d 3z 2 −r 2 by the indices 1-5, respectively, the corresponding diagonal elements in eV units are: ε 1 = 0, ε 2,3 = 0.028, ε 4 = 1.591, and ε 5 = 1.75. The nn hopping term, H t = i,e κλ t  Due to spatial symmetries, the structure of the nonzero matrix elements t (e) κλ in the orbital domain has also a compact form. In particular, for the nn hopping amplitudes in the xy plane, the diagonal elements are t      In the interaction term of the Hamiltonian (1), we employ the Slater parameterization of the intra-atomic electron-electron interaction matrix U κλµν (F 0 , F 2 , F 4 ) with the interaction strength F 0 = 2.1 eV, the Hund's couplingJ = (F 2 + F 4 )/14 = 0.66 eV, and the ratio of Slater integrals F 4 /F 2 = 0.625 [49]. These interaction parameters were found to provide a good quantitative agreement with the RIXS measurements in LaCoO 3 [21,22]. The employed amplitude of SOC ζ = 56 meV is based on the electronic-structure analysis in LaSrCoO 4 [50].