Superexchange Interactions in Orthorhombically Distorted Titanates RTiO3 (R= Y, Gd, Sm, and La)

Starting from the multiorbital Hubbard model for the t2g-bands of RTiO3 (R= Y, Gd, Sm, and La), where all parameters have been derived from the first-principles calculations, we construct an effective superexchange (SE) spin model, by treating transfer integrals as a perturbation. We consider four approximations for the SE interactions: (i) the canonical crystal-field (CF) theory, where the form of the the occupied t2g-orbitals is dictated by the CF splitting, and three extensions, namely (ii) the relativistic one, where occupied orbitals are confined within the lowest Kramers doublet obtained from the diagonalization of the crystal field and relativistic spin-orbit (SO) interactions; (iii) the finite-temperature extension, which consider the effect of thermal orbital fluctuations near the CF configuration; (iv) the many-electron extension, which is based on the diagonalization of the full Hamiltonian constructed in the basis of two-electron states separately for each bond of the system. The main results are summarized as follows. (i) Thermal fluctuations of the orbital degrees of freedom can substantially reduce the value of the magnetic transition temperature. (ii) The anisotropic and antisymmetric Dzyaloshinsky-Moriya interactions are rigorously derived and their implications to the magnetic properties are discussed. (iii) The CF theory, although applicable for YTiO3 and high-temperature structures of GdTiO3 and SmTiO3, breaks down in LaTiO3. Instead, the combination of the many-electron effects and SO interaction can be responsible for the AFM character of interatomic correlations in LaTiO3. (iv) The SE interactions in YTiO3 strongly depend on the details of the crystal structure. Distortions in the low-temperature structure tend to weaken the ferromagnetic interactions.


