format: text

text:
# Monte Carlo study of the scaling of universal correlation lengths in three-dimensional O(𝑛) spin models

## I Introduction

The concept of scaling, the observation that singular observables vary in a scale-free manner according to power laws when the driving parameter of a transition (temperature, magnetic field, …) is tuned towards a critical point, has since the first observations been a key ingredient of the theory of critical phenomena zinn-justin ; kadanoff:domb . Exploiting the symmetry of scale-invariance, forming the geometrical basis for the power-law behavior in the vicinity of a critical point, through the idea of real-space renormalization, scaling theory can be mapped on the behavior of finite systems near the transition point of the bulk system in the limit of diverging system sizes, the thermodynamic limit. This finite-size scaling (FSS) barber:72 ; fisher:74a ; barber:domb  occurs with scaling exponents generically linked to the exponents that govern scaling in the bulk system. Thus the apparent weakness of finite system size that hampers simulational approaches actually turns out to be their intrinsic strength, when exploring FSS means exploring thermal scaling binder:schlad ; bd:como .

The significance of scaling theory for the understanding of critical phenomena becomes quite exposed in the context of conformal field theory (CFT) for two-dimensional systems belavin:84a . In the course of exploiting the additional invariances of conformal symmetry one is able to split the critical point partition function of a lattice system into a sum over contributions from all the scaling variables present in a specific model. Consider a critical system on a $`L\times L^{}`$ lattice with toroidal boundary conditions; then the partition function decomposes as cardy:86a ; cardy:86b :

$$Z(L,L^{})=e^{fA+\pi c\delta /6}\underset{n}{}e^{2\pi x_n\delta },$$
(1)
where $`c`$ is the central charge of the considered theory, $`f`$ the bulk free energy per unit volume, $`\delta =L^{}/L`$, $`A=LL^{}`$, and the sum runs over the whole content of scaling operators with dimensions $`x_n`$. Thus, the knowledge of the operator content of a theory in connection with the corresponding scaling dimensions is equivalent to an “exact” solution of the model on finite lattices.

### Two Dimensions –

A particular example of a scaling relation in two dimensions that can be derived assuming conformal invariance of critical point entities concerns the two-point function in the limit of $`L^{}\mathrm{}`$. It it generally sufficient to assume translational, rotational, dilatational, and inversional invariance to imply conformal invariance henkel:book ; homogeneity, isotropy and scale invariance alone suffice to uniquely fix the critical, connected two-point function of an operator $`\varphi `$ in the infinite plane up to an overall normalization factor:

$$\varphi (z_1,\overline{z}_1)\varphi (z_2,\overline{z}_2)_c=(z_1z_2)^x(\overline{z}_1\overline{z}_2)^x,$$
(2)
where $`z_1`$, $`z_2`$ are complex co-ordinates parametrizing the plane. Then, one uses the logarithmic map

$$w=\frac{L}{2\pi }\mathrm{ln}z,z$$
(3)
to wrap the complex plane around an infinite length cylinder $`S^1\times `$ of circumference $`L`$ with co-ordinates $`w=u+iv`$, where $`v`$ measures the polar angle along $`S^1`$ and $`u`$ the longitudinal direction along $``$. Assuming conformally covariant transformation behavior of the (primary) operator $`\varphi `$, one arrives at an expression for the two-point function on the cylinder cardy:84a :

$`\varphi (w_1,\overline{w}_1)\varphi (w_2,\overline{w}_2)_c=\left({\displaystyle \frac{2\pi }{L}}\right)^{2x}\left({\displaystyle \frac{\left|z_1z_2\right|}{\left|z_1z_2\right|^2}}\right)^x=`$
$`\left({\displaystyle \frac{2\pi }{L}}\right)^{2x}\left(2\mathrm{cosh}{\displaystyle \frac{2\pi }{L}}(u_1u_2)2\mathrm{cos}{\displaystyle \frac{2\pi }{L}}(v_1v_2)\right)^x.`$ (4)
In the limit of large longitudinal distances $`|u_1u_2|L`$ and $`v_1=v_2`$, one is left with a purely exponential drop with a correlation length

$$\xi _{}=\frac{L}{2\pi x}.$$
(5)
Thus, utilization of conformal invariance yields a finite-size scaling relation including the amplitude, which is in contrast to renormalization group theory that usually gives the scaling exponents and only certain amplitude ratios, but not the amplitudes themselves. Since this result emerges from a field-theoretic description of statistical mechanics that does not take into account the microscopical details of the system, it is expected to be universalgriffiths:70a . Note, however, that this proposed universality goes beyond the usual notion of an universal quantity and comprises three different aspects: (i) the correlation length of a given operator should be the same within the associated universality class of models; (ii) when looking at different operators, on the other hand, the form of Eq. (5) should be left unchanged, all operator-dependent information being condensed in the scaling dimension $`x`$; (iii) finally, even when looking at models of different universality classes, all that should change are the scaling dimensions (and the definition of $`\varphi `$), the validity of Eq. (5) being untouched. Property (i) implies the “hyperuniversality” relation of Privman and Fisher privman:84a . In the following, we will refer to the whole extent of aspects (i)-(iii) exceeding the usual notion of universality with the term “hyperuniversal”. A corollary that is of importance for transfer matrix calculations that use an unnormalized (quantum) Hamiltonian results from taking the ratio of the correlation lengths of two primary operators, for example the densities of magnetization and energy which are usually primary for spin models:

$$\frac{\xi _\sigma }{\xi _ϵ}=\frac{x_ϵ}{x_\sigma }.$$
(6)
Because of the independence from the overall amplitude $`1/2\pi `$ of Eq. (5) this relation might still stay valid when changing the geometry in a way such that only this overall amplitude is altered. In terms of universality this constitutes a weaker form of the aspect (i) above, namely universality of amplitude ratios instead of amplitudes themselves; we will refer to this weaker property as (i’) in the following.

A suitable test-bed for these general field-theory results is, of course, given by the exactly solvable two-dimensional Ising model. Using Eq. (5) and the generic relations between scaling dimensions and the conventional critical exponents:

$$x_\sigma =\frac{\beta }{\nu },x_ϵ=\frac{1\alpha }{\nu },$$
(7)
giving $`x_\sigma =1/8`$ and $`x_ϵ=1`$ for the two-dimensional Ising model, one arrives at a ratio $`x_ϵ/x_\sigma =8`$. A direct evaluation of the spin-spin correlation length in the Onsager-Kaufman framework gives, as the leading term in the scaling series, $`\xi _\sigma =4L/\pi L/(2\pi \frac{1}{8})`$, in agreement with the CFT result domb:60 ; nightingale:76a ; derrida:82a . The same holds true for the leading scaling amplitude of the energy-energy correlation function bloete:83a , $`\xi _ϵ=L/2\pi `$. Both amplitudes have also been evaluated numerically to high precision in a Monte Carlo (MC) study diplom , resulting in perfect agreement with Eq. (5).

A possible alteration of the $`S^1\times `$ situation, namely changing the boundary conditions along the $`S^1`$-direction from periodic to antiperiodic has also been treated within the CFT framework, exploiting the fact that in the case of the ferromagnetic nearest-neighbor Ising model the antiperiodic boundary corresponds to the insertion of a seam of antiferromagnetic bonds along this boundary line. This calculation yields cardy:84b ; cardy:domb :

$$\xi _\sigma =\frac{4}{3\pi }L,\xi _ϵ=\frac{1}{4\pi }L,$$
(8)
again in good agreement with Monte Carlo data diplom . Note, however, that this last relation, in contrast to Eq. (5), is specific to the Ising model and the special choice of the densities of magnetization and energy as operators and thus is not “hyperuniversal” in the sense of properties (ii) and (iii) presented above.

The amplitude-exponent relation Eq. (5) for two-dimensional systems has been checked analytically or numerically and found valid for an impressive series of further models like the Potts model and its percolation limit derrida:82a , the XY model luck:82a , the symmetric eight-vertex model bloete:83a , and quantum spin models penson:84a  to name only the most prominent.

### Three Dimensions –

On leaving the domain of two-dimensional systems towards higher dimensions, the wealth of exact field theoretic calculations is instantly reduced to severe scarcity. The conformal group coincides with the set of holomorphic functions in the special case of spatial dimension $`d=2`$ and is thus infinite-dimensional as a group. For $`d3`$, unfortunately, it reduces to a simple Lie group with dimension $`D(d+1)(d+2)/2`$ for any Riemannian, connected manifold. As a consequence, only in two dimensions the postulate of conformal invariance is restrictive enough for a classification of the operator contents of the different universality classes and thus an exact solution of the critical theories within the limits of field-theory assumptions. For $`d3`$, on the other hand, the implications of the finite-dimensional conformal-group symmetry reach hardly beyond the consequences of plain renormalization group theory exploiting dilatational invariance. However, since inversional symmetry is still present, a transformation like Eq. (3) stays conformal in higher dimensions, now connecting the spaces $`^d`$ and $`S^{d1}\times `$. Applied to the two-point function one arrives at a scaling relation analogous to Eq. (5), namely $`\xi _{}=R/x`$, cp. Ref. cardy:85a , which contains the $`d=2`$ result as a special case assuming $`L=2\pi R`$, $`R`$ being the radius of $`S^{d1}`$. Since primarity of operators is a priori not well defined for $`d3`$, it is, however, unclear for which operators this relation should hold. A numerical analysis for this geometry, which has to cope with the fact that $`S^{d1}`$ for $`d3`$ is a truly curved space and thus hard to regularize by discrete lattices, will be presented in a separate publication prep .

On the other hand, the toroidal geometry $`S^1\times \mathrm{}\times S^1\times `$, which is much more convenient for numerical simulations, is not conformally flat and thus no CFT predictions exist for this case. In spite of this theoretically unfavorable situation a transfer matrix calculation for the Hamiltonian limit of the three-dimensional Ising model on the geometry $`S^1\times S^1\times T^2\times `$ by Henkel henkel:86a ; henkel:87a ; henkel:87b  rendered results still comparable to the situation for the $`S^{d1}\times `$ geometry. For the ratios of leading scaling amplitudes of correlation lengths for different boundary conditions (bc) he found

