The Two-Dimensional Infinite Heisenberg Classical Square Lattice: Zero-Field Partition Function and Correlation Length

The nonlinear σ-model has known a new interest for it allows to describe the properties of two-dimensional quantum antiferromagnets which, when properly doped, become superconductors up to a critical temperature notably high compared to other types of superconducting materials. This model has been conjectured to be equivalent at low temperatures to the two-dimensional Heisenberg model. In this article we rigorously examine 2d-square lattices composed of classical spins isotropically coupled between first-nearest neighbors (i.e., showing Heisenberg couplings). A general expression of the characteristic polynomial associated with the zero-field partition function is established for any lattice size. In the infinitelattice limit a numerical study allows to select the dominant term: it is written as a l-series of eigenvalues, each one being characterized by a unique index l whose origin is explained. Surprisingly the zero-field partition function shows a very simple exact closed-form expression valid for any temperature. The thermal study of the basic l-term allows to point out crossovers between land (l+1)-terms. Coming from high temperatures where the l=0-term is dominant and going to zero Kelvin, l-eigenvalues showing increasing l-values are more and more selected. At absolute zero l becomes infinite and all the successive dominant l-eigenvalues become equivalent. As the z-spin correlation is null for positive temperatures but equal to unity (in absolute value) at absolute zero the critical temperature is absolute zero. Using an analytical method similar to the one employed for the zero-field partition function we also give an exact expression valid for any temperature for the spin-spin correlations as well as for the correlation length. In the zero-temperature limit we obtain a diagram of magnetic phases which is similar to the one derived through a renormalization approach. By taking the low-temperature limit of the correlation length we obtain the same expressions as the corresponding ones derived through a renormalization process, for each zone of the magnetic phase diagram, thus bringing for the first time a strong validation to the full exact solution of the model valid for any temperature.