Introduction
Recently, the titanium perovskites RTiO 3 (where R is the rare-earth element or Y) have attracted considerable attention. These, formally isostructural and isoelectron materials, crystallize in the orthorhombic space group D 16 2h and have only one electron in the t 2g -shell of Ti-atoms. Nevertheless, the magnetic properties of RTiO 3 appear to be extremely sensitive to the magnitude and details of the lattice distortion. More distorted compounds with R= Yb -Gd (including YTiO 3 ) appear to be ferromagnetic (FM). The Curie temperature (T C ) typically varies from about 30 K (R= Y and Gd) till 60 K (R= Dy, Ho, and Tm) [1,2,3]. Less distorted compounds (R= Sm -La) form the so-called G-type antiferromagnetic (AFM) structure, where all nearest-neighbor (NN) spins are coupled antiferromagnetically. The Néel temperature (T N ) varies monotonously from 146 K (LaTiO 3 ) till 48 K (SmTiO 3 ). The complete phase diagram of RTiO 3 can be found in [2]. The origin of the FM-AFM transition and, more generally, the microscopic physical picture underlying the behavior of these Mott-Hubbard insulators with partially filled t 2g -shell is a matter of considerable debates. Different types of theories ranging from the degenerate orbital liquid [4] to the distortion-controlled crystal-field (CF) splitting of the atomic t 2g -levels [5,6,7] have been proposed. The main obstacle with the setting of the proper physical model is related to the fact that there are too many model parameters, which can drastically change the picture if they are chosen in a uncontrollable way. For example, the orbital liquid state is rapidly deteriorated by the lattice distortion away from the simple cubic structure. The conventional CF theories [5,7] crucially depend on the choice of the dielectric constant, which controls the magnitude of the t 2g -level splitting. Thus, in unbiased theories for RTiO 3 , it is very important to reduce to the minimum all arbitrariness related to the choice of the model parameters by using for these purposes first-principles electronic structure calculations. This direction is sometimes called "the realistic modeling of strongly correlated systems" [8].
Our previous works [8,9,10] were devoted to construction of the parameter-free lattice fermion model (more specifically -the multiorbital Hubbard model) for the lowenergy t 2g -bands of distorted perovskite oxides: The model itself is specified in the basis of six Wannier orbitals (the so-called spinorbitals), which are denoted by Greek symbols, each of which is a combination of two spin (s= ↑ or ↓) and three orbital (m= 1, 2, or 3) variables. The site-diagonal part ofĥ ij ≡ h αβ ij includes the CF Hamiltonianĥ CF i and the relativistic spin-orbit (SO) interactionĥ SO i . The off-diagonal (i =j) elements ofĥ ij stand for the transfer integralŝ t ij ≡ t αβ ij . Bothĥ CF i andt ij are diagonal with respect to the spin indices. U αβγδ are the matrix elements of screened Coulomb interactions. All parameters were determined in the ab initio fashion, in the frameworks of the density functional theory. Other details can be found in the original work [11] as well as in the review article [8].
In this paper we discuss how the lattice fermion model (1) for the t 1 2g compounds can be further mapped onto the the spin-1/2 model of the general form whereÂ ij is the three-dimensional tensor, describing various interactions between spins at the sites i and j. It can be further decomposed into isotropic Heisenberg (J ij ), antisymmetric Dzyaloshinsky-Moriya (d ij ), and symmetric anisotropic (τ ij ) interactions, 1 so that the model itself can be rewritten whereσ i = (σ x i ,σ y i ,σ z i ) is the vector of Pauli matrices in the spin subspace and the prefactor S 2 for S=1/2 is included in the definition of the model parametersÂ ij , J ij , d ij , andτ ij . Thus, J ij is the scalar, d ij = (d x ij , d y ij , d z ij ) is the vector, andτ ij = τ ab ij (where a and b stand to denote the x-, y-, or z-components in the orthorhombic coordinate frame) is the three-dimensional tensor satisfying the condition τ xx ij + τ yy ij + τ zz ij = 0. The notation (σ i ,σ j ) stands for the inner product and [σ i ×σ j ] -for the cross product of the vectorsσ i andσ j .
The goal of this work is twofold. On the one hand we will discuss different levels of approximations for the superexchange (SE) interactions underlying derivation of the spin Hamiltonian (2). On the other hand, we investigate the accuracy of realistic modeling for different types of materials. Particularly, how far one can go with the description of not only qualitative but also quantitative aspects of interatomic magnetic interactions in RTiO 3 . For these purposes we will consider four characteristic compounds: FM YTiO 3 and GdTiO 3 , and G-type AFM SmTiO 3 and LaTiO 3 . We will also present details of the parameters of the Hubbard Hamiltonian (1) and interatomic magnetic interactions (2)(3). After brief introduction of the crystal structure of RTiO 3 (Sec. 2) and description of parameters of the Hubbard Hamiltonian as derived from the firstprinciples electronic structure calculations (Sec. 3), we will turn directly to the analysis of the SE interactions. We will start with the one-electron approximation (referring to the form of the ground-state wavefunction in the atomic limit) and consider different types of the one-electron theories, such as the conventional CF theory at T =0 (Sec. 4.1), CF theory with the relativistic SO interactions, which explains the appearance of anisotropic and antisymmetric Dzyaloshinsky-Moriya interactions (Sec. 4.2). We will also estimate the effect of thermal fluctuations of the orbital degrees of freedom on the interatomic interactions between the spins (Sec. 5). In Sec. 4 we will consider a manyelectron extension of the SE theory, where the ground-state wavefunction is determined from the diagonalization of two-electron Hamiltonians constructed separately for each bond of the system. Particularly, we will argue that the behavior of LaTiO 3 is drastically differrent from other compounds. If in YTiO 3 , GdTiO 3 , and SmTiO 3 the type of the magnetic ground state is mainly determined by the CF splitting, which is clearly larger than the energy gain caused by the virtual hoppings underlying the SE processes and the relativistic SO interaction, in LaTiO 3 all three factors are at least comparable making the situation less straightforward. This revives the idea of our earlier work [9], 2 and now we will explicitly show how the combination of many-electron effects in the bonds and the relativistic SO interaction may explain the AFM character of interatomic correlations in LaTiO 3 . Finally, brief summary of the work and main conclusions will be presented in Sec. 6.

Crystal Structure
The titanates RTiO 3 (R= Y, Gd, Sm, and La) crystallize in the orthorhombic space group is D 16 2h in Schönflies notations (No. 62 in International Tables). An example of the crystal structure with the notations of four Ti-atoms composing the primitive cell is shown in Fig. 1 for the case of YTiO 3 [12]. All calculations have been performed in the experimental crystal structure. The low-temperature (LT) structure is currently available only for YTiO 3 and LaTiO 3 , which was measured at T = 2 K [2] and T = 8 K [13], respectively. Moreover, in order to make a direct comparison with our previous works [8,10], we also used another set of parameters for YTiO 3 , corresponding to the room temperature (RT) [12]. For SmTiO 3 and GdTiO 3 , the crystal structure was taken from [2], correspondingly for T = 100 K and 290 K.
Other details of the crystal structure and their implications to the properties of t 2g -compounds can be found in the previous publications [8,10].

