{'fraction_non_alphanumeric': 0.04347481053341622, 'fraction_numerical': 0.015077968971127013, 'mean_word_length': 4.19928483014716, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 10, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Low energy theory of the t-t'-t''-U Hubbard Model at half-filling: interaction strengths in cuprate superconductors and an effective spin-only description of La_2CuO_4

Spin-only descriptions of the half-filled one-band Hubbard model are relevant for a wide range of Mott insulators. In addition to the usual Heisenberg exchange, many new types of interactions, including ring exchange, appear in the effective Hamiltonian in the intermediate coupling regime. In order to improve on the quantitative description of magnetic excitations in the insulating antiferromagnetic phase of copper-oxide (cuprate) materials, and to be consistent with band structure calculations and photoemission experiments on these systems, we include second and third neighbor hopping parameters, t' and t'', into the Hubbard Hamiltonian. A unitary transformation method is used to find systematically the effective Hamiltonian and any operator in the spin-only representation. The results include all closed, four hop electronic pathways in the canonical transformation. The method generates many ring exchange terms that play an important role in the comparison with experiments on La2CuO4. Performing a spin wave analysis, we calculate the magnon dispersion as a function of U,t,t' and t''. The four parameters are estimated by fitting the magnon dispersion to the experimental results of Coldea et al. [Phys. Rev. Lett. 86, 5377, {2001}] for La2CuO4. The ring exchange terms are found essential, in particular to determine the relative sign of t' and t'', with the values found in good agreement with independent theoretical and experimental estimates for other members of the cuprate family. The zero temperature sublattice magnetization is calculated using these parameters and also found to be in good agreement with the experimental value estimated by Lee et al. [Phys. Rev. B 60, 3643 (1999)]. We find a value of the interaction strength U~8t consistent with Mott insulating behavior.

Spin-only descriptions of the half-filled one-band Hubbard model are relevant for a wide range of Mott insulators. In addition to the usual Heisenberg exchange, many new types of interactions, including ring exchange, appear in the effective Hamiltonian in the intermediate coupling regime. In order to improve on the quantitative description of magnetic excitations in the insulating antiferromagnetic phase of copper-oxide (cuprate) materials, and to be consistent with band structure calculations and photoemission experiments on these systems, we include second and third neighbor hopping parameters, t ′ and t ′′ , into the Hubbard Hamiltonian. A unitary transformation method is used to find systematically the effective Hamiltonian and any operator in the spin-only representation. The results include all closed, four hop electronic pathways in the canonical transformation. The method generates many ring exchange terms that play an important role in the comparison with experiments on La2CuO4. Performing a spin wave analysis, we calculate the magnon dispersion as a function of U, t, t ′ and t ′′ . The four parameters are estimated by fitting the magnon dispersion to the experimental results of Coldea et al. [Phys. Rev. Lett. 86, 5377, 2001] for La2CuO4. The ring exchange terms are found essential, in particular to determine the relative sign of t' and t ′′ , with the values found in good agreement with independent theoretical and experimental estimates for other members of the cuprate family. The zero temperature sublattice magnetization is calculated using these parameters and also found to be in good agreement with the experimental value estimated by Lee et al. [Phys. Rev. B 60, 3643 (1999)]. We find a value of the interaction strength U ≃ 8t consistent with Mott insulating behavior.

I. INTRODUCTION
High temperature superconductors have challenged almost every traditional concept and method of condensed matter theory. Ill understood issues concerning Fermi liquids and quantum critical behavior for example, may be intimately related to the problem of superconductivity itself. From the very beginning, it has therefore appeared important to develop a better understanding of quantum magnetism in general, given that the parent compounds of high-temperature superconductors are insulating antiferromagnets. For example, although these compounds are Néel ordered, it was proposed early on that, in two dimensions, the ground state of the Heisenberg antiferromagnet might be very close to being a resonating-valence bond state rather than a Néel ordered state 1 . This turns out not to be true, but it took some time to establish numerically that the ground state is indeed ordered 2,3 , in agreement with experiment 4,5 .
Issues relating to the quantum magnetism of the copper-oxide materials (cuprates) have resurfaced in the last few years. Detailed neutron scattering experiments on the parent high-temperature superconductor La 2 CuO 4 illustrate that the question of which model best describes these compounds 6 can now be addressed with precision. For example, while the magnetic properties of Copper Deuteroformate Tetradeuterate can be described in detail using the Heisenberg model 7 , La 2 CuO 4 exhibits clear deviations 6 from it.
Specifically, the observed spin wave dispersion along the antiferromagnetic Brillouin zone boundary is not predicted by the Heisenberg model. One can in principle account for these differences by introducing phenomenological two-spin and multi-spin exchange constants into the spin Hamiltonian. These spin Hamiltonians can be deduced, as shown by Dirac 8,9,10 , from general considerations of symmetry and permutation operators and the value of the exchange constants can be fixed by fitting to the experimental data. However, in our opinion, a more interesting approach is that taken by Coldea et al. 6 . The authors of Ref. [6] take the point of view, which we share, that the one-band Hubbard model is a more fundamental starting point for the description of the copper oxygen planes in these materials if one wants to connect the magnetism at the microscopic level with electronic correlation effects. Although in this paper we concentrate on half filling, it may also be valid on moving away from half-filling, into the region where the fermions eventually pair to give high-temperature superconductivity.
The one band Hubbard Hamiltonian H H is one of the simplest lattice models of interacting electron that admits a "rich" phenomenology resembling that of the cuprates. A key feature of this model is that it does describe both the insulating and metallic phases and the (Mott) transition between them. In interaction-induced insulators, the band is half-filled and insulating behavior occurs because of interactions, not because of band folding induced by long-range order. However, in twodimensions a half-filled band can be insulating because of a pseudogap induced by fluctuating antiferromagnetism, or because of the Mott phenomenon. In the former case (weak coupling), the pseudogap appears when the thermal de Broglie wavelength ξ th is less than the antiferromagnetic correlation length 11 ξ, while in the latter Mott insulating case (strong coupling), the system is insulating even when this condition is not satisfied. In La 2 CuO 4 neutron measurements have been done only in the regime where ξ > ξ th , so that antiferromagnetic correlations could be the source of the insulating behavior. In a very recent paper 12 Comanac et al. used singlesite dynamical mean-field theory (DMFT) to calculate the optical conductivity of La 2 CuO 4 . The authors of Ref. [12] found that La 2 CuO 4 , or at least its optical conductivity, is best parametrized by a U/t smaller than that necessary for the Mott transition in single-site DMFT. Consequently, the authors of Ref. [12] argue that correlations in the cuprates may not be as strong as generally believed and that these materials may not be Mott insulators but rather, that the development of antiferromagnetism is necessary to drive the the insulating state. On the other hand, cluster-DMFT calculations, which are more accurate in two dimensions, find that the inclusion of short-range antiferromagnetic correlations, which is possible in multi-site calculations, reduce considerably the critical (U/t) Mott necessary for the Mott transition compared with single-site DMFT 13,14,15,16 .
This leads to one of the key questions in this field. Do the strange properties observed in underdoped cuprates at finite temerature arise from interaction-induced localisation (Mott insulator) or from competing phases? If interactions are not strong enough to lead to a Mott insulator, competing phases arising at weak coupling are the only remaining possibility.
To find parameters appropriate for high-temperature superconductors, Coldea et al. used earlier theoretical results 17,18 relating the parameters of the Hubbard model, hopping energy t and on-site interaction U , to the exchange constants in an effective spin-Hamiltonian. The parameters in the Hubbard model were deduced by fitting the resulting effective low-energy spin theory to the spin-wave dispersion relation. This is quite a remarkable result. Indeed, had the Heisenberg limit accurately described the Hubbard model, only the value of t 2 /U could have been deduced from neutron experiments. The appearance of corrections to the Heisenberg model to higher order in powers of t/U allows one to obtain t and U separately, the former being an effective electronic band quantity.
However, to be consistent with a general message provided by Angle Resolved Photoemission Spectroscopy (ARPES) 19,20,21,22 experiments, optical experiments 23 and with results from band structure calculations 24,25 on FIG. 1: (Color online) Hopping processes and resulting spin interactions. a) shows the different hopping processes characterized by parameters t, t ′ and t ′′ in the Hubbard model considered in this paper. t, t ′ and t ′′ are hopping parameters between first, second and third nearest neighbors, respectively. b) At half-filling, and when t ′ = t ′′ = 0, canonical perturbation theory leads to order t 4 /U 3 to an effective spin-1/2 Hamiltonian, H (4) s , characterized by first (J1), second (J2) and third (J3) nearest neighbor exchange interactions as well as a 4-spin ring (cyclic) exchange interaction with coupling strength Jc. c) In a large S expansion, H (4) s can be recast to order 1/S as an effective spin Hamiltonian which only involves bilinear (pairwise) spin-spin exchange interactions between first (J eff 1 ), second (J eff 2 ) and third (J eff 3 ) nearest neighbors. To order 1/S, the effect of the ring exchange of strength Jc merely renormalizes the first and second nearest-neighbor exchange as discussed in Ref. [6]. d) illustrates an example of a four hops (ring exchange) electronic process that involves t ′ and t ′′ and contribute to order t 2 t ′ t ′′ /U 3 to the spin Hamiltonian H (4) s . e) shows that the ring exchange term illustrated in d) introduces in the 1/S approximation of H exchanges beyond those shown in c) and which are generated by all the hopping processes involving an allowed combination of t, t ′ and t ′′ to order 1/U 3 .
a wide variety of quasi two-dimensional cuprate materials, it is necessary to include in the Hubbard model both second (t ′ ) and third (t ′′ ) nearest-neighbor hopping, as they have been found to have sizable values in all the above studies (see Fig. 1). In this paper we thus address the following questions: What is the effect of these additional hopping constants on the magnetic properties of the cuprates, and on La 2 CuO 4 in particular? Can their "effective" values and that of U within the Néel or-der phase be determined by comparing experimental spin wave dispersion data with theoretical calculations? As a first step to address the question of the strength of t, t ′ , t ′′ and U in La 2 CuO 4 , we incorporate the effects of these extra hopping constants into a one-band Hubbard model from which we derive an effective spin-only Hamiltonian description of La 2 CuO 4 . As we shall see, the effects of t ′ and t ′′ are rather subtle and can even sometimes compensate each other. It is only when these hopping terms have been included in the derivation of the effective spin Hamiltonian, with all the subsequent new ring exchanges, that one can extract more meaningful values of t − t ′ − t ′′ − U from comparisons of theory with inelastic neutron scattering experiments. In particular, the values of t and U obtained by Coldea et al. 6 place the underlying Hubbard model in the intermediate coupling regime where the band width and interaction energies are comparable and hence in the regime where one might expect an insulator to metal transition. Our more detailed analysis leads to a larger values of U , repositioning material relevant model parameters more clearly within the Mott insulating regime, as we will discuss more thouroughly at the end of this paper.
Although we analyze experiments on La 2 CuO 4 , providing a specific materials context to our study, we believe that the parameters of the one-band Hubbard model that we extract below should be characteristic of the twodimensional copper-oxygen planes of cuprate superconductors in general: specifically, band structure calculations do reveal variations of band parameters from one compound to the next 25,26 , but overall, the band structure of CuO 2 planes is quite similar from one compound to the next 27 . Even away from half filling, whether for hole or electron doped materials, the one-band Hubbard model for CuO 2 planes seems to contain much of the Physics of high-temperature superconductors 28,29 . Our calculations should therefore have wider ranging applications than the detailed comparison with experiments on pure La 2 CuO 4 presented below.
In the next two subsections of this Introduction, we describe in more detail the central issues and key results of this paper. There are two main parts to our work. In the first part, Section II, we derive the effective Hamiltonian to order 1/U 3 , including t − t ′ − t ′′ . This derivation is exact and completely general. In the second part, Section III, we determine the spin wave energies, taking into account the quantum correction to the classical frequencies through the spin wave velocity renormalization factor, Z c (k), to lowest order in 1/S. We then determine a parameter set {t, t ′ , t ′′ , U } by fitting the experimental magnon excitation spectrum of Ref. [6]. The resulting set is in good agreement with values from ARPES experiments and band structure calculations on various cuprates. Further, we calculate the zero temperature staggered magnetization order parameter using the set {t, t ′ , t ′′ , U } and compare it with experimental values. A more complete solution of the spin dynamics for the effective spin Hamiltonian derived herein, including magnon-magnon interactions, would most likely require the use of numerical techniques 30 , since it would be an arduous task to include in the self-energy corrections the multitude of magnon creation and annihilation terms present beyond the non-interacting approximation.
A. Spin-only description of Hubbard model and La2CuO4 In many strongly correlated quantum mechanical systems, the separation of energy scales allows the development of a low energy effective theory, through the integration over the high energy degrees of freedom. In this context, spin-only descriptions of strongly-correlated electron systems are an excellent example. Here, in the limit of strong electron-electron interaction, elimination of the states with one or more doubly occupied sites reduces exponentially the dimension of the relevant low energy Hilbert space, facilitating the theoretical and numerical description of the problem. Through this procedure, the regime of particular interest, where energy scales actually compete, can be approached perturbatively, via the introduction of a small parameter, the ratio of the kinetic to potential energy scales. However, the price one pays for this dimensional reduction is the generation of longer range and many particle interactions, as one moves into this intermediate regime. The derivation of such a theory can be achieved using different methods, leading to distinct expressions for the effective Hamiltonian, which can be shown to be equivalent through the application of a unitary transformation 31 . One of these methods, the canonical transformation (CT) 18,32,33 , which we adopt below, is a convenient and systematic way of expanding a model like the Hubbard model to any order in the perturbation parameter.
In the regime where t/U remains a small parameter suitable for a perturbation theory, but where one moves away from the strictly Heisenberg limit (t/U → 0), the pairwise exchange interactions are joined by ring, or cyclic exchange terms that couple more than 2 spins 18 . Large enough ring exchange terms have been found to drive some model systems into exotic, intrinsically quantum mechanical ground states 34,35,36 . Indeed, there is currently a rapidly growing interest in the study of effective many-body non pairwise spin−spin interactions which lead to non-trivially correlated states in quantum spin systems 37,38 .
In most theoretical work, the effort has so far focused on models with only the nearest-neighbor Heisenberg exchange and ring exchange. However, when starting from a microscopic fermionic Hamiltonian the strength of the ring exchange coupling is determined from the small parameter t/U and is not a free parameter, nor is it the sole higher order term arising in the resulting spin only theory. To the same order in perturbation, t(t/U ) 3 , the canonical transformation also generates second, J 2 , and third, J 3 , neighbor spin-spin interactions through pro-cesses involving four electronic hops. The ring exchange, J c , is thus only one of a set of spin interactions generated to this order 18,31,33 and all should be taken into account. We remark that in the general approach of Dirac 8,9 , the coupling constants multiply spin permutation operators, rather than spin operators themselves and thus have a different definition from those here. One finds that permutation of four spins lead to both four spin ring exchange terms and pure two-spin exchange when written in terms of spin operators 10 .
Within the t − U Hubbard model, the inclusion of virtual hopping pathways of n hops thus introduces further neighbor spin couplings of order t(t/U ) n−1 . For a compound with a band structure characterized by several tight-binding parameters of comparable value, one should consider the role of direct further neighbor hops over length scales comparable with the nearest neighbor pathways of length n, as well as the multiple hop terms. Hence, to develop a theory accurate to order (1/U ) 3 one should include direct hopping to first, second and third nearest neighbors, characterized by energy scales t, t ′ and t ′′ , respectively. These hopping processes generate second and third neighbor Heisenberg exchange terms that compete with those generated by the four hop processes discussed above. There is evidence that this is the situation for the CuO 2 planes of parent superconducting materials: for example, in Sr 2 CuO 2 Cl 2 , which has a half filled band, the hopping constants t ′ and t ′′ are estimated to be t ′ /t ∼ −0.3, t ′′ /t ∼ 0.2, from comparisons of exact diagonalization with photoemission experiments 39 . Very similar values are found from comparisons with other types of calculations 40 and from band structure calculations for YBa 2 CuO 3 O 7 41 . In fact these values are typical for insulating materials with copper-oxygen planes, including both hole and electron doped cuprates 27,42 . For La 2 CuO 4 , band structure calculations 24,25 suggest smaller values of t ′ /t and t ′′ /t. ARPES experiments are not available for this compound but they are for the doped system La 2−x Sr x CuO 4 that becomes a high temperature superconductor. According to these experiments 19,20,21,22 , t ′ /t and t ′′ /t deduced from the shape of the Fermi surface are close to band structure values, although they are larger by about 30% for very lightly doped compound 19 . We will return to this issue later.
Given this range of values for t ′ and t ′′ one could expect the resulting exchange terms of order t ′ (t ′ /U ) or t ′′ (t ′′ /U ) to be of similar magnitude to those of order t(t/U ) 3 . This is one of the motivations for the work presented in this paper where we introduce all terms in a spin-only Hamiltonian that depend on t, t ′ and t ′′ and that are generated up to order 1/U 3 . We assess their importance, making particular reference to the magnetic properties of La 2 CuO 4 . We show that the canonical transformation method 18,31,33 can be easily generalized for the t − t ′ − t ′′ − U Hubbard model, giving a complete spin-only description for arbitrary first, second and third neighbor hopping, correct to order 1/U 3 . B. Spin-wave analysis and comparison with experiment.
It has been known for some time that a spin wave analysis of the spin 1/2 Heisenberg antiferromagnet on a square lattice, taken to leading order in 1/S, reproduces the zero point quantum fluctuations of the staggered moment to a good approximation 43 . This is because the second order terms renormalize the classical magnon frequencies but do not induce magnon-magnon interactions and make zero contribution to the magnetization. In this paper, we therefore take it as a reasonable first approximation to calculate quantum spin fluctuations to lowest order in 1/S and include the first order corrections to the classical frequencies. We have not addressed the question as to whether other physics, in particular magnon-magnon interaction is introduced to second order in 1/S. We compare our results with the inelastic neutron scattering data of Coldea et al. 6 and with the spin wave calculation therein where only nearest neighbor hopping t was considered. We find that including further neighbor hops allows for a more convincing description of the neutron data, with the set of parameters (t, t ′ , t ′′ , U ) that best describe the data being comparable to other estimates for cuprate materials. Our analysis also gives a new estimate for the total staggered magnetic moment which is in good quantitative agreement with the experimental estimate 44 .
We now describe in more detail the main features of the data that must be explained and how going beyond the nearest-neighbor hopping provides a better model. The inelastic neutron scattering data of Coldea et al. 6 provide an extensive description of the magnon dispersion in La 2 CuO 4 over the whole of the 1 st Brillouin zone. The analysis of Ref. [6] clearly shows that the nearest neighbor spin S = 1/2 Heisenberg Hamiltonian does not reproduce the magnon dispersion over the whole zone. Most notably, while a spin wave analysis of the Heisenberg model to order 1/S gives a flat magnon dispersion over the interval q BZ ∈ [(π, 0), (π/2, π/2)] (along the antiferromagnetic zone boundary), the experimental results show negative (downward) dispersion over this region (see Fig. 3). Going beyond 1/S for the Heisenberg model does not explain this negative dispersion feature 45,46 . Rather, introducing interactions between magnons gives rise to a small positive (downward) dispersion over this q BZ interval 45,46 . Including terms up to t(t/U ) 3 , which generates second and third neighbor and ring exchange interactions, does however provide the characteristic downward dispersion, and this is what allowed Coldea et al. to fit the experimental data within a 1/S spin wave analysis 47,48 . Recent analysis of triplon excitations in the copper oxide based ladder material La 4 Sr 10 Cu 24 O 41 also suggest that ring exchange interactions are present and are of a similar amplitude 49 to those found in Ref. [6] for It is important for the rest of the paper to expand briefly on this negative dispersion and to give an interpretation for it: the spin-only Hamiltonian, H s , obtained from the Hubbard Hamiltonian, H H , to order t 4 /U 3 , contains first, second and third nearest-neighbor exchanges, J 1 , J 2 and J 3 , respectively, as well as a fourspin ring exchange J c (see Eq. (16)). Within the Holstein-Primakoff spin-wave approximation, all terms in H (4) s , including the ring exchange, contribute to order 1/S quadratic magnon creation/annihilation terms to the resulting quadratic magnon Hamiltonian H and where the convention J eff n > 0 means antiferromagnetic exchange 51 . We note from Eq. 1 that J c leads to a "renormalization" of J 2 . In particular, for the t−U Hubbard model, J eff 2 = −16t 4 /U 3 (S = 1/2) is negative for all t/U . The negative value of J eff 2 favors ferromagnetic correlations across the diagonal of a plaquette and gives a downward dispersion along q BZ ∈ [(π, 0), (π/2, π/2)] 6 , as required. In the limit t/U → 0, one recovers the nearest-neighbor Heisenberg S = 1/2 Hamiltonian and the dispersion becomes flat along the zone boundary 52 .
Having provided an explanation for the origin of the negative (downward) dispersion along q BZ , we can now discuss the effect on the spin wave dispersion of including second nearest-neighbor hopping t ′ . When passing from the Hubbard model, including t ′ , to a spin only model, the leading order effect of t ′ is to generate an antiferromagnetic exchange between second-nearest neighbors that modify J 2 in H (4) s : which in turn modifies J eff 2 in the 1/S spin-wave approximation in Eq. (1): So, while for t ′ = 0 the dispersion is downward for all t/U , a nonzero t ′ competes with this trend since the second nearest-neighbor exchange generated to order 4(t ′ ) 2 /U , is antiferromagnetic and frustrating. This means that a sizable increase in t/U is required to achieve a good fit to the data of Ref. [6] when t ′ is included in the spin-only Hamiltonian to this order. The same observation has been made from a random phase approximation calculations 53,54 and also from quantum Monte Carlo calculations 55 , both on the half-filled one band Hubbard model. This evolution puts the best fit value for t/U into the intermediate coupling regime and perhaps very close, or maybe even beyond the critical value for a metal insulator transition. If we estimate this critical value from the value of U at which a finite gap in the density of states persists at finite temperature, then Fig. 5 of Ref. [56] suggests that for t ′ = t ′′ = 0 on the square lattice, the critical value is t/U ∼ 1/6 = 0.166. In the variational cluster approximation, which overestimates the effect of interactions, the critical value of t/U is larger but studies as a function of |t ′ | suggest that frustration leads to a decrease of this value 57 . The same trend as a function of frustration (t ′ ) occurs on the anisotropic triangular lattice 58 . Hence, with nonzero t ′ included in the model and presumably at play in the real material, it is ultimately important that successful fits to experiments on parent insulating compounds lead to smaller values of t/U than those where t ′ = 0, contrary to what is found in the references cited above where the third nearest neighbor hopping t ′′ is equal to zero. As discussed above, the third neighbor hopping term t ′′ is estimated to be of the same order as t ′ in a number of cuprate materials and should therefore also be taken into account 25 . However, the above discussion on the leading effect of t ′ might leave one wondering whether a perturbative spin-only Hamiltonian, starting from a half-filled one-band Hubbard model can provide a quantitative microscopic description of La 2 CuO 4 at all. It is this very question that has motivated us in pursuing the work reported here, and it is why we have derived the spin-only Hamiltonian H s from a t− t ′ − t ′′ − U Hubbard model, including all ring exchanges and hopping processes to order (1/U ) 3 . This procedure generates a large number of further neighbor and ring exchange paths, a small number of which make significant contributions. The main corrections come from terms of order (t ′ ) 2 /U and (t ′′ ) 2 /U but we also find that ring exchanges of order t ′ t ′′ t 2 /U 3 are also significant. The latter, in particular, determine the sign of t ′ /t ′′ . Of course, the resulting theory has now two more free parameters (i.e. t ′ and t ′′ ) than the original model 6 to fit the inelastic neutron scattering data and the reader may therefore not be surprised that we achieve a better fit than in Ref. [6]. However, as we show in Section III, the ensemble of parameters giving the best unbiased fit to the data is in good agreement with the values suggested from other sources for various cuprates 19,20,21,22,26,40 .
The rest of the paper is organized as follows. In Section II we present the method to obtain a consistent spin representation of the original Hubbard model up to third neighbor hops. As an application, in Section III, we investigate the consequences of applying the method to a general t − t ′ − t ′′ − U model by fitting the magnon dis-persion data of Ref. [6] for La 2 CuO 4 . Section III also discusses the procedure to calculate the staggered magnetization operator, and shows that our expansion of the t−t ′ −t ′′ −U model gives a value of the sublattice magnetization in good agreement with experiment 44 . We conclude the paper in Section IV. We have included a number of appendices to assist the reader with a few technical issues. Appendix A gives the various terms that contribute to the effective spin Hamiltonian with arbitrary t, t ′ and t ′′ up to order 1/U 3 . In Appendix B we show that the renormalization factor for the magnetization, coming from charge fluctuations 33 is the same in our approach as that found in an extension of the meanfield Hartree-Fock method of Ref. [59]. Appendix C gives the k dependence of the various terms coming in the 1/S spin-wave calculation. Appendix D discusses the spinwave velocity renormalization factor, Z c (k). Appendix E comments on the results of various constrained fits to the spin wave energies.
A derivation of the spin-only effective theory starting from the one band t − U Hubbard model has been given by several authors 17,18,32,33 . In this section we investigate the effects of including new parameters t ′ and t ′′ , for direct hops between second and third nearest neighbors, respectively, on the spin-only effective theory. We first recall the key steps in the unitary transformation method for the t − U Hubbard model. We then apply the method to the t − t ′ − t ′′ − U model and obtain the modified spin-only Hamiltonian generated by all possible virtual electronic paths up to order 1/U 3 , giving terms of order Derivation of the spin-only effective Hamiltonian We begin with a brief review of the derivation of the spin-only Hamiltonian of the one band nearest-neighbor Hubbard Hamiltonian, H H : The first term is the kinetic energy term that destroys an electron of spin σ at nearest neighbor site j and creates one on the nearest-neighbor site i. The second term is the on-site Coulomb energy U for two electrons with opposite spin to be on the same site i and where n i,σ = c † i,σ c i,σ is the occupation number operator at site i.
As introduced by Harris et al. 32 and developed further by MacDonald et al. 18 , the transformation relies on the separation of the kinetic part T into three terms that respectively increase by 1 (T 1 ), keep constant(T 0 ) or decrease by one (T −1 ) the number of doubly occupied sites.
We write : whereσ stands for up if σ is down and for down if σ is up. This separation comes from multiplying the kinetic term on the right by n i,σ + h i,σ = 1 and multiplying on the left by n j,σ + h j,σ = 1.
Applying the unitary transformation e iS to H H leads to an effective Hamiltonian,H s , where double occupancy is eliminated perturbatively, order by order in S. Using the relation: (10) and, defining as in Ref. [18] S is solved for, order by order. Hence, starting from a low energy vacuum of singly occupied sites,H s contains no terms that create or annihilate doubly occupied sites up to the order for which S has been determined. This leads to an expression forH s to order t(t/U ) 3 (see Refs. [18,33]): where the associated generator for the unitary transformation is: B. Derivation of the spin Hamiltonian.