Introduction
Since the middle of the eighties with the discovery of hightemperature superconductors [1], the nonlinear σ -model has known a new interest for it allows to describe the properties of two-dimensional quantum antiferromagnets such as La 2 CuO 4 [2][3][4][5][6][7][8][9][10]. These antiferromagnets, when properly doped, become superconductors up to a critical temperature T c notably high compared to other types of superconducting materials.
For studying the magnetic properties of such magnets Chakravarty et al. [6] have shown that it is necessary to consider the associated space-time which is composed of the crystallographic space of dimension d to which a time-like axis, namely called the iτ-axis, is added. The space-like axes are infinite but the time-like axis has a finite length called the "slab thickness" which is inversely proportional to the temperature T and hence goes to infinity as T goes to zero. As a result D = d + 1 is the space-time dimension: here d = 2 and D = 2+1. The nonlinear σ -model in 2+1 dimensions has been conjectured to be equivalent at low temperatures to the twodimensional Heisenberg model [11,12], which in turn can be derived from the Hubbard model in the large U-limit [13].
In their seminal paper, Chakravarty et al. [6] have studied this model using the method of one-loop renormalization group (RG) improved perturbation theory initially developed by Nelson and Pelkovits [14]. These authors have related the σ -model to the spin-1/2 Heisenberg model by simply considering [6,10]: (i) a nearest-neighbor s = 1/2 antiferromagnetic Heisenberg Hamiltonian on a square lattice characterized by a large ex-change energy; (ii) very small interplanar couplings and spin anisotropies.
In addition, they have pointed out that the long-wavelength, low-energy properties are well described by a mapping to a two-dimensional classical Heisenberg magnet because all the effects of quantum fluctuations can be resorbed by means of adapted renormalizations of the coupling constants. A low-temperature diagram of magnetic phases has been derived. It is characterized by three different magnetic regimes: the Renormalized Classical Regime (RCR), the Quantum Critical Regime (QCR) and the Quantum Disordered Regime (QDR). For each of these regimes Chakravarty et al. [6] have given a closed-form expression of the correlation length ξ exclusively valid near the critical point T c = 0 K. Finally these authors have shown that the associated critical exponent is ν = 1.
A little bit later Hasenfratz and Niedermayer published a more detailed low-T expression of ξ for the RCR case, exclusively [15]. Finally, also using a RG technique, Chubukov et al. [10] reconsidered the work of Chakravarty et al. [6] by detailing the static but also the dynamic low-T magnetic properties of antiferromagnets. They notably published exact expressions of ξ and the magnetic susceptibility χ (restricted to the case of compensated antiferromagnets), also exclusively valid near T c = 0 K, for each of the three zones of the magnetic diagram.
From an experimental point of view, at the end of the nineties, the first two-dimensional (2d) magnetic compounds appeared [16][17][18][19]. Some of them were composed of sheets of classical spins (i.e., manganese ions of spin S = 5/2) well separated from each others by nonmagnetic organic ligands. These 2d compounds were the first ones whose low-T magnetic properties were characterized by a quantum critical regime.
Thus the necessity of fitting experimental susceptibilities as well as the important theoretical conclusions of the respective works of Chakravarty et al. [6] and Chubukov et al. [10] motivated us to focus on the two-dimensional O(3) model developed on a square lattice composed of classical spins [20][21][22][23].
The mathematical framework common to our first series of articles was the following one: (i) we first considered the local exchange Hamiltonian where β = 1/k B T is the Boltzmann factor and J the exchange energy between consecutive spin neighbors. As a result the polynomial expansion describing the zerofield partition function Z N (0) directly appears as a characteristic l-polynomial, for the considered lattice.
We observed that, for most of the examined compounds showing a low-T QCR case, when fitting the corresponding experimental susceptibilities, the characteristic l-polynomial associated with the theoretical susceptibility χ could be restricted to the dominant term characterized by l = 0, as for Z N (0), in the physical case of an infinite lattice. In other words no mathematical study was necessary in spite of the fact that this assumption gave good results for the involved exchange energies J i.e., the exact corresponding tabulated experimental values, with a Landé factor value very close to the theoretical one G = 2 (in µ B /ℏ unit). However we also discovered that, for some compounds characterized by the same low-T quantum critical regime, it was necessary to take into account the terms l = 0 but also l = 1 (with m = 0) in the l-expansion of χ for obtaining a good fit of experimental susceptibilities.
Thus, from a theoretical point of view, the condition leading to choose the term l = 0 exclusively or the terms l = 0 and l = 1 (m = 0) in the common l-polynomial part shared by Z N (0) and χ remained a puzzling question. For all the experimental fits, the lowest possible value reached by temperature for ensuring a pure two-dimensional magnetic behavior was T = T 3d when the 3d-magnetic ordering appears.
We then observed that, if restricting the l-expansion of χ to the term l = 0, we had to fulfil the numerical condition k B T 3d /|J| ≥ 0.255. But, if compelled to consider the terms l = 0 and l = 1 (m = 0), we had k B T 3d /|J| ≥ 0.043. Finally the lowtemperature theoretical diagram of magnetic phases was restricted to a single phase, the quantum critical regime, in contradiction with the results derived near T c = 0 K, from a renormalization technique which points out three different magnetic regimes [6,10].
In order to solve these difficulties a full study of the characteristic l-polynomial associated with Z N (0) appeared as unavoidable. This is the aim of the present paper. Even if starting with the same mathematical considerations common to the first series of papers previously published and all restricted to the case l = 0 [20][21][22][23] this paper is intended as a new work because, in section 2 and for the first time, we establish the complete closed-form expression of the characteristic l-polynomial associated with Z N (0), valid for any lattice size, any temperature and any l.
The examination of the case of a finite lattice is out of the framework of the present article [24]. Then we exclusively consider the physical case of an infinite lattice (i.e., the thermodynamic limit) in section 2. We numerically show that, if studying the angular part of each l-term of the characteristic l-polynomial, the value m = 0 is selected. In addition we formally prove that the higher-degree term of the characteristic l-polynomial giving Z N (0) is such as all the l's are equal to a common value l 0 . Surprisingly we then obtain a very simple closed-form expression for Z N (0), valid for any temperature and any l.
Finally, in section 2, we report a further thermal numerical study of the l-higher-degree term. Thus and for the first time, this study allows to point out a new result i.e., thermal crossovers between two consecutive l-and (l+1)-eigenvalues. It means that the characteristic l-polynomial can be reduced to a single l-term within a given temperature range but, for the whole temperature range, all the l-eigenvalues must be kept. That allows to explain that the value l = 0 characterizes the dominant term for reduced temperatures such as k B T/|J| ≥ 0.255. For 0.255 ≥ k B T/|J| ≥ 0.043 we have l = 1 and so on. Finally l-eigenvalues ( ) l J λ β , with increasing l > 0, are successively dominant when temperature is decreasing down to 0 K. In the vicinity of absolute zero the dominant term is characterized by l → +∞.
As all the l-eigenvalues show a very close low-temperature behavior we then deal with a continuous spectrum of eigenvalues, confirming the fact that the critical temperature is T c = 0 K, in agreement with Mermin-Wagner's theorem [25].
From a mathematical point of view it means that Z N (0) is given by a series of continuous functions ( ) l J λ β − , the modified Bessel functions of the first kind, in the whole range of temperature so that, even if considering a specific temperature range inside which the Bessel function is dominant, no singularity can occur.
In section 3 we analytically show that the spin correlation is such as <S z > = 0 for T > 0 K whereas <S z > = ±1 for T = 0 K, again confirming the fact that T c = 0 K. Then, for the first time, in the thermodynamic limit, we obtain the exact closedform expression of the spin-spin correlation between any couple of lattice sites (0,0) and (k,k'), valid for any temperature.
For doing so we first show that all the correlation paths are confined within a closed domain called the " correlation domain" which is a rectangle whose sides are the bonds linking sites (0,0), (0,k'), (k,k') and (k,0) (theorem 1). Second we prove that open or closed loops are forbidden so that all the correlation paths show the same shortest possible length between any couple of lattice sites. All of them have the same weight i.e., they are composed of the same number of horizontal (respectively, vertical) bonds as the horizontal (respectively, vertical) sides of the correlation domain (theorem 2). This allows to derive an exact expression of the correlation length ξ also valid for any temperature.
In section 4 we examine the low-temperature behavior of the λ l (β|J|)'s. We retrieve the low-temperature magnetic phase diagram with 3 regimes. It is strictly similar to the one derived from a renormalization technique [6,10]. By taking the low-temperature limit of the correlation length ξ we obtain the same expressions as the corresponding ones derived through a renormalization process, for each zone of the magnetic phase diagram, thus bringing for the first time a strong validation to the full exact solution of the model valid for any temperature. At T c = 0 K we retrieve the critical exponent ν = 1, as previously shown [6,10].
Finally, near the critical point, the correlation length ξ x can be simply expressed owing to the absolute value of the renormalized spin-spin correlation | . 0,0 0,1 < > S S | between first-nearest neighbors i.e., sites (0,0) and (0,1). Section 5 summarizes our conclusions. The appendix gives all the detailed demonstrations necessary for understanding the main text, notably the lowtemperature study of key physical parameters.

Definitions
The general Hamiltonian describing a lattice characterized by a square unit cell composed of (2N + 1) 2 sites, each one being the carrier of a classical spin S i,j , is given by: , (1) with N → +∞ in the case of an infinite lattice on which we exclusively focus in this article and ex , 1 , 1 2 1, , where: In equation (2) J 1 and J 2 refer to the exchange interaction between first-nearest neighbors belonging to the horizontal lines and vertical rows of the lattice, respectively. J i > 0 (respectively, J i < 0, with i = 1, 2) denotes an antiferromagnetic (respectively, ferromagnetic) coupling. G i,j is the Landé factor characterizing each spin S i,j and expressed in µ B /ℏ unit.
Finally we consider that the classical spins S i,j are unit vectors so that the exchange energy JS(S + 1) ∼ JS 2 is written J. It means that we do not take into account the number of spin components in the normalization of S i,j 's so that S 2 = 1.
When =0 the zero-field partition function Z N (0) is defined as: where β = 1/k B T is the Boltzmann factor. In other words the zero-field partition function Z N (0) is simply obtained by integrating the operator ex exp( ) H −β over all the angular variables characterizing the states of all the classical spins belonging to the lattice.

Preliminaries
Due to the presence of classical spin momenta, all the operators ex , i j H commute and the exponential factor appearing in the integrand of equation (5) can be written: As a result, the particular nature of ex , i j H given by equation (2) allows one to separate the contributions corresponding to the exchange between first-nearest neighbor classical spins. In fact, for each of the four contributions (one per bond connected to the site (i,j) carrying the spin S i,j ), we have to expand a term such as exp (−AS 1 .S 2 ) where A is βJ 1 or βJ 2 (the classical spins S 1 and S 2 being considered as unit vectors). If we call Θ 1,2 the angle between vectors S 1 and S 2 , characterized by the couples of angular variables (θ 1 , ϕ 1 ) and (θ 2 , ϕ 2 ), it is possible to expand the operator exp(−AcosΘ 1,2 ) on the infinite basis of spherical harmonics which are eigenfunctions of the angular part of the Laplacian operator on the sphere of unit radius S 2 : In the previous equation the (π/2A) 1/2 I l+1/2 (−A)'s are modified Bessel functions of the first kind; S 1 and S 2 symbolically represent the couples (θ 1 , ϕ 1 ) and (θ 2 , ϕ 2 ). If we set: the λ l 's are nothing but the associated eigenvalues. Under these conditions, the zero-field partition function Z N (0) directly appears as a characteristic polynomial.
In the case of an infinite lattice edge effects are negligible so that it is equivalent to consider a lattice wrapped on a torus characterized by two infinite radii of curvature. Horizontal lines i = −N and i = N on the one hand and vertical lines j = −N and j = N on the other one are confused so that there are (2N) 2 sites and 2(2N) 2 bonds, with N → +∞. As a result Z N (0) can be written as: where F (i,j) is the current integral per site (i,j) (with one spherical harmonics per bond). Using the following decomposition of any product of two spherical harmonics appearing in the integrand of F (i,j) [26]

Principles of Construction of the Characteristic Polynomial Associated with the Zero-field Partition Function
The zero-field partition function given by equation (9) can be rewritten under the general form The examination of equation (13) giving the polynomial expansion of Z N (0) allows one to say that its writing is nothing but that one derived from the formalism of the transfer-matrix technique. Each current term appears as a product of two subterms: (i) a temperature-dependent radial factor containing a pro- Equation (13) can also be artificially shared into two parts labelled Part I and Part II of respective zero-field partition functions (0) so that the zero-field partition function Z N (0) can be written as with ( ) i.e., all the bonds are characterized by the same integer l but we can have a set of different relative integers m i,j ∈[−l,+l] and m' i,j ∈[−l,+ l] with m i,j = m' i,j or m i,j ≠ m' i,j . Part II appears as a product of "cluster" terms such as with n k < (2N) 2 and the condition n 1 + n 2 +... + n k = (2N) 2 . Thus, only n k bonds are characterized by the same integers l k , l' k and a collection of different relative integers m i,j ∈[−l k ,+l k ] and m' i,j ∈[−l' k ,+l' k ], with m i,j = m' i,j or m i,j ≠ m' i,j .

General Selection Rules for the Whole Lattice
The non-vanishing condition of each current integral F (i,j) due to that of C. G. coefficients allows one to derive two types of universal selection rules which are temperatureindependent.
The first selection rule concerns the coefficients m and m' appearing in equation (12). We have (2N) 2 equations (one per lattice site) such as: At this step we must note that, if each spherical harmonics the non-vanishing condition of the ϕ-part directly leads to equation (17). As a result, we can make two remarks: the SRm relation is unique; due to the fact that the ϕ-part of the The second selection rule is derived from the fact that the various coefficients l and l' appearing in equation (12) obey triangular inequalities as noted after this equation [24]. If M i,j ≠ 0 the determination of l i,j 's and l' i,j 's is exclusively numerical. If M i,j = 0 we have a more restrictive vanishing condition [26]: where K is a coefficient depending on l 1 , l 2 , l 3 and g [26]. In equation (12) , , , , we must have l i,j−1 + l' i+1,j + L i,j = 2A' i,j ≥ 0. Thus, if summing or substracting the two previous equations over l and l', we have (2N) 2 equations (one per lattice site) such as: (or equivalently where g' i,j or g'' i,j is a relative integer. We obtain two types of equations which are similar to equation (17) but now, in equation (19), instead of having a null right member like in equation (17), we can have a positive, null or negative but always even second member.

Zero-field Partition Function in the Thermodynamic Limit
The case of thermodynamic limit (N → +∞) is less restrictive than that of a finite lattice studied in previous articles [24] because we deal with a purely numerical problem. For simplifying the discussion when necessary we shall restrict the study to the case J = J 1 = J 2 without loss of generality. If not specified we shall refer to the general case As previously explained in Subsec. 2.3 the characteristic polynomial associated with Z N (0) for N finite or infinite is composed of two parts: Part I of current term and Part II whose current term is a product of "cluster" terms , so that, in both cases, the numerical study concerns the common term The m's and l's (respectively the m' 's and l' 's) appearing in F (i,j) (cf equation (12)) can vary or not from one site to another site.  (i) First let us consider the case of the m's. We have to solve a linear system of (2N) 2 equations (17) (one per site) but with 2(2N) 2 unknowns m i,j and m' i,j . As it remains 2(2N) 2 − (2N) 2 = (2N) 2 independent solutions over the set of relative integers m i,j and m' i,j (with here N → +∞) it means that there are (2N) 2 different expressions (i.e., here an infinity) for each local angular factor appearing in each term of the characteristic polynomial giving Z N (0) so that the statistical problem remains unsolved. Thus, at first sight, this result means that there is no unique expression for Z N (0).
In summary there are only (2N) 2 independent solutions i.e., (2N) 2 different sets of coefficients (m i,j ,m' i,j ) obeying equation (17). The simplest solution is given by the condition (ii) The study of l's is strictly similar to that of m's because we have to solve a linear system of (2N) 2 equations (19) (one per site) with 2(2N) 2 unknowns l i,j and l' i,j . We have (2N) 2 independent solutions over the set of integers l i,j and l' i,j and the particular solution l i,j = l' i,j ≠ 0 or l i,j = l' i,j = 0, for the whole lattice.
In other words, when N → +∞, a separate numerical study of integrals F (i,j) must allow one to select a unique m-value so that F (i,j) is maximum. First we restrict the set of integers l i,j , l i−1,j , l' i,j and l' i+1,j appearing in the integrand of F (i,j) to two different integers l i and l j . In that case, if setting In the infinite-lattice limit we expect that the highest eigenvalue must naturally arise in Part I of current term . If using equation (14) with this contribution. We make the assumption that it dominates all the other terms inside Part I as well as all the ones composing Part II defined by equations (15) and (16) [24]. This can occur in the whole temperature range or in a smaller temperature one if there exist thermal crossover phenomena among the set of eigenvalues.
In that case the dominant eigenvalue over a temperature range becomes subdominant when the temperature T is outside this range and a new eigenvalue previously subdominant becomes dominant and so on. In the case of 1d spin chains, we always have the same highest eigenvalue (within the framework described previously).
In this respect we have first studied the integral F (i,j) = F(l i ,l j ,m) given by equation (20) with l i = l j = l ≥ 0. In Figure  1 we have reported the ratio F(l,l,m)/F(0,0,0) vs l for various m-values such as m≤ l (with F(0,0,0) = 1/4π). We immediately observe that this ratio rapidly decreases for increasing m-values, for any l. However, we have zoomed the beginning of each curve corresponding to the case m= l. This trend is not followed but we always have F(l,l,m) < F(l,l,0) (see Figure 2). Second, in Figure 3, for m = 0, we observe that F(l i ,l j ,0) decreases for l j < l i . As a result, when N → +∞, the integral F(l,l,0) obtained if l = l i = l j appears as the dominant one i.e., so that the value m = 0 is selected. In addition this result shows that it is not necessary to consider 4 different integers l i,j and l' i,j in the integral F(l i ,l j ,m).   For sake of simplicity we now restrict to the case J = J 1 = J 2 . Under these conditions equation (13) can be rewritten in the thermodynamic limit: The notation In a first step we must wonder if all the current terms of the previous l-series must be kept in the first term of equation (22) i.e., if the series must be truncated, for a given range of temperature [ , where F L,L = F(L,L,0) is given by equation (10) reduced to m = m' = 0 as well as the following ratio: We have studied the thermal behavior of r l+1 /r l for various finite l-values. This work is reported in Figure 4. We observe that log 10 (r l+1 /r l ) shows a decreasing linear behavior with respect to k B T/J. We have zoomed Figure 4 in the very low-temperature domain ( Figure 5). If r l+1 /r l < 1 log 10 (r l+1 /r l ) < 0 and log 10 (r l+1 /r l ) > 0 if r l+1 /r l > 1. We can then point out a succession of crossovers, each crossover being characterized by a specific temperature called crossover temperature T CO . T CO is the solution of the equation: i.e., owing to equation (24): T . We have reported k B T CO /|J| vs l in Figure 6.
As expected we observe that T CO rapidly decreases when l increases. It means that, when the temperature tends to absolute zero, it appears a succession of closer and closer crossovers so that all the l-eigenvalues, characterized by an increasing l-value, successively play a role. But, when T ≈ 0 K, all these eigenvalues intervene due to the fact that the crossover temperatures are closer and closer. The discret eigenvalue spectrum tends to a continuum. As a result we can say that T = 0 K plays the role of critical temperature T c . This aspect will be more detailed later.
How interpreting this phenomena? In the 1d-case (infinite spin chain) we always have λ 0 (−βJ) as dominant eigenvalue in the whole range of temperature, the current integral F 0 being always equal to unity. In the 2d-case the situation is more complicated. The appearance of successive predominant eigenvalues is due to the presence of integrals F l,l ≠ 1, for any l > 0. A numerical fit shows that the ratio F l,l /F 0,0 increases with l according to a logarithmic law, more rapidly than the ratio |λ l (−βJ)/λ 0 (−βJ)| 2 which decreases when l increases for a given temperature. The particular case l → +∞ will be examined in a forthcoming section. Now, if taking into account the previous study, equation (22) can be rewritten as: where u max is given by equation (23) with l max = l i = l j = l and: As the ratios | , T are given by equation (14)) are positive and lower than unity S 1 (N,T) and S 2 (N,T) are absolutely convergent series.
In Appendix A.1 we have studied Z N (0) in the thermodynamic limit (N → +∞), for temperatures T > 0 K, in the whole range [0+ε,+∞[, with ε << 1. We show that, for a given range [ , (4 ) N N u π , with L= l max = l in equation (27). As the reasoning is valid for any [ , In the previous equation the special notation 0 t l ∞ + = ∑ recalls that the summation can be truncated due to the fact that each eigenvalue But, if considering the whole temperature range all the l-eigenvalues must be kept.
We must also note that, if dealing with a distribution of constant exchange energies J 1 and J 2 characterizing the horizontal and vertical lattice bonds, respectively, the infinite lattice can be described by the translation of these bonds along the horizontal and vertical axes of the lattice in the crystallographic space. As a result, if using a similar reasoning as the one used for expressing Z N (0), we can define a zero-field partition function per lattice site symbolically

Definitions
We first define the spin-spin correlation , , ' .
In the thermodynamic limit (29)). As a result the characteristic polynomial giving the numerator of , , ' .
The zero-field spin correlation .
The correlation function Γ k,k' is: if (i,j) is the site of reference. In this article we choose (0,0). In addition, due to the isotropic nature of couplings, we have , ' The general definition of the correlation length is: Along a horizontal lattice line k = 0 (x-crystallographic axis of the lattice) ξ = ξ x (respectively, k' = 0 and ξ = ξ y for a vertical lattice row, y-crystallographic axis of the lattice).
Using the general definition of the spin-spin correlation given by equation (31) and expanding the exponential part of the integrand on the infinite basis of spherical harmonics (cf equation (7)), we can write: for any site (i,j) (and a similar expression for site (i+k,j+k')). When , ( , ) (10)). Thus, if calculating 0,0 in the product of integrals appearing in equation (35).

Consequences
In this subsection we wish to calculate the numerator of the spin correlation with here m i,j = 0 and m' i,j = 0 for any site (i,j), in the thermodynamic limit (N → +∞). Then In the particular case , 0 i j l = which occurs at the beginning of each l-series expansion, we have 1 1 / 3 For the calculation of the spin correlation this transform can be equivalently applied to each of the four spherical harmonics appearing in ( , ) ' i j F given by equation ( with l i+1,j = l i,j−1 = l i,j = l' i,j = l. We immediately retrieve the calculation of integrals appearing in that of the zero-field partition function. As a result, if using the expansion of any product of two spherical harmonics given by equation (11) and their orthogonality condition in The non-vanishing condition of the current integral , which is due to that of the involved C. G. coefficients allows one to write down a universal temperature-independent selection rule concerning integers l (cf equation (19)). We This result rigorously proves that the critical temperature is absolute zero i.e., T c = 0 K.

Spin-spin Correlation Between Any Couple of Lattice Sites
In the thermodynamic limit (N → +∞) on which we exclusively focus, we are going to show that the spin-spin correlations can be derived from Z N (0). Notably the numerator show the same l-polynomial structure as Z N (0) which appears at the denominator. We have m k,k' = m' k,k' = 0, for any lattice site.
It means that, as for Z N (0), there also exist thermal crossovers between the new l-eigenvalues of the numerator. Their respective thermal domains of predominance are not necessary the same ones as those of eigenvalues appearing in Z N (0).  (31) and (35). We restrict the following study to k > 0 and k' > 0, without loss of generality. Due to the presence of cosθ 0,0 and cosθ k,k' appearing in the integrals F' (0,0) and F ' (k,k') which characterize the spin orientations at sites (0,0) and (k,k') we have to reconsider a new integration process. This process is similar to that one used for calculating Z N (0). It can be mainly carried out through two methods: It is useful to combine both methods. In a first step we choose method (i). In the dominant l-term the integrals F (i,j) involving sites located far from correlated sites (0,0) and (k,k') are characterized by a collection of integers l' i+1,j = l i,j−1 = l i,j = l' i,j = l for reasons explained in Subsec. 2.5. This part of the lattice constitutes the wing domain. When reaching the horizontal lattice lines i = 0 and k and the vertical ones j = 0 and k' whose respective intersections two by two are sites (0,0), (0,k'), (k,k') and (k,0) a special care must be brought. The inner domain defined by these two couples of lines is the correlation domain (see Figure 7).
A consequence of method (i) is that all the bonds located outside the correlation domain (or out-bonds) are characterized by the integer l, notably all the bonds linked to the frontier of the correlation domain. All the previous results are summarized in the following theorem: Theorem 1 (confinement theorem) In the thermodynamic limit, for calculating the numerator of the spin-spin correlation , respectively. We immediately retrieve the calculation of integrals appearing in that of the zero-field partition function. We express the products of pairs of spherical harmonics as C. G. series (cf equation (11)). For instance, if examining integrand 1, we can use the following combinations 1,0 1,0 1,0 Introducing these C. G. series in integral F' (0,0) given by equation (36) and using the orthogonality condition of spherical harmonics leads to L = L' i.e., notably L < = L' < and L > = L' > . Recalling that the characteristic polynomial associated with the numerator of where l' i+1,j = l i,j−1 = l i,j = l' i,j = l for any site (i,j) the unique solutions are l 0,0 = l ±1, l' 1,0 = l (case 1) and l 0,0 = l, l' 1,0 = l ±1 (case 2). All the other combinations between pairs of spherical harmonics lead to the same couple of solutions l 0,0 and l' 1,0 , for integrand 1 but also for integrands 2 and 3. From a mathematical point of view it also means that, for case 1 or 2, there are only two channels of integration leading to: if l > 0 a path beginning with a bond such as l 0,0 = l +1 and another one with l 0,0 = l −1 (case 1) or a path with l' 1,0 = l +1 and another one with l' 1,0 = l −1 (case 2).
(i) Case 1 (see Figure 7) We apply the decomposition law to the spherical harmonics    (48) and L > = min(2l, 2l + ε + ε'). In the previous equation As a result the corresponding contribution of site (0,0) to the numerator of As seen after equation (44) all the decompositions of products of spherical harmonics pairs as C.G. series only lead to two possible choices. If l' 1,1 = l the non-vanishing condition of F (0,1) (cf equation (48)) imposes l 0,1 = l 0,0 = l ±1 and the correlation path continues along the horizontal line i = 0 (case 1). If l 0,1 = l the non-vanishing condition of F (0,1) now imposes l' 1,1 = l ±1 and the correlation path continues along the vertical line j = 1: we then deal with a new correlation path called case 3 and detailed below (see Figure 7).
In summary, in case 1, we have two types of correlation path between sites (0,0) and (0,1): the bond is characterized by l+1 (F (0,1) = f l +1,l +1 ) or l−1 (F (0,1) = f l −1,l −1 ). This situation is similar for all the sites of line i = 0 i.e., between sites (0,1) and (0,k'−1). As a result the corresponding contribution to the Arriving at site (0,k') we have to determine l' 1,k' because l 0,k' = l' 0,k' = l due to the wing contribution and l 0,k'−1 = l ±1 due to the non-vanishing condition of integral F (0,k'−1) =f l ±1,l ±1 . That of integral F (0,k') gives l' 1,k' = l ±1 and F (0,k') = F (0,k'−1) = f l ±1,l ±1 . Then the work of integration is similar for the remaining sites of the vertical line j = k' between sites (1,k') and (k−1,k'), with for this later site l' k,k' = l ±1. The corresponding contribution to the numerator of If finally considering site (k,k') the integers l' k+1,k' and l k,k' have been already determined in the wing domain (l' k+1,k' = l k,k' = l) or along the correlation path (l' k,k' = l ±1). Concerning integral ( , ') ' k k ± F an independent study similar to that achieved at site (0,0) can be done. We have ( , ') is given by equations (45) In summary all the horizontal bonds of the correlation path are characterized by l 0,K' = l ±1 (0≤K'≤k'−1) between sites (0,0) and (0,k'−1) on the one hand and all the vertical bonds by l' K,k' = l±1 (1≤K≤k) between sites (0,k') and (k,k') on the other one. All the bonds of the correlation domain not involved in the correlation path are characterized by the integer l.
(iii) Case 3 (see Figure 7) This case is a mix of cases 1 and 2. At each current site (i,j) belonging to the correlation domain we can have l i,j = l ±1, l' i+1,j = l or l' i+1,j = l ±1, l i,j = l.
Due to the fact that it is impossible to go backwards when the integration process has been carried out any loop of the correlation path is forbidden. In addition the expression of the numerator of the spin-spin correlation is independent of the bond orientation chosen for the integration process i.e., between sites (0,0) and (k,k') or vice versa.
We conclude that (i) all the correlation paths contain the same number of hori zontal and vertical bonds; as a result all these paths show the same expression i.e., all the spin-spin correlations show a unique expression, as expected for this kind of lattice; (ii) these correlation paths correspond to the shortest possible length through the bonds involved between sites (0,0) and (k,k'); their total number is simply n = ' ' ; thus there are n l+1 = n paths whose bonds are characterized by the integer l+1 and n l−1 = n paths showing bonds characterized by l−1; they have the same weight w l+1 = w l−1 = n l±1 /2n = 1/2 (l > 0). In the particular case l = 0 w 1 = 1, w −1 = 0. Theorem 2 As loops are forbidden for all the correlation paths these paths have the same length inside the correlation domain. This length is the shortest possible one through the lattice bonds between any couple of correlated sites. Each path respectively involves the same number of horizontal and vertical bonds as the horizontal and vertical sides of the correlation rectangle, for a 2d-infinite square lattice.
We proceed as for the thermal study of the current term of the l-polynomial expansion of Z N (0) where we have defined a crossover temperature T CO through equation (26)  in the simplest case J = J 1 = J 2 without loss of generality. It can be extended to the case J 1 ≠ J 2 . As for equation (26) we have numerically solved this transcendental equation. We find that |P l−1 | ≤ 1 if T ≤ T CO,l−1 and |P l−1 | > 1 if T > T CO,l−1 on the one hand but |P l+1 | ≥ 1 if T ≤ T CO,l+1 and |P l+1 | < 1 if T > T CO,l+1 on the other one. Thus, this is the competition between the smooth decreasing l-law of the ratio f l+1,l+1 /F l,l > 1 (see Figure 8) which tends towards unity when l → +∞ (i.e., when T tends to T c = 0 K) and the T-law of the ratio λ l+1 (β|J|)/λ l (β|J|) involved in |P l+1 | which is responsible of such a crossover. For |P l−1 | the competition is between the increasing l-law of the ratio f l−1,l−1 /F l,l < 1 and the T-law of the ratio λ l−1 (β|J|)/λ l (β|J|) > 1.
Finally, owing to the numerical study reported in Figure 9, we have T CO < T CO,l+1 < T CO,l−1 . As a result it becomes possible to determine the new domains of thermal predominance [ 1, l T ± < , 1, l T ± > ] of the eigenvalues |P l±1 |. Their detailed classification is out of the framework of the present article. When l → +∞ i.e., as T approaches T c = 0 K, all the leigenvalues become equivalent but a common limit can be selected.

Correlation Length
The correlation length can be derived owing to equation (34). We have with on condition that , 1 i l z l , P i,l+ε (i = 1,2) K l+ε and the integral F l,l appearing in z l are respectively given by equations (50) and (20) in which l i = l j = l, m = 0. K l+ε → 1 as l → +∞.
In equation (54), if considering the spin-spin correlation between first-nearest neighbors derived from equation (49) and setting It means that the spin-spin correlation between firstnearest neighbors plays a fundamental role near the critical point. This is a hidden consequence of equation (51) itself derived from the classical character of spin momenta (cf equation (6)).

Preliminaries
For sake of simplicity we again reduce the study to the simplest case J = J 1 = J 2 without loss of generality. We examine the low-temperature behavior of the ratios P i,l+ε involved in the spin-spin correlation (cf equations (49), (50)), with here P i,l+ε = P l+ε (i = 1,2).
We first consider the ratio f l+ε,l+ε /F l,l (ε=±1) where integrals F l,l and f l+ε,l+ε are respectively given by equation (48) with ε = ε' = 0 for F l,l . The l-behavior of each of these ratios has been reported in Figure 8. f l+1,l+1 /F l,l > 1 and f l−1,l−1 /F l,l < 1 but f l±1,l±1 /F l,l → 1 as l→+∞ i.e., near the critical temperature T c =0 K. Indeed, if expressing the spherical harmonics involved in the definition of integral F l,l given by equation (10) in which m = 0, we have in the infinite l-limit [26] ( ) ,0 and the exact asymptotic result: Then, if taking into account equations (8) and (50), P l±1 behaves as  Intuitively, in the low-temperature limit, we must consider the three cases β|J| >> l, β|J| ∼ l and β|J| << l. The behavior of the Bessel function ( ) l I J −β in the double limit l → +∞ and β|J| → +∞ has been established by Olver [27]. In previous papers [24] we have extended this work to a large order l (but not necessarily infinite) and to any real argument β|J| varying from a finite value to infinity. The study of the Bessel differential equation in the large l-limit necessitates the introduction of the dimensionless auxiliary variables: Thus |ζ| plays the role of the inverse of an effective dimensionless coupling near T c . The correspondence with the writing of Chakravarty et al. is |ζ| = t −1 [6]. The numerical study of |ζ| is reported in Figure 10. We observe that there are two branches. |ζ| vanishes for a numerical value of |z 0 | −1 very close to π/2 so that there are 3 domains which will be physically interpreted in next subsection. Let T 0 be the corresponding temperature In the formalism of renormalization group T 0 is called a fixed point. In the present 2d case we have l → +∞. We then derive that T 0 → T c = 0 K as l → +∞ so that the critical temperarature can be seen as a non trivial fixed point. In other words it means that all the thermodynamic functions can be expanded as series of current term |T − T 0 | near T 0 ≈ T c = 0 K, in the infinite l-limit. Finally Figure 10 is nothing but the low-temperature diagram of magnetic phases and |ζ| defined by equation (62) gives the analytic expression of branches.
For convenience, we introduce the dimensionless coupling constant g at temperature T as well as its reduced value g : (65) g measures the strength of spin fluctuations. g is a universal parameter and is l-independent. At the critical point T 0 = T c we have g = 1. Owing to equation (64) the critical coupling g c can be written as: (66) Chubukov et al. have found that, at the critical temperature T c , the critical coupling is g c = 4π/Λ where Λ = 2π/a is a relativistic cutoff parameter (a being the lattice spacing) [10].
In the first case g c l = π/2 is obtained by a pure numerical method because this is the zero of the function |ζ| = f(1/|z|) given by equation (62) with 1/|z| = gl. The second case g c Λ = 4π supposes to take into account the volumic density of g c . Ac-cording to Chubukov et al., for D = 3, g c is such as where P = (ℏk,ℏω/c) is the relativistic momentum associated with the spin wave of wave vector k and energy ℏω. First we examine the problem of volumic density. As we focus on the static aspect i.e., the volume available to the spin momentum, the relativistic momentum P/ℏ reduces to S (in ℏ-unit This is precisely the result derived by Chakravarty et al. as the mathematical solution of the recursion relations established between g and t = |ζ| −1 through a one-loop renormalization process [6]. Thus g c = 0 when d = 1 as expected and g c = 4π when d = 2. Second, for obtaining a relation between l and Λ, we must also consider the volume of the lattice unit cell of spacing a. For sake of simplicity we restrict to the case D = 3. In equation (66) we introduce the density of exchange energy |J|S(S+1)/a 2 ∼ |J|S 2 /a 2 = |J|/a 2 as S = 1 in our conventional writing. In the associated D-space-time (D = 3) the volume of the phase space is V φ = V S a 3 = a 3 /4π. It means that the value of g per phase space volume is g/ V φ = (4π/a 3 )g.
As a result g = g c l given by equation (66)  Thus the l-index of the Bessel functions appears as the ratio of two different scales of reference Λ and Λ * , respectively associated with the lattices of spacing a and a * = a/2l << a.
We introduce: (i) the thermal de Broglie wavelength λ DB ; (ii) the low-temperature spin wave celerity By definition we must have λ DB >> a (71) i.e., Λ = 2π/a >> 2π/λ DB or equivalently with Λ * = 2lΛ Thus Λ −1 (or Λ *−1 ) appears as a short distance cutoff. No such intrinsic cutoff exists for the imaginary variable τ. At the critical point T c = 0 K λ DB → +∞ (as well as L τ ) and spins are strongly correlated. For T > T c λ DB and L τ become finite and diminish as the spin-spin correlation magnitude when T increases. The adequate tool for estimating this correlation between any couple of spins is the correlation length ξ. As a result λ DB (or L τ ) appears as the good unit length for measuring ξ.
Under these conditions we generalize the application of the cutoff parameter Λ. |z| defined by the second of equation (62) can be rewritten if using equations (65) and (69) spacing a (i.e., in the Λ-scale) is ξ Λ = ξ 2 2 a (or ξ Λ Λ=

Low-temperature Behaviors of the Spin-Spin Correlation
As previously seen (cf equations (49), (50)) the spin-spin correlation is expressed owing to the ratios At first sight all these quantities seem to strongly depend on Λ * whereas they must show a universal behavior near the critical point. As a result scaling parameters i.e., parameters which are Λ * -independent must be introduced so that the spin-spin correlation as well as the correlation length are scale-independent near T c = 0 K, as expected.
Chakravarty et al. [6] have introduced the physical parameters ρ s and ∆ defined as: In the 2d-case ρ s and ∆ have the dimension of an energy JS 2 (in our case J). ρ s is the spin stiffness of the ordered ground state (Néel state for an antiferromagnet) and ∆ is the T=0-energy gap between the ground state and the first excited state. In the framework of the classical spin approximation the spectrum is quasi continuous. In our case it means that ∆ is very small.
At the critical point T c = 0 K g = 1: ρ s and ∆ vanish and, near critically, we have ρ s << |J| and ∆ << |J| where |J| finally appears as the bare value of ρ s and ∆ i.e., their value at 0 K.
For all the previous reasons we are then led to introduce the following parameters: where the factor 4π appears in ∆ for notational convenience.
As a result we can define the following scaling parameters: where the factor 2π also appears for notational convenience.
As ρ s and ∆ vanish at T 0 = T c , x 1 and x 2 become infinite at this fixed point. They are scaling parameters as well as * * c | | / z z . From a physical point of view and as noted by Chakravarty et al. [6] as well as by Chubukov et al. [10], these parameters control the scaling properties of the magnetic system.
There is an analytical continuity between x 1 and x 2 when T 0 = T c . As a result there are only 3 domains of predominance: x 1 << 1 (T < T c and |ζ|/4π < 1− g , Zone 1) i.e., ρ s >> k B T, x 2 << 1 (T > T c and |ζ|/4π < g − 1, Zone 2) i.e., ∆ >> k B T; finally x 1 >> 1 (T < T c and |ζ|/4π > 1− g , Zone 3) i.e., ρ s << k B T and x 2 >> 1 (T > T c and |ζ|/4π > g − 1, Zone 4) i.e., ∆ << k B T. Along the line T = T c ( g = 1), we directly reach the Néel line (see Figure 11). Each of these domains previously described corresponds to a particular magnetic regime. The physical meaning of each regime can be derived from the low-temperature study of the ratio The corresponding behaviors are reported in Figure 12. The asymptotic expansions of * * | | Λ ζ can then be derived for the four zones of the magnetic diagram. We have  ( ) The ratio α = (1 5 ) / 2 + is the golden mean. As a result, starting from a closed expression of |ζ| given by equation (62) We tend towards an assembly of noncorrelated spins. In Zones 3 (x 1 >> 1) and 4 (x 2 >> 1) . 0,0 0,1 < > S S → 1 as in Zone 1. As a result this is the low-temperature behavior of the correlation length ξ which is going to allow the characterization of the magnetic order nature.

Low-temperature Behaviors of the Correlation Length
If using the expression of the correlation length given by equation (53) as well as the scale invariance property near T c = 0 K given by equation (77) the measure ξ of ξ x and ξ τ is such as where a is the lattice spacing as well as a * =a/2l; |ζ * |Λ * is defined by equations (81) We exactly retrieve the result first obtained by Hasenfratz and Niedermayer [15] and confirmed by Chubukov et al. [10]. This characterizes the Renormalized Classical Regime We deal with the Quantum Disordered Regime (QDR) characterized by ∆ >> k B T. Owing to equation (80) we have ∆ = k B T/x 2 so that ξ τ ≈ L τ x 2 << L τ as x 2 << 1. Equivalently, due to the fact that L τ = B / c k T ℏ and a = ℏc/2|J| we have ξ x ≈ 2ax 2 << 2a: we then pass from no T=0-order to a short-range order when T increases. The magnetic structure is made of spin dimers or aggregates of spin dimers organized in Kadanoff blocks of small size ξ τ that we can assimilate to blobs weakly interacting between each others. We deal with a spin fluid. The detailed study is out of the framework of the present article. From a formal point of view it is often phrased in term of Resonating Valence Bonds (RVB) between pairs of quantum spins (considered here in the classical spin approximation) [29,30]. In Zones 3 (x 1 >> 1) and 4 (x 2 >> 1) we have have ξ QDR < ξ QCR < ξ RCR (cf Figure 14). Thus Kadanoff blocks show a smaller size when passing from Zone 1 to Zone 2 through Zones 3 and 4. As a result each behavior of the correlation length characterizes a magnetic regime. All the predominance domains of these regimes are summarized in Figure 13. At the frontier between Zones 3 (ρ s << k B T) and 4 (∆ << k B T) i.e., along the vertical line reaching the Néel line at T c , x 1 and x 2 become infinite so that: i.e., τ ξ ≈ L τ as C −1 is close to unity, as predicted by the renormalization group analysis [6,10]. As τ ξ diverges according to a T −1 -law the critical exponent is: in the D-space-time. The low-temperature behaviors of ξ have been reported in Figure 14.
Owing to previous results the correlation length can also be written as where . 0,0 0,1 < > S S is the renormalized spin-spin correlation between first-nearest neighbors (0,0) and (0,1) as J = J 1 = J 2 , expressed near T c = 0 K. We retrieve the result predicted in: This result is exclusively valid for 2d magnetic systems characterized by isotropic couplings because the correlation function reduces to the spin-spin correlation.
As a result it becomes possible to characterize the nature of magnetic ordering owing to the T-decreasing law derived from equation (98) 1 . 0,0 0,1 where u recalls the nature of the magnetic regime: u = RCR (Zone 1), QCR (Zones 3 and 4) and QDR (Zone 2). As previously seen we have ξ QDR < ξ QCR < ξ RCR so that Thus, in Zone 1 (Renormalized Classical Regime), we have a strong long range order in the critical domain whereas in Zones 3 and 4 (Quantum Critical Regime) the magnitude of magnetic order is less strong. In Zone 2 (Quantum Disordered Regime) we deal with a very short magnitude characteristic of a spin fluid.

Conclusion
In this paper, if restricting the study of the two-dimensional Heisenberg square lattice composed of classical spins to the physical case of the thermodynamic limit, we have obtained the exact closed-form expression of the zero-field partition function Z N (0) valid for any temperature.
The thermal study of the basic l-term of Z N (0) has allowed to point out a new phenomenon: thermal crossovers between two consecutive eigenvalues. When T → 0, l → +∞. As a result all the successive dominant eigenvalues become equivalent so that T c = 0 K.
In addition we have exactly retrieved the low-temperature diagram of magnetic phases already obtained through a renormalization group approach [6,10].
If using a similar method employed for expressing Z N (0) we have derived an exact expression for the spin-spin correlations and the correlation length ξ valid for any temperature.
The T=0-limit of ξ shows the same expression as the corresponding one obtained through a renormalization process but exclusively valid near the critical point T c = 0 K, for each low-temperature regime, with the good critical exponent ν = 1, thus validating the closed-form expressions obtained for Z N (0), the spin-spin correlations and the correlation length, respectively. Appendix