Parameters of Hubbard Model
In this section we summarize parameters of the Hubbard model (1), as obtained from the first-principles electronic structure calculations. The procedure used for the construction of the low-energy model (1) as well as for the calculation of the model parameters was described in many details in [8,11]. It is true that the procedure itself relies on a number of approximations. Nevertheless, we would like to emphasize that apart from these approximations, we did not use any adjustable parameters. In this sense, the procedure is parameter-free.
Basically, there are three sets of the model parameters describing The behavior of parameters of the crystal field is explained in table 1, which shows the eigenvalues and eigenvectors obtained after the diagonalization ofĥ CF i . One can Table 1. Eigenenergies (measured in meV from the lowest energy) and eigenvectors obtained after the diagonalization of the crystal-field Hamiltonianĥ CF i . The eigenvectors are expanded over the basis of the xy, yz, z 2 , zx, and x 2 -y 2 orbitals, in the orthorhombic coordinate frame. clearly see that the CF splitting tends to quench the orbital degrees of freedom for all considered compounds, except LaTiO 3 : the lowest level, split off by the crystal field, is nondegenerate and separated from the next level by an energy gap of at least 100 meV. Thus, in the atomic limit, the single electron occupies the lowest t 2g -level, and the degeneracy of the ground state is lifted by the crystal field. The CF splitting decreases in the direction YTiO 3 →GdTiO 3 →SmTiO 3 →LaTiO 3 . Especially, in LaTiO 3 the CF splitting is exceptionally small. In YTiO 3 , there is a clear dependence of the CF splitting on the crystal structure. Basically, there are two effects: (1) the order of the middle and highest t 2g -levels is reversed in the RT structure (in comparison with the LT one); (2) the CF splitting itself is substantially smaller in the RT structure. The structure of t 2g -orbitals, which are split by the crystal field, is rather similar in YTiO 3 , GdTiO 3 , and SmTiO 3 : apart from some numerical differences, the coefficients of the expansion over the basis function have similar values and similar phases. The same tendencies are clearly seen in the RT structure of YTiO 3 after the interchanging of the middle and highest t 2g -levels. Nevertheless, the structure of the t 2g -levels in LaTiO 3 appears to be different. The behavior of transfer integrals in the NN bonds is explained in table 2. The transfer integrals are presented in the local coordinate frame of orbitals, which Table 2. Transfer integrals in the bonds 1-2 and 1-3 (measured in meV) in the local coordinate frame corresponding to the diagonal representation of the crystal-field splitting. The positions of the atomic sites are explained in Fig. 1. Transfer integrals in other bonds can be obtained fromt 12 andt 13 using the symmetry operations of the space group D 16 2h . compoundt 12t13  table 1).
The screened on-site Coulomb interactions for the t 2g -band are expressed through the 3×3×3×3 matrices, in the basis of three Wannier orbitals [8,10,11]. For the distorted compounds, the structure of these matrices is rather complex and may involve many independent parameters. These matrices are directly used in the next section for the analysis of SE interactions without any additional approximations or simplifications. Nevertheless, just for explanatory purposes in this section, the full matrix U αβγδ was fitted in terms of two Kanamori parameters [14]: the Coulomb interaction U between the same t 2g -orbitals and the exchange interaction J . The fitting implies that U αβγδ has the same symmetry as in the spherical environment of isolated atoms. For example, in the process of fitting, it was assumed that the Coulomb interaction U ′ between different t 2g -orbitals is related to U and J by the identity U ′ = U − 2J [14], and the parameter U is the same for all three t 2g -orbitals. The results of such fitting are shown in table 3. The Table 3. Results of fitting of the effective Coulomb interactions in terms of two Kanamori parameters: the intraorbital Coulomb interaction U and the exchange interaction J . All energies are measured in eV. screening of the Coulomb interaction U is sensitive to the local environment in solids: generally, U is larger for the more distorted YTiO 3 and smaller for the least distorted LaTiO 3 [8]. On the contrary, the exchange interaction J is less sensitive to the details of the screening and close to the atomic limit [15].