The effective HamiltonianH
(4) s in Eq. (12) is still defined in terms of fermion operators entering the T (k) operators. When focusing on the magnetic properties of the half-filled Hubbard model, it is convenient to recast the effective HamiltonianH s in Eq. (12) in a spin-only notation. For this, we use a mapping 18 between the subspace of the Hubbard model with singly occupied sites and the Hilbert space of a spin S = 1/2 system. The mapping is: The spin Hamiltonian H (k) s acting on this space is derived from the Hamiltonian acting on the occupation number subspace (see Ref. 18 for more details). Once it is written in the spin 1/2 basis it can be transformed into an explicitly SU(2) invariant form using: αp is the Pauli matrix for site p, with electron in spin state α p . From this one finds the spin-only Hamiltonian evaluated to third order in t/U for the t−U Hubbard model 18,33 : where < i, j >, << i, j >>, and <<< i, j >>> denote sums over first, second and third nearest neighbors, respectively and where < i, j, k, l > identifies a square plaquette with sites i, j, k, l defining the four corners of an elementary square plaquette.
C. Inclusion of t ′ and t ′′ Having reviewed the procedure for deriving H (k) s for nearest-neighbor hopping t only, we now discuss how to include second and third nearest neighbor hopping t ′ and t ′′ . We start with an extended Hubbard Hamiltonian: where t ′ is the hopping constant to the 2 nd nearest neighbor and t ′′ to the 3 rd nearest neighbor (see Fig. 1); j α is the α nearest neighbor of i. As above, we define the operators T m , taking into account t ′ and t ′′ : where t i,j = t if i and j are nearest neighbors, t i,j = t ′ if i and j are second nearest neighbors, and t i,j = t ′′ if i and j are third nearest neighbors and 0 otherwise. The commutation relations in (10) do not change and the expression for the unitary transformation in (13) remains formally the same. On the other hand, the resulting spin Hamiltonian gets drastically modified by the inclusion of t ′ or t ′′ . In particular, many new plaquette or ring exchange terms are generated. The complete expression for the spin Hamiltonian is reported in Appendix (A). As in previous work 18,33 , and even more so here, given the much increased complexity of the Hamiltonian, a computer program was used to collect together all the terms in the projection of Eq. 15 which ultimately lead to a globally SU(2) invariant effective spin-only Hamiltonian H