$$\begin{array}{ccc}\hfill \xi _\sigma /\xi _ϵ& =& 3.62(7)\text{for periodic bc,}\hfill \\ \hfill \xi _\sigma /\xi _ϵ& =& 2.76(4)\text{for antiperiodic bc.}\hfill \end{array}$$
(9)
A comparison with the (inverse) ratio of the corresponding scaling dimensions,

$$x_ϵ/x_\sigma =\frac{(1\alpha )/\nu }{\beta /\nu }=\frac{2(\nu d1)}{\nu d\gamma }=2.7264(13),$$
(10)
(cp. Table 1 and Eq. (7)) showed that even though the original expectation to possibly find agreement in the case of periodic boundary conditions as in the two-dimensional case was not met, the data are consistent with the relation Eq. (6) for the unorthodox case of antiperiodic boundary conditions. Note that one has to compare ratios in this case, because the quantum Hamiltonian used in the calculation is defined only up to on overall normalization constant. This result is in qualitative agreement with a Metropolis MC simulation by Weston weston:90 , who found ratios $`\xi _\sigma /\xi _ϵ`$ of about $`3.7`$ for periodic and $`2.6`$ for antiperiodic boundary conditions, respectively. Considering these striking observations it seems interesting to check whether this behavior is just a coincidence or special feature of the Ising model or instead indicates a general property of critical models on this special three-dimensional geometry.

The rest of the paper is organized as follows. In Sec. II we introduce the general class of models we want to examine and present the way we are going to discretize the three-dimensional geometry $`T^2\times `$. We discuss simulation methods, observables, estimators for measurements and parameters of the simulations. In Sec. III we outline the statistical tools used for the data analysis. It is quite hard to extract high-precision information about correlation lengths from MC simulation data; we will thus discuss the special path of data analysis we are going to proceed along and present details of the statistical tools used there for. This tool-set is “calibrated” with simulations of the two-dimensional Ising model, where exact results for comparison are available. In Sec. IV we discuss the results for the correlation lengths ratios of our simulations for the Ising, XY and (generalized) Heisenberg models. Our results, already briefly announced in Ref. prl:99a , confirm Henkel’s findings on a high level of accuracy. Furthermore this behavior seems to carry through for the whole class of O($`n`$) spin models and is thus far from being a “numerical accident”. In Sec. V we try to rank our numerical findings in the context of the classification of universality presented above. The type of the model considered enters not only via a variation of the scaling dimensions, but also influences the overall prefactor of Eq. (5). Sec. VI is devoted to the discussion of the relation of our finite-$`n`$ results to the spherical model, which is connected to the limit $`n\mathrm{}`$ of the class of O($`n`$) spin models. The classic identification of both models seems to break down as soon as (multi-point) correlation functions are considered. The final Sec. VII contains our conclusions.

## II Models and Simulation

Throughout this paper we consider classical, ferromagnetic, zero-field, nearest-neighbor, O($`n`$) symmetric spin models with Hamiltonian

$$=J\underset{ij}{}𝝈_i𝝈_j,𝝈_iS^{n1}.$$
(11)
The underlying lattice is taken to be simple cubic with dimensions $`L_x\times L_y\times L_z`$. Special cases of this class of models include the Ising ($`n=1`$), XY ($`n=2`$), and Heisenberg ($`n=3`$) models. This Hamiltonian has the advantage of representing a whole class of models with critical points in three dimensions, tuned by the single parameter $`n`$. According to the $`T^2\times `$ geometry we set $`L_x=L_y`$ and apply periodic or antiperiodic boundary conditions in the $`x`$ and $`y`$ directions. In both cases we use periodic boundary conditions in the $`z`$-direction to eliminate surface effects that are also absent in the $`L_z\mathrm{}`$ case assumed in Eq. (4). To reduce effects of finite size in $`z`$-direction one has to ensure that $`L_z\xi _{}`$, a concrete rule will be given below.

In view of the problem of critical slowing down, we use the Wolff single cluster update algorithmwolff:89a  for all O($`n`$) model simulations, cp. prl:99a . The adaption of this update procedure to the case of antiperiodic boundary conditions along the torus directions is straightforward if one exploits the above mentioned equivalence of an antiperiodic boundary to the insertion of a seam of antiferromagnetic bonds along the boundary line for the case of nearest-neighbor interactions. Considering the Ising model or, alternatively, embedded Ising spins for $`n>1`$ models wolff:90a , this means: adjacent spins interacting antiferromagnetically are connected with a bond obeying the Swendsen-Wang probability $`p=1\mathrm{exp}(2\beta J)`$ in case of opposite orientation and are left unbonded in case of identical orientation. This rule exactly reflects the change in energy compared to the ferromagnetic case and thus trivially satisfies detailed balance.

The main observables of our simulations are the connected correlation functions of the densities of magnetization and energy:

$$\begin{array}{ccc}\hfill G_\sigma ^c(𝐱_1,𝐱_2)& =& 𝝈(𝐱_1)𝝈(𝐱_2)𝝈𝝈,\hfill \\ \hfill G_ϵ^c(𝐱_1,𝐱_2)& =& ϵ(𝐱_1)ϵ(𝐱_2)ϵϵ.\hfill \end{array}$$
(12)
We define the energy density as a local sum over the nearest neighborhood $`𝐱^{}`$ of a spin $`𝐱`$ ($`𝐱^{}\mathrm{nn}𝐱`$):

$$ϵ(𝐱)=\frac{J}{2}\underset{𝐱^{}\mathrm{nn}𝐱}{}𝝈(𝐱)𝝈(𝐱^{}),$$
(13)
the factor $`1/2`$ ensuring that $`E=_𝐱ϵ(𝐱)`$. It is straightforward to construct a bias-reduced estimator for the case of $`(𝐱_2𝐱_1)\widehat{e}_z`$, corresponding to the correlation length $`\xi =\xi _{}`$: first, taking advantage of the translation invariance of the systems along the $`z`$-axis established by a periodic boundary, one can average over the “layers” $`i|z_2z_1|=\text{const}`$. To improve on that consider a “zero-mode projection” wj:93a , i.e. define layered variables

$$\overline{𝒪}_t(z)=\frac{1}{L_xL_y}\underset{𝐱^{},z^{}=z}{}𝒪_t(𝐱^{}),$$
(14)
where $`𝒪_t=𝝈_t`$ or $`ϵ_t`$ denotes the times series of MC measurements, and consider the estimator

$`\widehat{G}_𝒪^{c,}(i)`$ $`=`$ $`{\displaystyle \frac{1}{T}}{\displaystyle \underset{t=1}{\overset{T}{}}}{\displaystyle \frac{1}{L_z}}{\displaystyle \underset{|z_2z_1|=i}{}}\overline{𝒪}_t(z_1)\overline{𝒪}_t(z_2)`$ (15)
$`\left({\displaystyle \frac{1}{TL_z}}{\displaystyle \underset{t=1}{\overset{T}{}}}{\displaystyle \underset{z}{}}\overline{𝒪}_t(z)\right)^2,`$
where $`T`$ denotes the length of the MC time series. This estimator obviously does not directly measure $`G^{c,}`$, but inspecting the continuum form Eq. (4) reveals that the deviation stemming from transversal cross-correlations entering the estimator declines exponentially with increasing longitudinal distance $`i`$ and thus becomes irrelevant for the long-distance behavior we are interested in. Numerical investigations confirm that these considerations stay correct when passing to three dimensions diplom . In the large-distance regime zero-mode projection reduces the variance of correlation function estimates by a factor inversely proportional to the layer volume $`L_xL_y`$. Note that the given estimator for the disconnected part $`𝒪^2`$ has a bias that vanishes as $`1/T`$ in the large-$`T`$ limit.

As mentioned above, periodic boundary conditions in $`z`$-direction eliminate surface effects associated with this direction, but still effects of finite $`L_z`$ will trigger deviations from the $`L_z\mathrm{}`$ limit assumed in Eq. (5). Inspecting the form of Eq. (4) in the limit of distances $`i\xi _{}`$ one expects longitudinal correlations according to:

$$G^{c,}(i)e^{i/\xi _{}}+e^{(L_zi)/\xi _{}},$$
(16)
i.e. the exponential decay is superimposed by an exponentially increasing part. Thus, using too small values of $`L_z`$ results in an effective underestimation of correlation lengths. In order to satisfy $`L_z\xi _{}`$ in a systematic way, i.e. to keep this effect away from the region of clear signal for measuring the correlation lengths, and assuming linear scaling of correlation lengths according to $`\xi _{}=AL_x`$, one has to keep the ratio $`L_z/\xi _{}=L_z/AL_x`$ fixed and therefore has to scale $`L_z`$ proportionally to $`L_x`$. Simulations for the case of the two-dimensional Ising model show that these finite-size effects are negligible compared to the statistical errors for $`L_z/\xi _{}10`$ and lengths of time series of about $`10^6`$ to $`10^7`$ measurements diplom . Adding a safety margin the longitudinal system sizes for the simulations in three dimensions where chosen such that $`L_z/\xi _{}15`$, the scaling amplitude $`A`$ being estimated from a simulation of an “oversized” system. Since $`\xi _\sigma >\xi _ϵ`$ for all models under consideration, the amplitude $`A_\sigma `$ of the spin-spin correlation length scaling is significant for the satisfaction of this condition. Note that from Eq. (15) increasing $`L_z`$ also has the positive side effect of improving the statistics of the correlation function estimation.

In order to judge the efficiency of the used cluster update algorithm and to ensure reasonable usage of computer time we evaluated integrated autocorrelation times $`\tau _{\mathrm{int}}`$ using a binning technique flyvbjerg:89a . The strong asymmetry of the model lattices reduces the average size of clusters and thus Wolff’s cluster update algorithm does not perform as good as on (hyper-)cubic lattices, resulting in increased autocorrelation times. Since measurements of $`\widehat{G}^{c,}`$ are computationally expensive compared to update steps, but the statistical gain vanishes with increasing $`\tau _{\mathrm{int}}`$, measurements were done with frequencies of about $`1/\tau _{\mathrm{int}}`$. Approaching the low-temperature phase, antiperiodic boundary conditions in the torus directions produce a spatially stable boundary of the geometric clusters along the antiferromagnetic seam, which in turn enforces a second boundary along the $`z`$ direction. This results in a further reduction of the average cluster size compared to the periodic boundary case. Fig. 1 shows typical configurations for the case of the (two-dimensional) Ising model.