Superexchange Interactions in the One-Electron Approximation
The superexchange interaction in the bond i-j is related to the gain of the kinetic energy, which is acquired by an electron residing at the atomic site i in the process of virtual hoppings in the subspace of unoccupied orbitals at the atomic site j, and vice versa [16,17]. Let us consider first the approximation adopted in [10,18], where the energy gain caused by virtual hoppings in the bond i-j was computed in the following way: Namely, we start with the lattice of isolated atoms, each of which accommodates one electron, and describe the ground-state wavefunction of such a reference system by the single Slater determinant G. Since in the atomic limit there is only one t 2g -electron per one Ti-site, this is essentially one-electron problem and all many-electron effects emerge only in the process of virtual hoppings. Typically, this justifies the use of the single-determinant approximation for G in the case of titanates [10]. By denoting the occupied one-electron orbitals at the sites i and j as ϕ i and ϕ j , respectively, the Slater determinant G for the bond i-j takes the following form: where the symbols "1" and "2" stand for the coordinates of two electrons, associated with the sites i and j. After that we treat the transfer integralst ij as a perturbation.
Since the Coulomb repulsion U (and U ′ ) is large, it is sufficient to consider the excited configurations accommodating only two t 2g -electrons, which contribute to the second order perturbation theory with respect tot ij . The eigenvalues and eigenvectors of such excited configurations at the site j are denotes as E jM and |jM , respectively. They are obtained after the diagonalization of the Coulomb interactions, in the basis of all possible two-electron Slater determinants at the site j. For six t 2g spinorbitals, there are 6(6−1)/2=15 such determinants [10]. Thus, each |jM is constructed from several Slater determinants and takes into account the correct multiplet structure in the atomic limit. Examples of the atomic multiplet structures corresponding to the exciting two-electron configurations can be found in [10].P j in (4) is a projector operator, which enforces the Pauli principle and suppresses any transfer of an electron into the subspace of occupied orbitals at the site j. In practice,P j projects |jM into the subspace spanned by the Slater determinants of the form where ϕ j is the occupied orbital at the site j and ψ j is any orbital residing at the site j, which can be reached from ϕ i by the transfer integralst ij . Thus, the expression (4) for the energy gain combines elements of the one-electron approximation for |G with the exact many-electron treatment for the excited twoelectron states. On the one hand, since |jM is a combination of several Slater determinants, the method takes into account some many-electron effects in the atomic limit, which may have interesting consequences on the magnetic properties of RTiO 3 . 3 On the other hand, the form of |G is restricted by the single Slater determinant. In this sense, this is an one-electron approach. That is why, in the following we will refer to the expression (4) as to an one-electron approximation for the SE interactions. The multi-determinant analog of (4) will be considered in Sec. 5.

Crystal-Field Theory for the Superexchange Interactions at T =0
In the crystal-field theory, it is assumed that the form of the occupied orbitals {ϕ i } is totally controlled by the crystal field, which is much larger than the energy gain caused by the virtual hoppings (4) as well as the energies of thermal fluctuations. Then, in the one-electron approximation, the isotropic spin couplings J ij is related to the energy gains (4) by the following expression , and the indices ↑ and ↓ explicitly show the spin states of the occupied orbitals ϕ i and ϕ j . The values of J ij between the nearest and some next nearest neighbors are listed in table 4. It also shows the corresponding type of the magnetic ground state expected in the RTiO 3 compounds. The magnetic transition Table 4. Results of the crystal-field theory for RTiO 3 : isotropic superexchange interactions (J ij , measured in meV), corresponding magnetic transition temperatures (T C,N , measured in K) obtained in the random phase and mean-field approximation (results of the mean-field approximation are shown parenthesis), and the type of the magnetic the ground state (GS, where "F" stands for the ferromagnetic state and "A"for the A-type antiferromagnetic state). The positions of the atomic sites are explained in Fig. 1. Depending on the magnetic ground state, the notations T C and T N stand for the Curie and Néel temperature, respectively. temperature is estimated in the random phase approximation (RPA, Appendix A). Thus, in the CF theory, the ground state appears to be FM in the case of GdTiO 3 and SmTiO 3 , and A-type AFM in the case of LaTiO 3 . The ground state of YTiO 3 depends on the crystal structure: FM and A-type AFM for the RT and LT structure, respectively. Nevertheless, as we will see in the next section, both structures yield the same type of the noncollinear magnetic ground state, which combines elements of the FM and A-type AFM ordering within one magnetic structure, when the relativistic SO interaction is taken into account. To certain extent the FM interactions in RTiO 3 can be also stabilized by considering explicitly the e g -band of these compounds [5]. Nevertheless, at the present stage the situation is not completely clear because the same effects would act against the G-type AFM ground state in the case of SmTiO 3 and LaTiO 3 .
The magnetic interactions in YTiO 3 (RT) are in reasonable agreement with the ones obtained in the Hartree-Fock approximation for the Hubbard model (1) with the same parameters of the crystal structure (J 12 = 3.2 ÷ 3.9 meV and J 13 = 1.0 ÷ 1.2 meV, depending on the magnetic state [10]). To certain extend, the same is true for LaTiO 3 , although since the CF splitting is small and allows for the additional change of the orbital ordering in each magnetic state [17], the parameters of interatomic magnetic interactions obtained in the Hartree-Fock approximation for LaTiO 3 exhibit a strong dependence on the magnetic state in which they are calculated (J 12 = 1.0 ÷ 4.5 meV and J 13 = −1.2 ÷ −4.9 meV [10]). Nevertheless, results of the CF theory (table 4) are within the parameters range obtained previously in the Hartree-Fock approximation for LaTiO 3 [10].
RPA considerably reduces the magnetic transition temperature in comparison with the mean-field approach. This effect is particularly strong in YTiO 3 due to the quasi-two-dimensional character of the SE interactions, which according to the Mermin-Wagner theorem suppress the long-range magnetic order at finite T [19]. 4 Nevertheless, J 13 is finite and the reduction of the magnetic transition temperature has only logarithmic dependence of the anisotropy g = |J 13 /J 12 | of NN interactions: [20]. 5 Therefore, the transition temperature remains finite. Moreover, the Curie temperature is overestimated by factor two even in the case of YTiO 3 . The correct G-type AFM ground state is reproduced neither in SmTiO 3 nor in LaTiO 3 .