III. PHYSICAL RESULTS
A. Spin wave dispersion

Spin wave calculation
In order to obtain the magnetic excitation spectra presented in Fig. 3 and Fig. 4 we perform a 1/S spin wave calculation. The spin operators are written in terms of boson operators through a Holstein-Primakoff 60,61 1/S expansion for a bipartite Néel ordered square lattice: In reciprocal space, one obtains to order S the following general expression of the spin Hamiltonian H (23) Note that here and elsewhere in the paper, the sums in reciprocal space are over the magnetic Brillouin zone, which for a square lattice with N spins contains N/2 sites. H (4) s,quad can be diagonalized through a Bogoliubov transformation giving H and where The magnons for a given wave vector k are thus two fold degenerate with eigenfrequencies 60 : In this expression, the functions A k and B k take into account all the bilinear and ring exchange interactions entering H (4) s and given in Appendix A. Detailed expressions for A k and B k can be found in Appendix C.
The staggered magnetization operator, M † s , is, in reference to sublattice a, defined conventionally for a spinonly Hamiltonian as Rewriting the previous equation in k space and introducing the Bogoliubov transformation above, we arrive at a standard expression 60 : At zero temperature, all α † k α k and β † k β k are zero, and one has for the ground state expectation value of M † s , M † s : give a contribution to the diagonal part of the magnon self energy 43 . This gives rise to a 1/S correction to the spinwave energies ǫ k in Eq. (27). This renormalization of the magnon energy scale is important for the quantitative comparison between our calculation and the experimental results of Coldea et al. 6 and for the subsequent determination of the parameters {t, t ′ , t ′′ , U } and hence we evaluate it for our effective spin Hamiltonian.
Expanding the spin Hamiltonian to fourth order in the boson operators of the Holstein-Primakoff transformation, one finds: H s,quart is now treated as a first order perturbation to the quadratic spin-wave Hamiltonian H (4) s,quad . To this order the magnon energy is shifted but the eigenvectors are unchanged and magnon-magnon interactions are not generated. The shift is found by diagonalizing the set of 2 × 2 matrices with elements 0|γ k H We hence find a correction δǫ k to the magnon energy, and a renormalized magnon energỹ ǫ(k): where we refer to ∆ k as the magnon energy correction and Z c (k) the magnon energy renormalization factor. For the Heisenberg model, the product ǫ k ∆ k corresponds to the leading contribution to the magnon self energy calculated in Ref. [43]. To order t 2 /U and with t ′ and t ′′ set to zero, ∆ k is uniform over the whole Brillouin Zone and we find: as obtained in Ref. [43]. Details of this calculation are given in Appendix D.
When non-zero values of t ′ and t ′′ are included, ∆ k is no longer constant over the Brillouin zone, but the twofold degeneracy for each wave vector k remains. The renormalizing factors Z c (k) are calculated numerically. For a given set of t ′ and t ′′ values, a finite size scaling analysis was used to extrapolate the sums over k ′ in Appendix D over the Brillouin zone to the thermodynamic limit. In Fig. 2 we show the dispersion of Z c (k) obtained for the parameter set (t, t ′ , t ′′ , U ) (see Eq. 38) giving the best fit to the magnon dispersion data 6 from Coldea et al., as described in the next section. For these values the mean value of ∆ k over the Brillouin zone is found to be We also calculate the dispersion of ∆ k over the Brillouin zone for this set of parameters, finding: where σ is the standard deviation of the distribution of values. Details of the calculation can be found in Appendix D. From our results here one can see that the magnon energy renormalization Z c (k), calculated to order 1/S, is indeed an important element of the quantitative comparison between our microscopic theory and the experimental results of Coldea et al. 6 . The values obtained for Z c (k) change the energy scale of the magnon dispersion by around 20% of what would have been obtained if the a bare Z c (k) = 1 had been used. The Z c factors therefore have proportionate consequences for the values of U and t deduced by comparisons with experiment.