## III Data Analysis

Having sampled correlation functions according to Eq. (15) and assuming the functional form $`G^{c,}(i)=a\mathrm{exp}(i/\xi _{})+b`$, we refrain from using instrinsically unstable non-linear three-parameter fits and resort to the following estimator instead,

$$\widehat{\xi }_𝒪(i)=\mathrm{\Delta }\left[\mathrm{ln}\frac{\widehat{G}_𝒪^{c,}(i)\widehat{G}_𝒪^{c,}(i\mathrm{\Delta })}{\widehat{G}_𝒪^{c,}(i+\mathrm{\Delta })\widehat{G}_𝒪^{c,}(i)}\right]^1,$$
(17)
which eliminates the additive and multiplicative constants $`a`$ and $`b`$ above. Note that it is not allowed to assume $`b=0`$ a priori for time series of finite length, cp. Ref. prl:99a . Apart from stability considerations this approach allows for computational simplifications, since correlation functions can be sampled irrespective of normalization and the biased estimation of the disconnected part $`𝒪^2`$ can be dropped. In addition, Eq. (17) simplifies the distinction of the long-distance part of the correlation function from the short-distance region: as the explicit two-dimensional expression Eq. (4) implies, exponential decay will only occur asymptotically, but with deviations decaying themselves exponentially; apart from that, lattice artefacts that are not reflected in the continuum expression Eq. (4) additionally distort the short-distance behavior. Fig. 2 shows an example plot of the spin-spin correlation length estimates $`\widehat{\xi }_\sigma (i)`$ for the Ising model. The transition from the short-distance region that should not be used for the final estimate to the purely exponential long-distance behavior is clearly visible. The parameter $`\mathrm{\Delta }`$ in Eq. (17) can be used to tune the signal-noise ratio for the correlation length estimate; increasing $`\mathrm{\Delta }`$ dramatically reduces the apparent statistical fluctuations in $`\widehat{\xi }(i)`$, cp. Fig. 2. Note, however, that the reduction of variances for individual distances $`i`$ is accompanied by an increase of cross-correlations between estimates for adjacent estimates, so that the error of an average over a region of distances becomes minimal for a value $`\mathrm{\Delta }`$ clearly below its allowed maximum. As a compromise, we use $`\mathrm{\Delta }2\xi _ϵ`$ for both estimators $`\widehat{\xi }_\sigma (i)`$ and $`\widehat{\xi }_ϵ(i)`$.

Naive estimates for the statistical errors (variances) of complex, non-linear combinations of observable measurements like the estimator Eq. (17) are extremely biased due to two effects: even for quite sparse measurements with frequencies around $`1/\tau _{\mathrm{int}}`$ successive elements of the time series are still correlated, generically leading to systematic underestimation of variances. This effect is being eliminated by the grouping together of measurements to sub-averages of length $`\mu `$ (“binning”) flyvbjerg:89a , which leads to an asymptotically uncorrelated time series of length $`T^{}=T/\mu `$ used in the further process of error estimation. For the production-run time-series the bin size was chosen to regularly include several thousand measurements, which is far in the asymptotical regime. Secondly, the strong non-linearity of estimators like Eq. (17) forbids the use of the usual formula for the standard deviation of a set of measurements. A common solution to this problem is the use of the Gaussian error propagation formula, which, however, only uses a lowest order Taylor series approximation to the functions and assumes Gaussian distribution of the mean values, i.e. long enough time series for all observables. A far more general ansatz is given by resampling techniques such as the “jackknife” efron:82  that apply to a quite general set of probability distributions and capture function non-linearities exactly. The jackknife variance and bias estimators mimic the brute force error estimation method of comparing $`k`$ completely independent MC time series of lengths $`T^{}`$ and applying the naive estimates: removing single elements (i.e. bins) of a single time series of length $`T^{}`$ one by one results in $`T^{}`$ time series of length $`T^{}1`$, e.g. for the correlation function estimates:

$$\widehat{G}_{(s)}(i)=\frac{1}{T^{}1}\underset{ts}{}\widehat{G}_t(i),$$
(18)
resulting in jackknife-block estimates for the correlation length of:

$$\begin{array}{ccc}\hfill \widehat{\xi }_{(s)}(i)& =& \mathrm{\Delta }\left[\mathrm{ln}\frac{\widehat{G}_{(s)}(i)\widehat{G}_{(s)}(i\mathrm{\Delta })}{\widehat{G}_{(s)}(i+\mathrm{\Delta })\widehat{G}_s(i)}\right]^1,\hfill \\ \hfill \widehat{\xi }_{()}(i)& =& \frac{1}{T^{}}\underset{s}{}\widehat{\xi }_{(s)}(i).\hfill \end{array}$$
(19)
Then the jackknife estimate of variance is given by:

$$\widehat{\mathrm{VAR}}(\widehat{\xi }(i))=\frac{T^{}1}{T^{}}\underset{s=1}{\overset{T^{}}{}}\left(\widehat{\xi }_{(s)}(i)\widehat{\xi }_{()}(i)\right)^2.$$
(20)
Note the missing factor of $`1/(T^{}1)^2`$ as compared to the usual variance estimate which accounts for the trivial correlation between the $`T^{}`$ jackknife-block estimates. One can show that this estimator, apart from the reweighting prefactor $`(T^{}1)/T^{}`$, is strictly conservative, i.e. deviations from the true variance are always positive efron:82 . Similarly, the resampling scheme provides an estimate for the bias of estimators, namely:

$$\widehat{\mathrm{BIAS}}(\widehat{\xi }(i))=(T^{}1)(\widehat{\xi }_{()}(i)\widehat{\xi }(i)),$$
(21)
and thus offers a bias corrected correlation length estimate as $`\stackrel{~}{\xi }(i)=T^{}\widehat{\xi }(i)(T^{}1)\widehat{\xi }_{()}(i)`$. Since in non-pathological cases the bias of an estimator vanishes with increasing length of the time series, the jackknife bias estimate provides a good check for whether the considered series are long enough to neglect bias. A jackknife error estimate for these bias-corrected estimators is possible iterating the jackknife resampling scheme to second order berg:92a .

Since Eq. (17) gives a vector of estimators for the correlation length instead of only a single one, an improved final estimate can be achieved by an average over the $`\widehat{\xi }(i)`$. However, as for example Fig. 2 reveals, only a certain range of distances $`i=i_{\mathrm{min}},\mathrm{},i_{\mathrm{max}}`$ is suited for this purpose, where the lower bound $`i_{\mathrm{min}}`$ results mainly from small-distance deviations as reflected by Eq. (4), whereas the large distance bound $`i_{\mathrm{max}}`$ cuts off the region where the signal of exponential fall-off drops below the size of statistical fluctuations, so that error estimates become inaccurate and eventually the estimator Eq. (17) becomes maldefined due to negative arguments of the logarithm. Conventionally, averaging over the estimates $`\widehat{\xi }(i)`$ for $`i=i_{\mathrm{min}},\mathrm{},i_{\mathrm{max}}`$ would be done with weights $`\alpha _i1/\sigma ^2(\widehat{\xi }(i))`$ that minimize the theoretical variance of the mean value. This prescription, however, neglects correlations between the individual estimates. Note that cross-correlations between adjacent estimates $`\widehat{\xi }(i)`$ are quite large, not only because large scale fluctuations of the correlation functions are dominant, but also since the used estimator Eq. (17) explicitly introduces such correlations increasing in range with increasing $`\mathrm{\Delta }`$. As a simple variational calculation shows, for the case of correlated variables to be averaged over, one has to choose the weights according as

$$\alpha _k=\frac{_i(\mathrm{\Gamma }^1)_{ik}}{_{i,j}(\mathrm{\Gamma }^1)_{ij}},$$
(22)
in order to minimize the variance of the mean value. Here, $`\mathrm{\Gamma }_{p\times p},p=i_{\mathrm{max}}i_{\mathrm{min}}+1`$, denotes the covariance matrix of the $`\widehat{\xi }(i)`$. $`\mathrm{\Gamma }`$ itself can be estimated within the jackknife resampling scheme as:

$`\widehat{\mathrm{CORR}}_{ij}\widehat{\mathrm{CORR}}(\widehat{\xi }(i),\widehat{\xi }(j))=`$
$`={\displaystyle \frac{T^{}1}{T^{}}}{\displaystyle \underset{s=1}{\overset{T^{}}{}}}\left(\widehat{\xi }_{(s)}(i)\widehat{\xi }_{()}(i)\right)\left(\widehat{\xi }_{(s)}(j)\widehat{\xi }_{()}(j)\right).`$ (23)
The fact that, considering Eq. (22), variance and covariance estimates directly influence the final results for the correlation lengths, gave the motivation for the quite careful statistical treatment presented above.

Finally, the selection of the regime $`i=i_{\mathrm{min}},\mathrm{},i_{\mathrm{max}}`$ can, besides the obvious eyeball method, also be done in a way based on statistical criteria. Interpreting the average over the $`\widehat{\xi }(i)`$ as a fit of the estimated $`\widehat{\xi }(i)`$ values to the trivial function $`f(\widehat{\xi })=\overline{\xi }=\mathrm{const}`$, the systematic deviations from the plateau regime for very small and very large distances $`i`$ should be clearly reflected in quality-of-fit parameters. Thus, looking at the $`\chi ^2`$-distribution,