Relativistic Spin-Orbit Interaction and the Effective Spin Hamiltonian
In the relativistic generalization of the crystal-field theory, the occupied orbital ϕ i at each Ti-site is obtained after the diagonalization of the one-electron Hamiltonianĥ CF i +ĥ SO i , which combines the crystal field and the relativistic SO interactionĥ SO where ξ i ≈ 20 meV is related to the spherical part of the one-electron potential [8]. The Hamiltonian is constructed in the basis of six t 2g spin-orbitals. The lowest eigenstate obtained after the diagonalization is the Kramers doublet, whose eigenvectors can be formally denotes as ϕ 1 i and ϕ 2 i . Then, we consider a liner combination and find the coefficients c 1 i and c 2 i from the condition that the averaged spin moment, e ±a i = ϕ ±a i |σ i |ϕ ±a i corresponding to ϕ ±a i , has the maximal projection along the ±a = ±x, ±y, and ±z axes in the orthorhombic coordinate frame. 6 Then, we use these orbitals in the expression (4) for the energy gain, and calculate 36 parameters T (ϕ ±a i , ϕ ±b j ) corresponding to all possible combinations of ±a and ±b at the sites i and j. In the one-electron (mean-field) approximation, the magnetic part of T (ϕ ±a i , ϕ ±b j ) should 4 Strictly speaking, this behavior does not seem to be consistent with the experimental inelastic neutron scattering data, which are typically interpreted in terms of the three-dimensional isotropic Heisenberg model [21]. However, at present there is no clear consensus on this matter, neither on theoretical nor on experimental side. For example, the experimental orbital ordering pattern determined in [22] is more consistent with the anisotropic (quasi-two-dimensional) structure of interatomic magnetic interactions [10,23]. Thus, the problem requires additional study, both on theoretical and experimental sides. 5 Moreover, the dependence T C,N ∼ −(ln g) −1 can be further modified by the longer range interactions J 23 and J 23 ′ between neighboring ab-planes. 6 The actual procedure was based on the numerical maximization of the function correspond to the energies, which are obtained from (2) by replacing the Pauli matriceŝ σ i by the unit vectors e ± i , describing different directions of spin at the site i. In total, the magnetic part of T (ϕ ±a i , ϕ ±b j ) is characterized by 9 independent parameters, which constitute the tensorÂ ij of interatomic magnetic interactions. The values of A ij for two inequivalent NN bonds are given in table 5 and the decomposition of A ij into J ij , d ij , andτ ij is presented in table 6.
One can clearly see that the Table 5.  (see table 1). Nevertheless, it is interesting to note that the relativistic SO interaction also contributes to the isotropic interactions J ij , and all of them decrease (with some exception of LaTiO 3 ) after taking into account the SO interaction in comparison with results of the pure CF theory in table 4. The directions of the magnetic moments in the ground state and the values of the magnetic transition temperature in the mean-field approximation (Appendix B) are summarized in table 7. As expected, the magnetic ground state is noncollinear. The type of the magnetic ground state for the space group D 16 2h can be formally denoted as X-Y-Z, where X, Y, and Z is the magnetic structure (F, A, C, or G) formed by the projections of the magnetic moments onto the orthorhombic axes a, b, and c, respectively. By comparing the values of the magnetic transition temperature with the ones reported in table 4, one can clearly see that the SO interaction tends to additionally increase T C , Table 6. Isotropic Heisenberg interactions (J ij ), antisymmetric Dzyaloshinsky-Moriya interactions (d ij ), and symmetric anisotropic (τ ij ) interactions associated with the nearest-neighbor bonds 1-2 and 1-3 as obtained in the crystal-field theory with the relativistic spin-orbit interaction. All interactions are measured in meV. The positions of the atomic sites are explained in Fig. 1.  Table 7. Magnetic ground state, direction of magnetization in the ground state and the magnetic transition temperature (T C,N , measured in K) in the mean-field approximation as obtained in the crystal-field theory with the relativistic spin-orbit interaction. The vector of magnetization is referred to the site 1 in Fig. 1.

