A general procedure to evaluate many-body spin operator amplitudes from periodic calculations: application to cuprates

A general procedure is presented which permits the form of an extended spin Hamiltonian to be established for a given magnetic solid and the magnitude of its terms to be evaluated from spin polarized, Hartree–Fock or density functional calculations carried out for periodic models. The computational strategy makes use of a general mapping between the energy of pertinent broken-symmetry solutions and the diagonal terms of the spin Hamiltonian in a local representation. From this mapping it is possible to determine not only the amplitude of the well-known two-body magnetic coupling constants between near-neighbor sites, but also the amplitudes of four-body cyclic exchange terms. A scrutiny of the on-site spin densities provides additional information and control of the many broken-symmetry solutions which can be found. The procedure is applied to the La2CuO4, Sr2CuO2F2, Sr2CuO2Cl2 and Ca2CuO2Cl2 square lattices and the SrCu2O3 ladder compound. It is shown that a proper description of the magnetic structure of these compounds requires that two- and four-body terms are explicitly included in the spin Hamiltonian. The implications for the interpretation of recent experiments are discussed.


Introduction
The discovery of the anomalous properties of high-T c superconducting cuprates (HTCSs) in the late eighties has triggered a considerable interest in the crystal and electronic structure of these compounds from both experimental and theoretical points of view [1]- [5]. An enormous research effort on ceramic materials, mostly based on copper oxides, has been carried out continuously to try to improve the properties of known structures and synthetic pathways and has resulted in the synthesis of a wide variety of cuprates. The impressive richness of lowdimensional magnetic behavior of the different copper compounds can be, to a large extent, traced back to the stacking of the distorted CuX 6 octahedra (or CuX 5 pyramidal or CuX 4 planar units) in the lattice [4,6]. Most of the cuprate based materials are formed by almost independent CuO 4 units with distant apical ligands to complete the strongly distorted CuO 4 X 2 or CuO 4 X units and, hence, the structure is dominated by the link between CuO 4 units. Depending on the nature of the counter ions and on the number of links between the different CuO 4 units, different structures can be formed ranging from the typical lamellar two-dimensional (2D) structure of the HTCSs to many lower dimensional structures by different combinations of edgesharing and corner-sharing CuO 4 plaquettes (or CuO 4 units) that can give rise to spin ladders (e.g. Sr n−1 Cu n O 2n−1 series with n 2 [7]- [9]), zigzag spin chains (e.g. SrCuO 2 [10]) and quasi-1D systems (e.g. A 2 CuO 2 (A = Ca,Sr) [11,12] or Li 2 CuO 2 [13]) formed by edge-sharing CuO 4 units.
The electronic ground state of this kind of materials is usually described by the open shell nature of the Cu 2+ ions arranged in the CuO 4 units in which the Cu(3d 9 ) atomic configuration gives rise to a d x 2 -y 2 type hole with the lobes pointing towards the O ions. The resulting Cu-O-Cu pathways range from ∼ 90 • to 180 • and they are responsible for the rich variety of low-dimensional magnetic structures dominated by moderate ferromagnetic (FM) to strong antiferromagnetic (AFM) interactions. From the theoretical point of view, these systems are strongly correlated in nature, making standard band theory techniques based on density functional theory (DFT) unable to accurately describe either their valence or low energy spectrum [14]. However, it has been shown that hybrid exchange-correlation functionals can provide reliable descriptions of strongly correlated transition metal magnetic systems ( [15] and references therein). 3 Most of the HTCS materials have a lamellar structure in which strong AFM interactions take place along 180 • Cu-O-Cu bonds in edge sharing Cu 4 O 4 plaquettes leading to a 2D network of effective spin S = 1/2 particles. These strong magnetic interactions observed in the HTCS are thought to be fundamental ingredients of the high-T c superconductivity microscopic mechanism [5,16]. Since these compounds may be regarded as effective S = 1/2 spin lattices, their low energy spectrum and collective properties are assumed to be governed by a simplified Heisenberg Hamiltonian as in equation (1) which only accounts for the magnetic coupling J i j between nearest-neighbor (NN) centers i and j. This is in agreement with the widely accepted general picture for HTC superconductivity involving a 'Heisenberg sea' where holes or electrons are introduced by doping the perfect structures. However, it has been claimed that to fully understand the magnetic excitations, infrared and neutron scattering spectra of 2D [17]- [23] and spin ladder cuprates [24,25] it is necessary to extend the spin Hamiltonian as in equation (2) H In this expression the constants 1/4 and 1/16 have been introduced to define the zero of energy as that of the FM solution. In this spin model, the signs and amplitudes of the local intersite magnetic interactions govern the collective properties of a spin lattice. They appear as the basic ingredients of the effective spin Hamiltonian which in full generality involves not only the two-body exchange J i j but also other interactions, such as those represented by the four-body cyclic term J ring , or even higher-order terms. The largest two-body couplings are expected to occur between NN sites although next-NN (NNN) interactions may be non-negligible or even of the same order of magnitude (cf CuGeO 3 system [26,27]). Regarding the four-body operator terms, they may be important in Cu 4 O 4 plaquettes since their origin lies in the cyclic circulation of electrons around the ring. Their importance and that of analogous cyclic six-body effects have been pointed out in other types of half-filled band systems such as the π system of conjugated organic molecules [28]. Similar four-body terms are crucial to describe the ground state properties of 3 He [29].
The direct determination of the amplitude of the many-body terms of an extended spin Hamiltonian such as the one in equation (2) from experiment is, in general, impossible. Spin ladders with a variety of intersite distances represent an especially difficult case. In general, the experimental determination of the coupling constants is based on a series of hypotheses about the negligibility of interactions between 'remote' sites. From these hypotheses, a given spin model is assumed and validated only from a numerical fit of the thermodynamic or spectral properties. It is customary to consider NN interactions only although this may be an excessive simplification and eventually can lead to contradictory estimates of the dominant couplings. This is precisely the case of the cuprate spin ladders for which conflicting values of the J rung /J leg ratios ranging from 0.5 to 1.0 were proposed. Interactions involving sites at a longer distance or 4 involving various sites had to be invoked to rationalize the different experimental results arising from different techniques. One may think for instance in inter-ladder interactions, diagonal terms in the plaquettes or four-body cyclic effects. Indeed, several theoretical studies seem to consistently indicate that four-body terms (J ring ) are crucial [21]- [24]. Here, it is important to point out that very recently Toader et al [30] provided strong experimental evidence of the importance of the J ring term in La 2 CuO 4 with J ring ≈ 0.5 J . This value is comparable to the pairing energies and strongly suggests that the resulting circulating currents could have an important role in the mechanism of superconductivity. This is in contrast with previous estimates of J ring for 2D and spin ladder cuprates, obtained either from indirect measurements or from numerical simulations with an extended Heisenberg model, which propose substantially smaller amplitudes with J ring ∼ 0.3 J [19,21,22,24,31]. The origin of these discrepancies relies on the choice made by different authors for the magnitude of the other coupling constants in the spin Hamiltonian of equation (2). Hence, the J ring term as evaluated by Toader et al relies on NN and NNN coupling constants of J = 111.8 meV and J d = −11.4 meV, respectively, extracted from one of the various fittings of the magnon spectrum [22]. However, it is important to point out that this value for J is smaller than another experimental estimate of 135 ± 6 meV obtained with a NN Heisenberg Hamiltonian [32]. Moreover, one should note that the present estimate of a FM J d differs from previous theoretical predictions [33] and from fitting to experimental measurements on materials with similar exchange paths [34]. However, one should also recognize that the theoretical study by Annett et al [33], takes only into account twoelectron processes to evaluate this J d term while, as clearly explained by Toader et al [30] and references therein, the FM character of J d arises precisely from three-electron exchange processes around a plaquette, which are of the same order of magnitude as four-electron exchange processes around the same plaquette and which are not taken into account in [33,34]. The discussion above illustrates the difficulties faced by experimentalists when attempting to extract the magnitude of the important terms and the need for independent accurate and unbiased theoretical predictions. For the SrCu 2 O 3 ladder compound a similar situation is encountered; the recent Raman response experiments by Schmidt et al suggest J rung = 140 meV, J ring /J rung = 0.2 and J leg /J rung = 1.5 [35], in contrast with previous work indicating a more isotropic behavior between rungs and legs [36].
From the preceding discussion it is clear that an accurate prediction of the various magnetic interactions entering in the spin Hamiltonian as in equation (2) is not only highly desirable to understand the ground state properties of this kind of systems but urgently needed. One could, for instance, start from some approximate electronic Hamiltonian such as a single-band model involving only the magnetic centers or a two-band model involving also the electrons of the bridging ligands. However, this is likely to introduce even more problems since it is difficult to assess the accuracy of such a multiparametric approach. An alternative and straightforward way to an unbiased estimate is the direct evaluation of the amplitude of the relevant magnetic interactions on a realistic model system treating all the electrons in a large enough basis set. If practicable, this pragmatic approach will present the advantage of providing an independent, unbiased and consistent prediction of the amplitude of relevant exchange parameters that can solve the difficulties encountered by the usual fitting techniques used by experimentalists.
Unfortunately, the direct calculation of the important terms is not a simple task since it also requires the use of a model for the solid. In a first approach, one may neglect translational symmetry and define a properly embedded finite cluster, a fragment of the periodic lattice with two, three or four magnetic sites with their coordination ligands and perform the best ab initio (using the exact Hamiltonian) explicitly correlated calculations. This strategy, in particular when using the so-called difference dedicated configuration interaction (DDCI), has provided very consistent and reliable evaluations of the two-body interactions in a large series of perovskites, oxides and 2D cuprates [37]- [41]. However, this approach rapidly faces computational limits due to the need for very large embedded cluster models. This is especially the case when attempting to extract four-body terms in ladder compounds [42]- [45]. Alternatively, one may exploit the periodic symmetry of the crystal but in this case explicitly correlated calculations are not feasible and one must rely on spin-polarized mean-field type approaches. In such a case, only the FM solution may be considered as a valid approximation of an eigenstate of the problem with properly defined spin quantum numbers [46,48]. The solutions with lower values of the square of the total spin operator S 2 (lower magnetization) cannot be properly treated. Fixing a priori different localized periodic spin distributions one may obtain a set of distinct selfconsistent solutions of different energies. However, these solutions are not eigenfunctions of the S 2 operator; they are (spin) broken-symmetry solutions. We will show that assigning the energy of the broken-symmetry solutions to the expectation value of the corresponding spin distribution of a Heisenberg Hamiltonian enables one to obtain estimates of the magnetic coupling constants. Such approaches have been rather extensively used in the field of molecular magnetism and may also be employed on embedded clusters to extract the two-body terms [48,49]. The mean-field calculations may use the exact Hamiltonian, and the corresponding information comes either from unrestricted broken-symmetry Hartree-Fock (UHF) solutions or from the similar ones obtained by means of DFT approaches. In the first case, the resulting evaluations of the AFM couplings are severely underestimated [15], [50]- [52]. In the second one, the introduction of electron correlation effects improves the result but the numerical values are highly sensitive to the chosen exchange-correlation potentials. Notice that the local density approximation (LDA) and the different generalized gradient approximations (GGA) in DFT frequently lead to metallic solutions, closed shell in nature [15], which miss the main physical features of the low energy states and cannot provide any information about the magnetic couplings. Actually, the best results are obtained with hybrid functionals initially proposed by Becke [53,54], for which a percentage of Fock exchange is mixed with a given DFT exchange-correlation functional. However, the results are again strongly dependent on the amount of Fock exchange included in the exchange potential [55,56]. Previous calculations have shown that the best percentage, at least for NiO [15] and cuprates [57,58], is around 35% Fock. Here, it is worth pointing out that hybrid DFT calculations need to include the non-local Fock contribution which can be accurately computed when the total density is expressed in terms of Gaussian type orbitals (GTO) basis sets. Recently, however, it has been shown that a reasonable accuracy can be reached employing plane waves as well [59].
In the present work, we report a generalization of the periodic symmetry-broken (UHF or DFT) approach which permits one to obtain not only the various NN and NNN most important two-body terms, as initially done for MF 2 (M = Cu, Ni, Co, Fe, Mn) compounds [52,60], but to predict higher-order terms such as the four-body cyclic terms appearing in equation (2). To this end broken-symmetry hybrid DF computations are carried out for the La 2 CuO 4 , Sr 2 CuO 2 F 2 , Sr 2 CuO 2 Cl 2 and Ca 2 CuO 2 Cl 2 2D square lattices and the SrCu 2 O 3 two-leg spin ladder. A general theoretical framework is presented that permits one to extract accurate values for the many-body terms of extended spin Hamiltonians from these periodic first-principle calculations. In addition, it will be shown that these periodic calculations confirm the importance of the FM interladder exchange and of the four-body intra-plaquette operator. 6 This paper is organized as follows: in section 2, we develop the relevant methodological aspects to extract the amplitude of the parameters of a general spin Hamiltonian from periodic all-electron calculations of the electronic structure of the systems. First, the logics leading to a mapping between spin eigenfunctions and symmetry-broken solutions is developed in a general scheme starting from the Hubbard model for a plaquette with few remarks on the calculated solutions. Section 3 describes the computational details of the calculations. The results obtained for the planar cuprates and the ladder system are reported in section 4 and, finally, in section 5, we present the conclusions of this investigation and further extensions of the theory.