$$\widehat{\chi }^2=\underset{i,j=i_{\mathrm{min}}}{\overset{i_{\mathrm{max}}}{}}(\widehat{\xi }(i)\overline{\xi })(\widehat{\mathrm{\Gamma }}^1)_{ij}(\widehat{\xi }(j)\overline{\xi }),$$
(24)
will be a good criterion for judging the “flatness” of the plateau regime $`i_{\mathrm{min}},\mathrm{},i_{\mathrm{max}}`$ included in the average. Again, as an estimator $`\widehat{\mathrm{\Gamma }}_{ij}`$ for the covariance matrix one can use the jackknife expression $`\widehat{\mathrm{CORR}}_{ij}`$. Then finding the optimal region of distances for the average is equivalent to the optimization problem $`|\widehat{\chi }^2/g1|\mathrm{min}`$, with $`g=i_{\mathrm{max}}i_{\mathrm{min}}=p1`$ denoting the number of degrees of freedom of the fit. However, this ansatz of optimization bears some uncertainties: minimizing the distance of $`\widehat{\chi }^2/g`$ from $`1`$ supposes that the optimal choice includes estimates $`\widehat{\xi }(i)`$ whose dispersion around $`\overline{\xi }`$ is exactly reflected by the estimated variances. In view of the jackknife’s tendency to overestimate errors it might be more favorable to minimize $`|\widehat{\chi }^2/g|`$ itself. Furthermore, considering the statistical nature of the data, the absolute minimum of $`|\widehat{\chi }^2/g1|`$ or $`|\widehat{\chi }^2/g|`$ sometimes happens to be an isolated fluctuation, far apart from the bulk of next-to-optimal solutions. Finally, this optimization procedure tends to result in minimal values for very small regime sizes $`p`$ since the fit becomes trivial for very small numbers of points; this, however, conflicts with another possible goal of optimization, namely the minimization of the overall variance of the final result. To circumvent these problems we resort to considering the whole two-dimensional distribution $`\widehat{\chi }^2/g(i_{\mathrm{min}},i_{\mathrm{max}})`$. It is characterized by a rather flat plateau regime for intermediate values of $`i_{\mathrm{min}}`$ and $`i_{\mathrm{max}}`$ and steep increases at the boundaries, cp. Fig. 3. A good recipe for the determination of bounds is then given by first choosing a preliminary $`i_{\mathrm{min}}`$ well above the steep ascent for small $`i`$; then a plot like Fig. 3(a) allows to determine the upper bound $`i_{\mathrm{max}}`$. Finally, a plot of $`\{\widehat{\chi }^2/g|i_{\mathrm{max}}=\mathrm{const}\}`$ determines the final lower bound $`i_{\mathrm{min}}`$, cp. Fig. 3(b).

To test the methods of data analysis described in this section we performed simulations of the two-dimensional Ising model. Using a series of systems with $`L_x=5,\mathrm{},20`$ and finite-size scaling fits including an effective higher-order correction term of the form $`\xi (L_x)=AL_x+BL_x^\kappa `$, we find for the leading correlation lengths scaling amplitudes $`A_{\sigma /ϵ}`$ final estimates for the case of periodic boundary conditions of $`A_\sigma =1.27374(81)`$ and $`A_ϵ=0.1583(17)`$, in excellent agreement with the exact results $`A_\sigma =4/\pi 1.27324`$ and $`A_ϵ=1/2\pi 0.15915`$, cp. Eq. (5). For the case of antiperiodic boundary conditions we arrive at $`A_\sigma =0.42410(30)`$ and $`A_ϵ=0.07984(38)`$, compared to CFT results of $`A_\sigma =4/3\pi 0.42441`$ and $`A_ϵ=1/4\pi 0.07958`$, cp. Eq. (8).

## IV Results: amplitude ratios

Let us now turn to the three-dimensional geometry $`T^2\times `$ and the determination of amplitude ratios according to Eq. (6). We report the results of simulations for the O($`n`$) spin models for $`n=1`$, $`2`$, $`3`$, and $`10`$.

### Ising Model –

Simulations of the Ising model were done at an inverse temperature given by a high-precision MC estimate of the bulk critical coupling in three dimensionstalapov:96 , $`\beta _c=0.2216544(3)`$. We use a temperature reweighting technique to check for the influence of the uncertainty of $`\beta _c`$ on the final results ferrenberg:88a ; ferrenberg:88ae . We find it completely negligible compared to the statistical errors for the case of the Ising model. To enable a proper FSS analysis including sub-leading terms we performed simulations for transverse system sizes $`L_x=4,\mathrm{\hspace{0.17em}5},\mathrm{},\mathrm{\hspace{0.17em}20},\mathrm{\hspace{0.17em}25},`$ and $`30`$, scaling $`L_z`$ accordingly. Adapting the frequency of measurements to the autocorrelation times, about $`2\times 10^6`$ and $`8\times 10^6`$ nearly independent measurements were recorded for the systems with periodic and with antiperiodic boundary conditions, respectively. Collecting the final estimates $`\overline{\xi }`$ for the correlation lengths one ends up with a scaling plot like that shown in Fig. 4. The scaling behavior is quite linear, however, as plots of the amplitudes $`\overline{\xi }/L_x`$ reveal, corrections to the purely linear scaling behavior are clearly resolvable, cp. Fig. 5. As an aside, Fig. 5(b) additionally shows jackknife bias corrected estimators according to Eq. (21); for the given length of time series bias effects of our estimator Eq. (17) can clearly be neglected.

Returning to the two-dimensional case for a moment, it is easy to see the source for the leading correction term in the correlation length scaling. In the framework of conformal field theory the effect of lattice discretization as well as the influence of non-linearity of scaling fields that increase with the distance from criticality (i.e. the thermodynamic limit in our case) can be included in considerations using conformal perturbation theory henkel:book . A formal perturbation expression for the spin-spin correlation length including the effect of a perturbing operator coupled with strength $`a_k`$ is to first order given by

$$\xi _{\sigma }^{}{}_{}{}^{1}=\frac{2\pi }{L}\left[x+2\pi a_k(𝐂_{1k1}𝐂_{0k0})\left(\frac{2\pi }{L}\right)^{x_k2}\right],$$
(25)
where the perturbing operator has dimension $`x_k`$ and the coefficients $`𝐂_{nkn}`$ result from the operator product expansion (OPE). One finds reinicke:87a  that to lowest order the only non-vanishing amplitude belongs to an operator that corresponds to the breaking of rotational symmetry by the square lattice as compared to the continuum solution. It has dimension $`x_k=4`$ leading to $`1/L^2`$ corrections, in agreement with the first-order expansion of the exact result nightingale:76a :

$$\xi _{\sigma }^{}{}_{}{}^{1}(L)=\frac{2\pi }{L}\left[\frac{1}{8}2\pi \frac{1}{768\pi }\left(\frac{2\pi }{L}\right)^2\right].$$
(26)
A similar effect will be present in the three-dimensional systems, but the correction exponent can no longer be evaluated analytically. Fig. 5 shows that the sign of the leading correction term is unchanged in three dimensions for the systems with periodic boundary conditions, whereas it is reversed for the systems with antiperiodic boundary. This stays true for the other O($`n`$) spin models discussed below. To account for corrections to scaling we fit the correlation lengths data to the functional form

$$\xi (L_x)=AL_x+BL_x^\kappa ,$$
(27)
treating the correction exponent $`\kappa `$ as an additional fit parameter. Due to the presence of higher-order corrections, however, the resulting values of $`\kappa `$ have to be taken as effective exponents, that will in general differ from Wegner’s correction exponent $`\omega `$. Therefore we decided to keep $`\kappa `$ as a parameter, despite of existing field-theory estimates for $`\omega `$, cp. zinn-justin . We take into account the effect of neglecting higher-order correction terms by successively dropping points from the small $`L_x`$ end while monitoring the quality-of-fit parameters $`\chi ^2/g`$ or $`Q`$ to find a compromise between fit stability and precision of the final amplitudes $`A`$. The range of sizes $`L_x`$ used is indicated by the range of the solid lines in Fig. 5. Our results for the scaling amplitudes and their ratios as listed in Table 3 and the ratio of scaling dimensions according to Eq. (10) show precise agreement in the sense of Eq. (6) for the case of antiperiodic boundary conditions and clear mismatch for a periodic boundary. This is in agreement with the results of Henkel henkel:87a  and Weston weston:90 , but at a level of accuracy that makes a casual coincidence very unlikely.

### XY Model –

The XY model is, as well as the Heisenberg models, accessible to cluster update methods using the embedded cluster representation wj:chem , which we made use of. The simulations were performed at the coupling $`\beta _c=0.454167(3)`$, which is an average of recent literature estimates, cp. Table 2. Using the same transverse system sizes $`L_x=L_y`$ as for the Ising model, but adjusting the lengths $`L_z`$ according to the different correlation length amplitudes, we took between $`1\times 10^6`$ and $`16\times 10^6`$ measurements, using measurement frequencies around $`1/\tau _{\mathrm{int}}`$ as above. Fig. 6 shows the amplitude scaling plot of the spin-spin correlation length for periodic boundary conditions. The additional curves are results of a temperature reweighting analysis, trying to judge the effect of critical coupling uncertainties. The precision of the data is well illustrated by the fact that, reweighting our results to the minimum and maximum estimated critical couplings, respectively, cited in Table 2, results in a variation of the scaling curves far beyond the range covered by the remaining statistical errors. Nevertheless, reweighting to the $`1\sigma `$-range inverse temperatures $`\beta _c\mathrm{\Delta }\beta `$ and $`\beta _c+\mathrm{\Delta }\beta `$ as given above triggers deviations at most comparable to the error estimates of the statistical analysis. The intermediate maximum of the curve for $`\beta _{\mathrm{min}}`$, however, might be an artefact indicating that $`\beta _{\mathrm{min}}`$ is already too far away from the simulation temperature to allow for reliable reweighting. The effect of temperature variation is generally observed to be smaller for the antiperiodic boundary systems; furthermore, it is more important for the case of the spin-spin correlation length since here statistical errors are clearly smaller than for the energy-energy correlation length estimates. Thus, Fig. 6 shows the largest effect observed.

Fitting the final correlation length results $`\overline{\xi }_{\sigma /ϵ}`$ to the functional form Eq. (27), we arrive at the final estimates for the leading amplitudes given in Table 3. Comparing these to the ratio of scaling dimensions resulting from the averaged critical exponent estimates of Table 2 and Eq. (10), we again find Eq. (6) confirmed for antiperiodic boundary conditions only; this behavior is obviously not specific to the Ising model.

### Heisenberg Model –

The $`n=3`$ Heisenberg model case is treated analogously to the XY model. Table 2 gives the critical parameter estimates used for the simulations and comparison. With statistics similar to that for the $`n=1`$ and $`n=2`$ cases, the simulations confirm the findings of the Ising and XY models, cp. Table 3 for details. For the case of the energy-energy correlation length of the systems with periodic boundary conditions the gathered statistics did not suffice for a stable non-linear fit including corrections according to Eq. (27). We thus performed a simple linear fit dropping the correction term. This, however, results in an error estimate which is not quite realistic and, furthermore, induces a systematic underestimation of the amplitude since one expects $`B_ϵ<0`$, cp. Fig. 5(b). From the results of the other models this effect is estimated to be about $`2\sigma `$-$`3\sigma `$ in magnitude.