Thermal Fluctuations of the Orbital Degrees of Freedoms
In this section we consider the effect of thermal fluctuations of the orbital degrees of freedom on interatomic magnetic interactions between the spins. For these purposes, it is convenient to work in the local coordinate frame corresponding to the diagonal representation of the CF Hamiltonianĥ CF i at each atomic site. Then, we assume that for each projection of spin, the three-component occupied t 2g -orbital at the site i can be presented in the form of the real vector in the basis of three CF orbitals listed in table 1, where 0 ≤ θ i ≤ π/2 and 0 ≤ φ i ≤ π due to the invariance of both crystal field and pair interactions (4) with respect to the inversion v i → −v i . In these notations, the point θ i =φ i =0 corresponds to the lowest CF orbital. This situation was already considered in Sec. 4.1 in the limit of large CF splitting. However, since the CF splitting is finite, other orbital configurations can contribute to the thermodynamic averages of physical quantities at elevated temperatures. Thus, similar to Sec. 4.1, we assume that the CF splitting is much larger that the energy gain (4) caused by the virtual hoppings, but can become comparable with the energies of thermal fluctuations, and consider the finitetemperature extension of the CF theory. Then, the thermal average of (4) is given by where and T (θ i , φ i , θ j , φ j ) is the pair interaction (4) constructed from the occupied t 2g -orbitals of the form (5) for either ferromagnetic (↑↑) or antiferromagnetic (↑↓) configurations of spins (in the following denoted as T ↑↑ ij and T ↑↓ ij , respectively). The numerical integration in (6) was performed by using the the Metropolis algorithm [24,25]. An example of the temperature dependencies of T ↑↑ ij (T ) and T ↑↓ ij (T ) for YTiO 3 and LaTiO 3 is shown in Fig. 2. 7 Then, using T ↑↑ ij (T ) and T ↑↓ ij (T ), one can evaluate the temperature dependence of interatomic interactions between the spins, 2J ij (T )=T ↑↑ ij (T )−T ↑↓ ij (T ), averaged over all orbital configurations, and self-consistently solve the equation for the magnetic transition temperature in RPA (Appendix A), where the value of T C,N is used as the argument of J ij (T ), and the procedure is repeated until reaching the self-consistency with respect to T C,N . Then, we obtain the following values of the Curie temperature T C = 60, 68, and 37 K for YTiO 3 (RT), GdTiO 3 , and SmTiO 3 , respectively. YTiO 3 (LT) and LaTiO 3 are expected to develop the A-type AFM order with the Néel temperature T N = 64 and 52 K, respectively. Thus, the thermal fluctuation systematically decreases the values of the magnetic transition temperature (by 10-44 %). As expected, the largest change is observed in LaTiO 3 , which has the smallest CF-splitting, and in YTiO 3 , due to the quasi-two-dimensional character of interatomic interactions. Nevertheless, SmTiO 3 remains FM and LaTiO 3 -A-type AFM, contrary to the experimental data.