Experimental data and previous results
The magnon energy spectrum for La 2 CuO 4 obtained through inelastic neutron scattering 6 is shown in Fig. 3. The data follow a trajectory through the Brillouin zone covering, in particular, the line along the magnetic zone boundary with wave vector in the interval q BZ = [(π, 0), (π/2, π/2)], where the previously discussed downward dispersion is manifest. Also included in this figure is the fit of Ref. [6] using the Hamiltonian (16). The parameters used there to make the fit at a temperature of 10 K are 6 One can see that the energy scale, set by U and the ratio t/U in Eq. (37) are quite different from that found by Coldea et al. 6 . Also, while the effects of t ′ and t ′′ are neglected in Ref. [6], the ARPES results unambiguously show that t ′ and t" have non-negligible values within the CuO 2 planes 22,40 , although band structure calculations 24 In other words, based on various experiments on a number of cuprates, one would have "anticipated" a t/U value for La 2 CuO 4 somewhat smaller than the value in Eq. (36) from Ref. [6], although the experimental uncertainty allows for some overlap between the two. Most importantly, from a theoretical point of view, t/U = 0.135 corresponds to U = 7.4t, which is smaller than the tight binding band width of 8t. This value may therefore be a little small 56   Energy (meV) (π/2, π/2) (π, π) (π, 0) (π/2, π/2) (0, 0) (π, 0) Energy (meV) (π/2, π/2) (π, π) (π, 0) (π/2, π/2) (0, 0) (π, 0)