### O(10) Model –

To gain additional evidence and in order to facilitate considerations concerning the $`n\mathrm{}`$ limit, giving a clear picture of systematic dependencies on the parameter $`n`$, we also simulated the $`n=10`$ generalized Heisenberg model. Since, of course, in the past much less effort has gone into the investigation of the O($`n`$) model with $`n>4`$, there are quite few estimates of the critical coupling. We thus here use a single high-temperature series estimate of $`\beta _c=2.42792(8)`$ butera:97a . The implementation of the Wolff cluster update algorithm has to cope with the technical intricacy of generating pseudo-random numbers equally distributed on a hyper-sphere, see Appendix A for details. Due to this complication we only simulated systems up to a transversal size of $`L_x=20`$ and reduced the number of measurements to $`2\times 10^6`$. The critical exponents for comparison, given by a plain average over some recent estimates sokolov:95a ; butera:97a ; kleinert:99a , are:

$$\nu =0.8713(75),\gamma =1.721(14).$$
(28)
Table 3 shows again agreement between amplitude and exponent ratios only for the case of antiperiodic boundaries. Note that, as critical exponent estimates become rare with increasing $`n`$, the correlation length ratio estimate already reaches the precision of the scaling dimension ratio estimate. Checking the influence of the critical coupling uncertainty we find it only important compared to statistical errors in the case of the spin-spin correlation length for periodic boundary systems; the results reweighted to $`\beta _\pm =\beta _c\pm \mathrm{\Delta }\beta `$ are $`A_\sigma ^{}=0.670805(56)`$ and $`A_\sigma ^+=0.671432(65)`$, respectively. This, however, does not noticeably influence the error of the ratio estimate, since here the error of the estimate of $`A_ϵ`$, which is much larger, is dominant.

We thus find the linear amplitude-exponent relation Eq. (6) confirmed for several spin models in three dimensions with the peculiarity that one has to insert a seam of antiferromagnetic bonds along the $`T^2`$-directions to restore the two-dimensional situation.

## V Results: “Meta” Amplitudes

Comparing our results for the three-dimensional geometry $`S^1\times S^1\times `$ to the CFT conjecture for the case of two dimensions, we are interested in the respective ranges of validity in terms of the classification of universality aspects given above in the Introduction. The fact that our simulations of the isotropic lattice representation of the O($`n`$) universality classes give results in agreement with the strongly anisotropic quantum Hamiltonian representation used by Henkel in his transfer matrix calculations for the case of the Ising model henkel:86a ; henkel:87a ; henkel:87b , indicates that the considered amplitude ratios are universal, i.e. (i’) holds. Apart from that, Henkel henkel:86a  explicitly checked for universality of amplitude ratios by the insertion of an irrelevant perturbing operator and found it confirmed for both cases of boundary conditions. However, strictly speaking, there is no proof of universality for the cases $`n>1`$. The universality aspect (i) above, i.e. universality of the amplitudes themselves, could not be checked in Henkel’s calculations, because the quantum Hamiltonian is only defined up to an overall normalization constant. Yurishchev yuri:94a ; yuri:97a  considered the behavior of an anisotropic Ising model and found varying correlation lengths amplitudes on variation of the ratios of couplings in the different directions. This, however, is no argument against amplitude universality since anisotropy is represented by marginal instead of irrelevant operators. On the other hand, amplitude ratios stay universal even with respect to those marginal perturbations, in consistency with Henkel’s strongly anisotropic Hamiltonian limit calculations. In fact it has been argued that for all systems below their upper critical dimension correlation length scaling amplitudes are universal quantities privman:84a .

Having found very good agreement in three dimensions between ratios of correlation lengths and scaling dimensions according to Eq. (6) for the case of antiperiodic boundary conditions, it is interesting to see what the overall, operator-independent, “meta” amplitude $`𝒜`$ according to:

$$\xi _{\sigma /ϵ}=A_{\sigma /ϵ}L_x=\frac{𝒜}{x_{\sigma /ϵ}}L_x,$$
(29)
that was $`𝒜=1/2\pi `$ for two-dimensional periodic systems, cp. Eq. (5), becomes in three dimensions, in particular whether it is again model independent. Since our results for the spin-spin correlation lengths are always more precise than those for energy-energy correlation lengths, we use $`\overline{\xi }_\sigma `$ to determine $`𝒜`$. The estimates for the spin-spin scaling dimension $`x_\sigma `$ resulting from the corresponding estimates of bulk critical exponents $`\nu `$ and $`\gamma `$ listed in Tables 1 and 2 and Eq. (28) are $`x_\sigma =0.5182(4)`$ (Ising), $`x_\sigma =0.5188(9)`$ (XY), $`x_\sigma =0.5160(17)`$ (Heisenberg), and $`x_\sigma =0.512(12)`$ ($`n=10`$), respectively. Thus, inserting our results for $`A_\sigma `$ listed in Table 3, we obtain for the “meta” amplitudes $`𝒜(n)`$:

$$𝒜=A_\sigma x_\sigma =\{\begin{array}{cc}0.12278(43)\hfill & \text{Ising}\hfill \\ 0.12510(37)\hfill & \text{XY}\hfill \\ 0.12622(49)\hfill & \text{Heisenberg}\hfill \\ 0.1325(30)\hfill & n=10\hfill \end{array}.$$
(30)
These values can additionally be compared with an analytical result that is available for the case of the spherical model, which is commonly believed to be identical to the $`n\mathrm{}`$ limit of the O($`n`$) spin model stanley:68a . Again using the Hamiltonian formulation, Henkel and Weston henkel:88a ; henkel:92a  found that the amplitude exponent relation Eq. (6) is valid for the spherical model on $`S^1\times S^1\times `$ for both kinds of boundary conditions, periodic and antiperiodic. This is due to the fact that the quantum Hamiltonian factorizes into a set of uncoupled harmonic oscillators. The amplitude $`𝒜`$ for the case of antiperiodic boundary conditions was found to be $`𝒜0.13624`$ henkel:92a ; allen:93a . Plotting this value together with the finite-$`n`$ results of Eq. (30) shows an apparently smooth variation of the “meta” amplitudes with the order parameter dimension $`n`$, the eyeball extrapolation of the finite-$`n`$ values to $`1/n0`$ matching the spherical model result, cp. Fig. 7(a). Facing this variation, the hypothesis of a “hyperuniversal” amplitude $`𝒜(n)=𝒜`$ that does not depend on $`n`$, as was the case for the two-dimensional systems, can be clearly ruled out. Thus, type (iii) universality of the classification above gets broken when passing from two to three dimensions. The matching of the finite-$`n`$ values with the universal spherical model amplitude, on the other hand, indicates universality also of the finite-$`n`$ amplitudes and thus universality of type (i) above.

Even without a scaling law of the type Eq. (6) being valid for the case of periodic boundary conditions, one can nevertheless plot the corresponding combination $`A_\sigma x_\sigma `$ for this case also, as is illustrated in Fig. 7(b). The values are:

$$A_\sigma x_\sigma =\{\begin{array}{cc}0.4240(17)\hfill & \text{Ising}\hfill \\ 0.3912(7)\hfill & \text{XY}\hfill \\ 0.3719(12)\hfill & \text{Heisenberg}\hfill \\ 0.3439(78)\hfill & n=10\hfill \end{array}.$$
(31)
The corresponding value for the spherical model is given by $`A_\sigma x_\sigma 0.3307`$, cp. brezin:82a ; henkel:92a . The finite-$`n`$ values again run smoothly into the spherical model limit.

## VI The limit of infinite spin dimensionality

While the finite-$`n`$ amplitudes of Fig. 7 fit well to the spherical model result, this is not the case for the correlation lengths ratios themselves. From inspection of Fig. 8 the smooth variation of correlation length ratios for finite $`n`$ does not fit at all to the spherical model result of Henkel and Weston henkel:88a ; henkel:92a  that gives a ratio $`A_\sigma /A_ϵ=2`$ for both, periodic and antiperiodic boundary conditions. By eyeball extrapolation one would instead expect the amplitude ratios to reach values around $`4`$ for antiperiodic and around $`5\frac{1}{3}`$ for periodic boundary conditions in the limit $`n\mathrm{}`$. And indeed, accepting the validity of a linear amplitude-exponent relation according to Eq. (6) for the case of antiperiodic boundary conditions and using the usual relations for the connection between scaling dimensions and bulk critical exponents, namely Eq. (7), one would expect $`x_\sigma =1/2`$ and $`x_ϵ=2`$ since $`\beta =1/2`$, $`\nu =1`$ and $`\alpha =1`$ for the spherical model. The resulting ratio $`x_ϵ/x_\sigma =4`$ perfectly agrees with the eyeball extrapolation of our finite-$`n`$ data. However, by inspection of the energy-energy correlation function in the Hamiltonian limit and using factorization arguments, Henkel henkel:88a  conjectured $`x_ϵ=1`$ instead, resulting in the ratio $`A_\sigma /A_ϵ=2`$, in contrast to the relation Eq. (7). Taking into account the obvious agreement of eyeball extrapolation and spherical model calculation for the amplitudes $`𝒜(n)`$ that were calculated from the spin-spin correlation length amplitude as $`𝒜(n)=A_\sigma x_\sigma `$, cp. Fig. 7, it becomes obvious that the mismatch is entirely due to the behavior of the energy-energy correlations. Note also that, since the specific heat is constant in the low-temperature phase of the spherical model in three dimensions, interpreting this as an effectively vanishing specific-heat exponent $`\alpha ^{}=0`$ leads to an effective energetic scaling dimension $`x_ϵ^{}=1`$. This, in fact, implies a violation of the scaling relation Eq. (7), which is of the hyper-scaling type, for the case of the spherical model.

Puzzled by this striking mismatch, we performed a roughly explorative MC simulation directly in the spherical model, which rendered results in qualitative agreement with an amplitude ratio of $`A_\sigma /A_ϵ=2`$ as suggested by the analytical calculation. Then, it is natural to ask whether there is a contradiction with Stanley’s result on the equivalence of the $`n\mathrm{}`$ limit of the O($`n`$) model and the spherical model berlin:52 , which has been, after some debate over mathematical subtleties helfand:69a , rigorously proven kac:71a . The precise statement that can be proven is the identity of the partition functions or, equivalently, free energies of the two models in the thermodynamic limit for the whole temperature range, even independent of the order of taking the limits $`n\mathrm{}`$ and $`N\mathrm{}`$ (the thermodynamic limit). Since multi-point correlation functions do not follow from the (source-free) partition function, this does not say anything about the behavior of these functions in those two models. A direct calculation in the spherical model, cf. Appendix B, results in a simple factorization property of the long distance behavior of the connected energy-energy correlation function for all temperatures in one and two dimensions and in the high-temperature phase down to $`T_c`$ in three dimensions. If the four-point function of the spherical-model spins is denoted by $`C_{ijkl}`$, one has:

$`C_{ii+1jj+1}C_{ii+1}^2`$ $`=C_{ij}C_{i+1j+1}+C_{ij+1}C_{i+1j}`$ (32)
$`2C_{ij}^2,|ji|\mathrm{},`$
where $`C_{ij}`$ are the corresponding two-point functions. This confirms Henkel’s results for the Hamiltonian formulation henkel:88a  on more general grounds.

Considering the $`n\mathrm{}`$ limit of the O($`n`$) model, on the other hand, reveals that the connected part of the energy-energy correlation function vanishes in the first-order saddle-point approximation that is being used for the comparison of the two models, cf. Appendix C. This is in agreement with general considerations for the large $`n`$ model by Brézin brezin:82a . For the case of the one-dimensional spin chain, the connected energy-energy correlation function even vanishes exactly for all $`n`$, so that one can rule out an agreement of the two limits to higher order of the steepest-descent expansion in this case. Thus the mismatch of finite-$`n`$ extrapolations and spherical model results of Fig. 8 has some well-defined mathematical reason.

Starting from the observation that the curves of Fig. 8 for the amplitude ratios seem to be quite parallel as a function of (finite) $`n`$ for the both kinds of boundary conditions, we also plotted the collapsed ratio $`\frac{A_\sigma /A_ϵ}{x_ϵ/x_\sigma }`$ that should be unity if the amplitude-exponent relation Eq. (6) holds true. Inspecting Fig. 9, this is, according to our above results, of course the case for antiperiodic boundary conditions. Moreover, and a priori somewhat unexpected, this ratio seems to be also quite constant for the case of a periodic boundary, stabilizing around a value compatible with $`4/3`$ within statistical errors. Note that the exceptionally small error of the value for $`n=3`$ (the Heisenberg model) and its apparent deviation towards a larger ratio is due to the impossibility to fit the $`n=3`$ energy-energy correlation lengths to a scaling law including a correction term as mentioned in Sec. IV. Statistically, the data are consistent with a fit to a constant ($`Q=0.08`$), and perfectly so when dropping the $`n=3`$ point ($`Q=0.4`$).

In view of this observation one might argue that the asymptotic scaling relation Eq. (6) in three dimensions has to be replaced by a generalized ansatz of the form

$$\frac{\xi _\sigma }{\xi _ϵ}=R\frac{x_ϵ}{x_\sigma },$$
(33)
with an overall, model independent factor $`R`$ that depends only on the boundary conditions and happens to be just $`1`$ for the case of an antiperiodic boundary. For the amplitude scaling law this would lead to an asymptotic form

$$\xi _{\sigma /ϵ}(n)=R\frac{𝒜(n)}{x_{\sigma /ϵ}}L_x,$$
(34)
cp. Eq. (5). Accepting such a generalized ansatz, a least-squares fit of the collapsed ratios of Fig. 9 to a constant $`R`$ gives $`R=1.0037(45)`$ for antiperiodic boundary conditions, underlining the validity of the original amplitude-exponent relation Eq. (6), or alternatively Eq. (33) with $`R=1`$, for this case. For the periodic-boundary systems, on the other hand, we arrive at $`R=1.3546(76)`$ (omitting the $`n=3`$ point), indeed statistically consistent with the conjectured value of $`4/3`$.

This somewhat diminishes the at first sight apparently exceptional importance of choosing antiperiodic boundary conditions in three dimensions. Taking into account the smooth amplitude variation of Fig. 7(b) the same universality statements hold for periodic and for antiperiodic boundary conditions.

## VII Conclusions

We performed extensive MC simulations for several representatives of the class of O($`n`$) spin models. Concentrating on the geometry of three-dimensional slabs $`S^1\times S^1\times `$ we found a simple inversely linear relation between the leading scaling amplitudes of the correlation lengths of the magnetization and energy densities and the corresponding scaling dimensions valid to high accuracy for the Ising ($`n=1`$), XY ($`n=2`$), Heisenberg ($`n=3`$), and $`n=10`$ generalized Heisenberg models for antiperiodic boundary conditions along the torus directions. This is the analogue of the CFT result in two dimensions with periodic boundary conditions applied. There is evidence for the universality not only of amplitude ratios (type (i’) of our classification in the Introduction), but also of scaling amplitudes themselves (type (i)). To definitely decide the question whether universality in the sense (ii) above, i.e. condensation of all operator dependent information in the scaling dimensions, is present, further operators would have to be considered. Independence, apart from changes in the scaling dimension, of the scaling amplitudes from the model under consideration, i.e. type (iii) universality, is explicitly broken for three dimensions as compared to the two-dimensional case: we find a smooth variation of the overall “meta” amplitudes $`𝒜(n)=A_\sigma (n)x_\sigma (n)`$, depending on the order-parameter dimension $`n`$. It might be interesting to consider further classes of models, such as for example Potts models, to see whether any of these properties are specific to the O($`n`$) spin model class.

Considering the deviation of the periodic boundary correlation lengths ratios from the corresponding inverse scaling dimension ratios, the validity of a scaling law of the form Eq. (6) can be definitely ruled out for this case. Generalizing this ansatz with an overall factor $`R`$ depending on boundary conditions as in Eq. (33), however, we find it fulfilled also for the case of periodic boundaries with a factor $`R`$ independent from $`n`$ and taking a value compatible with $`4/3`$. In view of that, the fact that $`R=1`$ for the case of antiperiodic boundary conditions might be rather a coincidence than a “deep” physical property. Taking into account that in two dimensions the corresponding prefactors are specific to the operators considered, cp. Eq. (8), makes it probable that a similar behavior occurs in three dimensions, destroying type (ii) universality. It might be interesting to analyze the behavior of correlation lengths in the four-dimensional geometry $`S^1\times S^1\times S^1\times `$ to check whether a scaling law of the generalized form Eq. (33) can be retained and if so, how the factor $`R`$ depends on the dimensionality of the lattice.

Trying to match our finite $`n`$ results with analytical calculations for the spherical model we found a striking mismatch of the data concerning energy-energy correlations. Inspecting the four-point functions directly in the spherical model and the O($`n\mathrm{}`$) model limit we find that both results do not match to first order of the saddle-point approximation in general dimensions and to all orders in one dimension. Thus, the idea of equivalence of the two models has to be limited to its original extent, namely the identity of partition functions in the thermodynamic limit. Quantities not directly related to the partition function, like multi-point correlation functions, do not necessarily have to coincide. Further work has to be done to possibly evaluate exactly the correlation lengths ratios in the $`n\mathrm{}`$ limit for both sorts of boundary conditions.

Since, still, there is no explanation of the findings concerning the correlation lengths ratios for finite $`n`$ in terms of a field-theoretic or otherwise exact approach, we would like to encourage further research in this direction.

###### Acknowledgements.

The authors thank K. Binder, J. L. Cardy, and M. Henkel for useful discussions. MW gratefully acknowledges support by the “Deutsche Forschungsgemeinschaft” through the “Graduiertenkolleg Quantenfeldtheorie”.

## Appendix A Equal distribution of random numbers on a hyper-sphere

Consider a probability density in polar co-ordinates $`f(\varphi ,\theta )`$ equally distributed on the $`2`$-sphere $`S^2`$, i.e.:

$$\frac{f(\varphi ,\theta )\mathrm{d}\varphi \mathrm{d}\theta }{\mathrm{sin}\theta \mathrm{d}\varphi \mathrm{d}\theta }=\mathrm{const}.$$
(35)
Factorizing $`f(\varphi ,\theta )=p(\varphi )q(\theta )=\mathrm{const}\mathrm{sin}\theta `$, and taking into account the normalization condition $`d\mathrm{\Omega }f(\varphi ,\theta )=1`$, one finds:

$$f(\varphi ,\theta )=p(\varphi )q(\theta )=\frac{1}{2\pi }\frac{1}{2}\mathrm{sin}\theta .$$
(36)
Pseudo-random number generators usually generate numbers equally distributed in the unit interval $`[0,1]`$. How does this transform to an arbitrary distribution? Let a random variable $`z`$ be distributed with a density $`g(z)`$ and transform according to $`z^{}=\omega (z)`$; the density $`h(z^{})`$ then follows from the equation

$$g(z)\mathrm{d}z=h(z^{})\mathrm{d}z^{}=h(\omega (z))\omega ^{}(z)\mathrm{d}z.$$
(37)
Thus, for random numbers $`z`$ equally distributed in $`[0,1]`$ the transformation $`\theta =\mathrm{arccos}(12z)`$ gives the desired distribution $`q(\theta )=\frac{1}{2}\mathrm{sin}\theta `$. This form is being used for the simulations of the $`n=3`$ Heisenberg model. For general polar co-ordinates in $`^n`$, $`x_1=r\mathrm{cos}\theta _1`$, $`x_2=r\mathrm{sin}\theta _1\mathrm{cos}\theta _2`$, up to $`x_n=r\mathrm{sin}\theta _1\mathrm{}\mathrm{sin}\theta _{n1}`$, where $`0\theta _i\pi `$, $`0\theta _{n1}<2\pi `$ is understood, the volume element is given by:

$$\mathrm{d}V=r^{n1}\mathrm{sin}^{n2}\theta _1\mathrm{sin}^{n3}\theta _2\mathrm{}\mathrm{sin}\theta _{n2}\mathrm{d}r\underset{i}{}\mathrm{d}\theta _i,$$
(38)
so that one has for the factors $`f^{(i)}(\theta _i)`$ of an equally distributed density $`f(\theta _1,\mathrm{},\theta _{n1})=_if^{(i)}(\theta _i)`$:

$`f^{(i)}(\theta _i)={\displaystyle \frac{1}{\gamma (ni1)}}\mathrm{sin}^{ni1}\theta _i,i<n1,`$
$`f^{(n1)}(\theta _{n1})={\displaystyle \frac{1}{2\pi }},`$ (39)
with normalization factors $`\gamma (k)=\sqrt{\pi }\mathrm{\Gamma }(\frac{k+1}{2})/\mathrm{\Gamma }(\frac{k}{2}+1)`$. Thus, for $`z_i`$ equally distributed in $`[0,1]`$ the transformations $`z_i(\theta _i)`$ are given by:

$$z_i(\theta _i)\mathrm{int}(\theta _i)=\frac{1}{\gamma (ni1)}d\theta _i\mathrm{sin}^{ni1}\theta _i,$$
(40)
for $`i<n1`$. The integrals can be evaluated analytically for each $`\theta _i`$. There is, however, no closed form expression for the inverse transformation $`\theta _i(z_i)`$ that is needed to generate random vectors equally distributed on the hyper-sphere $`S^{n1}`$. The trivial workaround solution of sampling equally distributed in the hyper-cube $`L^n=[1,1]\times \mathrm{}\times [1,1]`$, discarding the complement $`L^n\backslash B^n`$ and projecting the remaining points on the sphere $`S^{n1}`$, suffers from asymptotically vanishing efficiency, since the ratio of used to discarded volumes vanishes with increasing $`n`$ exponentially as $`\pi ^{n/2}/2^n\mathrm{\Gamma }(\frac{n}{2}+1)`$. We thus resorted to a numerical inversion of $`z_i(\theta _i)`$ using interpolation between the pre-calculated points of a binary tree.

## Appendix B Energy-energy correlation function in the spherical model

Consider the spherical model of Berlin and Kac berlin:52  consisting of “spins” $`ϵ_i`$ with the constraint:

$$\underset{i=1}{\overset{N}{}}ϵ_i^2=N,$$
(41)
where $`N`$ denotes the number of lattice sites. For ease of reference we use the notation of the original paper here; thus, the $`ϵ_i`$ are not to be confused with the local energy densities defined above in Eq. (13). The Hamiltonian is:

$$=J\underset{ij}{}ϵ_iϵ_j.$$
(42)
Using the Fourier representation of the $`\delta `$-constraint Eq. (41) the partition function can be written as:

$`Z_N`$ $`=`$ $`{\displaystyle \frac{A_N^1}{2\pi i}}{\displaystyle _{\alpha _0i\mathrm{}}^{\alpha _0+i\mathrm{}}}dse^{Ns}{\displaystyle \mathrm{}dϵ_1\mathrm{}dϵ_N}`$ (43)
$`\times \mathrm{exp}(s{\displaystyle \underset{i}{}}ϵ_i^2+K{\displaystyle \underset{ij}{}}ϵ_iϵ_j),`$
choosing $`\alpha _0`$ such that the singularities in $`s`$ of the integrand are excluded from the integration volume. $`A_N`$ ensures the correct normalization of the integral measure and $`K=\beta J`$ denotes the coupling. Diagonalizing the quadratic form $`_{ij}ϵ_iϵ_j`$ with eigenvalues $`\lambda _j`$ via an orthogonal transformation $`ϵ_i=_jV_{ij}y_j`$, the Gaussian integration over the $`ϵ_i`$ can be performed:

$`{\displaystyle \mathrm{}dy_1\mathrm{}dy_N\mathrm{exp}[\underset{j}{}(sK\lambda _j)y_j^2]}`$
$`=\pi ^{N/2}\mathrm{exp}[{\displaystyle \frac{1}{2}}{\displaystyle \underset{j}{}}\mathrm{ln}(sK\lambda _j)]`$ (44)
so that,

$`Z_N`$ $`=`$ $`A_N^1\pi ^{N/2}2Ke^{\frac{1}{2}N\mathrm{ln}2K}{\displaystyle \frac{1}{2\pi i}}{\displaystyle _{z_0i\mathrm{}}^{z_0+i\mathrm{}}}dz`$ (45)
$`\times \mathrm{exp}[N2Kz{\displaystyle \frac{1}{2}}{\displaystyle \underset{j=1}{\overset{N}{}}}\mathrm{ln}(z{\displaystyle \frac{1}{2}}\lambda _j)],`$
where $`s=2Kz`$. This expression can be evaluated in the saddle point limit $`N\mathrm{}`$ depending on the distribution of the eigenvalues $`\lambda _i`$ for a given lattice. Now consider the two-point function,

$$C_{ij}ϵ_iϵ_j=\underset{r,s}{}V_{ir}V_{js}y_ry_s=\underset{r}{}V_{ir}V_{jr}y_r^2,$$
(46)
where the last equality follows from the symmetry of the partition function Eq. (43). Compared to the Gaussian integration Eq. (44) the insertion of a factor $`y_r^2`$ in the integrand gives an additional factor of:

$$\frac{1}{2(sK\lambda _r)}=\frac{1}{4K(z\frac{1}{2}\lambda _r)}.$$
(47)
The corresponding integral over $`z`$ can also be evaluated in the saddle point approximation berlin:52 . Now, analogously, consider the four-point function:

$$C_{ijkl}ϵ_iϵ_jϵ_kϵ_l=\underset{r,s,t,u}{}V_{ir}V_{js}V_{kt}V_{lu}y_ry_sy_ty_u.$$
(48)
Here, again, only paired occurrences of the $`y_m`$ survive due to the inversion symmetry:

$`C_{ijkl}`$ $`={\displaystyle \underset{r}{}}V_{ir}V_{jr}V_{kr}V_{lr}y_r^4+{\displaystyle \underset{rs}{}}V_{ir}V_{jr}V_{ks}V_{ls}y_r^2y_s^2+`$ (49)
$`{\displaystyle \underset{rs}{}}`$ $`V_{ir}V_{js}V_{kr}V_{ls}y_r^2y_s^2+{\displaystyle \underset{rs}{}}V_{ir}V_{js}V_{ks}V_{lr}y_r^2y_s^2.`$
The insertion of $`y_r^4`$ under the Gaussian integral gives an additional factor of $`3/[4(sK\lambda _r)^2]=3/[16K^2(z\lambda _r/2)^2]`$, whereas $`y_r^2y_s^2`$ gives $`1/[16K^2(z\lambda _r/2)(z\lambda _s/2)]`$, so that the diagonal terms left out in Eq. (49) are reinserted:

$`C_{ijkl}=`$ $`{\displaystyle \underset{r,s}{}}`$ $`(V_{ir}V_{jr}V_{ks}V_{ls}+V_{ir}V_{js}V_{kr}V_{ls}`$ (50)
$`+V_{ir}V_{js}V_{ks}V_{lr})y_r^2y_s^2`$
Now performing the $`z`$-integration of Eq. (45) in the saddle point limit $`N\mathrm{}`$ is equivalent to just inserting the saddle point value $`z=z_s`$ for the factors given above, whenever a normal saddle point exists. As Berlin and Kac have shown, this is the case for all finite temperatures in one and two dimensions and for $`TT_c`$ in three dimensions; in the low-temperature phase, the saddle point “sticks” to its value for $`T=T_c`$. Then, the four-point function simply factorizes, so that, comparing Eq. (50) to the expression Eq. (46) for the two-point function it is clear that:

$$C_{ijkl}=C_{ij}C_{kl}+C_{ik}C_{jl}+C_{il}C_{jk},$$
(51)
and, finally, considering the connected energy-energy correlation function, one has:

$`C_{ii+1jj+1}C_{ii+1}^2`$ $`=C_{ij}C_{i+1j+1}+C_{ij+1}C_{i+1j}`$ (52)
$`2C_{ij}^2,|ji|\mathrm{},`$
so that the energy-energy correlation function is trivially related to the spin-spin correlation function. Note that Eq. (51) would follow from Wicks’s Lemma for the Gaussian model. This especially confirms the factor-two relation $`x_ϵ/x_\sigma =2`$ between the corresponding scaling dimensions derived by Henkel using transfer matrices henkel:88a . The factorization property can also be seen in the grand-canonical formulation of the spherical model, the “mean” spherical model lewis:52a , where the hard constraint Eq. (41) is being replaced by its thermodynamical average, so that one can leave out the problematic $`z`$-integration above. There has been some debate over the coincidence of the thermodynamic limit of the two models, which is now believed to be settled yan:65a .

## Appendix C Energy-energy correlation function in the limit of infinite spin dimensionality

The treatment of the partition function of the O($`n`$) model in the $`n\mathrm{}`$ limit is quite analogous to that of the spherical model, cp. stanley:68a . For the comparison of the $`n\mathrm{}`$ limit with the spherical model the constraint $`𝝈_i𝝈_i=1`$ of Eq. (11) has to be replaced by $`𝝈_i𝝈_i=n`$. We write the partition function of the model as:

$`Z_N^{(n)}(K)=A_{N}^{(n)}{}_{}{}^{1}{\displaystyle \mathrm{}d\sigma _1^{(1)}\mathrm{}d\sigma _N^{(n)}\underset{j}{}\delta (n𝝈_j^2)}`$
$`\times \mathrm{exp}[K{\displaystyle \underset{ij}{}}{\displaystyle \underset{\nu }{}}\sigma _i^{(\nu )}\sigma _j^{(\nu )}],`$ (53)
where $`A_N^{(n)}`$ ensures the correct normalization. Rewriting the $`\delta `$-constraints to the Fourier representation, one now has to introduce $`N`$ variables $`\left\{t_i\right\}`$, arriving at:

$`Z_N^{(n)}(K)=A_{N}^{(n)}{}_{}{}^{1}\left({\displaystyle \frac{K}{2\pi i}}\right)^N{\displaystyle _{\mathrm{}}^+\mathrm{}}\mathrm{}{\displaystyle _{\mathrm{}}^+\mathrm{}}d\sigma _1^{(1)}\mathrm{}d\sigma _N^{(n)}`$
$`\times {\displaystyle _i\mathrm{}^{+i\mathrm{}}}\mathrm{}{\displaystyle _i\mathrm{}^{+i\mathrm{}}}\mathrm{d}t_1\mathrm{}\mathrm{d}t_N\mathrm{exp}(Kn{\displaystyle \underset{i}{}}t_i)`$
$`\times {\displaystyle \underset{\nu =1}{\overset{n}{}}}\mathrm{exp}(K{\displaystyle \underset{i}{}}t_i\sigma _{i}^{(\nu )}{}_{}{}^{2}+K{\displaystyle \underset{ij}{}}\sigma _i^{(\nu )}\sigma _j^{(\nu )}).`$ (54)
Interchanging the order of integrations one is again left with integrals of Gaussian type that are easily solved transforming the spin variables orthogonally according to $`\sigma _i^{(\nu )}=_jV_{ij}y_j^{(\nu )}`$. Note that the transformation is symmetric in the component index $`\nu `$ of the spins. The calculation is given in more detail for the case of a one-dimensional chain below. Here, we again consider the relation between two-point and four-point correlation functions. We take the two-point function to be

$$C_{ij}\frac{1}{n}𝝈_i𝝈_j=\sigma _i^{(\nu )}\sigma _j^{(\nu )},$$
(55)
where the last equation for any $`\nu =1,\mathrm{},n`$ follows from the O($`n`$) symmetry of the model in the unbroken, high-temperature phase. Using the same arguments of Gaussian integration as for the case of the spherical model, the four-point function:

$$C_{ijkl}\frac{1}{n^2}(𝝈_i𝝈_j)(𝝈_k𝝈_l),$$
(56)
again decomposes in terms of the diagonal variables $`y_i^{(\nu )}`$ as:

$`C_{ijkl}`$ $`=`$ $`{\displaystyle \frac{1}{n^2}}{\displaystyle \underset{r,t,\mu ,\nu }{}}V_{ri}V_{rj}V_{tk}V_{tl}y_{r}^{(\mu )}{}_{}{}^{2}y_{t}^{(\nu )}{}_{}{}^{2}`$ (57)
$`+{\displaystyle \frac{1}{n^2}}{\displaystyle \underset{r,s,\mu }{}}V_{ri}V_{sj}V_{rk}V_{sl}y_{r}^{(\mu )}{}_{}{}^{2}y_{s}^{(\mu )}{}_{}{}^{2}`$
$`+{\displaystyle \frac{1}{n^2}}{\displaystyle \underset{r,s,\mu }{}}V_{ri}V_{sj}V_{sk}V_{rl}y_{r}^{(\mu )}{}_{}{}^{2}y_{s}^{(\mu )}{}_{}{}^{2}.`$
In the saddle point limit, which now corresponds to $`n\mathrm{}`$, this expression factorizes in terms of two-point functions as:

$$C_{ijkl}=C_{ij}C_{kl}+\frac{1}{n}C_{ik}C_{jl}+\frac{1}{n}C_{il}C_{jk},$$
(58)
so that the “mixed” terms are suppressed with $`1/n`$. This asymmetry results from the preset pairing of the spin component indices $`\mu `$ and $`\nu `$ in the four-point function. As a consequence, the connected part of the energy-energy correlation function:

$$C_{ii+1jj+1}C_{ii+1}^2=\frac{1}{n}C_{ij}C_{i+1j+1}+\frac{1}{n}C_{ij+1}C_{ji+1},$$
(59)
vanishes in the first-order saddle-point approximation. Thus, any non-vanishing contributions that are to be expected from our numerical results, have to come from sub-leading terms in the steepest-descent expansion. The correspondence of the $`n\mathrm{}`$ limit to the spherical model seems only to hold to leading order of the saddle-point approximation.

In the broken, low-temperature phase Eq. (55) has to be replaced by

$$C_{ij}=\frac{1}{n}𝝈_i𝝈_j\underset{\nu }{\mathrm{max}}\sigma _i^{(\nu )}\sigma _j^{(\nu )}C_{ij}^{\mathrm{max}},$$
(60)
so that the factorization property of the four-point function Eq. (58) becomes

$$C_{ijkl}C_{ij}C_{kl}+\frac{1}{n}C_{ik}^{\mathrm{max}}C_{jl}^{\mathrm{max}}+\frac{1}{n}C_{il}^{\mathrm{max}}C_{jk}^{\mathrm{max}},$$
(61)
and again the connected part of the energy-energy correlation function is $`O(1/n)`$, vanishing in the first-order saddle-point limit.

For the case of an one-dimensional lattice the first-order saddle-point approximation is exact as can be checked by explicit calculation. Consider an open chain of O($`n`$) spins<sup>1</sup><sup>1</sup>1Considering a closed chain is technically much more intricate, cp. thompson:domb , but, of course, gives the same results in the thermodynamic limit $`N\mathrm{}`$.. The partition function is given by the general expression Eq. (53) with the nearest-neighbor sum $`_{ij}𝝈_i𝝈_j`$ replaced by the one-dimensional expression $`_i𝝈_i𝝈_{i+1}`$. Following Stanley stanley:69a , we factor out the integration over the last spin $`𝝈_N`$, which has the form:

$`𝒵^{(n)}(K)={\displaystyle \frac{K}{2\pi i}}{\displaystyle \mathrm{}d\sigma ^{(1)}\mathrm{}d\sigma ^{(n)}_i\mathrm{}^{+i\mathrm{}}du}`$
$`\times \mathrm{exp}[uK(n{\displaystyle \underset{\nu }{}}\sigma _{}^{(\nu )}{}_{}{}^{2})]\mathrm{exp}[K{\displaystyle \underset{\nu }{}}c_\nu \sigma ^{(\nu )}],`$ (62)
where $`c_\nu \sigma _{N1}^{(\nu )}`$. Inserting the unity factor $`\mathrm{exp}[K\alpha _0(n_\nu \sigma _{}^{(\nu )}{}_{}{}^{2})]`$ and choosing $`\alpha _0`$ sufficiently large to exclude the singularities, one has:

$`𝒵^{(n)}(K)`$ $`=`$ $`{\displaystyle \frac{K}{2\pi i}}{\displaystyle _{\alpha _0i\mathrm{}}^{\alpha _0+i\mathrm{}}}dve^{vKn}{\displaystyle \underset{\nu }{}}{\displaystyle d\sigma ^{(\nu )}}`$ (63)
$`\times \mathrm{exp}[K(v\sigma _{}^{(\nu )}{}_{}{}^{2}c_\nu \sigma ^{(\nu )})],`$
where $`vu+\alpha _0`$. Square completion and a change of variables $`w=2v`$ gives:

$`𝒵^{(n)}(K)`$ $`=`$ $`\left({\displaystyle \frac{2\pi }{K}}\right)^{n/2}{\displaystyle \frac{K}{4\pi i}}{\displaystyle _{2\alpha _0i\mathrm{}}^{2\alpha _0+i\mathrm{}}}dw`$ (64)
$`\times \mathrm{exp}[{\displaystyle \frac{1}{2}}nK(w+1/w)]w^{n/2}`$
$`=`$ $`{\displaystyle \frac{1}{2}}K(2\pi /K)^{n/2}I_{n/21}(nK),`$
which is an integral representation of the modified Bessel function of the first kind. Thus, the spin integrations can be done successively, the full partition function being given by:

$$Z_N^{(n)}(K)=[(nK/2)^{1n/2}\mathrm{\Gamma }(n/2)I_{n/21}(nK)]^{N1},$$
(65)
where the $`\mathrm{\Gamma }`$ function enters through the normalization factor $`A_{N}^{(n)}{}_{}{}^{1}`$ and the last integration, which corresponds to $`𝒵^{(n)}(0)`$. Considering the two-point function, an additional factor $`𝝈_i𝝈_j`$, $`i<j`$, is inserted in the integrand of Eq. (53). Again starting the integration with the last spin $`𝝈_N`$, the first $`Nj`$ integrations are unaltered. The integration over $`𝝈_j`$ gives additional factors of $`c_\nu /2v`$ from the Gaussian integration Eq. (64), where now $`c_\nu \sigma _{j1}^{(\nu )}`$, so that one is left with

$$\stackrel{~}{𝒵}^{(n)}(K)=\frac{1}{2}K\left(\frac{2\pi }{K}\right)^{n/2}I_{n/2}(nK)\underset{\nu }{}\sigma _i^{(\nu )}c_\nu ,$$
(66)
and the form of the integrand for the next integrations is unchanged. The integration over $`𝝈_i`$ adds a factor of $`n`$ since $`c_\nu `$ above becomes $`\sigma _i^{(\nu )}`$ and $`_\nu \sigma _i^{(\nu )}\sigma _i^{(\nu )}=n`$, followed by another $`i1`$ integrations of the partition-function type. With $`uu(nK)=I_{n/2}(nK)`$ and $`vv(nK)=I_{n/21}(nK)`$ one arrives at:

$$\frac{1}{n}𝝈_i𝝈_j=\frac{v^{N1+ij}u^{ji}}{v^{N1}}=(u/v)^{ji}.$$
(67)
From this it is straightforward to derive the form of the four-point function by analogy:

$`{\displaystyle \frac{1}{n^2}}(𝝈_i𝝈_j)(𝝈_k𝝈_l)`$ $`=`$ $`v^{Nl}u^{lk}v^{kj}u^{ji}v^{i1}v^{1N}`$ (68)
$`=`$ $`(u/v)^{(lk)+(ji)},`$
where $`i<j<k<l`$ is understood. For the special case of energy-energy correlations one has:

$$\frac{1}{n^2}(𝝈_i𝝈_{i+1})(𝝈_j𝝈_{j+1})=(u/v)^2,$$
(69)
which does not depend on the distance $`|ji|`$. Hence the connected energy-energy correlation function vanishes exactly even for finite $`n`$ in one dimension. The $`n\mathrm{}`$ limit of this expression is given by:

$$\frac{1}{n^2}(𝝈_i𝝈_{i+1})(𝝈_j𝝈_{j+1})=\frac{4K^2}{[1+\sqrt{1+(2K)^2}]^2}.$$
(70)