Multi-Determinant Approach for Superexchange Interactions
By summarizing results of the previous sections, we note the following shortcomings of the one-electron approach: • In the case of LaTiO 3 and SmTiO 3 , it fails to predict the correct G-type AFM ground state. Although small G-type AFM component (along the a-axis) is expected in the theoretical magnetic ground state of SmTiO 3 after including the relativistic SO interaction, that of LaTiO 3 does not involve the G-type AFM arrangement. • Even for the FM compounds YTiO 3 and GdTiO 3 , the magnetic transition temperature is typically overestimated by factor two. To certain extent, it can be reduced by thermal fluctuations of the orbital degrees of freedom. On the other hand, the relativistic SO interaction acts in the opposite direction and additionally increases T C (at least on the level of mean-field approximation).
In this section, we investigate whether these problems can be resolved by considering the many-electron effects.
As was pointed out in the beginning of Sec. 4, the main drawback of all previous considerations was that, although the expression (4) for the energy gain takes into account the many-electron effects in the intermediate configurations t 2 2g , which can be reached in the process of virtual hoppings, it was combined with the one-electron approximation for the ground-state wavefunction G in the atomic limit. The purpose of this section is to go beyond this one-electron approximation and to clarify the role played by the many-electron effects in the problem of SE interactions. Namely, for each bond i-j, we construct the complete basis of two-electron Slater determinants of the form, where ψ i denotes the one-electron spin-orbital residing at the site i. In practice, these orbitals were obtained from the diagonalization of the CF Hamiltonianĥ CF i (table 1). In total, there are 6 such spin-orbitals and 36 Slater determinants of the form (7). Then, for each bond, we calculate matrix elements of the pair interactions in the basis of these Slater determinants, and combine them with the matrix elements of one-electron operators of the crystal fieldĥ CF i +ĥ CF j and (optionally) the SO interactionĥ SO i +ĥ SO j . The corresponding 36×36 Hamiltonian matrix is denoted aŝ Then, the spectrum of two-electron states in each bond is obtained after the diagonalization ofĤ ij . Thus, the basic idea (and approximation) behind this treatment is that the entire lattice is divided into "molecules", and the electronic structure of each "molecule" can be considered independently from other "molecules". We will mainly focus on the origin of AFM correlations in LaTiO 3 .
First, let us discuss results without relativistic SO interaction. The spectrum of two-electron states, obtained after the diagonalization ofĤ ij for the bonds 1-2 and 1-3 is shown in Fig. 3. There is a drastic difference of LaTiO 3 from other compounds. In YTiO 3 , SmTiO 3 , and GdTiO 3 , four low-energy levels, including one spin-singlet and three degenerate spin-triplet states, are clearly separated from the next states by an energy gap of at least 100 meV, corresponding to the value of the CF splitting (table 1). However, in LaTiO 3 , the singlet-triplet splitting in the low-energy part of the spectrum is at least comparable with energy gap, separating these levels from the next two-electron states. For example, in the bond 1-3, the singlet-triplet splitting is about 20 meV. However, the next level is located only 4 meV higher than the triplet state. Thus, the crystal-field theory, although applicable for YTiO 3 and high-temperature structures of SmTiO 3 and GdTiO 3 , definitely breaks down in the case of LaTiO 3 , where the splitting of the two-electron levels caused by the SE effects in the low-energy part of the spectrum is at least comparable with the CF splitting.
The parameter of the isotropic exchange coupling in the bond i-j can be expressed through the energies of the low-lying spin-singlet (E S ) and spin-triplet (E T ) states, where the prefactor 1/4 stands for S 2 , which according to (3)  decrease the magnetic transition temperature. This behavior is solely related to the quasi-two-dimensional character of interatomic magnetic interactions. Thus, it seems unlikely that the two-electron effects alone will stabilize the G-type AFM ground state in the case of LaTiO 3 . Therefore, we investigate the last possibility related to the relativistic SO interaction. For the most of the considered systems, the SO interaction is small in comparison with the CF splitting and does not significantly change the distribution of the two-electron states. The typical example for the LT structure of YTiO 3 is shown in Fig. 4, which is practically identical to the spectrum without the SO interaction (Fig. 3). For the FM systems, the splitting of the low-lying spin-triplet states by the SO interaction is very small. Nevertheless, below we will see that this splitting may have rather interesting consequences on the magnetic properties. Again, one clear exception is LaTiO 3 , where the strength of the SO interaction (ξ = 21 meV) is Table 8. Parameters of isotropic exchange coupling (J ij , measured in meV) derived from the singlet-triplet splitting of the low-energy levels, which were obtained from the diagonalization of pair and crystal-field interactions in the complete basis of twoelectron Slater determinants separately for each bond of the system; corresponding magnetic transition temperatures (T C,N , measured in K) obtained in the random phase approximation; and the type of the magnetic ground state. Depending on the magnetic ground state, the notations T C and T N stand for the Curie and Néel temperature, respectively.
compound  Table 9. Parameters of interatomic spin correlations in the ground state calculated separately for each bond of the system after the diagonalization of the two-electron Hamiltonian including pair, crystal field, and relativistic spin-orbit interactions. The atomic positions are explained in Fig. 1. For the FM systems, this picture may substantially differ from the one without the SO interaction. However, it should be noted that the details of the ground state in the FM case are determined by the small splitting of the low-lying spin-triplet levels by the SO interaction. For the most of the systems this splitting is small and typically varies from few hundredths of meV to one meV. For example, for the FM bond 1-2 in YTiO 3 (LT), the second and third levels are split from the lowest one by 0.56 meV and 0.59 meV, respectively. Therefore, when the temperature exceeds T ∼ 0.6 meV/k B ≈ 7 K, the proper picture for the magnetic interactions in this compound would correspond to the (thermal) average σ a Similar results can be obtained for other systems forming the FM bonds. After the averaging, the longitudinal spin correlations become isotropic, while all transversal correlations are considerably smaller. This behavior is consistent with the form of the interatomic magnetic interactions tensorÂ ij obtained on the level of one-electron approximation (table 5). Another effect, which can mix the two-electron states within the lowest manifold of certain bond is related to the virtual hopping in neighboring bonds, which involve one of the atomic sites of the original bond. Note that the manyelectron effects will mix the lowest CF orbital at certain atomic site with other states.
The symmetry of these states as well as the magnitude of the mixing will be generally different for different bonds. Since each atom participate in several bonds, this effect will lead to the addition mixing of the low-lying (spin-triplet) states obtained from the diagonalization of the two-electron Hamiltonian separately for each bond of the system. Nevertheless, the quantitative estimate of this effect requires a more rigorous solution of the many-electron problem, which is beyond the scopes of the present work.
In the case of LaTiO 3 , all longitudinal correlations in the bond 1-3 are antiferromagnetic, while the transversal correlations are small. The correlations in the bond 1-2 are strongly anisotropic: if yyand zz-components favor the FM coupling, the xx-correlations are antiferromagnetic. The splitting between the lowest and the next two two-electron levels in the bond 1-2 is about 4 meV, which is considerably larger than the splitting obtained in other (more distorted) compounds. Therefore, the effect is expected to be more robust against the thermal fluctuations mixing the spin-triplet states. Strong transversal correlations are also expected. Thus, the xxcorrelations appears to be antiferromagnetic simultaneously in the bonds 1-2 and 1-3, being consistent with the G-type AFM structure. On the basis of this analysis, we expect that if the G-type AFM order took place in LaTiO 3 , it would be more likely developed by the x-(a-) projections of the magnetic moments in the orthorhombic coordinate frame.