FIG. 4: (Color online) Magnon energy in La2CuO4
, as a function of wave vector k, across the Brillouin zone at a temperature of 10 K. The experimental data (red squares) are from reference 6 . The fit is made using the parameters given by (38).
From this discussion it is clear that a more detailed analysis of the experimental magnon dispersion that accounts for t ′ and t ′′ is needed.

Fitting procedure and fit to the experimental data
Our theory, H s (t, t ′ , t ′′ , U ) in Appendix A, now contains four independent parameters, t, t ′ , t ′′ , U which we fit to the ensemble of experimental data points from Ref. [6]. As we show in Appendix D, Z c (k) is a function of t/U, t ′ /t, t ′′ /t and k. It thus should be calculated over the ensemble of data points for each iteration of the fitting algorithm. As the expression for Z c (k) contains inner sums ( k ′ ) over the Brillouin zone, the fitting procedure is time consuming. A finite size scaling analysis of the convergence of these sums towards the thermodynamic limit is discussed in Appendix D. The fitting procedure is as follows. We choose a first point at wave vector k with experimental magnon energy ǫ exp (k) and minimize the quantity δ(k) = |ǫ(k) − ǫ exp (k)| with respect to t, t ′ , t ′′ , U . From this first fit we extract a fac-tor 4t 2 /U that fixes the energy scale for the ensemble of points and allows us to writeǫ(k) as (4t 2 /U ) times a function of t/U, t ′ /t, t ′′ /t and k. This energy scaling step can be achieved for many values of t/U, t ′ /t and t ′′ /t and the constraints on these variables are not very high at this stage. The best fit parameter set is now found by minimizing the ensemble of variables {δ(k)(U/4t 2 )} with respect to t/U, t ′ /U, t ′′ /U , using least squares fit. Once the final value of t/U is established, the value of t, and hence U , t ′ and t", are deduced. The values found are: The error bars, obtained through the calculation of the least-squares function χ 2 , represent a deviation of less than 5% from the best fit values. If the dispersion in Z c (k) is ignored, it only needs to be calculated once for each iteration of the ensemble of points, thus considerably speeding up the procedure. Ignoring the dispersion and taking Z c (k) = Z c (k) = constant, the best fit parameter set for this case gives the values: Comparing the two sets, taking the size of the error bars into account and considering the fact that Z c (k) represents only the leading correction to the magnon energies, we conclude that the dispersion of ∆ k does not have a significant effect on the determination of the parameters.
The best fit to the experimental data using parameter set (38) is shown in Fig. (4). Comparing with Fig. (3), one observes that introducing t ′ and t" leads to a clear improvement in the quality of the fit through the experimental points. One might have expected this, given the increased number of free parameters. However, comparing the set (38) with (37), one can also see that the parameters correspond more realistically with those thought to describe other cuprates such as, for example, SrCuO 2 Cl 2 . This is particularly the case for the energy scales t and U which are found to be larger than those in Eq. (36). The ratios t ′ /t and t"/t are also in fair agreement with values for other cuprates and, perhaps most interesting, the fitting procedure even finds the correct signs for these ratios. The ratio t/U is slightly reduced compared with (36) and more in line with a value that one would expect for a system in the Mott insulating regime, although the change is less spectacular than in the case of the absolute energy scales t and U , with the difference in the t/U value remaining within the experimental error bars. However, considering the overall comparison between calculations and the experiments of Ref. [6], it would appear that adding all terms in t, t ′ and t ′′ up to order 1/U 3 does allow a quantitative description of the magnon dispersion in La 2 CuO 4 , which is an improvement over the equivalent procedure which only includes the t and U parameters.