Theoretical background and methodology
Let us assume that one may calculate a mean field single determinant description of the FM state in a given unit cell. This wavefunction is usually expressed in terms of symmetry-adapted (delocalized) Bloch functions. An appropriate localizing transformation of the singly occupied functions will define atom-centered Wannier magnetic orbitals [61,62]. The FM solution for the system may be written as where core stands for the closed-shell part of the wavefunction and a, b, c, . . . , n represent the magnetic orbitals localized on the different magnetic centers. In this picture there is one electron per site and one can refer to the corresponding determinant as a valence bond (VB) neutral form (see also the definition in [63]). In the following, the expectation energy Φ Sz,max |Ĥ |Φ Sz,max = E 0 will be taken as energy origin (E 0 = 0). The n electrons in n orbitals define a half-filled band problem. Each of the possible distributions of the magnetic electrons in the magnetic orbitals defines a determinant and the set of these determinants spans an orthogonal VB basis set. The matrix elements of the exact Hamiltonian in this basis define a valence configuration interaction matrix. In order to understand the physical effects arising from the interaction between the α and β (or spin up and spin down) electrons, we shall consider first a finite set of neighbor sites, the periodicity being introduced in a subsequent step. One must first recall that the lowest energy determinants are the neutral ones (in the VB sense), for which any magnetic orbital is singly occupied. The ionic determinants set, presenting positively and negatively charged sites, are of higher energy. Let us consider a given VB neutral determinant as built with the local magnetic orbitals of the FM solution and where the bars indicate, as usual, beta spin electrons (hence, i is synonymous of iα andī of iβ). The expectation energy (within the exact non relativistic Hamiltonian) of such a determinant is given by equation (5) where α(I ) and β(I ) are the set of spin orbitals of α and β spin in Φ I and iα and jβ the corresponding spin orbitals. This expectation energy is higher than that of Φ Sz,max since it misses the direct (positive) exchange integrals between the α and β spin orbitals. It is reasonable to assume that the nonzero K i j integrals concern NN and NNN site pairs only. These integrals 7 have been estimated in cuprates [42]- [45], [57,64,65] and play an important role in defining the final value of the magnetic coupling constant. However, for the perturbative development given below we consider that the difference among the different K i j integrals is small enough to be neglected in front of the excitation energies to the VB ionic determinant and a common zeroorder energy E 0 will be kept for the VB neutral determinants in the perturbative developments. The strategy will now consist of establishing a relation between the diagonal elements of the Heisenberg Hamiltonian and the energies of symmetry-broken UHF solutions. This is accomplished by means of suitable perturbation expansions.
The Heisenberg Hamiltonian is obtained from the exact one by considering the space of all VB determinants Φ I as a model space [66] in the frame of the effective Hamiltonian theory of Bloch or des Cloizeaux [67,68]. Let us callP the projector on this subspace of dimension N the effective HamiltonianĤ eff is defined by the N eigen equationŝ where the vectors | m are the N eigenvectors ofĤ having the largest projections on the model space, and E m the corresponding eigenvalues. Since all VB neutral determinants have the same spatial part, the effective Hamiltonian may be written in terms of spin operators and, in principle, it will involve many-spin operators. Due to the large energy gap between the neutral and VB ionic determinants the effective Hamiltonian may be approached through the quasi-degenerate perturbation theory (QDPT). Let us consider the diagonal matrix element associated with a determinant Φ I . To the second-order where we have considered that the direct exchange integrals K i j are small in front of the energy differences between the neutral and the VB ionic configurations in the evaluation of the secondorder corrections. The determinant Φ I interacts with all the determinants obtained by charge transfers between α and β sites leading to ionic determinants of the form and The amplitude of the hopping integral t i j rapidly decreases with the distance between sites i and j. Comparing Φ I |Ĥ eff |Φ I with equation (2) one obtains 8 which is a rather well-known relation. Let us suppose now that one tries to optimize the energy of a symmetry-broken determinant Φ I obtained from Φ I by relaxing the spin orbitals. To the first-order this optimization is equivalent to a CI between Φ I and the singly excited determinants leading to the spinorbital rotations giving delocalization tails to the α spin orbital on site i on to the orbitals of the neighbor sites occupied by β spin electrons (or equivalently, which only differs from I by second-order effects and its energy only differs from the second-order corrected energy associated with I by thirdorder terms Hence, Equation (19) relates the calculated energy of the I' broken symmetry solution for the exact Hamiltonian, relative to the FM state, to the corresponding energy expression of the Heisenberg Hamiltonian in equation (1) and it provides a first useful relationship to extract the relevant twobody magnetic coupling constants from the different broken-symmetry solutions. Since this relation is also valid when Φ I represents a periodic spin distribution, equation (19) provides a practical way to extract these effective parameters from periodic calculations. This is precisely the procedure followed in previous works [52,60]. However, it is possible to use the perturbation developments above to generalize the procedure in such a way that extraction of higher-order terms in the spin Hamiltonian as in equation (2) becomes straightforward. Before continuing the development a caveat is necessary since the Wannier functions as in equation (3) or (4) refer to the whole crystal. In practice, one has to limit the number of spin distributions to a finite number which is given by all possible spin permutations in a large enough supercell. Convergence of the results with respect to the model chosen can be indeed verified by extending the supercell. However, since the magnetic coupling constants are local parameters [37], the results are not expected to vary significantly upon enlarging the supercell used to extract the relevant parameters. Therefore, starting from n non equivalent spin distributions Φ I , Φ J , . . . , Φ M , one may, in principle, determine m quantities J i j and, of course, discriminate the leading terms from 9 the negligible ones. As will be shown later on, it may happen that several spin distributions lead to identical relationships between the calculated broken-symmetry energies and the set of J i j values. However, if one performs a large enough number of broken-symmetry calculations it will be always possible to obtain a sufficiently large number m of independent equations from different guesses, and since one only looks for n < m spin couplings, it becomes also possible to evaluate the consistence of the calculated set of magnetic coupling constants. The procedure is not only general but allows an auto check of the obtained results.
In the forthcoming discussion, we will show how to obtain the higher order terms by exploiting the perturbation development introduced above. However, before proceeding with such a generalization, it is convenient to consider the on-site spin densities. In principle, an approximate estimation of the spin density on site i in spin mean-field solution I, ρ i (I ), can be obtained from the first-order corrected wavefunction. Notice, however, that for simplicity we are heretofore using I to denote Φ I . Let us recall that this perturbation treatment-using the Hubbard model Hamiltonian as zero-order Hamiltonian-gives, in an approximate way, the effect of relaxing the orbitals from the FM reference to the self consistent orbitals for the particular broken symmetry I solution, assuming that the relaxation comes from up to first-order mixing neutral and VB singly ionic determinants. It can be shown that ρ i (I ) is thence given by The second term in equation (20) reflects the impact of the VB singly ionic determinants in the final spin density obtained with the relaxed orbitals. If the largest J amplitudes concern NN pairs of sites and they are of equal or of similar magnitude, one may expect that where λ = 2(t i j /U ) 2 and n I NNβ is the number of β spin NN atoms of a given α spin in the distribution I. In some lattices the FM solution, which has been introduced above as the starting vector defining the set of orthogonal magnetic orbitals, may have a very high energy, due to the spatial vicinity of the magnetic centers. In such a case the energy minimization of the FM determinant may produce irrelevant magnetic orbitals, with exceedingly large delocalization on the ligands. The magnetic orbitals which are appropriate for the description of the lowenergy states keep a stronger metal character. A signature of this problem is observed when the atomic spin densities at the highly frustrated, high-energy unrestricted solutions no longer obey equation (21). Therefore, the obedience to this equation may be used as a criterion to retain the guess distribution I as a valid member of the set. In practice, when the FM solution becomes exceedingly high, it is convenient to express all energy differences with respect to the lowest AF solution.
The mapping procedure above described can be extended to provide information concerning the high-order terms of the spin Hamiltonian in equation (2) such as four-body cyclic operator amplitudes. In a rigorous QDPT development starting from the Hubbard Hamiltonian and using the set of VB neutral and singly ionic forms as model and outer spaces, respectively, one may find that these four-body spin operators appear at fourth-order and imply a cyclic circulation of the electrons around a plaquette. Such a circulation of the electrons is impossible in a plaquette with all spins parallel ( S z = ±2), but it generates off-diagonal and diagonal corrections when the total S z value for the plaquette is 0 or ±1. It has been shown [28,42] that for a plaquette with S z = ±1 the fourth-order diagonal corrections due to the cyclic circulation of electrons is It reflects processes of two types, involving the major spin electrons or the minor spin electrons through clockwise or anticlockwise movements; hence the factor 4 in equation (22). It is now possible to show that the variational symmetry-broken solutions Φ I associated to Φ I effectively incorporate these effects. To the second-order the orbital i may be expressed as with c (1) i j = t i j /U and as represented by the corresponding diagrams, where Φ I is taken as the vacuum state and the unrestricted self consistent solution involves the fourth-order diagrams which correspond to the minor spin circulation appearing in the last scheme of the plaquette. A similar identification can be performed for the major spin electron movements, (25) and the corresponding wavefunction diagrams and the corresponding fourth-order energy diagrams For a S z = 0 plaquette it has been shown from a QDPT expansion [28] that the fourth-order corrections due to the cyclic circulations of spins result in a diagonal energy shift equal to It is also possible to identify term by term the contributions such as which involves the circulation of the α spin electrons with a fourth-order correction included in the variational treatment leading to Φ I .
The factor two with respect to the S z = 1 case comes from the time-ordering degree of freedom between the second and third interactions and the sign change is due to the hole-hole interaction. For the spin alternant distribution ijkl in the plaquette the mechanisms involve doubly ionic intermediate states for ijkl as There are 46 of such processes, each of them contributing by t i j t jk t kl t li /2U 3 = g i jkl /2. The corresponding energy diagrams are of the type and they involve the components of the UHF function on the doubly ionic determinants of the form jjll as obtained at the second-order through the diagrams Hence, the cyclic contributions to the energies of the symmetry-broken single determinantal solutions are at the fourth-order at least, identical to the cyclic QDPT diagonal corrections Φ I |Ĥ eff(4) |Φ I . It has been shown by some of us [42,43,45] that these corrections may be expressed in terms of four-spin operators (equation 2) as and, hence, the four-body corrections can be identified as the diagonal terms of the biquadratic operator with amplitudes given by the average values of equation (28). In the following section examples of the applications of this general formula will be given and discussed in detail.

Computational details
The calculation of the energy of each relevant magnetic solution has been carried out at the experimental geometry 4 and using different exchange-correlation potentials as implemented in the CRYSTAL03 code [73]. These are the well known B3LYP hybrid method [74,75] and the Fock-35 scheme proposed by Moreira et al to reach a balanced description of many electronic structure properties of NiO [15], which gives satisfactory values of NN magnetic coupling constants in molecular magnets and solid state magnetic cuprates [76,77]. Crystalline orbitals are built as linear combinations of Bloch functions which in turn are built from atomic basis sets (AOs) optimized for the crystal environment. The AOs are contracted real spherical harmonic Gaussian type functions (GTFs). Extended all-electron basis sets have been used to describe the Cu and O atoms whereas effective core pseudopotentials have been used to represent inner electrons of the remaining ions. These standard all-electron and pseudopotential basis sets have been previously used by Moreira and Dovesi [58,78] and Su et al [79]. They correspond to the optimized basis sets for the corresponding ions and the outer isolated sp or sp and d exponents are re-optimized for a given environment. Also, the basis sets can be obtained from the CRYSTAL site (http://www.chimifm.unito.it/). Strict convergence criteria and a set of 105 points in the irreducible Brillouin zone have been used to ensure a numerical accuracy of 10 −7 Hartree per formula unit.

Results
In this section, we will discuss the application of the methodology and general procedure described above to the case of 2D square lattices which are common to many cuprates which may become superconducting under doping and to the case of ladder compounds which have a more complex structure.

Spin Hamiltonian parameters for the 2D square lattices
For the cuprate planar square lattices which may be seen as constructed of adjacent plaquettes, three interactions are expected to be important. These are J between NN, J d between NNN through the main diagonal of a plaquette, and, finally, J ring which is the four-body ring interaction in the plaquette. The amplitudes of J, J d and J ring have been estimated from the energy differences corresponding to the magnetic solutions schematically shown in figure 1. One of these solutions corresponds to a FM alignment of all magnetic moments whereas the other three are symmetrybroken AFM arrangements. Table 1 reports the relationship between the above-described effective parameters and the calculated energy for UHF and each of the two DF methods. These energy differences univocally determine the amplitudes of the three parameters since, in this case, the corresponding set of equations is linearly independent. The calculated values of J, J d and J ring for 2D cuprates are reported in table 2. Concerning the leading NN spin coupling amplitudes (J) it is worth noting that UHF largely underestimates them, as found in many previous works [57,76,80]. Density functional calculated values generally represent an improvement over UHF results although they exhibit the typical strong dependence on the exchange-correlation functional already pointed out by Martin and Illas [55,56]. Hence, B3LYP overestimates the J values by a factor close to 1.5 whereas the values predicted by the Fock-35 potential are close to experiment as also expected from previous work on similar systems [15,76,80]. Therefore, one can take the Fock-35 values as a reliable prediction. In all cases, the NNN diagonal interaction is predicted to be small and of AFM character, in contrast to what is suggested by some fits to experimental data [81]. However, the most important result of the present work concerns the cyclic exchange J ring amplitude. The general procedure described in detail in section 2 has permitted the first direct estimate of the amplitude of the J ring terms from ab initio periodic calculations [82]. For La 2 CuO 4 , the outcome of the present periodic approach using the Fock-35 parameterization is in agreement with the other only available theoretical estimates for this compound arising from cluster calculations performed by either explicitly correlated wavefunctions [42,43] or DFT based calculations [45,76] and experimental estimations of J [22,32]. In addition, the present Fock-35 estimate of J ring ∼ 35.8 meV is in excellent agreement with the indirect evaluations of Coldea et al (J ring ∼ 38 ± 8 meV) [22] and the simulations of Mizuno et al (J ring ∼ 40 meV) [31]. Notice, however, that the J ring /J = 0.25 ratio predicted here is smaller than that reported by Toader et al [30], which is of 0.5, but in excellent agreement with a more generally accepted ratio J ring /J = 0.3 [21,22,31] for this compound. Finally, we would like to point out that the experimental determination of this term is still under discussion and even the existence and importance of four-body terms in La 2 CuO 4 have raised some controversy [83]. The present work provides unbiased first principles results that fully  [84]. It is also worth pointing out that the Fock-35 calculations predicts an AFM J d ∼ 8.8 meV again consistent with the values provided by embedded cluster calculations [42,44,77]. Hence, it will be of great interest to repeat the fit in Toader et al experiments [30] by using the present estimate of both J and of J d or, in a bottom up approach, use the present values for the three effective parameter to check consistency with respect to experiment. For Ca 2 CuO 2 Cl 2 , Sr 2 CuO 2 F 2 and Sr 2 CuO 2 Cl 2 , there are not previous theoretical or experimental values for the amplitude of the ring exchange term. Present calculations predict values for Ca 2 CuO 2 Cl 2 and Sr 2 CuO 2 Cl 2 that are slightly smaller than the corresponding ones in La 2 CuO 4 . This is consistent with larger NN and NNN distances in Ca 2 CuO 2 Cl 2 and Sr 2 CuO 2 Cl 2 compared to La 2 CuO 4 (3.87 and 3.97 versus 3.81 Å, for NN and ∼ 5.47 and ∼ 5.62 versus 5.38 Å for NNN distances). For Sr 2 CuO 2 Cl 2 the agreement between the estimate of the NN interaction (J ∼ 130 meV) and the experimental value (J = 125 ± 6 meV [85,86]) is excellent, the J d coupling is also predicted to be AFM and the estimated J ring /J ratio ∼ 0.20 is somewhat smaller than that for La 2 CuO 4 but still in a realistic range. For Sr 2 CuO 2 F 2 the results are closer to La 2 CuO 4 ones despite the differences in crystal structure (3.86 versus 3.81 Å, for NN and ∼ 5.45 and ∼ 5.62 versus 5.38 Å for NNN distances).
Before closing this section it is interesting to note that the magnitude of J alone does not discriminate Sr 2 CuO 2 F 2 and Ca 2 CuO 2 Cl 2 from La 2 CuO 4 (stoichiometric Sr 2 CuO 2 Cl 2 corresponds to a extremely stable structure and attempts to synthesize doped phases by cation substitution of interstitial anion excess have been unsuccessful). Taking the reliable Fock- We would like to point out that this is also in line with previous results indicating that a linear relationship exists between J/t and T c [41,76]. This is because, a correct value of J/t also indicates a correct value of t/U and hence of J ring /J (∼ (t/U ) 2 , cf equation (22) and [45]).

Spin Hamiltonian parameters for the SrCu 2 O 3 spin ladder compound
For the SrCu 2 O 3 spin ladder, ten different magnetic solutions have been obtained from the spin distributions shown in figure 2. Due to an excessive spin frustration between the ladders, the FM solution becomes exceedingly high in energy and indeed physically meaningless. The shortrange intra-ladder repulsion forces the unpaired electron to be excessively delocalized on the oxygen atoms. Therefore, the AFM gs solution is used as energy origin. For the remaining nine broken-symmetry solutions with different magnetic orders, we derived a set of nine equations which are shown in table 3. These equations have been derived assuming that only five magnetic interactions may have non-negligible amplitudes: these are NN interactions along the legs J l , or the rungs J r , the NNN J d interactions in the plaquettes, the interaction between legs of different ladders J i , and the four-body operator J ring amplitude. At variance from the abovediscussed cuprates, the resulting equations are linearly dependent, the system of equations is overdetermined and a least-square procedure has been used to find the optimum set of values for the five magnetic interactions considered. The optimum values are reported in table 4. The consistency of the set of equations is almost complete since the standard deviation between the DF calculated energies and those obtained from the equations in table 3 and the so-obtained set of magnetic interactions is of 1.5 meV for the B3LYP results and even smaller for the Fock-35 set. This consistency indicates that all non-negligible interactions have been included in the spin model Hamiltonian.
For all the methods the dominant interactions are the couplings across the rungs (J r ) and along the legs (J l ) which all methods predict to be nearly equal in magnitude-i.e. (J r /J l ) ≈ 1as expected from the similarity of the Cu-O-Cu exchange paths in agreement with previous two magnetic sites cluster model calculations [36] although a moderate cluster size dependence on the magnetic interaction along the leg has been observed [77]. Nevertheless, enlarging the cluster model by introducing a third copper atom belonging to the very close neighbor ladder results again in J r /J l ≈ 1 [71]. In any case, this dependence evidences the difficulty of embedding finite clusters for certain types of materials and supports the present periodic approach that does not depend on the cluster design used to represent the real material.
The accuracy of the present calculations can be judged from the amplitude predicted for J r through the most realistic Fock-35 potential (∼ 155 meV) which is consistent with the recent Raman response experiments [35]. The estimated J d value is AFM as in the 2D cuprates and a non-negligible interladder FM exchange (J i ≈ −0.22J l ) is found and, consequently, one should not consider this system as formed by non-interacting ladders. Regarding J ring , both its amplitude and the J ring /J r ratio are larger than for 2D cuprates (table 2), but still close to 0.3 in agreement with various sets of experimental informations [21]- [24]. The agreement is particularly good with the values reported by Schmidt et al [35] although these authors do not consider J i . This neglect is probably the reason for the strong anisotropy in their fitting which results in a J l /J r = 1.5 ratio, a value which is difficult to accept from the Cu-O-Cu distances. Such large anisotropy between J l and J r is not supported from the present periodic calculations. We also found that repeating the least square fitting of the energies in table 2 but neglecting the terms involving J ring results in a larger anisotropy in the thus J r and J l calculated values. These results strongly suggest revising previous fittings using the present estimates as starting points.  (12)). Previous works have shown that t values are not very sensitive to the exchange potential [57,77]. This suggests that the exceedingly large J value predicted by the B3LYP method effectively implies a too small U value and hence an insufficient on-site two-electron repulsion, as discussed elsewhere [57,80]. This incorrect description of electron-electron correlations in the B3LYP functional has a much more dramatic consequence. In fact, it has been shown by Malrieu and Maynau [28] that J ring scales as t 4 /U 3 (cf equation (22)) and, hence, a too small U value will result in an even larger overestimation of J ring as clearly shown in tables 2 and 4. This interpretation is supported by the results obtained from the UHF method. In this case the values   (21), is also shown.
of J are largely underestimated due to the lack of electronic correlation effects, a well-known effect [39,55,56,90,91]. Therefore, one may conclude that UHF largely overestimates the on site two-electron repulsion and, consequently, J ring obtained at this level of theory must be very small. This is indeed the case, the UHF estimate of J ring is vanishingly small. These results are manifestations of the sensitivity of J ring /J with respect to the |t/U | ratio; this is J ring /J scaling as |t/U | 2 . If J (calc)/J (exp) = λ, then J ring (calc)/J ring (exp) ≈ λ 3 . To conclude, the DFT based approaches effectively include electronic correlation effects that lead to a decrease of the effective U parameter which was severely overestimated by the mean-field Hartree-Fock method. However, the B3LYP potential appears to lead to unphysical too low U values with important consequences in the resulting physical description. Scrutiny of the atomic spin densities, defined according to Mulliken population analysis, brings consistent additional information. In order to check the validity of equation (21), we report in figure 3 the dependence of the spin densities on the metal centers as a function of the number of opposed spins on NN sites. Notice that this information is extracted from the various magnetic 'phases' for each compound. As suggested by equation (20), a good linear dependence of the local spin density on the number of spin alternations is observed for each of the Hamiltonians. This result evidences the different strengths of the electronic delocalization in these three Hamiltonians and the slope of the correlation can be used to estimate the t/U ratio. For La 2 CuO 4 , the t/U estimated ratios are 0.044, 0.093 and 0.110 for UHF, Fock-35 and B3LYP, respectively. These results nicely compare with previous estimates of this ratio from finite cluster calculations: 0.080 and 0.120 for Fock-35 and B3LYP [57]. Also, these results confirm the general trend discussed above, namely an excessive localization in UHF, resulting in a too large unscreened effective U, and a slightly too strong delocalization in B3LYP with a concomitant underestimation of the effective U. In the ladder compound, the analysis of the correlation with respect to the number of spin alternations is not univocal since the bonds are not identical and, therefore, is not carried out.
As mentioned in section 2, it is observed that some highly frustrated and high energy selfconsistent solutions provide local spin densities which deviate significantly from the behaviour expected from equation (21). This is due to the physically meaningless strong metal-to-ligand delocalization of the magnetic electrons. Hence, the spin density analysis provides a criterion to eliminate these spurious solutions. This problem has been essentially observed in the ladder where the J values are larger as a result of the shorter Cu-Cu distances and, also, due to the FM interladder interaction. As a consequence, the total spin frustration becomes exceedingly costly. Accordingly, we strongly recommend using a sufficiently large set of low energy solutions in the fitting procedure.

Summary, conclusions and possible extensions
The extraction of NN spin couplings in magnetic lattices from symmetry-broken mean-field calculations has been frequently performed on molecular systems and solids (see [48] and references therein). The present work uses the benefit of the large multiplicity of symmetrybroken periodic solutions and proposes a new, general and unbiased scheme to predict the amplitude of the parameters defining a general spin Hamiltonian from DFT periodic calculations. An important property of the present procedure is that one does not need to make any assumption on the relative amplitude of these terms. Instead, it relies on a mapping approach between the energy of pertinent magnetic solutions and the diagonal terms of the spin Hamiltonian in a local representation. In addition, the overdetermined set of equations provides a test for the completeness of the interactions considered in the spin model. It has also been shown that the spin densities of the symmetry-broken solutions may be rationalized and they furnish additional information on the t/U ratio through equation (20). Hence, the analysis of the spin density distributions of these various solutions provides useful and consistent information on the t/U (or J/t) ratio and, hence, on the electronic delocalization. As a corollary one can point out that the observed correlation between T c at optimal doping and the J/t ratio [41,76] can be formulated as well as a correlation with t/U ratio that brings in a more physically intuitive description. Finally, this work also points out the extreme sensitivity of the four-body operators amplitudes to the choice of the density functional and especially to the ratio of Fock exchange. For the time being the use of a 35% ratio has led to consistent results for cuprates and other transition metal compounds. The question of the transferability of this parametric quantity to other magnetic ions remains open.
A very important outcome of the present work concerns the prediction of the four-body cyclic operators in the general spin Hamiltonian in equation (2). The application of the general procedure discussed at length in section 2 to a series of compounds provides an independent confirmation of the importance of four-body terms in these materials. In particular, for La 2 CuO 4 the present results provide further support to the conclusions by Toader et al [30] based on the analysis of neutron diffraction experiments. Moreover, we supply reliable values for spin Hamiltonian parameters of other key compounds such as Ca 2 CuO 2 Cl 2 , Sr 2 CuO 2 F 2 and Sr 2 CuO 2 Cl 2 . These permit us to conclude that the importance of four-body terms is likely to be similar for most of the HTCS related cuprates and even more important in ladder compounds, in agreement with the recent experiments of Schmidt et al [35]. Likewise, it is also suggested that the fit to neutron scattering data should be revised by considering alternative values for both J and J d magnetic coupling terms. In addition, the present study provides further evidence that the four-body term and the FM inter-ladder exchange introducing spin frustration between legs in the SrCu 2 O 3 ladder should not be neglected.
To conclude this study, it is worth mentioning some possible generalizations of the procedure presented in section 2. One may, for instance, wonder whether the six-body spin operators would have significant amplitudes in such lattices. The cyclic circulation of electrons in six-member rings-related to the so-called chemical aromaticity-is responsible for the appearance of such operators at the sixth-order of perturbation theory. It has been demonstrated [28] that while the four-body operator results in a coupling = 40 t 4 /U 3 Ĥ the six-body matrix element resulting in a rotation of six spins is given by = 504 t 6 /U . 5

Ĥ
The large prefactor is due to the huge number of processes leading to a full permutation of spins in a six-membered ring. Now, notice that a set of two fused plaquettes defines such a sixmembered ring. It would be worth checking whether, despite the smallness of the t/U ratio, such many-body operators are negligible in cuprate lattices. From the perturbative arguments above, and in view of the calculated t/U ratio, one expects the ratio J six-ring /J four-ring ∼ 0.1. Consequently, this question has not been addressed in the present study.
Yet another generalization would concern lattices with S > 1/2 magnetic sites such as Ni(d 8 ) ions, NiO being the paradigm of these materials. The same strategy is applicable to evaluate not only the NN interactions but four-body exchange amplitudes as well. This will be considered in further work which will also show that one may, in principle, use symmetry-broken solutions to evaluate the possible importance of biquadratic spin exchange terms.