Summary and Conclusions
Starting from the multiorbital Hubbard Hamiltonian for the t 2g -bands of orthorhombically distorted titanates RTiO 3 (R= Y, Gd, Sm, and La), where all the parameters were derived from the first-principles electronic structure calculations, we considered different levels of approximations for the construction of the spin-only superexchange model. Namely, after considering the conventional crystal-field theory, where all interatomic magnetic interactions are solely determined by occupied t 2g -orbitals, that are split off by the crystal distortion, we consecutively incorporate the effects of the relativistic spin-orbit interaction, thermal fluctuations of the orbital degrees of freedom, and many-electron effects related to the virtual electron hoppings in the bonds.
Even in the conventional CF theory, the interatomic magnetic interactions appear to be extremely sensitive to the details of the orbital structure, and small change of the orbital structure, caused by either crystal distortions or thermal fluctuations, can have a profound effect on the magnetic properties.
Particularly, the use of the room-and low-temperature structure for YTiO 3 provides rather different sets of parameters of the spin Hamiltonian. The additional distortion in the low-temperature structure tends to weaken the ferromagnetic interactions: it substantially reduces the FM coupling J 12 in the orthorhombic ab-plane and makes the coupling along the c-axis weakly AFM. Similar tendency was reported for LaTiO 3 , where the use of the low-temperature structure [13] could help to stabilize the G-type AFM state [5,7,26]. Apparently, some of the problems encountered in the present work, such as the incorrect magnetic ground state in SmTiO 3 and the overestimation of the Curie temperature in GdTiO 3 , may be resolved by using the low-temperature structural data for these compounds, which are not available today.
Moreover, the thermal fluctuations of the of the orbital degrees of freedom near the CF configuration can substantially reduce the value of the magnetic transition temperature (up to 40%).
The relativistic spin-orbit interaction, which is responsible for the appearance of anisotropic and antisymmetric Dzyaloshinsky-Morita interactions, leads to the noncollinear magnetic alignment. The role of the relativistic effects increases in the direction YTiO 3 →GdTiO 3 →SmTiO 3 →LaTiO 3 , when the crystal distortion decreases.
The crystal-field theory, although applicable for YTiO 3 and high-temperature structures of GdTiO 3 and SmTiO 3 , definitely breaks down in the case of LaTiO 3 , which has the smallest CF-splitting and where other factors, such as the relativist SO interaction and many-electron effects in the bonds, start to play an important role. Particularly, we found that the combination of the latter two factors could explain the G-type AFM character of interatomic correlations in LaTiO 3 . where J(q) is the Fourier image of J ij and Q is the vector describing the magnetic structure in the ground state, to the case of the complex lattices containing n magnetic atoms in the primitive cell. First, we transform the parameters of the Heisenberg model to the local coordinate frame: J ij → J ij = s ij J ij , where s ij = +1 and −1 corresponds to the FM and AFM arrangement of spins in the bond i-j in the magnetic ground state. For a collinear magnetic configuration, this procedure is equivalent to the shift of the origin of the Brillouin zone in the right-hand side of (A.1). Then, we construct the dynamical matrix: is the Fourier image of J ij acting between the atomic sublattices ℓ and ℓ ′ , N is the number of the primitive cells, and R i is the radius-vector of the site i. Then, by diagonalizing J(q), one can find its eigenvalues ω ℓ (q) (the "magnon energies") and by repeating the RPA arguments [27] derive the following expression for the magnetic transition temperature: The mean-field approximation for the relative magnetization σ i (T ), which is the temperature average ofσ i , is formulated in the following way. The mean field (or the molecular field), corresponding to the spin Hamiltonian (2), is given by Then, the temperature average ofσ i in the molecular field h i has the following form