Comments on the influence of t ′ and t ′′
Referring to Appendix A, one can see that the inclusion of t ′ and t ′′ and all four hop processes involving t, t ′ and t ′′ introduces a large number of new ring exchange terms. From this observation one might ask which of the many terms make the most difference and could we have got away with a less complex effective spin Hamiltonian? In this section we address these questions and justify our choice to do a complete calculation of all four hop processes.
The first question is whether or not we can limit ourselves to just one extra parameter t ′ or t ′′ . This point has already been discussed in the Introduction (Section I B). We argued there, (Eq. (3)), that t ′ increases the strength of the antiferromagnetic second neighbour exchange, J 2 → t 4 /U 3 + 4(t ′ ) 2 /U . This in turn modifies the effective second neighbour coupling in the spin wave analysis, J eff 2 → 4(t ′ ) 2 /U − 16t 4 /U 3 making it less strongly ferromagnetic, or even antiferromagnetic if t ′ is sufficiently large. As the negative magnon dispersion along the zone boundary 6 is a consequence of a ferromagnetic (i.e. negative) J eff 2 , one expects that fitting the calculated magnon dispersion to the experimental data will require a larger value of t/U if t ′ = 0 but t ′′ is ignored. This is indeed the case; introducing t ′ to order (t ′ ) 2 /U and setting t ′′ = 0, as well as neglecting all terms of order (t ′ t) 2 /U 3 and (t ′ ) 4 /U 3 , the following best fit parameters are found: The results are very close to those found from a random phase approximation 53,54 and also from quantum Monte Carlo calculations 55 . It is clear that the addition of t ′ only has taken the parameter set in the wrong direction: t/U is now further from most estimates than the case when t ′ = 0. Also, with U = 6t the system could be very near the metallic phase 56 , and not deep in the Mott insulating phase, as is the case for La 2 CuO 4 . As we discussed in the introduction, for U = 6t the t ′ = 0 model is right at the metal-insulator transition and the presence of a finite t ′ increases even further the value of U necessary to position the system in the insulating phase 29,57 . The leading effect of third-neighbor hopping t ′′ is also to increase an exchange interaction, this time J 3 , through the relation : However, in this case antiferromagnetic (positive) exchange also contributes, as does the ring exchange term, to a negative (downward) dispersion along the magnetic zone boundary. Including just t ′′ , one would therefore expect the best fit parameters for the data to include a particularly small ratio of t/U , which is indeed the case.
Hence we see that, as far as this characteristic dispersion is concerned, t ′ and t ′′ play opposing roles. It is therefore clear that one needs to include both hopping constants to get a good fit to the experimental data as we have illustrated in Fig. 4 with parameters (Eq. (37)) similar to those in other cuprates 23,27,41 . Given that individually, terms of order t 2 t ′2 /U 3 are only one tenth of the magnitude of terms of order t ′2 /U or of the terms of order t 4 /U 3 , one might expect the good fit reported in Section ( 38) to be also obtained when taking into account terms of order (t ′ ) 2 /U and (t ′′ ) 2 /U only, and ignoring the host of more complex four hop terms generated by the further neighbor hopping. Not so! If one neglects them and only considers t ′ and t ′′ to order 4(t ′ ) 2 /U and 4(t ′′ ) 2 /U ), one obtains the following set of best fit parameters: Within numerical error this is the same data set as Coldea et al. ( 36), with t ′ /t = t ′′ /t = 0. The rather surprising and interesting conclusion is therefore that, on their own, the two-hop further neighbor terms make no net contribution. Rather, it is the four-hop processes that make the difference. It seems that their small value is compensated for by their multiplicity; that is, by the fact that a given site i appears in a large number of diagrams of order 1/U 3 . Hence, their contribution is increased by an order of magnitude and in this intermediate coupling regime, the global effect of the near-neighbor four hop terms, the majority of which are ring exchanges, is as important as that for the two hop further neighbor processes (Note that as t ′ and t" are both zero in Eq. (42), Z c (k) is once again given by (33) over the whole zone). An alternative procedure would be to impose some partial constraints on the parameters (t, t ′ , t ′′ , U ), using values obtained from, say, band structure calculations. Results of such fits are briefly discussed in Appendix E.
Finally, we remark that in our fitting procedure, the symmetry is broken between data sets with positive and negative values for the ratio t ′ /t ′′ . The best fit occurs for a negative ratio t ′ /t ′′ (see Eq. (37)). Recall that a canonical particle-hole transformation at fixed chemical potential changes the sign of all hopping integrals and leaves the filling invariant, since 2 − n = n when n = 1 3,62 . An additional change of phase of the creation-annihilation operators on one of the sublattices by π restores the sign of nearest-neighbor hopping t, but not that of t ′ and t ′′ . This shows that t ′ /t and t ′′ /t can have arbitrary sign but that the sign of t ′ /t ′′ is physically relevant. The symmetry breaking between positive and negative signs for t ′ /t ′′ can only be accessed by including four hop processes. As one can see from Appendix A, most of the spin interaction terms are even in powers of t, t ′ and t ′′ . However, a small number of terms have odd powers of t ′ and t ′′ . These are the terms that determine the sign of the hopping ratios. This can be seen in Fig. 5 where we show the χ 2 parameter from the fitting algorithm as a function of t ′ /t and t ′′ /t (the value of t/U being set to the best fit value). If the terms proportional to t 2 t ′ t ′′ /U 3 are deleted from H (4) s (see Fig. 5(a)), the χ 2 parameter has four fold symmetry and there are four points of best fit independently of the sign of t ′ /t and t ′′ /t. On the other hand, when these terms are included (see Fig. 5(b)), the symmetry is clearly broken. The minima for positive ratio t ′ /t ′′ become broader and shallower, while those for negative sign are narrower and deeper. The fitting procedure hence favors the two points with a negative sign for the ratio t ′ /t ′′ . To lift the remaining degeneracy we have to appeal to band structure calculations 24,25,41 or ARPES experiments 19,20,21,22,40 , both of which find t ′ /t negative, hence t ′′ /t positive.

B. Magnetization and spin Hamiltonian
Having found a set of parameters (t, t ′ , t ′′ , U ) that suitably describes the experimental spin-wave dispersion data, we now explore how this new set and the derivation of the effective spin Hamiltonian, H (4) s (t, t ′ , t ′′ , U ) affects the value of the zero temperature Néel order parameter.

Magnetization operator
In the construction of effective theories all operators, O H , defined in the original model must be canonically transformed before they can be exploited in a calculation within the spin only theory. That is, within the effective theory, O H becomes O s = e iS O H e −iS and the expectation value in the ground state is defined by Here |0 H and |0 s = e iS |0 H are the ground state wave vectors in the original Hubbard (H) and spin-only (S) models. We have recently shown 33 that incorporating this transformation has important consequences for the ground state magnetization as one moves into the intermediate coupling regime. On application of this procedure on the magnetization operator we find 33 new quantum fluctuations arising from the charge delocalization over closed virtual loops of electronic hops. These spin independent fluctuations allow us to reconcile an apparent paradox concerning the behavior of the ground state magnetization, as a function of t/U . As discussed above, the leading effect of including processes to order t(t/U ) 3 is to introduce effective second and third neighbor spin interactions (see Eq. (1)). Within the lowest 1/S order spin wave approximation, the resulting effective ferromagnetic second neighbor interaction reduces the transverse quantum spin fluctuations and stabilizes the Néel order of the spin model. This would naively seem to imply that the staggered magnetization should be an increasing function of t/U as the system departs from the Heisenberg limit, t/U → 0; a result which is difficult to justify on physical grounds. Indeed, one might expect that as t/U is increased, the enhanced electron mobility would lead to a progressive return to the nonmagnetic metallic state, with the staggered moment decreasing as t/U increases, as found in Ref. [59], rather than increasing. In Ref. [33] we showed that this is indeed the case: the new quantum fluctuations, or "charge fluc-tuation" channel, arising from virtual doubly occupied states, counter the effect of generating a ferromagnetic second neighbor interaction. Consequently, the ground state magnetization for the Hubbard model, when calculated using the spin-only description, is indeed a decreasing function of t/U . The main steps of the calculation are reviewed below and then extended to include the hopping parameters t ′ and t ′′ .
Within the Hubbard Hilbert space, the staggered magnetization operator M † H is defined by: where: The magnetization in the ground state |0 H of the Hubbard model is defined as: Working within the effective spin theory, we have to express all the operators in the spin language. Mathematically this means that the unitary transformation in (13) has to be applied to all operators including M † H : The commutation relations between M † H and the different operators in (13) are 77 Expressed in terms of these operators, the staggered magnetization operator reads 33 : (51) The symmetry associated with half filling ensures that contributions to this expression corresponding to odd numbers of hops, such as the second term in the righthand side of Eq. (51), are zero. Writing the expression in terms of S = 1/2 spin operators 33 we find: In the case of periodic boundaries, or in the thermodynamic limit where boundaries can be neglected, all sites become equivalent and the correction to the magnetization operator becomes a multiplicative factor 33 : whereM † s = 1 N i S z i e iQ·ri is the naive expression for the magnetization operator in the spin language 33 .

Role of t ′ and t ′′
Inclusion of the higher order hopping constants t ′ and t ′′ introduces extra paths for charge delocalization and we might expect new charge fluctuation terms to appear. The hopping operators T m in Eqs. (19,20,21), can be generalized to include t ′ and t ′′ and the commutation relations in Eqs. (48,49,50) are redefined: for hops involving t,t ′ and t ′′ , respectively. The staggered magnetization operator (51) is thus also generalized. Dropping the first order terms in 1/U that do not contribute to the expectation value in a singly-occupied state, M † s takes the form: (57) However, it turns out that the operatorsT ′ m andT ′′ m defined in Eqs. (55,56) are identically equal to0 78 . From this result, it follows that the first corrections to Eq. (53) that depend on t ′ and t ′′ appear beyond the second order in the perturbation scheme. Hence we find that to order (t µ /U ) 2 the expression for M † s (t, t ′ , t ′′ ) is unchanged and given by Eq. (53). However, this does not mean that introducing t ′ and t ′′ has no effect on the expectation value for the staggered magnetization. Firstly, the inclusion of t ′ and t ′′ in the magnon dispersion modifies the ensemble of parameters (t, t ′ , t ′′ , U ) that best fit the data. Secondly, the zero point fluctuations are controlled by (t/U, t ′ /t, t ′′ /t), through the k dependence of u k , v k , A k and B k in Section III A 1 and this, as we discuss below, directly affectsM † s . The charge renormalization factor has also been calculated using a Hartree-Fock mean-field method as in Ref. [59], in which transverse spin fluctuations are neglected (see Appendix B). This calculation on a square lattice gives the same factor, 1 − 8(t/U ) 2 , independently of t ′ and t ′′ . Physically this comes about because, in the spin-density wave ground state (i.e. Néel order, Q = (π, π)), electron hops between two sites that belong to the same Néel ordered sublattice are prohibited by the Pauli principle. Hence, same-sublattice hops by t ′ and t ′′ do not contribute to charge renormalization to order (t ′ ) 2 /U , (t ′′ ) 2 /U .

Consequences
An experimental estimate of the zero temperature sublattice ordered moment in La 2 CuO 4 has been determined from neutron scattering. Lee et al. 44 find where g is the Cu 2+ Landé g-factor. The g-factor for Cu 2+ in layered cuprates such as La 2 CuO 4 has apparently not been determined experimentally 63 . However, a typical value for a distorted octahedral environment is given by Abragam and Bleaney 64 , g ≃ 2.2. We use this value to compare theoretical estimates of the ordered moment with the experimental value. A spin wave calculation for the nearest-neighbor Heisenberg antiferromagnet gives 2,65,66 M † s nn ≃ 0.304, which is in good agreement with the value obtained from quantum Monte Carlo. This gives for the moment, M sw = g M † sw nn , M sw ≈ 0.67µ B , well in excess of the above 0.55µ B value determined experimentally 44 . Using the parameter set of Coldea et al. 6 but neglecting charge renormalization leads to an increased expectation value: M † s ≃ 0.32, taking the moment to 0.70µ B , hence further from the experimental estimate of Ref. [44]. Including the charge renormalization factor 33 1 − 8t 2 /U 2 of Eq. (53) gives M † s ≃ 0.27 and M (t) ≃ 0.59µ B in better agreement with experiment. Including both the charge renormalization and the effects of t ′ and t ′′ we find from which we obtain which is in even better agreement with experiment. Given the uncertainty in the experimental measurement value in Eq. (37) and in the theoretical value that should be taken for g, it is difficult to make in-depth comment on the difference between M tt ′ t ′′ and M t . However, the results do serve to illustrate the importance of the charge renormalization at this level of approximation. Here, it is an essential element in the quantitative agreement with experiment, and it is only through the inclusion of spin independent quantum corrections (e.g. virtual double site-occupancy) that quantitative agreement can be achieved.

IV. DISCUSSION
It is known from ARPES measurements and band structure calculations that t ′ and t ′′ are not negligible compared to t in a number of copper oxide materials. We have thus introduced second and third neighbor hopping parameters, t ′ and t ′′ , into the Hubbard Hamiltonian and included all closed four hop virtual electronic pathways into the canonical transformation to derive an effective spin-only Hamiltonian for the case of a half filled band. This Hamiltonian contains many ring exchange terms. These kinds of terms have been the subject of much discussion lately 33,34,35,36,37,38 . Through this calculation we are able to test the capacity of the one band Hubbard model to describe the magnetic properties of the antiferromagnetic parent high-temperature superconductor La 2 CuO 4 . We find in general good quantitative agreement between the predictions of the t − t ′ − t ′′ − U model and experimental measurements of the magnetic properties of the compound 67 .
An exact solution of the spin-only Hamiltonian we have derived is clearly far out of reach. In fact, because of the frustration and ensuing sign problem introduced by the various frustrating bilinear spin-spin couplings and the many ring exchange terms, even a quantum Monte Carlo attack on the low temperature properties of this model seems difficult to envisage in the near future. However, it has long been known that a spin-wave analysis up to leading order in a 1/S expansion accurately reproduces the zero temperature staggered moment in the Heisenberg model on a square lattice 2 . This is because of the cancellation of the magnon-magnon interaction terms to second order in 1/S for this particular lattice 43 . This suggested that a similar spin wave analysis can be a good starting point to describe the observed magnetic excitation spectrum in La 2 CuO 4 and this is indeed what we have found in this paper.
In our spin wave calculations, we calculated the magnetization to leading order in 1/S and included the first order correction to the classical spin wave frequencies.
The magnon energies are thus renormalized by a spinwave renormalization factor Z c (k). Since the classical spin-wave frequencies are already of order 1/U 3 , we only retained terms of order t 2 /U , (t ′ ) 2 /U and (t ′′ ) 2 /U in cal-culating Z c (k). Hence, to this order, it does not depend on the ring exchange terms which are of order 1/U 3 . In the end, Z c (k) raises the energy scale by about 20% over the whole of the Brillouin zone, allowing good agreement between experiment and theory without the introduction of any arbitrary scale factors.
To make contact with experiments on La 2 CuO 4 the spin wave calculation was set up starting from the classical antiferromagnetic Néel ground state. Within this framework, the general effect of the ring exchange is to reduce the transverse quantum spin fluctuations rather than increase them (see however the discussion of charge renormalization in the following paragraph). With the parameters coming out of the calculation, we therefore seem far from entering an exotic quantum phase driven from large quantum spin fluctuations about the Néel ordered state. Yet, the ring exchanges are found to play a very important role. In particular, it seems that the large numbers of such terms and their large dimensionless prefactors (see Appendix A) outweighs their small O(1/U 3 ) amplitude, although a more extensive study is required here to understand this point in detail. Further, ring exchange plays a crucial role in determining the relative sign of the ratios t ′ /t and t ′′ /t. There is experimental evidence from ARPES measurements and band structure calculations that t ′ /t is negative while t ′′ /t is positive in a number of cuprates. In our calculation there are terms of order t 2 t ′ t ′′ /U 3 which break the t ′′ /t ′ → −t ′′ /t ′ symmetry in the effective spin-only Hamiltonian H (4) s , allowing us to deduce the correct sign for the ratio of t ′ /t ′′ compared with experiment. In our procedure, we determined unconstrained values of U, t, t ′ and t ′′ by comparing our magnon dispersion data with that from experiments on La 2 CuO 4 6 . The best fit is found with one positive (t ′′ /t) and one negative (t ′′ /t ′ ) ratio, exactly as in ARPES experiments and band calculations and, for example, our best fit data set is in good agreement with those from the ARPES experiments on Sr 2 CuO 2 Cl 2 (compare Eq. (38) with Eq. (37)).
At the level of approximation at which the canonical transformation is performed 18,32,33 , finite charge mobility renormalizes the magnetic moment calculated from the spin only Hamiltonian by a factor (1 − 8t 2 /U ) 33 . We found that t ′ and t ′′ do not contribute further to this factor. This is because direct second and third neighbor hops take electrons from one site to another on the same sublattice, hence virtual double occupancy is excluded by the Pauli exclusion principle. The expectation value for the zero temperature staggered moment coming from the spin only Hamiltonian is thus scaled by the same (1−8t 2 /U ) factor as in the case where t ′ = t ′′ = 0 33 . The relevance of this term can be estimated by comparing theoretical and experimental estimates of the total magnetic moment. Here, we also find good agreement between the estimates of Lee et al in Ref. [44] and our calculated value. The charge renormalization is important here, as it scales the magnetization obtained by considering solely transverse spin fluctuations, by an amount clearly in excess To close, and reiterate, one of the noteworthy results of this work is that the presence of the ring exchange terms in the effective spin-only theory allow us to determine the relative sign of t ′ and t ′′ . More generally, the theory presented herein, restricted to the reduced Hilbert space of spin degrees of freedom only, gives results for the electronic band parameters that are compatible with ARPES experiments and theoretical band calculations for a variety of cuprate materials. In addition, the value of U ≃ 8t that we find is consistent with the origin of the insulating behavior in parent high-temperature superconductors being mostly due to Mott Physics. As a final word, we propose that our results provide further evidence that the t−t ′ −t ′′ −U one-band Hubbard model in the intermediate coupling regime gives a consistent unified description of high-temperature cuprate superconductors and of their parent insulating phases. The purpose of this appendix is to give the new spinspin interaction terms appearing in the spin-only Hamiltonian including t ′ and t ′′ up to order 1/U 3 in the canonical transformation. The sole effect of introducing t ′ and t ′′ to order (t ′ ) 2 /U and (t ′′ ) 2 /U is to renormalize the coupling constants already present in the spin Hamiltonian (Eq. 16) generated from the t − U Hubbard model, up to order t 4 /U 3 . Including t ′ and t ′′ up to order, 1/U 3 generates many more terms that further renormalize these coupling constants. However, a number of terms with new "topologies" appear and these play a key role on the results discussed in this paper. For reasons of compactness we introduce the following notation: Figure 6 illustrates the sites involved in t ′ and t ′′ hopping i j m k l FIG. 6: Label of the different sites involved in the many-spin terms induced by t ′ and t ′′ to order 1/U 3 (see also Fig. (1)). The symbol means that the corresponding sites participate in the expression of the coupling interaction, whereas the symbol means that these sites are transparent in the electronic process. That is the electron 'hops over' the site in question. As there is at present considerable interest in diluted Mott-Hubbard systems, such as in the La 2 Cu x Zn 1−x O 4 / La 2 Cu x Mg 1−x O 4 / 73,74,75,76 , we keep track, in the derivation of the spin Hamiltonian, of the occupation of the sites visited by the electrons, so that ǫ i = 1 if a site i is occupied by a spin and ǫ i = 0 if the site is not occupied. Here we only consider the case of 100% Cu 2+ site occupation (x = 1).
The sublattice magnetization for nesting wave vector Q is defined by: where we consider the case where the magnetization is polarized along theẑ direction. |Ω is the spin density wave ground state and Q. The spin density operator S i q for arbitrary q wave vector is defined by: N being the number of sites. We define new operators γ c k,α and γ v k,α through the Bogoliubov transformation: and inject (B7) into the Hartree-Fock factorization 59 of Hamiltonian (B3), finding: where the density wave gap ∆ is given by: The Bogoliubov coefficients for this transformation u k and v k are: With the valence band filled and the conduction band empty, since the calculation here is done at zero temperature, it follows from the above relationships that hence: Putting (B9) into (B12) gives the following self-consistent equation for the sublattice magnetization: We remark that equation (B13) is independent of ǫ ′ k or ǫ ′′ k . Considering t/U ≪ 1, we can expand the square roots of (B13) in ǫ k /U , after which we find: (B14) Further simplifications lead to: To first order in t/U we obtain: N k (cos(k x ) + cos(k y )) 2 , and 1 N k (cos(k x ) + cos(k y )) 2 = 1 2 , we conclude that: This result is equivalent to that obtained through the unitary transformation method in Eq. (53) if transverse spin fluctuations are neglected, as they are in the Hartree-Fock approach. Mathematically, it is the property ǫ k = −ǫ k+Q that makes nearest-neighbor hopping different from the other two hops which both obey, instead, ǫ ′ k = ǫ ′ k+Q , ǫ ′′ k = ǫ ′′ k+Q . Note that ǫ k , ǫ ′ k , ǫ ′′ k in this appendix are the quasi-particle energies, which should not be confused with the spin-wave (magnon) energies in Section III A 1.
Physically, since t ′ and t ′′ are hops within the same sublattice (parallel spins) of the classical Néel antiferromagnetic solution, in the strong coupling limit the Pauli principle prohibits hopping between these sites. Hence, to leading order, t ′ and t ′′ cannot change double occupancy or, in other words, produce charge fluctuations (or "charge renormalization").
from which the magnon energy correction ∆ k ≡ ∆ǫ k /ǫ k can be computed. ∆ k converges numerically to give uniformly over the Brillouin zone, as found previously in Ref. [43]. When t ′ or t ′′ are brought into the picture the leading effects are to create second and third nearest neighbor interactions J 2 (t ′ ) = 4(t ′ ) 2 /U and J 3 (t ′′ ) = 4(t ′′ ) 2 /U producing new contributions to the quartic Hamiltonian: The contribution coming from t ′ reads: while that from t ′′ reads: (D6) We now find that the total contribution of the terms coming from t, t ′ and t ′′ is no longer independent of wave vector, but the two-fold degeneracy for each k remains. The evolution of Z c (k) over the Brillouin zone is shown in Fig. 2 for the parameter set (38).
Dividing the contributions to Z c (k) into three parts coming from t, t ′ and t ′′ and denoting an average over the Brillouin zone by · k we find    ∆ k = 0.187, ∆ ′ k = 0.019, ∆ ′′ k = 0.013, using the parameter set (38). Note that here ∆ ′ k and ∆ ′′ k do not factorize analytically in Eq. (D5) or (D6), rather, for each wave vector, we calculate ∆ ′ k ≡ δǫ ′ k /ǫ ′ k and similarly for ∆ ′′ k . Note also that as t ′ and t ′′ become non-zero, even the first term evolves. This is because ǫ k is itself a function of t ′ and t" and so it changes as these parameters are switched on.
As discussed in Section III A 4, incorporating the kdependence of Z c (k) makes the fitting procedure somewhat computationally cumbersome and slow from a CPU speed point of view. There are two reasons for this. Firstly, it is due to the need to constantly redo the "internal" sum over k ′ in Eqs. (D2), (D5) and (D6) whenever the t, t ′ , t ′′ and U parameters are readjusted. Secondly, to compound this problem the CPU time is further increased by noting that the sums over k ′ in Eqs. (D2), (D5) and (D6) converge somewhat slowly with system size. k ′ in equations (D2,D5,D6) This is illustrated in Fig. 7, for a magnetic Brillouin zone of linear dimension L, for the parameter set (38) (here the lattice parameter is taken to be unity) Re-plotting the data in Fig. 8, we show that ∆ k scales with system size as where α is a constant, from which we find: The width of the dispersion of ∆(k) over the Brillouin zone, σ(∆(k)), scales in a similar way and we eventually find Finite size effects lead to negligible corrections compared to the least squares fitting procedure, for L ∼ 100. A magnetic Brillouin zone of this size was used to make the fits described in the main text. (BSCCO, Bi-2212) and YBa 2 Cu 3 O 6+δ (YBCO), values of t ′ ∼ −0.3 eV, t ′′ ∼ 0.2 eV, and U > 10 eV are commonly used 22 . However, these values are not truly universal among the cuprates. For example, in Ref. [25], local density approximation (LDA) calculations do predict variations from one compound to the other. In particular, for La 2 CuO 4 , t ′ /t = −0.17 is obtained, and which corresponds to the value used in ARPES data analysis 21 .
A recent ARPES experiment by Yoshida et al. 19 on doped La 2−x Sr x CuO 4 finds that t ′ varies slightly between compositions. For x = 0.03, they obtain t ′ /t ≈ −0.2 and t ′′ /t ′ ≈ 0.5. To the best of our knowledge, there exist no ARPES measurements on undoped La 2 CuO 4 , but one may extrapolate to t ′ /t = −0.21 for this composition. In this context, we have performed a constrained fit to the spin-wave dispersion data of Ref. [6], imposing t ′ /t = −0.2, t ′′ /t ′ = 0.5 and t = 0.25 eV and only allowing U to vary. We were unable to find a reasonable value of U with such constraints. Indeed, we found U ∼ 1.6 ev with a poor quality of fit.
We then allowed both t and U to vary. The motivation being that the effective ratios t ′ /t = −0.2 and t ′′ /t ′ = 0.5 could be well determined by ARPES, while the value of t determined from the electronic bandwidth could be strongly renormalized. The best fit then obtained was t = 0.3430 eV and U = 2.55 eV, hence U/t = 7.42. Compared to the fit with free t, t ′ , t ′′ and U values, the overall χ 2 , with where E i is n'the magnon energy data point, is increased by roughly 25% for the fit with t ′ /t and t ′′ /t ′ constrained compared with the fit with the values in Eq. (38).
This value of t ≈ 0.34 eV agrees roughly with those found by Coldea et al. 6 , and the value U ≈ 2.55 is inbetween that of Ref.
[6 and the one found in this work and reported in Eq. (38). In other words, for the constrained values of t ′ and t ′′ used in such a fit, the values are sufficiently "small" that the results for t ′ = t ′′ = 0 of Ref. [6] are more or less recovered. In other words, the ratio t/U is a crucial parameter in the fit. Working with a fit with t ′ /t and t ′′ /t fully constrained leads to a (reduced) ratio t/U compared to the unconstrained fit of Eq. (38). Yet, this reduced U = 2.55 eV also gives a reduced t = 0.34 ev, with the result that while the unconstrained fit has t/U = 0.126, the constrained fit has t/U = 0.34/2.55 = 0.133 (I just calculated) and Coldea et al. finds t/U = 0.135. Given that t ′ /t ′′ is essentially the same in the constrained and unconstrained fits, one can think of the constrained fit as an intermediate point in switching on t ′ and t ′′ at constant t ′ /t ′′ . Indeed, we find the following progression of values, moving from the Coldea et al. results, the constrained fit results and the unconstrained fit results: U = 2.3 eV, 2.55 eV and 3.34 eV, t = 0.3 eV, 0.34 eV and 0.42 eV, and t/U = 0.135, 0.133, and 0.126. This progression of parameters appears monotonous. Given that we expect t ′ and t ′′ to frustrate Néel ordering, we also expect the metal insu-lator transition also to appear for smaller values of t/U with increasing t ′ and t ′′ . Our results are reassuring in this respect because we find the fitted t/U value getting progressively small as t ′ and t ′′ are turned on, consistent with the system remaining in the insulating Néel ordered state.