format: text

text:
# Magnetic fluctuations in the classical XY model: the origin of an exponential tail in a complex system.

## I Introduction

### A Motivation for the Present Work

The fluctuations in a global measure of a many body system are often assumed to be of Gaussian form about the mean value . This assumption is nearly always true if the system in question can be divided into statistically independent microscopic or mesoscopic elements , as dictated by the central limit theorem (see App. A). However, in correlated systems, where this is not the case, there is no universal reason to expect the central limit theorem to apply. The fluctuations can then take on a multitude of different mathematical forms, including those of other, well-defined, limit distributions.

In this context, the most studied correlated systems are critical systems. At the critical point of a second order phase transition, a correlation length, $`\xi `$, diverges from the microscopic scale (taken as unity through out the paper). It is only cut off, in an ideal world, by the macroscopic, or integral scale $`L`$. The probability density function (PDF) for the order parameter $`m`$ associated with the diverging correlation length is essentially the exponential of the free energy $`P(m)\mathrm{exp}(F(m)/k_\mathrm{B}T)`$ and takes on an approximately Gaussian form as long as the Landau approximation, $`F(m)a+bm^2+\mathrm{}`$, is valid. Close to the critical point the Landau approximation breaks down and the PDF becomes non-Gaussian. The key assumption of the renormalisation group theory of critical phenomena is that the critical PDF remains scale invariant in the thermodynamic limit and can be obtained from the fixed point of a renormalisation group transformation  (see App. A). Thus, renormalisation group theory can be regarded as a generalisation of the central limit theorem to systems that are correlated over all length scales. The critical PDFs can be termed “universal”, in that, when properly normalised, they depend on at most a few basic symmetries that define the universality class of the system. A non-Gaussian and universal PDF is therefore a direct signature of the fluctuation driven critical phenomena that have revolutionised modern statistical mechanics . Analytical and numerical work  on Ising, Potts, and XY models has shown that a generic feature of such systems is a skewness, with large fluctuations below the mean, towards small order parameter values.

Correlations that are both strong and long range are a feature not only of critical phenomena but also of systems driven far from equilibrium. However, in the case of driven systems, the absence of a microscopic theory means that one has to rely heavily on empirical observations from experiment and numerical simulation. Labbé et al.  have shown that the PDF for the energy injected into a closed turbulent flow at constant Reynolds number is also non-Gaussian and universal. In this case “universal” means that the PDF, when suitably normalised, does not depend on Reynolds number or several other parameters (for example, the type of fluid). The PDF again has a marked skewness with an apparent exponential tail for fluctuations towards low energies.

The present work is motivated by our empirical observation  that the universal PDF of energy fluctuations in the turbulence experiment  is, within experimental error, of the same functional form as that of the universal $`P(m)`$ for the critical system that we have studied. The latter is the spin wave limit to the low temperature phase of the $`2D`$-XY model  that is known to capture the critical behaviour of the full $`2D`$-XY model . The distribution is shown in Fig. 2: it is asymmetric, with fluctuations below the mean approaching an exponential asymptote, while those above the mean approach a double exponential. This observation led us to the proposition that many systems, both equilibrium and non-equilibrium, sharing the property of long range correlations and multi-scale fluctuations, should share the same features, at least to a good approximation . The proposition appears to be strikingly confirmed in ref.  where, from numerical simulation, similar behaviour is observed in a number of different systems: for order parameter fluctuations in the two-dimensional Ising model and in the two-dimensional percolation problem, as well as as for fluctuations in global quantities for models of forest fires and avalanches, driven into a self-organised critical state. This appears to contradict the idea that the PDF should depend on the particular universality class of the model in hand. One possible way of accounting for our observations is that many universality classes share common features, with the differences between them appearing either outside the range of physical observation, or being hidden by experimental error. There are therefore many open questions regarding a possible and much desired connection between critical phenomena and non-equilibrium systems as well as regarding the details of the PDF in critical systems. It is these questions that we address in the current paper, via an analytic study of the PDF for order parameter fluctuations in a finite XY model in arbitrary dimension.

### B Normalisation of the Order Parameter

We discuss order parameter fluctuations of finite systems in terms of distributions that are calculated in the thermodynamic limit, $`N\mathrm{}`$. As discussed in App. A, it is essential to normalise the order parameter by an appropriate power of $`N=L^D`$ in order to obtain a distribution of finite width, or, equivalently, a form for $`P(m)`$ that is independent of system size. By extending the scaling hypothesis to include finite systems , the following form of $`P(m)`$ has been proposed:

$`P(m,L)L^{\beta /\nu }P_L(mL^{\beta /\nu },\xi /L).`$ (1)
Here $`\beta `$ and $`\nu `$ are the conventionally defined critical exponents for the magnetisation and correlation length $`\xi `$ respectively . The appropriate normalisation of the order parameter is provided by the factor $`L^{\beta /\nu }`$ while fixing different ratios $`\xi /L`$ will in principle result in an infinity of different limit distributions as $`L\mathrm{}`$. We concentrate on the case of a truly critical system with correlations over all length scales, which should result in maximum deviation from the Gaussian form. Here the dependence on $`\xi `$ can be dropped from eqn. (1), and $`P_L(m)`$ should closely approximate a single universal function of the variable $`mL^{\beta /\nu }`$ for all values of $`L`$. In this form it is independent of the microscopic details of the system, although it could indeed depend on the universality class of the transition through the critical exponents.

Equation (1) is demystified somewhat by recognising that the normalising factor $`L^{\beta /\nu }`$ is, in such an ideal system, proportional to the mean value of the order parameter, $`m`$. Further, one of the key properties of a critical system is that the standard deviation of the distribution, $`\sigma `$, scales with system size in the same way as the mean value. This property, which is a direct result of the hyperscaling relation and which we refer to as the hyperscaling condition, means that eqn. (1) can alternatively be written in the form $`P(m)=1/\sigma P_L(m/\sigma )`$. Thus $`\sigma `$ provides, as might be intuitively expected, the correct normalisation of the order parameter, such that a reasonable PDF of finite width is obtained in the limit $`N\mathrm{}`$. In this paper, in addition, we shift the distribution with respect to the mean value and define

$`\sigma P(m)`$ $`=`$ $`\mathrm{\Pi }(\theta ),`$ (2)
where

$`\theta ={\displaystyle \frac{mm}{\sigma }}.`$ (3)
In this representation one expects the PDF to fall, in the thermodynamic limit, onto a single universal curve. Provided that finite-size corrections to scaling are negligible one should observe data collapse onto this universal curve for large but finite system sizes.

### C The Two Dimensional XY Model

The model that we study, the harmonic spin wave limit to the XY model, is defined in section II for the case of two dimensions, $`D=2`$. This is the dimension of most interest in the present context, as the system is at its lower critical dimension. At low temperature the coupling $`J/k_\mathrm{B}T`$ is an exactly marginal variable that characterises a line of critical points in zero applied field . The critical line is separated from the paramagnetic phase by the Kosterlitz-Thouless-Berezinskii phase transition at $`T_{\mathrm{KTB}}`$ . The critical phase that exists below this temperature is an attractive subject of investigation from both an analytic and a numerical point of view. Its physics is entirely captured by the harmonic Hamiltonian  with the result that many calculations can be performed microscopically, without the need to use renormalisation techniques, or the scaling hypothesis. From a numerical point of view simulation results near a single, isolated, critical point are often complicated by a shift in the effective critical temperature by an amount scaling to zero as $`L^{1/\nu }`$ , making it unclear exactly which temperature should be studied. Indeed numerical studies of Ising and Potts models  do find distributions whose form depends on temperature in the critical region (see App. A). In the $`2D`$-XY model, as the system is critical over a range of temperatures there are no such technical problems and data for $`P(m)`$ can be collected at all points below $`T_{\mathrm{KTB}}`$. These factors make the $`2D`$-XY model an ideal system with which to study the effects of critical correlations.

The finite-size scaling for the $`2D`$-XY model has been discussed in our previous publications . In this work, we began, following Berezinskii  and Rácz and Plischke , an exact calculation for the PDF of order parameter fluctuations. This calculation is completed and presented in detail in the current paper (section II). It shows explicitly that the non-Gaussian behaviour in the $`2D`$-XY model stems from the influence of all length scales from the microscopic to the macroscopic scale. We propose that the same is true for other complex systems including those driven far from equilibrium, which provides a basis for understanding the apparent overlap of their PDFs and provides an unexpected experimental motivation for studying a system as simple as the $`2D`$-XY model.

Two results coming out of our calculation are worthy of note at this stage. The first is an exact analytic result that is rather surprising given the previous discussion and the general belief concerning the dependence of the PDF on universality class: shifting the curve with respect to the mean, eqn. (2), gives us universal data collapse, not only for all system size but also for all temperatures for which the harmonic Hamiltonian is valid. The ratio of exponents $`\beta /\nu `$ depends linearly on temperature, from which we deduce that the PDF is independent of the value of the exponents along the line of critical points. One should note however that these points are rather special and the result cannot necessarily be generalised to all critical points: not all the usual critical exponents are defined. For example the exponents $`\beta `$ and $`\nu `$ are not individually defined, but their ratio is . The usual scaling relations are valid in terms of the ratio $`\beta /\nu `$ only and this “weak scaling”  means that there is only one independent critical exponent, $`\eta =2\beta /\nu `$ , compared with two for a regular critical point. This is all that is required for the analysis leading to eqn. (1), but is not sufficient to ensure a unique functional form for the general problem with two exponents. However, it does seem consistent with the idea that only small differences separate results for different universality classes.

The second result that it is relevant to mention at this stage concerns the finite size scaling data collapse of eqn. (2). We find that the hyperscaling property $`\sigma m`$ is not a necessary condition for data collapse onto a non-Gaussian function. With our definition (2), the first two moments of $`P(m)`$ fall trivially out of the calculation and all that is required for data collapse is that the moments $`\theta ^p`$ for $`p>\mathrm{\hspace{0.33em}2}`$, are independent of the system size. This is the most general condition for non-Gaussian data collapse, while the PDF only satisfies the scaling hypothesis in the form of eqn. (1), if the hyperscaling condition is satisfied. We give, in section II.A, an explicit example where data collapse onto the universal curve of the $`2D`$-XY model occurs, but where the hyperscaling condition is not satisfied. If we make an expansion of the order parameter about a perfectly ordered state ($`m=1`$) in powers of temperature, keeping only the linear term, then $`m`$ diverges logarithmically with system size , while the standard deviation is a constant. The ratio $`\sigma /m`$ is actually an increasing function of system size throughout the physical domain. It is only when the order parameter is correctly defined on the interval $`[0,1]`$ that the hyperscaling relation is re-established, but written in the form (2) the two distinct variables have the same universal PDF, even outside the range of temperature and system size for which the development gives an accurate representation of the true order parameter. This result is more than a mathematical curiosity; the harmonic approximation for the $`2D`$-XY model maps directly onto the Edwards-Wilkinson (EW) model  for interface growth and the linearised order parameter is related to the square of the interface width, $`w`$: $`m=1w^2`$. Our PDF therefore corresponds precisely to that for interface width fluctuations and for which, in two dimensions, the hyperscaling condition for the observable $`w^2`$ is explicitly violated.

### D Organisation of the Paper

The rest of the paper is organised as follows. In section II we present detail of the calculation for the PDF in the $`2D`$-XY model *** For convenience, throughout this paper we use the term “XY model” to refer to either the model defined by the spin wave Hamiltonian or the full XY model over a temperature range in which the spin wave approximation is valid. This should not cause any confusion in reading the present paper, but our choice of terminology should be born in mind when comparing to other work on the XY model. We show explicitly that it is a universal function of system size and of temperature and find an exact expression for the characteristic function (section II A). Transforming the distribution numerically, we compare it in detail with extensive Monte Carlo and molecular dynamics simulations of the full XY model and show that it is clearly the complete solution of the problem (section II B). We calculate the asymptotic values of the distribution for large deviations below and above the mean, which we find to be exponential and double exponential respectively (sections II C, and II D).

In section III we try to fit the computed PDF to standard functions by comparing the moment expansion of the generating function with those of the Fourier transform of the test function. Three functions are considered:

$`\mathrm{\Pi }(\theta )`$ $``$ $`\{\begin{array}{cc}\mathrm{exp}a\left[\theta s\mathrm{exp}(\theta s)\right],\hfill & \\ & \\ \frac{1}{s\theta }\mathrm{exp}\left\{\left[\mathrm{log}(s\theta )a\right]^2\right\},\hfill & \\ & \\ (s\theta )^{\nu /21}\mathrm{e}^{a(s\theta )}.\hfill & \end{array}`$ (9)
The PDF is fitted to an excellent approximation over the physical range by the first two functions, while the third give a reasonable but slightly inferior fit. The first function, with $`a`$ an integer, comes from extremal statistics (section III A). It is Gumbel’s third asymptote, corresponding to the PDF for the $`a^{th}`$ largest value from ensembles of $`N`$ random numbers . The interpretation with $`a`$ non-integer (we find $`a=\pi /2`$) is not clear, but a connection between critical phenomena and extremal statistics is a very appealing concept . The second function is a generalised log-normal distribution (section III B). Unlike the first curve, it does not have the correct asymptotic forms but despite this it fits just as well over the physical domain. The third function is a $`\chi ^2`$ distribution describing identical and statistically independent degrees of freedom (section III C). It gives reasonable qualitative agreement indicating that a good, zeroth order description of a correlated system is in terms of a reduced number of statistically independent variables. However, this description has its limits, as shown by the fact that this function fits the exact PDF slightly less well than the other two. This variety of different fits suggests that one should treat the physical interpretations that they offer with caution; however even with this caveat in mind they still represent useful mathematical tools. To investigate this point further, in section III D we derive an approximate functional form for the curve using an analysis due to Pearson which reconstructs the PDF from the four principle moments, which in this case have been calculated analytically. The Pearson analysis gives a quite different function which also gives a good description of the exact PDF over a physical range of $`\theta `$. This serves to emphasise that, given zero mean and unit variance, the shape of the curve over a typical experimental range is essentially defined by its skewness, $`\gamma `$, and kurtosis, $`\kappa `$. Therefore, an alternative way of summarising the observed universality , is that $`\gamma `$ and $`\kappa `$, for several different systems, have the same scale-invariant values as they do for the XY model.

In section IV we extend our calculation to $`D`$ dimensions, which apart from $`D=2`$ are all non-critical. Despite this, we find evidence of the integral scale for all dimensions $`D<4`$. For $`D=1`$, the PDF for the linearised order parameter shows an exponential tail. However, we show numerically that the PDF for the correctly defined order parameter is quite different and is just what one would expect for a paramagnetic system without correlations (section IV A). The case $`D=3`$ holds a final surprise (section IV B): despite the long range order of the low temperature phase the PDF is still not a Gaussian function. The temperature is a dangerously irrelevant variable in the ordered phase of the $`3D`$-XY model with the result that the susceptibility remains weakly divergent at low temperature . The result of this divergence is that the asymptotes of the PDF for large fluctuations are exponential below the mean and $`\mathrm{exp}(\mathrm{Cst}\theta ^3)`$ above the mean. The hyperscaling relation, in this case, is again violated. The divergence disappears at the upper critical dimension and we find a truly Gaussian PDF for $`D4`$.

In section V we conclude by returning to the physical reasons for the exponential tail in the PDF. The XY model is diagonisable in reciprocal space reducing it to a model of statistically independent degrees of freedom: spin wave amplitudes at wave vector $`𝐪`$, $`\varphi _𝐪`$. The amplitudes $`\varphi _𝐪^2`$ diverge at small $`q`$ and are these modes that give the non-Gaussian fluctuations. In one-dimension they completely destroy magnetic order, in two-dimensions they give critical behaviour and between two and four dimensions they give remnant critical behaviour in the form of a dangerously irrelevant variable.

## II Probability density function for the order parameter in the $`2D`$-XY model

### A Analytic Expression

The $`2D`$-XY model is defined by the Hamiltonian

$`H=J{\displaystyle \underset{i,j}{}}\mathrm{cos}(\theta _i\theta _j),`$ (10)
where the angles $`\theta _i`$ refer to the orientation of classical spins $`𝐒_i`$ confined in a plane and where the sum is over nearest neighbour spins. In the following we consider a square lattice of side $`L`$, with periodic boundaries. The order parameter is a two dimensional vector $`𝐦`$ which, in zero field is free to point in any direction. We define the instantaneous magnetisation as the scalar $`m=|𝐦|`$

$`m={\displaystyle \frac{1}{N}}{\displaystyle \underset{i=1}{\overset{N}{}}}\mathrm{cos}(\theta _i\overline{\theta }),`$ (11)
where $`\overline{\theta }=\mathrm{tan}^1(_i\mathrm{sin}\theta _i/_i\mathrm{cos}\theta _i)`$ is the instantaneous magnetisation direction. Within small corrections, which disappear in the thermodynamic limit, this corresponds to the more conventional definition

$$m=\frac{1}{N}\sqrt{\left(\underset{i=1}{\overset{N}{}}𝐒_i\right)^2.}$$

For all temperatures below $`T_{\mathrm{KTB}}`$ the renormalisation group trajectories flow, at large length scale, towards a regime where only spin-wave excitations are relevant . The physics of the low temperature phase is therefore completely captured by the quadratic Hamiltonian

$`H={\displaystyle \frac{J}{2}}{\displaystyle \underset{i,j}{}}\left(\theta _i\theta _j\right)^2.`$ (12)
We therefore restrict ourselves, in the following calculation to this Hamiltonian and neglect the periodicity of the variables $`\theta _i`$. Our calculation can not therefore take into account the presence of vortex pairs. Close to $`T_{\mathrm{KTB}}`$ in two dimensions and also in one dimension, where free vortices are relevant variables we would expect a deviation from the behaviour shown in Fig. 2. This point is discussed further below.

We now calculate the PDF $`P(m)`$ that the system be in a state with magnetisation $`m`$, using the standard property that a probability density function is entirely defined by the value of its moments . Indeed, $`P(m)`$ can be expressed in terms of its characteristic function, $`\stackrel{~}{P}(x)`$:

$`P(m)={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{2\pi }}\mathrm{e}^{imx}\stackrel{~}{P}(x),`$ (13)
which can in turn be expanded in a Taylor series whose coefficients are the moments $`m^p`$ :

$`\stackrel{~}{P}(x)={\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{x^p}{p!}}{\displaystyle \frac{^p\stackrel{~}{P}}{x^p}}|_{x=0}={\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(ix)^p}{p!}}m^p,`$ (14)
so that

$`P(m)={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{2\pi }}\mathrm{e}^{imx}{\displaystyle \underset{p=0}{\overset{\mathrm{}}{}}}{\displaystyle \frac{(ix)^p}{p!}}m^p.`$ (15)
Eqn. (15) assumes that the series converges and that all the moments exists. Note that this last feature demands that $`P(m)`$ falls off faster than any power-law of $`m`$.

The program for calculating $`P(m)`$ is therefore to calculate the moments $`m^p`$, sum the series and transform the final result. To this end it is useful to define the Green function in Fourier space,

$`G(𝐪)`$ $`=`$ $`{\displaystyle \frac{1}{42\mathrm{cos}q_x2\mathrm{cos}q_y}},`$ (16)
where $`q_x`$ and $`q_y`$ take the discrete values $`\frac{2\pi }{L}n`$ of the Brillouin zone with $`n=0,\mathrm{}L1`$. We also define the set of constants $`g_k=_𝐪G(𝐪)^k/N^k`$. The value of $`g_1`$ diverges logarithmically with system size, illustrating the critical nature of the low temperature phase : $`g_1=\frac{1}{4\pi }\mathrm{log}(CN),C=1.8456`$ . The values of $`g_k,k2`$ are independent of $`N`$ in the thermodynamic limit. We find: $`g_2\mathrm{3.8667\hspace{0.33em}10}^3`$, $`g_3\mathrm{7.5719\hspace{0.33em}10}^5`$, $`g_4=\mathrm{1.7626\hspace{0.33em}10}^6`$ and that for large $`k`$, $`g_k`$ behaves like $`(2\pi )^{12k}/2(k1)`$, see App. B.

The first moment is easily calculated within this approximation (see App. B and refs. ). One finds that $`m`$ decreases algebraically with the size, as one would expect from finite-size scaling

$`m`$ $`=`$ $`(NC)^{k_\mathrm{B}T/8\pi J}.`$ (17)
As discussed above, while the critical exponents $`\beta `$ and $`\nu `$ are not individually defined for the $`2D`$-XY model, their ratio is  and the system obeys what Kosterlitz referres to as weak scaling . Through eqn. (17) the ratio of exponents is defined: $`\beta /\nu =\eta /2=T/4\pi J`$.

For higher moments we need a more systematic approach. A specific property of the quadratic Hamiltonian (12) is that the moments can be calculated using the tools of Gaussian integration . In particular, by the application of Wick’s theorem, propagators of order $`2p`$ in reciprocal space can be exactly expressed in terms of quadratic propagators so that the $`p^{th}`$ moment is proportional to $`m^p`$. One finds

$`m^p`$ $`=`$ $`m^p{\displaystyle \frac{1}{(2N)^p}}{\displaystyle \underset{𝐫_1,\mathrm{},𝐫_p}{}}{\displaystyle \underset{\sigma _1,\mathrm{},\sigma _p=\pm 1}{}}\mathrm{exp}\left[{\displaystyle \frac{\tau }{2}}{\displaystyle \underset{ij}{}}\sigma _iG_R(𝐫_i𝐫_j)\sigma _j\right],`$ (18)
where $`\tau `$ is the reduced temperature $`k_\mathrm{B}T/J`$ and $`G_R(𝐫)`$ the regularised Green function $`_{q0}G(𝐪)\mathrm{exp}(i𝐪𝐫)/N`$. In order to compute each moment of order $`p`$, we have to evaluate the sums over the positions and operators $`\sigma _i`$. The idea is to expand the exponential term (18) and introduce a diagrammatic representation of each quantity computed. For example, we represent $`\sigma _iG_R(𝐫_i𝐫_j)\sigma _j`$ by a line between $`i`$ and $`j`$ on a lattice of $`p`$ sites. The general term of the expansion is then a set of graphs with a combinatorial factor for the symmetries. Since $`\sigma _i^2=1`$, only closed diagrams are relevant, the factor $`2^p`$ being cancelled by the sum over all the $`\sigma _i`$. The factor of $`\tau ^k`$, is common to all graphs with $`k`$ lines connected together, with an even connectivity at each vertex. For example, up to the second order term in $`\tau `$, we have

$`m^p`$ $`=`$ $`m^p\left[1+\left({\displaystyle \frac{\tau }{2}}\right)^2{\displaystyle \frac{1}{2!}}2p(p1){\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐫}{}}G_R^2(𝐫)+\mathrm{}\right].`$ (19)
The term $`_𝐫G_R^2(𝐫)/N^2=_{𝐪\mathrm{𝟎}}G(𝐪)^2/N^2=g_2`$ is the value of the one loop graph with two lines, as shown in Fig. 1a. There is an additional symmetry factor $`2\times p(p1)`$, which is the number of possible positions for such diagrams connecting 2 lines on a closed graph on a lattice of $`p`$ points. For the third order term in $`\tau `$, we have only one diagram with 3 vertices, of value $`g_3`$, Fig. 1b. The symmetry factor is equal to $`p(p1)(p2)\times 4\times 2`$. The factor $`4\times 2`$ comes from the number of possible ways of connecting 3 lines together. For the $`4^{th}`$ order term there are three different graphs, two of which are shown in Fig. 1c. The first has 3 loops and 2 vertices, the second, of value $`g_4`$, has one loop and 4 vertices. The third graph, not shown, consists of two disconnected one loop graphs of the type shown in Fig. 1a. In general, at each order in $`\tau `$, we have the product of different closed diagrams, with one or many loops. It appears that the values of multiple loop graphs, like the first one in Fig. 1c), are zero in the thermodynamic limit. We therefore find that only the one loop diagrams are relevant and the value for such a diagram, with $`k`$ lines and $`k`$ vertices is $`g_k`$. We can now express the $`p^{th}`$ moment of the magnetisation as

$`{\displaystyle \frac{m^p}{m^p}}`$ $`=`$ $`1+{\displaystyle \underset{k2}{}}\left({\displaystyle \frac{\tau }{2}}\right)^k{\displaystyle \frac{1}{k!}}{\displaystyle \underset{r1}{}}{\displaystyle \underset{\stackrel{k_1+\mathrm{}+k_r=k}{k_i2}}{}}g_{k_1}\mathrm{}g_{k_r}C(k_1,\mathrm{},k_r)\times p(p1)\mathrm{}(pk+1),`$ (20)
with $`C(k_1,\mathrm{},k_r)`$ a combinatorial factor which takes into account the possible ways of putting together $`k`$ lines on $`r`$ graphs, the first with $`k_1`$ lines, the second with $`k_2`$ lines, etc…, including the symmetries. For example, the factor associated with one triangle is $`C(3)=4\times 2`$. It is then relatively easy to show that

$`C(k_1,\mathrm{},k_r)={\displaystyle \frac{2^{kr}k!}{(k_1+\mathrm{}+k_r)(k_2+\mathrm{}+k_r)\mathrm{}k_r}}.`$ (21)
Next, we can use the fact that every diagram is invariant by the action of the group $`𝒮_r`$ of permutations of its $`r`$ single elements, so that, instead of (21), one can use a more convenient form for the combinatorial factor:

$`{\displaystyle \frac{1}{r!}}{\displaystyle \underset{\sigma 𝒮_r}{}}C(k_{\sigma (1)},\mathrm{},k_{\sigma (r)})={\displaystyle \frac{1}{r!}}{\displaystyle \frac{2^{kr}k!}{k_1\mathrm{}k_r}}.`$ (22)
Setting $`f_{k_i}=g_{k_i}(\tau )^{k_i}/2k_i`$ we arrive at the result

$`{\displaystyle \frac{m^p}{m^p}}`$ $`=`$ $`1+{\displaystyle \underset{k2}{}}{\displaystyle \underset{r=2}{\overset{k}{}}}{\displaystyle \frac{1}{r!}}{\displaystyle \underset{k_1+\mathrm{}+k_r=k}{}}f_{k_1}\mathrm{}f_{k_r}\times p(p1)\mathrm{}(pk+1)`$ (23)
$`=`$ $`1+{\displaystyle \underset{k2}{}}{\displaystyle \underset{r=2}{\overset{k}{}}}{\displaystyle \frac{1}{r!}}{\displaystyle \underset{k_1+\mathrm{}+k_r=k}{}}f_{k_1}\mathrm{}f_{k_r}{\displaystyle \frac{^k}{z^k}}z^p|_{z=1}`$ (24)
$`=`$ $`\mathrm{exp}\left[{\displaystyle \underset{k=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{g_k}{2k}}(\tau )^k_z^k\right]z^p|_{z=1}.`$ (25)
For $`p=2`$ we find $`m^2/m^2=1+g_2\tau ^2/2`$ and defining $`\sigma =\sqrt{m^2m^2}`$ we thus arrive at the hyperscaling condition that the ratio $`\sigma /m`$ is independent of the system size. Hence

$`\sigma `$ $`=`$ $`\sqrt{{\displaystyle \frac{g_2}{2}}}{\displaystyle \frac{k_\mathrm{B}T}{J}}m.`$ (26)
One can now substitute for $`m^p`$ in (15) using (23) and after re-arranging the summations the distribution can finally be expressed as an integral, depending on the values of the one loop diagrams $`g_k`$ only

$`P(m)={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{2\pi }}\mathrm{exp}\left[ix\left(mm\right)+{\displaystyle \underset{k=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{g_k}{2k}}(i\tau mx)^k\right].`$ (27)
Changing variables, $`xx/\sigma `$ and using (26) we find

$`P(m)={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{2\pi \sigma }}\mathrm{exp}\left[ix{\displaystyle \frac{mm}{\sigma }}+{\displaystyle \underset{k=2}{\overset{\mathrm{}}{}}}{\displaystyle \frac{g_k}{2k}}\left(ix\sqrt{{\displaystyle \frac{2}{g_2}}}\right)^k\right],`$ (28)
which is the principal result of ref. . Defining $`\sigma P(m)=\mathrm{\Pi }(\theta )`$, we see that the function $`\mathrm{\Pi }`$ depends uniquely on the variable $`\theta =(mm)/\sigma `$ and the $`g_k,k2`$. As the $`g_k`$ are constants in the thermodynamic limit, $`\mathrm{\Pi }(\theta )`$ is a universal function, independent of both system size and temperature. The asymmetry comes from the fact that the ratios $`g_k/(g_2/2)^{k/2},k3`$ are non zero and this constitutes the abnormal influence of the integral scale. If, in the thermodynamic limit, $`k=2`$ were the only non-zero term one would arrive at a Gaussian PDF centered on $`m`$. Departure from a Gaussian function is typically characterised by the skewness, $`\gamma =\theta ^3`$ and kurtosis, $`\kappa =\theta ^4`$ . We find

$`\gamma `$ $`=`$ $`{\displaystyle \frac{g_3}{(g_2/2)^{3/2}}}=0.8907,`$ (29)
$`\kappa `$ $`=`$ $`3+3{\displaystyle \frac{g_4}{(g_2/2)^2}}=4.415.`$ (30)
Although we can calculate the asymptotic behaviour of $`g_k`$ for large $`k`$, we are not able to compute the constants analytically and so we cannot sum the series (28). However we can transform it into a very much more useful form by keeping $`N`$ large but finite and inverting the sums over $`𝐪`$ and $`k`$. The even and odd terms are separated and summed independently and we eventually find:

$`\mathrm{\Pi }(\theta )={\displaystyle _{\mathrm{}}^{\mathrm{}}}\sqrt{{\displaystyle \frac{g_2}{2}}}{\displaystyle \frac{dx}{2\pi }}\mathrm{exp}\left\{ix\theta \sqrt{{\displaystyle \frac{g_2}{2}}}{\displaystyle \underset{𝐪\mathrm{𝟎}}{}}\left[{\displaystyle \frac{i}{2}}xG(𝐪)/N{\displaystyle \frac{i}{2}}\mathrm{arctan}\left(xG(𝐪)/N\right)+{\displaystyle \frac{1}{4}}\mathrm{log}\left(1+x^2G(𝐪)^2/N^2\right)\right]\right\}.`$ (31)
The sum over $`𝐪`$ and the integral over $`x`$, in (31) can now be performed numerically, allowing the evaluation of $`\mathrm{\Pi }(\theta )`$.

### B Comparison with Simulation

To test the above calculation and to verify its scaling properties we have carried out extensive numerical simulations of the $`2D`$-XY model with full cosine interaction, eqn. (10), for different values of temperature and system size. In addition, we have also done microcanonical molecular dynamics (MD) simulations to check the possible dependence of the PDF for fluctuations on the statistical ensemble.

The Monte Carlo simulations were performed with $`10^8`$ Monte Carlo steps per spin, with $`10^6`$ steps used for equilibration. The MD simulation was carried out for systems of $`N`$ classical rotators, , with Hamiltonian

$`H_R={\displaystyle \underset{i=1}{\overset{N}{}}}{\displaystyle \frac{\dot{\theta }_i^2}{2}}+J{\displaystyle \underset{ij}{}}\left[1\mathrm{cos}(\theta _i\theta _j)\right].`$ (32)
The equations of motion were integrated numerically, using a Verlet algorithm. In order to explore the low-temperature fluctuations regime the initial configuration of the system was chosen with the spins pointing in the same direction and with a Gaussian distribution of momenta. The system was then equilibrated for a time of $`10^610^7`$ sweeps and data collected over a time span of $`10^810^9`$ sweeps according the size of the system. Note that one cannot use the harmonic interaction (12) to study deterministic dynamics in the microcanonical ensemble, as this would allow no coupling between the spin wave modes and no evolution would be possible. The non-linearity of the cosine interaction allows mixing between the normal modes and the sampling of equilibrium states. Here we do not report work at high enough energies to allow vortex formation  with any significant probability. Rather, the non-linearity plays the role of the heat bath in the canonical ensemble, while the physics is still correctly described by the harmonic part of the interaction.

The numerical integration of eqn. (31), performed with a fast Fourier transform (FFT) algorithm , is shown in Fig. 2, where it is compared with Monte Carlo results for $`T/J=0.1`$ and $`N=32^2`$. The theoretical curve is clearly in extremely good agreement with the numerical data. The curve is asymmetric, with what appears to be an exponential tail for fluctuations below the mean, with a much more rapid fall off in amplitude, for fluctuations above the mean.

In Fig. 3 we show the PDF for fluctuations in $`m`$ obtained from MC simulation for fixed system size and varying temperature, as well as MD for fixed temperature and different system sizes. The result of ref.  and section II of this paper is that, for the harmonic Hamiltonian, eqn (12), $`\mathrm{\Pi }(\theta )`$ is independent of both system size and temperature, while we have explicitly tested this result against the PDF generated for the full Hamiltonian, eqn. (10). Qualitative agreement is clearly excellent, independently of the ensemble used, but there are small systematic deviations in the tails, when observed on a logarithmic scales , as shown in Figs. 4 and 5. We can only expect agreement between the analytic result and simulation in the range of temperature sufficiently below $`T_{\mathrm{KTB}}`$ such that vortex pairs do not influence the PDF . Even in the absence of vortices one must expect small variations from our theoretical result for small system sizes that stem from the utilisation of (10) rather than (12). In a renormalisation group treatment the non-linearities of Hamiltonian (10) scale away on changing the length scale and the Hamiltonian is replaced by an effective harmonic Hamiltonian at higher temperature . For example, at $`T/J=0.7`$, for $`L=32`$, we find $`m=0.76`$ from simulation, while eqn. (17) gives $`m=0.81`$. The effective coupling constant can be calculated by expanding the cosine and approximating the nonlinear terms using a Hartree approximation . Renormalisation of the non-linearities introduces a microscopic length scale $`a^{}`$ which gives small corrections when compared with the calculated PDF. However, this length scale is fixed by the temperature and the corrections should scale away as the ratio $`a^{}/L0`$. This scenario is confirmed in Fig. 5, where data are shown at $`T/J=0.7`$, for $`L=8,\mathrm{\hspace{0.33em}16},\mathrm{\hspace{0.33em}32},\mathrm{\hspace{0.33em}64}`$ and compared with the theoretical curve. Deviations from the theoretical result are observed for $`L=8`$ and $`L=16`$ but the PDF clearly approaches the predicted scale independence for the larger system sizes.

Near $`T_{\mathrm{KTB}}`$ vortices influence the PDF, however the vortex population decreases exponentially moving away from $`T_{\mathrm{KTB}}`$  and they only make their presence felt within the physical domain in a small band of temperatures near the transition. In this regime the data do not fit on the universal curve  but a detailed discussion of this point is outside the scope of this paper.

### C $`P(m)`$ for the Linearised Order Parameter

As eqn. (31) is independent of temperature, one should be able to obtain it at low temperature where the magnetisation is approximately

$`m`$ $`=`$ $`1{\displaystyle \frac{1}{2N}}{\displaystyle \underset{i}{}}(\theta _i\overline{\theta })^2.`$ (33)
In fact, using this expression one can arrive at (31) in a more straightforward manner. What is perhaps surprising is that the calculation, using (33) is valid for all temperatures below $`T_{\mathrm{KTB}}`$, even for temperatures where (11) and (33) represent different physical quantities.

Using the Hamiltonian (12), we have

$`P(m)={\displaystyle \frac{1}{Z}}{\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{2\pi }}{\displaystyle \underset{i}{}d\theta _i\mathrm{exp}\left\{ix\left[m1+\frac{1}{2N}\underset{i}{}\theta _i^2\right]\frac{1}{2\tau }\underset{i,j}{}\theta _iG_{ij}^1\theta _j\right\}}`$ (34)
where $`G_{ij}^1`$ is the inverse Green’s function operator connecting sites $`i`$ and $`j`$ with non-zero elements for $`i`$ and $`j`$ nearest neighbours , and $`Z=\left(detG^1/2\pi \tau \right)`$ is the partition function.

It is easy to integrate the Gaussian integral by transforming into reciprocal space. Defining the trace Tr of any function of $`G`$ as the sum for $`𝐪0`$ of the same function of $`G(𝐪)`$ and and using $`m=1\tau \mathrm{Tr}G/2N`$ we find

$`P(m)={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{2\pi }}\mathrm{exp}\left[ix\left(mm\right)ix{\displaystyle \frac{\tau }{2}}\mathrm{Tr}G/N{\displaystyle \frac{1}{2}}\mathrm{Tr}\mathrm{log}\left(\mathrm{𝟏}ix\tau G/N\right)\right].`$ (35)
We can now use the fact that $`\sigma =\sqrt{g_2/2}\tau `$ in this approximation, to transform (35) into a dimensionless and universal form

$`\mathrm{\Pi }(\theta )`$ $`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\sqrt{{\displaystyle \frac{g_2}{2}}}{\displaystyle \frac{dx}{2\pi }}\mathrm{exp}\left[ix\theta \sqrt{{\displaystyle \frac{g_2}{2}}}i{\displaystyle \frac{x}{2}}\mathrm{Tr}G/N{\displaystyle \frac{1}{2}}\mathrm{Tr}\mathrm{log}\left(\mathrm{𝟏}ixG/N\right)\right]`$ (36)
$`=`$ $`{\displaystyle _{\mathrm{}}^{\mathrm{}}}\sqrt{{\displaystyle \frac{g_2}{2}}}{\displaystyle \frac{dx}{2\pi }}\mathrm{exp}\left[i\mathrm{\Phi }(x)\right]`$ (37)
which is the same expression as (28,31), once we separate the real and imaginary parts of the integrand. This demonstration proves that the only relevant graphs are those with only one loop, the others being zero in the thermodynamic limit.

Within this linear approximation the mean magnetisation $`m`$ and the standard deviation $`\sigma `$ do not scale in the same way with system size: while $`m=1(T/8\pi J)\mathrm{log}(CN)`$, $`\sigma =\sqrt{g_2/2}\tau `$ is a temperature dependent constant. This exact result can be verified by applying eqn. (14) to eqn. (35) and calculating $`m`$ and $`m^2`$ directly. The fact that we find the same universal function for the two calculations, when written in the form (2) shows explicitly that the hyperscaling result, $`\sigma /mO(1)`$ is not a necessary condition for non-Gaussian data collapse. Rather, it seems that hyperscaling is a consequence, in these circumstances, of the correct definition of $`m`$ as an order parameter on the interval $`[0,1]`$.

The Gaussian limit of the $`2D`$-XY model is identical to the Edwards-Wilkinson model of interface growth and the linear approximation for the order parameter is related to the square of the interface width $`m=1w^2`$. The PDF for $`w^2`$ has been studied in one  and two  dimensions together with extensions to the EW model, including non-linearity . All models give non-Gaussian PDF’s with the same qualitative features as Fig. 2. These models provide an important microscopic link between equilibrium and non-equilibrium systems and suggest that a formalism could exist that incorporates the statistical features that we have observed to be shared, at a global level, between such different systems.

### D Asymptotes of $`\mathrm{\Pi }(\theta )`$ for Large Fluctuations

As a first step towards an analytic form for $`\mathrm{\Pi }(\theta )`$ one can approximate (31) beyond the Gaussian approximation by retaining only the elements $`(g_2,g_3)`$. In this case, the solution is proportional to the Airy function

$`\mathrm{\Pi }(\theta )\mathrm{exp}\left({\displaystyle \frac{1}{6\alpha }}\theta \right)\mathrm{Ai}\left[{\displaystyle \frac{1}{(3\alpha )^{1/3}12\alpha }}{\displaystyle \frac{1}{(3\alpha )^{1/3}}}\theta \right],`$ (38)
where $`\alpha =2^{3/2}g_3/3g_2^{3/2}0.296876`$. The $`g_3`$ term assures that it is not symmetric on reversing the sign of $`\theta `$. We find that the approximation reproduces qualitatively the apparent exponential behaviour for $`\theta 1`$:

$`\mathrm{\Pi }(\theta )\left\{2\sqrt{\pi }\left[{\displaystyle \frac{1}{(3\alpha )^{1/3}12\alpha }}{\displaystyle \frac{1}{(3\alpha )^{1/3}}}\theta \right]^{1/4}\right\}^1\mathrm{exp}\left\{{\displaystyle \frac{1}{6\alpha }}\theta {\displaystyle \frac{2}{3}}\left[{\displaystyle \frac{1}{(3\alpha )^{1/3}12\alpha }}{\displaystyle \frac{1}{(3\alpha )^{1/3}}}\theta \right]^{3/2}\right\}`$ (39)
However the approximation does not allow us to extract the asymptote above the mean, as for $`\theta >0`$ the Airy function develops oscillations.

A more fruitful approach is to look at the saddle points of the integrand (31), from which one can extract both asymptotes. If $`\theta 1`$, an expansion near $`x=0`$ is not very satisfactory and we must rather seek the solution for the extrema of the whole integrand, $`\mathrm{\Phi }(x)/x=0`$. We find:

$`\sqrt{{\displaystyle \frac{g_2}{2}}}\theta ={\displaystyle \frac{1}{2}}\mathrm{Tr}{\displaystyle \frac{G^3}{N^3}}{\displaystyle \frac{x^2}{1+x^2G^2/N^2}}{\displaystyle \frac{i}{2}}\mathrm{Tr}{\displaystyle \frac{G^2}{N^2}}{\displaystyle \frac{x}{1+x^2G^2/N^2}}.`$ (40)
If $`\theta `$ is negative and $`x`$ real, the real part of the second term is always positive and there is no solution to this equation. We therefore seek a solution for $`x`$ pure complex, $`x=iy`$. In this case, eqn. (40) becomes

$`\sqrt{{\displaystyle \frac{g_2}{2}}}\theta ={\displaystyle \frac{1}{2}}\mathrm{Tr}{\displaystyle \frac{G^2}{N^2}}{\displaystyle \frac{y}{1+yG/N}}=\phi (y)`$ (41)
The function $`\phi `$ has simple poles at $`y=4\pi ^2,8\pi ^2,32\pi ^2,\mathrm{}`$ and its asymptotic value near the first pole $`y_0=4\pi ^2`$ is $`\phi (y)2/(yy_0)`$. The extremum of the integrand satisfies the condition $`y^{}y_02\sqrt{2/g_2}/\theta >y_0`$, for $`|\theta |`$ large and we can deform the real path of the integration so that it passes through the extremum on the imaginary axis. Near the extremum, we can expand the integrand up to second order in $`yy^{}`$ and perform a Gaussian integration:

$`\mathrm{\Pi }(\theta ){\displaystyle _{\mathrm{}}^{\mathrm{}}}\sqrt{{\displaystyle \frac{g_2}{2}}}{\displaystyle \frac{dx}{2\pi }}\mathrm{exp}\left[i\mathrm{\Phi }(iy^{})+i{\displaystyle \frac{1}{2}}(xiy^{})^2\mathrm{\Phi }^{\prime \prime }(iy^{})\right]`$ (42)
We finally find that the asymptotic value of the distribution varies as

$`\mathrm{\Pi }(\theta )|\theta |\mathrm{exp}\left(4\pi ^2\sqrt{{\displaystyle \frac{g_2}{2}}}\theta \right)`$ (43)
We have superimposed the asymptotic result (43) and the full numerical integration for $`N=101^2`$ of (31) in Fig. 6. The amplitude of eqn. (43) is chosen so that the curves are slightly displaced, to allow comparison of the slopes. The asymptotic solution is in excellent agreement even for $`\theta `$ values where the PDF shows a distinct deviation away from exponential behaviour and only fails for $`\theta >2`$. Further out in the tail, in the range $`10<\theta <4`$, $`\mathrm{log}(\mathrm{\Pi })`$ is approximately linear. However the value of the slope is not the argument of the exponential in (43), $`4\pi ^2\sqrt{g_2/2}1.736`$. The logarithmic corrections given by the the term $`|\theta |`$ are significant over the whole of this range, but the curvature is so small that the data can be fitted to an effective exponential $`\mathrm{\Pi }(\theta )\mathrm{exp}(\alpha \theta )`$, with $`\alpha =1.56867\mathrm{}`$. The data only approaches true exponential behaviour for $`\theta <30`$, which is completely outside any imaginable physical range. Strictly speaking it is therefore more correct to speak of pseudo-exponential, $`x\mathrm{exp}(\alpha x)`$, for the asymptote below the mean.

For large and positive $`\theta `$ a solution of eqn. (41) exists for large and positive $`y`$. A reasonable approximation is to replace $`G`$ by $`1/q^2`$ and perform the integration

$`\phi (y)`$ $``$ $`{\displaystyle \frac{1}{2}}{\displaystyle _{q=2\pi /\sqrt{N}}^{2\pi }}{\displaystyle \frac{Nd^2q}{4\pi ^2}}{\displaystyle \frac{1}{N^2q^4}}{\displaystyle \frac{y}{1+y/Nq^2}}`$ (44)
$``$ $`{\displaystyle \frac{1}{4\pi }}{\displaystyle _{2\pi /\sqrt{y}}^{\mathrm{}}}{\displaystyle \frac{dq}{q(1+q^2)}}{\displaystyle \frac{1}{8\pi }}\mathrm{log}y`$ (45)
A more precise computation gives $`\phi =\mathrm{log}(y)/8\pi +\widehat{a}+1/2y+\mathrm{}`$ where $`\widehat{a}`$ is a numerical constant which can be computed exactly. An analytical study (see App. C) gives

$`\widehat{a}={\displaystyle \frac{1}{24}}+{\displaystyle \frac{\gamma }{4\pi }}{\displaystyle \frac{1}{4\pi }}\mathrm{log}(4\pi ){\displaystyle \frac{1}{2\pi }}\mathrm{log}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}\left[1\mathrm{exp}(2\pi k)\right]=0.11351444337\mathrm{}`$ (46)
For large $`\theta `$, the saddle point of the integrand is therefore located at $`y^{}=\mathrm{exp}8\pi (\widehat{a}+\sqrt{\frac{g_2}{2}}\theta )`$, and the asymptotic value for $`\mathrm{\Pi }`$ follows from a Gaussian integration of (42):

$`\mathrm{\Pi }(\theta )\mathrm{exp}\left[{\displaystyle \frac{1}{8\pi }}\mathrm{e}^{8\pi \left(\sqrt{\frac{g_2}{2}}\theta \widehat{a}\right)}+8\pi \sqrt{{\displaystyle \frac{g_2}{2}}}\theta \right]`$ (47)
Comparing the asymptote with the full curve we again find that the true asymptote only fits accurately outside the physical domain, although the data is clearly consistent with a very rapid fall off in the PDF for $`\theta `$ above the mean.

## III Fitting to known functional forms

The obvious question now arises: is the PDF generated by the characteristic function (31) of known functional form? We do not have a definitive answer to this question, as we are not able to transform (31) analytically. In the absence of an answer, we test the PDF against three skewed functions, shown in eqn. (9), which describe statistics in different physical situations. These are: a modified Gumbel function, characteristic of problems where extreme values dominate the sum over many contributions; a log-normal distribution, characteristic of statistically independent multiplicative processes and a $`\chi ^2`$ distribution which describes the PDF of a quantity made up of a finite number of positive definite microscopic variables. The analysis is the same in all three cases, but is only shown in detail for the modified Gumbel function: each curve has 4 parameters, but once the value of the first is chosen the others are fixed by normalisation and the constraints $`\theta =0`$, $`\theta ^2=1`$. The family of one parameter curves are Fourier transformed and the first four terms in a Taylor expansion are set equal to those for the generating function, which fixes the value of the free parameter. The method takes into account the skewness of the curve but not the kurtosis and its accuracy is ultimately limited. The goodness of fit can be measured by comparing the ratio of higher order terms of the expansion of the test and generating functions. For an exact solution all higher ratios would be equal to unity, while for a poor fit they diverge rapidly from this value. Other functions could be tested in the same way and an exact solution may well exist in the statistics literature, unknown to us.

The method described above is quite similar to that due to Pearson , who realised a century ago that, in all practical situations, knowledge of the first four moments of a distribution is sufficient to generate a curve, fitting any set of data points . Pearson developed a phenomenological differential equation containing the numerical values of the moments, whose solution gives the fitting function. A Pearson analysis is performed on the calculated PDF at the end of the section.

### A The Generalised Gumbel Distribution

The asymptotes (43) and 47) are of the same general form as those for Gumbel’s first asymptote distribution from the theory of extremal statistics : defining $`z`$ to be the $`a^{th}`$ largest value from a set of $`z_i,i=1,N`$ random numbers taken from a generator $`f(z)`$; the PDF for $`z`$ is

$$g_a(z)=\frac{a^a\alpha _a}{\mathrm{\Gamma }(a)}\mathrm{exp}\left\{a\left[\alpha _a(zu_a)+\mathrm{e}^{\alpha _a(zu_a)}\right]\right\}.$$
(48)
$`\mathrm{\Gamma }(a)`$ is the gamma function; $`u_a`$ is the value of $`z`$ such that $`a`$ of the $`N`$ random numbers are greater than $`z`$. $`F(z)`$ is the probability of having $`a`$ of the values less than $`z`$, such that $`F(u_a)=1a/N`$. $`\alpha _a`$ is referred to as the intensity: $`\alpha _a=(N/a)f(u_a)`$. In conventional statistics, $`a`$ would of course be an integer. However, in what follows we are going to see an irrational number appearing.

The function (48) has an exponential tail for fluctuations towards large values of $`z`$, the opposite of the PDF, in Fig. 2. We therefore make a change of variables $`m_z=1z,\theta _z=(m_zm_z)/\sigma _z`$ which makes a mirror reflection of eqn. (48). Within the linear approximation for the order parameter this corresponds to the relevant variable being the sum of the spin wave amplitudes $`z(1/2N)_i(\theta _i\overline{\theta })^2`$ . Changing variables we find

$`\sigma _z\mathrm{\Pi }_G(\theta _z)`$ $`=`$ $`w\mathrm{e}^{ab(\theta _zs)a\mathrm{e}^{b(\theta _zs)}}`$ (49)
$`b`$ $`=`$ $`\alpha _a\sigma _z`$ (50)
$`s`$ $`=`$ $`(1m_zu_a)/\sigma _z`$ (51)
$`w`$ $`=`$ $`{\displaystyle \frac{a^a\alpha _a}{\mathrm{\Gamma }(a)}}\sigma _z.`$ (52)
Eqn. (49) is also the distribution for the $`a^{th}`$ smallest random number from the set $`z_i`$. After some algebra one can show that

$`b`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{\mathrm{\Gamma }(a)}}{\displaystyle \frac{^2\mathrm{\Gamma }(a)}{a^2}}\left[{\displaystyle \frac{1}{\mathrm{\Gamma }(a)}}{\displaystyle \frac{\mathrm{\Gamma }(a)}{a}}\right]^2}`$ (53)
$`s`$ $`=`$ $`{\displaystyle \frac{1}{b}}\left[\mathrm{log}(a){\displaystyle \frac{1}{\mathrm{\Gamma }(a)}}{\displaystyle \frac{\mathrm{\Gamma }(a)}{a}}\right].`$ (54)
Now re-writing (43) and (47), one finds

$`\mathrm{\Pi }(\theta )`$ $``$ $`\{\begin{array}{ccc}|\theta |\mathrm{exp}\left(\frac{\pi }{2}b\theta \right),\hfill & & \theta 0\hfill \\ & & \\ \mathrm{exp}\left[\frac{\pi }{2}\mathrm{e}^{b(\theta s)}+c\theta \right],\hfill & & \theta 0,\hfill \end{array}`$ (58)
with $`b=8\pi \sqrt{g_2/2}1.105`$, $`s=0.745`$ and $`c=b`$. These asymptotes differ only slightly from those for a generalised Gumbel function with $`a=\pi /2`$: firstly through the term $`|\theta |`$ for fluctuations below the mean and secondly through the term $`\mathrm{exp}(c\theta )`$ above the mean: the coefficient $`c=(\pi /2)b`$ for the modified Gumbel function, while $`c=b`$ for the true asymptote. These differences are enough to ensure that the modified Gumbel’s equation is not an exact solution to eqn. (28), however the comparison is so close that it is tempting try to get a good fit to (31) by solving for the constants $`a`$, $`b`$, $`s`$ and $`w`$.

Fourier transforming (49) gives

$`\mathrm{\Pi }_G(\theta )={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{2\pi }}{\displaystyle \frac{w}{b}}\mathrm{exp}\left(ix\theta isx+i{\displaystyle \frac{x}{b}}\mathrm{log}aa\mathrm{log}a\right)\mathrm{\Gamma }\left(ai{\displaystyle \frac{x}{b}}\right)={\displaystyle _{\mathrm{}}^{\mathrm{}}}{\displaystyle \frac{dx}{2\pi }}\mathrm{exp}\left[i\mathrm{\Phi }_G(x)\right].`$ (59)
We can compare $`\mathrm{\Phi }_G(x)`$ with $`\mathrm{\Phi }(x\sqrt{2/g_2})`$, assuming that the two Fourier transforms are nearly equal. The four constants should be calculated by minimising the difference between the two functions. To do this we can set the first four coefficients of the Taylor expansion of these functions equal. For $`\mathrm{\Phi }_G(x)`$ we have:

$`\mathrm{\Phi }_G(x)`$ $`=`$ $`ia\mathrm{log}ai\mathrm{log}(w/b)i\mathrm{log}\mathrm{\Gamma }(a)+\left[s\mathrm{\Psi }(a)/b+\mathrm{log}(a)/b\right]x+{\displaystyle \frac{i}{2b^2}}\mathrm{\Psi }^{}(a)x^2`$ (60)
$`+`$ $`{\displaystyle \frac{1}{6b^3}}\mathrm{\Psi }^{\prime \prime }(a)x^3{\displaystyle \frac{i}{24b^4}}\mathrm{\Psi }^{\prime \prime \prime }(a)x^4{\displaystyle \frac{1}{120b^5}}\mathrm{\Psi }^{(4)}(a)x^5+\mathrm{}`$ (61)
where $`\mathrm{\Psi }(z)`$ is the digamma function $`\mathrm{\Gamma }^{}(z)/\mathrm{\Gamma }(z)`$. For $`\mathrm{\Phi }`$ we have:

$`\mathrm{\Phi }(x\sqrt{2/g_2})={\displaystyle \frac{i}{2}}x^2{\displaystyle \frac{\sqrt{2}g_3}{3g_2^{3/2}}}x^3i{\displaystyle \frac{g_4}{2g_2^2}}x^4+{\displaystyle \frac{2\sqrt{2}g_5}{5g_2^{5/2}}}x^5+\mathrm{}`$ (62)
We therefore find that the four constants satisfy the relations:

$`{\displaystyle \frac{b}{w}}`$ $`=`$ $`{\displaystyle \frac{\mathrm{\Gamma }(a)}{a^a}},sb=\mathrm{log}a\mathrm{\Psi }(a),`$ (63)
$`b^2`$ $`=`$ $`\mathrm{\Psi }^{}(a),b^3g_3\left({\displaystyle \frac{2}{g_2}}\right)^{3/2}=\mathrm{\Psi }^{\prime \prime }(a)`$ (64)
The first three equations arise from the constraints of normalisation of the distribution, while the last expresses these constraints in terms of $`g_2`$ and $`g_3`$. The equations can be solved numerically. We find

$`a`$ $`=`$ $`1.5806801,b=0.9339355`$ (65)
$`s`$ $`=`$ $`0.3731792,w=2.1602858`$ (66)
The constants $`b`$ and $`s`$ calculated in this way are shifted slightly from the values extracted from the asymptotes, but $`a`$ is close to our very appealing first try $`\pi /2`$. Taking this value and calculating the constants $`b`$, $`s`$ and $`w`$ from normalisation one finds:

$`a`$ $`=`$ $`\pi /2,b=0.938`$ (67)
$`s`$ $`=`$ $`0.374,w=2.14,`$ (68)
in very satisfying agreement with the first method of calculation.

Given this solution, we can compute the coefficient ratio for the higher order terms in (61) and (62):

$`{\displaystyle \frac{1}{\mathrm{\Phi }_G^{(4)}(0)}}{\displaystyle \frac{^4\mathrm{\Phi }(x\sqrt{2/g_2})}{x^4}}|_{x=0}`$ $`=`$ $`{\displaystyle \frac{12g_4b^4}{g_2^2\mathrm{\Psi }^{\prime \prime \prime }(a)}}=0.9265029,`$ (69)
$`{\displaystyle \frac{1}{\mathrm{\Phi }_G^{(5)}(0)}}{\displaystyle \frac{^5\mathrm{\Phi }(x\sqrt{2/g_2})}{x^5}}|_{x=0}`$ $`=`$ $`{\displaystyle \frac{48\sqrt{2}g_5b^5}{g_2^{5/2}\mathrm{\Psi }^{(4)}(a)}}=0.8267429`$ (70)
The ratio of coefficients clearly does diverges from unity, but is does so slowly, indicating that the modified Gumbel function should be a good fit to the curve over the physical range. This is confirmed in Fig. 7 where we compare (49), using the values (65), with the exact result, from eqn. (31). On a natural scale the agreement is remarkably good over the entire range, with the only visible deviation coming around the maximum of the PDF, where the Gumbel curve is very slightly lower. On a logarithmic scale there is excellent general agreement over the whole of the plotted range, but a slight deviation can be observed for probabilities below $`10^3`$. For fluctuations below the mean the deviation is because the true asymptotic behaviour is quasi-exponential, $`x\mathrm{exp}(\alpha x)`$ and has a slight curvature, as discussed in the previous section. The results therefore confirm that, although the generalised Gumbel function is an excellent approximation for the PDF (31), it is not an exact solution.

¿From these results it is very tempting to take the generalised Gumbel function, with $`a`$ exactly $`\pi /2`$ as a working analytic expression for the PDF. However the connection with extremal statistics remains an open question . As discussed in section V, the spin wave Hamiltonian (12) is diagonalised in reciprocal space and the problem can be formulated in terms of a set of statistically independent variables. The PDF for extreme values of statistically independent variables can only follow three different asymptotic , or limit functions as the thermodynamic limit is taken. The only possible limit functions from extremal statistics of the Gumbel form discussed here are for $`a`$ integer, with $`a=1`$ for the biggest or smallest values.

Chapman et al. , have recently argued that the PDF for global quantities in any system with identifiable excitations on scales up to the system size should be dominated by extreme values. They shown that the PDF of extreme values among $`10^5`$ Gaussian random number generators approximates to a Gumbel function with $`a=\pi /2`$. This is not one of the predicted asymptotes , and we suggest that the deviation must be due to a very slow approach to the limit function with system size. It therefore does not seem to be a correct description of the $`2D`$-XY data as we do have a limit function which is well represented by eqn. (49) with $`a=\pi /2`$. However, if the results of  are relevant for non-equilibrium phenomena such as turbulence and self-organised criticality it would suggest the interesting property that corrections to the asymptotic forms, or limit functions,are a generic feature of these systems.

### B Generalised log-normal Distribution

The generalised log-normal distribution has the form

$$\mathrm{\Pi }_L(\theta )=\frac{w}{\sqrt{2\pi \sigma _L^2}(s\theta )}\mathrm{exp}\left\{\frac{1}{2\sigma _L^2}\left[\mathrm{ln}(s\theta )a\right]^2\right\},$$
(71)
with $`w=1`$. Following the same procedure as before, the generating function $`\mathrm{\Phi }_L(x)`$ can be developed as a power series

$$\mathrm{\Phi }_L(x)=x(\theta s+\mathrm{e}^{a+\sigma _L^2/2})+i\frac{x^2}{2}(\mathrm{e}^{2a+2\sigma _L^2}\mathrm{e}^{2a+\sigma _L^2})x^3(\frac{1}{6}\mathrm{e}^{3a+9\sigma _L^2/2}+\frac{1}{3}\mathrm{e}^{3a+3\sigma _L^2/2}\frac{1}{2}\mathrm{e}^{3a+5\sigma _L^2/2}).$$
(72)
Comparing (72) with (62) one finds the following expressions for $`s,a`$ and $`\sigma _L`$:

$`s`$ $`=`$ $`\mathrm{e}^{a+\sigma _L^2/2}`$ (73)
$`a`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}\left(\mathrm{e}^{2\sigma _L^2}\mathrm{e}^{\sigma _L^2}\right)`$ (74)
$`{\displaystyle \frac{\sqrt{2}}{3}}{\displaystyle \frac{g_3}{g_2^{3/2}}}`$ $`=`$ $`{\displaystyle \frac{1}{6}}\mathrm{e}^{3a}\left(\mathrm{e}^{9\sigma _L^2/2}+2\mathrm{e}^{3\sigma _L^2/2}3\mathrm{e}^{5\sigma _L^2/2}\right).`$ (75)
Eliminating $`a`$ and $`\sigma _L`$ leads to a cubic equation for $`s`$ in terms of $`\alpha =\frac{(g_2/2)^{3/2}}{g_3}=1/|\gamma |`$:

$$s^33\alpha s^2\alpha =0,$$
(76)
which could be solved exactly. We have solved it numerically, verifying that there exists one real and two complex roots. We find

$`s=3.45981`$, $`a=1.20109`$, $`\sigma _L=0.28325`$.

The function, with these parameters, is compared with the calculated PDF in Fig. 7. The general quality of fit is again excellent over the plotted range, with very small systematic deviations occurring in the wings of the distribution. It does not have the correct asymptotes; either exponential on the left or double exponential on the right, but as we have shown in the previous section, the true asymptotic behaviour is only reached outside the plotted regime, which explains why such a good fit can be achieved.

We have not, for the moment been able to develop a physical reasoning associated with the log-normal function and the origin, $`s=3.4\mathrm{}`$ although related to $`\gamma `$, seems rather arbitrary, but we do not exclude an explanation in terms of random multiplicative processes.

Note that log-normal distribution does appear in surface dynamics. Namely, starting with a flat interface as an initial condition, the short-time limit of the $`D=1`$ Edwards-Wilkinson dynamics yields a log-normal distribution for the interface width .

### C Generalised $`\chi ^2`$ Distribution

The $`\chi ^2`$ distribution for $`\nu `$ statistically independent degrees of freedom has the form

$$\mathrm{\Pi }_\chi (\theta )=w(s\theta )^{\nu /21}\mathrm{e}^{a(s\theta )},$$
(77)
with

$`w`$ $`=`$ $`{\displaystyle \frac{a^{\nu /2}}{\mathrm{\Gamma }(\nu /2)}}`$ (78)
$`\nu `$ $`=`$ $`2a^2.`$ (79)
As in the case of the Gumbel function, the generating function can be found in closed form:

$$\mathrm{\Phi }_\chi (x)=x(\theta s)+i\frac{\nu }{2}\mathrm{ln}(1ix/a),$$
(80)
whose development up to $`4^{th}`$ order in $`x`$ leads to

$$\mathrm{\Phi }_\chi (x)=x(\theta s)+\frac{\nu }{2a}x+i\frac{\nu }{4a^2}x^2\frac{\nu }{6a^3}x^3i\frac{\nu }{8a^4}x^4+O(x^5)+\mathrm{}$$
(81)
This series can again can be compared with (62) to give

$`s`$ $`=`$ $`{\displaystyle \frac{\nu }{2a}}`$ (82)
$`a`$ $`=`$ $`\sqrt{{\displaystyle \frac{\nu }{2}}}=s`$ (83)
$`\nu `$ $`=`$ $`{\displaystyle \frac{g_2^3}{g_3^2}}={\displaystyle \frac{8}{\gamma ^2}},`$ (84)
with numerical values

$`\nu =10.07155`$, $`a=2.24405`$ , $`s=a`$ , $`w=2.31233`$.

Comparing the function, shown in Fig. 7 with these parameters with the calculated curve, there is reasonably good agreement but this time deviation can be seen when plotted both on real and logarithmic scale. On the logarithmic scale the deviation is stronger than for the other fitting functions.

One can see that describing the correlated system as a finite number of degrees of freedom is a reasonably good approximation. It is an appealing concept and the calculation yields a system size independent number which depends uniquely on the skewness: $`\nu =g_2^3/g_3^2=8/\gamma ^2`$. If $`\gamma `$ developed towards zero, then $`\nu `$ would diverge and the $`\chi ^2`$ interpretation would be consistent with a Gaussian distribution. However, quantitatively it is not correct and the true description is a many body one . The difference between the two curves can be quantified by considering the ratio of the $`4^{th}`$ order terms:

$`\mathrm{\Phi }_\chi (x)^{(4)}`$ $`=`$ $`i{\displaystyle \frac{1}{2\nu }}`$ (85)
$`\mathrm{\Phi }(x)^{(4)}`$ $`=`$ $`i{\displaystyle \frac{g_4}{2g_2}},`$ (86)
so that $`\frac{\mathrm{\Phi }(x)^{(4)}}{\mathrm{\Phi }_\chi (x)^{(4)}}0.0238`$ which is very far from $`1`$.

### D Pearson’s Curve

Pearson  described an ingenious method of deriving a functional form for a PDF to fit experimental data, given the first four moments of the latter. He considered the differential equation

$$\frac{d\mathrm{ln}y}{dx}=\frac{x+b}{b_0+b_1x+bx^2}$$
(87)
and showed that if $`y`$ is a distribution then the parameters $`b`$, $`b_0`$, $`b_1`$ are specific functions of the first four moments. The expression can then be integrated to give (within a normalisation factor) an approximate functional form for the PDF, which by definition has the same principal moments as the data to be fitted. The success of Pearson’s approach relies on the observation that PDFs with the same moments are approximately coincident over the range of a few standard deviations, which is exactly the range of experimental interest. In the present case the mean is zero and the standard deviation is set to unity, so the shape of the curve depends only on the skewness, $`\gamma `$, and kurtosis, $`\kappa `$.

We find $`\gamma =g_3(2/g_2)^{3/2}=0.8907`$, and $`\kappa =3+3g_4(2/g_2)^2=4.415`$, which gives the following solution:

$$y=y_0\frac{(\beta \xi )^q}{(\alpha \xi )^p}$$
(88)
in which $`\xi =x0.39723`$, $`\beta =2.4787`$, $`\alpha =11.430`$, $`q=10.249`$, $`p=47.267`$, $`y_0=\mathrm{exp}(105.02)`$. Equating $`y(x)=\mathrm{\Pi }(\theta )`$, the fit to the exact expression (Fig. 8) is good between $`x=6`$ and $`x=2`$, but the very large numbers involved in eqn. (88) suggest that this functional form has no physical significance. From this analysis we can conclude that data collapse observed in refs.  should be interpreted as meaning that the third and fourth moments scale with $`\sigma `$ and $`L`$ in the same way as they do in the critical $`2D`$-XY model.

## IV Distribution in the $`D`$-dimensional Gaussian model

In this section, we study the asymptotics of the distribution function in general dimension $`D`$. It is straightforward to generalise the development from eqn. (28) to eqn. (31) for arbitrary dimension by redefining $`G(𝐪)`$ for dimension $`D`$ and summing over a $`D`$ dimensional Brillouin zone. The generalised expression (31) can then be numerically transformed to give $`\mathrm{\Pi }(\theta )`$. The results for $`D=1`$ and $`D=3`$ are shown in Figs. 10 and 11, where they are compared with data from Monte Carlo and molecular dynamics simulations. There is again excellent agreement showing that eqn. (31) is essentially exact in the low temperature regime, where the Hamiltonian (12) is valid. At higher temperatures the full Hamiltonian (10) generates vortex structures, the eqn. (18) is no longer valid and the derivation of eqn. (31) breaks down. Within the low temperature approximation there are three regimes: $`D<2`$, $`2<D<4`$, and $`D4`$, in addition to the special case $`D=2`$. The different regimes can be seen from a dimensional analysis of $`g_1`$ and $`g_2`$. As deviation from Gaussian behaviour is due to the abnormal influence of the integral scale in the form of infrared divergences, we approximate replacing $`G`$ by $`1/q^2`$ and re-calculate the $`g_k`$ by performing integrals over the Brillouin zone between $`2\pi /N^{1/D}`$ and $`2\pi `$. This procedure gives the correct qualitative behaviour, but there is a difference between the discrete sums and the integrals over the Brillouin zone, even in the thermodynamic limit (see App. C). The correct qualitative behaviour is

$`g_1`$ $``$ $`\{\begin{array}{ccc}C_{1,D}N^{(2D)/D},\hfill & D<2\hfill & \\ & & \\ A_1\mathrm{log}N+B_1,\hfill & D=2\hfill & \\ & & \\ C_{1,D},\hfill & D>2\hfill & \end{array}`$ (94)
and

$`g_2`$ $``$ $`\{\begin{array}{ccc}C_{2,D}N^{2(2D)/D},\hfill & D<4\hfill & \\ & & \\ \left(A_2\mathrm{log}N+B_2\right)/N,\hfill & D=4\hfill & \\ & & \\ C_{2,D}/N,\hfill & D>4.\hfill & \end{array}`$ (100)
The lower and upper critical dimensions, $`D=2`$ and $`D=4`$, are marked by the logarithmic divergence of $`g_1`$ and $`g_2`$ respectively.

Using the linearised order parameter (33) we find for $`D<2`$ that $`g_1`$ diverges as a power of $`N`$ giving $`m=1\tau C_{1,D}N^{(2D)/D}`$, which is a poor approximation for a thermodynamic quantity bounded on the interval $`[0,1]`$. Once outside this restricted low temperature region, $`\tau 1/\left[C_{1,D}N^{(2D)/D}\right]`$, both the linear approximation for the order parameter and the quadratic approximation for the Hamiltonian break down and there is a divergence in the behaviour of the PDF, as calculated from (31) and as simulated numerically. The system is, of course, disordered at all temperatures, so that the correct $`m`$ and $`\sigma `$ both vary as $`1/\sqrt{N}`$ and the PDF for the vector order parameter is a two dimensional Gaussian function centered on $`𝐦=0`$. The PDF for $`m`$, as defined in (11), is $`P(m)m\mathrm{exp}(m^2/2\sigma ^2)`$, analogous to a Maxwellian distribution of molecular speeds, and the thermodynamic system satisfies the central limit theorem (see App. A). As we have already seen, for $`D=2`$ the situation is different, as there is a large region of temperature where the quadratic Hamiltonian correctly describes the physics even though eqn. (33) is not a good approximation. In this regime of temperature, the PDF $`\mathrm{\Pi }(\theta )`$, for parameters (33) and (11) are however identical.

For dimension $`D>2`$, the low temperature expansion for the order parameter gives consistent results for all $`N`$, as long range order is stable and $`m1C_{1,D}\tau `$ is well defined. Above $`D=4`$, our results agree with mean field theory ($`D=\mathrm{}`$) where all sites are connected together. Here, $`m1\tau /4`$ and $`\sigma \tau /2\sqrt{2N}`$, and for large but finite $`N`$, the universal function $`\mathrm{\Pi }`$ is simply a Gaussian

$`\mathrm{\Pi }(\theta )={\displaystyle \frac{1}{\sqrt{2\pi }}}\mathrm{exp}\left({\displaystyle \frac{1}{2}}\theta ^2\right),`$ (101)
which corresponds to the central limit theorem for a collection of $`N`$ independent oscillators, each of expectation value $`m`$ and standard deviation $`\tau /2\sqrt{2N}`$.

### A Low Temperature Calculation in $`D=1`$

If the low temperature calculation for $`D<2`$ is not terribly pertinent for the thermodynamic system, it it highly relevant for the interface problem in the context of the EW model  and is exactly solvable in $`D=1`$ . In this case, computing the different $`g_k`$, we find

$`g_1=N/12,g_2=N^2/720,\mathrm{},g_p={\displaystyle \frac{2\zeta (2p)}{(2\pi )^{2p}}}N^pp1,`$ (102)
where $`\zeta (k)=_{i=1}^{\mathrm{}}i^k`$ is the Riemann zeta function . The expectation value of the magnetisation and standard deviation are

$`m`$ $`=`$ $`\mathrm{exp}\left({\displaystyle \frac{\tau N}{24}}\right)1{\displaystyle \frac{\tau N}{24}},`$ (103)
$`\sigma ^2=\left({\displaystyle \frac{1}{N}}{\displaystyle \underset{𝐫}{}}\mathrm{cosh}\tau G_R(𝐫)1\right)m^2`$ $`=`$ $`\left({\displaystyle _0^1}\mathrm{cosh}\tau N(x^2x+1/6)𝑑x1\right)m^2\tau ^2N^2,`$ (104)
which means that the ratio $`m/\sigma `$ scales as $`1/N`$, although for the parameters of the interface problem $`w^2/\sigma _{w^2}O(1)`$. We evaluate the universal distribution $`\mathrm{\Pi }`$, using the general eqn. (36) with $`G`$ defined for the $`D=1`$. After some algebra, we find for $`\mathrm{\Pi }(\theta )`$

$`\mathrm{\Pi }(\theta )`$ $`=`$ $`{\displaystyle \frac{dx}{2\pi }\mathrm{exp}\left[ix\left(\theta \frac{\sqrt{360}}{12}\right)\underset{k=1}{\overset{\mathrm{}}{}}\mathrm{log}\left(1\frac{ix\sqrt{360}}{2\pi ^2k^2}\right)\right]}`$ (105)
$`=`$ $`{\displaystyle \frac{dx}{2\pi }\mathrm{exp}\left[ix\left(\theta \frac{\sqrt{360}}{12}\right)\mathrm{log}\left(\frac{\mathrm{sin}\sqrt{ix\sqrt{360}/2}}{\sqrt{ix\sqrt{360}/2}}\right)\right]}`$ (106)
$`=`$ $`{\displaystyle \frac{dx}{2\pi }\mathrm{exp}\left[i\mathrm{\Phi }(x)\right]}`$ (107)
This expression is related directly to the function $`\stackrel{~}{\mathrm{\Phi }}`$ of eqn. (11) in ref. : $`\mathrm{\Pi }(\theta )=\stackrel{~}{\mathrm{\Phi }}(224\theta /\sqrt{360})`$. The method used in  is based on path integration, but the results are the same as our saddle point method, used to compute the asymptotics. Setting $`x=iy`$, the extrema of $`\mathrm{\Phi }`$ satisfy the equation

$`\theta {\displaystyle \frac{\sqrt{360}}{12}}`$ $`=`$ $`{\displaystyle \frac{\sqrt{360}}{2\pi ^2}}{\displaystyle \underset{k=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{k^2+y\sqrt{360}/2\pi ^2}}`$ (108)
$`=`$ $`{\displaystyle \frac{\sqrt{360}}{2\pi ^2}}\left({\displaystyle \frac{1}{2y\sqrt{360}/2\pi ^2}}+{\displaystyle \frac{\pi }{2\sqrt{y\sqrt{360}/2\pi ^2}}}\mathrm{coth}\pi \sqrt{y\sqrt{360}/2\pi ^2}\right).`$ (109)
For $`\theta 1`$, $`y`$ is close to the first pole $`2\pi ^2/\sqrt{360}`$ of the right hand side of 108), which is similar to the $`2D`$ case (41) except that the $`1D`$ extrema function is easier to evaluate. Performing the saddle point computation, we find that $`\mathrm{\Pi }`$ behaves asymptotically as

$`\mathrm{\Pi }(\theta )\mathrm{exp}\left[2\pi ^2\theta /\sqrt{360}\right]`$ (110)
which is the same as . The asymptote, (110), is drawn on Fig. (10 b), where it can be compared with the full calculation and with simulation. The exponential tail is extremely well defined and the predicted slope is clearly correct.

In the regime of fluctuations above the mean, for $`y1`$, $`\theta `$ is close to the constant $`\sqrt{360}/12`$, and no extrema exist for $`\theta `$ beyond this value. In this case, $`y\sqrt{360}/[8(\sqrt{360}/12\theta )^2]`$, and the saddle point approximation leads to the following asymptotic value for $`\mathrm{\Pi }`$ near this upper limit

$`\mathrm{\Pi }\left(\theta \right)\left(\sqrt{360}/12\theta \right)^{5/2}\mathrm{exp}\left({\displaystyle \frac{1}{8}}\sqrt{360}/(\sqrt{360}/12\theta )\right),`$ (111)
which is the same result as . We refer the reader to ref.s  for the precise coefficients in both asymptotic limits.

In conclusion, we find that for $`\theta 1`$, the universal distribution again has an exponential tail, while for large fluctuations above the mean the PDF shoots to zero as $`\mathrm{exp}[3\theta _0/2(\theta \theta _0)]`$, with $`\theta _0=\sqrt{360}/12`$. This upper limit corresponds to the constraint that $`m1`$.

It is worth pointing out in some detail here that the exponential tail in the one-dimensional problem is not the result of critical fluctuations. The small deviations in angle $`(\theta _i\theta _j)`$ constitute a random walk with $`w\sqrt{1m}`$ being the radius of gyration, which scales correctly as the square root of the walk length, $`L`$. The $`1D`$ linear order parameter, or interface problem is therefore nothing more than a simple random walk , but despite this the PDF is, as shown in Fig. 10: a standard result for such a walk is that the mean radius of gyration is proportional to the mean end to end distance, $`S`$ of the walk. It is easily shown that the PDF $`P(S)`$ is Gaussian . Changing variable from $`S`$ to $`X=S^2`$ one finds $`P(X)X^{1/2}\mathrm{exp}(X/X_0)`$; a trivial distribution with an exponential tail. The PDF for $`w^2`$ has the same exponential tail, but does not show the essential singularity at $`w^2=0`$ ($`m=1`$) and we conclude the rather surprising property of a random walk, that the PDF for the radius of gyration and for the end to end distance are not the same. The origin of this difference is that the average angle $`\overline{\theta }`$, corresponding to the center of mass of an equivalent random walk, fluctuates with $`L`$ in the same way as the radius of gyration itself and this lack of self averaging removes the essential singularity from the PDF at $`w^2=0`$.

### B Asymptotic Solutions in General Dimension

We first evaluate the asymptotic value of $`\mathrm{\Pi }`$ for positive $`\theta `$ by solving the saddle point of (36), rescaling the variable $`x\sqrt{g_2/2}x`$ for convenience. For $`D<2`$, the ratio $`g_1/\sqrt{g_2}`$ is independent of the system size and, with $`x=iy`$, the equation to solve is

$`\theta {\displaystyle \frac{g_1}{\sqrt{2g_2}}}{\displaystyle _{\mathrm{Cst}/N^{1/D}}^{\mathrm{Cst}}}{\displaystyle \frac{Nq^{D1}}{q^2N+y\sqrt{2/g_2}}}𝑑q`$ (112)
where Cst is a constant. By setting $`N^{1/D}q/\sqrt{y}q`$, we find that, for large and positive $`y`$

$`{\displaystyle \frac{g_1}{\sqrt{2g_2}}}\theta y^{(D2)/2}{\displaystyle _{\mathrm{Cst}/\sqrt{y}}^{\mathrm{Cst}N^{1/D}/\sqrt{y}}}{\displaystyle \frac{q^{D1}}{1+q^2}}𝑑qy^{(D2)/2}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{q^{D1}}{1+q^2}}𝑑q`$ (113)
which means that $`\theta `$ is close to the upper bound $`g_1/\sqrt{2g_2}`$. Replacing the asymptotic value of $`y`$ for the extrema in the function $`\mathrm{\Phi }`$ (36), we find that

$`\mathrm{log}\mathrm{\Pi }(\theta )\mathrm{Cst}\left({\displaystyle \frac{g_1}{\sqrt{2g_2}}}\theta \right)^{D/(D2)}+\mathrm{Logarithm}\mathrm{corrections},\theta {\displaystyle \frac{g_1}{\sqrt{2g_2}}},D<2.`$ (114)
The logarithmic corrections come partly from the Gaussian integration around the saddle point and partly from other terms in (113) which are not accurately evaluated within our approximation. Note again that $`D=2`$ is a special case as, instead of (113) we have a logarithmic divergence (see eqn. (44)) and subsequently a double exponential fall in $`\mathrm{\Pi }`$ for large $`\theta `$. For the interval $`2<D<4`$ the ratio $`g_1/\sqrt{g_2}`$ and the integral (113) are no longer finite and so we look to eqn. (41) for the asymptotic behaviour:

$`\theta {\displaystyle \frac{1}{g_2}}{\displaystyle _{\mathrm{Cst}/N^{1/D}}^{\mathrm{Cst}}}{\displaystyle \frac{d^Dq}{Nq^4}}{\displaystyle \frac{y}{1+y\sqrt{2}/(\sqrt{g_2}Nq^2)}}.`$ (115)
By again setting $`N^{1/D}q/\sqrt{y}q`$ and using the fact that $`g_2N^{2(D2)/D}`$ is finite (94), we arrive at

$`\theta y^{(D2)/2}{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{q^{D3}dq}{1+q^2}},y1`$ (116)
The integral is convergent for $`2<D<4`$ and by replacing the value for $`y`$ in the saddle point approximation, we get the asymptotic form for $`\mathrm{\Pi }`$, in the limit of large and positive $`\theta `$:

$`\mathrm{log}\mathrm{\Pi }(\theta )\mathrm{Cst}\theta ^{D/(D2)}+\mathrm{Logarithm}\mathrm{corrections},\theta 1,\mathrm{\hspace{0.33em}\hspace{0.33em}2}<D<4.`$ (117)
In 3 dimensions, we therefore expect that the logarithm of the distribution falls off like $`\theta ^3`$, well above the mean. We have not tested this in detail, but the PDF does fall off more slowly for $`D=3`$ than $`D=2`$, in qualitative agreement with the predictions here. Finally, we note that throughout the range $`2<D<4`$ the universal PDF is non-Gaussian, but the hyperscaling relation is invalid: $`m/\sigma g_1/\sqrt{g_2}N^{(D2)/D}`$.

For $`D>4`$, $`g_2`$ decreases as $`1/N`$, consequently, (116) has to be modified. We find, instead of (116), that

$`\theta y^{(D2)/2}N^{(4D)/4}{\displaystyle _{N^{(D4)/4D}/\sqrt{y}}^{N^{1/4}/\sqrt{y}}}{\displaystyle \frac{q^{D3}dq}{1+q^2}}y,N1.`$ (118)
We can, in fact, replace the integrand inside the integral by $`q^{D5}dq`$ since the integration domain is large, from which we find that the saddle point is proportional to $`\theta 1`$ and deduce that $`\mathrm{\Pi }`$ is Gaussian on the right hand side of the curve. The same is true for $`D=4`$ despite the logarithmic divergence of $`g_2`$.

In the opposite limit $`\theta 1`$, for both $`D=1`$ and $`D=2`$ the asymptotic value of the distribution falls down exponentially (43,110). We would now like to evaluate this limit in general dimensions. In both cases the coefficient of $`\theta `$ is related to the value of $`g_2`$, i.e. $`C_{2,D}`$. Rewriting the eqn. (41) with discrete sums (see also App. C), we have

$`\theta ={\displaystyle \frac{N^{2(2D)/D}}{16\pi ^4g_2}}{\displaystyle \underset{m_i0}{}}{}_{}{}^{^{}}{\displaystyle \frac{1}{(m_1^2+\mathrm{}+m_D^2)}}{\displaystyle \frac{y}{(m_1^2+\mathrm{}+m_D^2)+y\sqrt{2/g_2}N^{(2D)/D}/4\pi ^2}}`$ (119)
where the sum excludes $`m_i=0`$, $`i=1,\mathrm{},D`$. The saddle point equation has a solution $`y`$ which is the pole nearest the origin, $`y=4\pi ^2\sqrt{g_2/2}N^{(D2)/D}`$, i.e. for sets of $`\{m_i\}`$ with one element equal to 1, the others being zero. For $`D<4`$ and large $`N`$, this pole is finite since $`g_2`$ compensates $`N^{2(D2)/D}`$, so that its value is simply $`y=4\pi ^2\sqrt{C_{2,D}/2}`$. Applying the saddle point integration, we find that the dominant term in the logarithm of $`\mathrm{\Pi }`$ is, below the mean

$`\mathrm{log}\mathrm{\Pi }(\theta )4\pi ^2\sqrt{{\displaystyle \frac{C_{2,D}}{2}}}\theta ,\theta 1,D<4,`$ (120)
and is linear in $`\theta `$ for every dimension below 4. Included in Fig. 11 for $`D=3`$ is a fit, on the left hand side of the form (120), with $`C_{2,3}`$ calculated numerically. There is again excellent agreement, which convincingly confirms the presence of the exponential tail. In fact, true exponential behaviour is reached for smaller values of $`\theta `$ than for $`D=2`$.

For $`D>4`$, the value of this pole diverges like $`N^{(D4)/2D}`$, and the previous solution fails. In fact, the solution (118) for positive $`\theta `$ and $`y`$ is also valid for negative values, if $`q^{D1}dq/(q^2+1)`$ is replaced by $`q^{D1}dq/(q^21)`$. Since the integration domain is far from the pole of the denominator, we can approximate the integrand in both cases by $`q^{D5}dq`$, and we get the same result as (118). We therefore, finally conclude that $`\mathrm{\Pi }`$ is also Gaussian on the left hand side of the curve and the central limit theorem applies for $`D>4`$.

## V Conclusion

Probability functions with exponential, rather than Gaussian behaviour are a common feature of complex systems . For example, the PDF for velocity differences at microscopic scales, in fully developed turbulence show exponential tails . This appears to be true in turbulence, not only for microscopic quantities, but also for global quantities; the energy injected into a closed turbulent flow being a very well controlled and documented example . Following these observations we have proposed that this is also a generic feature of complex systems . In this paper we have shown that, for the low temperature phase of the XY model, a critical system at equilibrium, analogous behaviour occurs when a few long wavelength and large amplitude modes make their presence felt in the global measure, which is typically a sum over $`O(N)`$ degrees of freedom. The exponential tail can occur in three physically different situations. The first is in two-dimensions, when the system is critical and fluctuations occur over all length scales. The second is in one-dimension, when the system is not critical, but an exponential tail occurs for a particular global measure, relevant to problems of interface growth, whose moments are completely dominated by the integral scale. The third is in three-dimensions, also non-critical, where despite stable long range order, the large amplitude long wavelength modes continue to make their presence felt. The detailed form of the PDF in these three cases are quite different and easily discernible in experiment. In table 1 we show the evolution of the skewness and the kurtosis with spatial dimension. The deviations from the Gaussian limit are largest in one dimension, and decrease continually to zero at $`D=4`$. We propose that the difference in the form of the PDF could be used as an experimental signature of the underlying physics.

From the general evolution shown in table 1, one might expect a dependence on shape, with dimensional crossover as the length scale in one direction changes from microscopic to macroscopic. This is indeed the case and for example, in two-dimensions, the skewness and kurtosis of the PDF calculated from eqn. (31) increases towards the values for $`D=1`$ if the ratio of lengths in the $`x`$ and $`y`$ directions, $`L_x`$ and $`L_y`$, are varied continuously from unity. It would be extremely interesting to establish if the same is true when the length scales are varied in turbulence experiments and numerically, in the models of self-organised criticality.

To see how the anisotropy of the PDF comes from the long wavelength excitations, we give an analysis in reciprocal space: the Hamiltonian (12) is diagonalised

$`H={\displaystyle \frac{J}{2}}{\displaystyle \underset{𝐪>0}{}}G(𝐪)^1\mathrm{}\{\varphi _𝐪\}^2,`$ (121)
where $`\varphi _𝐪`$ is the discrete Fourier transform of $`\theta _i`$ and the sum is over the Brillouin zone , with the thermodynamic variable for each $`𝐪`$ taken as the real part of $`\varphi _𝐪`$. Defining $`m_𝐪=(1/2N)\mathrm{}\{\varphi _𝐪\}^2`$ the linear order parameter can be written $`m=1_𝐪m_𝐪`$, where the $`m_𝐪`$ are statistically independent variables with PDF

$`P(m_𝐪)`$ $`=`$ $`\sqrt{{\displaystyle \frac{\beta Jq^2N}{4\pi }}}m_𝐪^{1/2}\mathrm{e}^{\beta JNq^2m_𝐪}.`$ (122)
Here, as we are principally interested in the modes at small $`q=|𝐪|`$ we have, without loss of generality, approximated $`G(𝐪)^1q^2`$. The PDF for $`m`$ is thus nothing more than the composite PDF for a set of independent spin wave modes or an “ideal gas” of particles, whose only peculiarity is that the mass term varies as $`q^2`$. The Goldstone modes have wave vector $`q=2\pi /L`$ and hence make contributions of $`O(1)`$ to $`m`$, while the modes on the zone edge, with $`q=\pi `$ have only microscopic amplitude. This dispersion in amplitudes is the key to the unusual behaviour for $`D=2`$, as it violates one of the conditions for the central limit theorem to apply to a sum of statistically independent variables: that the individual amplitudes do not differ by too much. However it is not true that the Goldstone modes, by themselves give the complete PDF. The mean value $`_qm_q_{2\pi /L}^\pi q^2n(q)𝑑q`$ where $`n(q)q^{D1}`$ is the density of states. For $`D=2`$ both limits of the integral are required and a detailed calculation gives $`_qm_q=(\eta /4)\mathrm{log}(CN)`$, with $`C=1.87`$ and with critical exponent $`\eta =T/2\pi J`$. The anomalous term $`\mathrm{log}N`$ therefore reflects the fact that modes from all over the Brillouin zone are relevant for $`m`$ and through eqn. (23), for the higher moments $`m^p`$.

For $`D=1`$ only the lower limit of integration is required, the upper limit can be set to $`\mathrm{}`$ and the constants $`g_p`$ are proportional to $`N^p`$. As a result the linear development of the order parameter in small angles, (33), is a very poor approximation for the thermodynamic quantity defined on the interval $`[0,1]`$. The two expressions, (11) and (33), describe different physical quantities. The former is directly related to the interface width in the Edwards-Wilkinson model of interface growth. The PDF for the full order parameter is consistent with an uncorrelated system, that is, a paramagnet with two-dimensional order parameter. For the linear order parameter the PDF, shown in Fig. 10, does have an exponential tail, but this is not the result of critical fluctuations, it is the property of a simple random walk. We remark further that dependence on a macroscopic length scale does not, in itself, indicate critical behaviour. Rather, critical behaviour is exemplified by the case $`D=2`$, where all length scale are important between the microscopic and macroscopic cut off.

$`D=3`$ represents the opposite of the one dimensional case: $`m`$ is controlled by the upper limit of integration and the result is unchanged by setting the lower limit to zero. However, despite long range order being stable and the system not being critical in the low temperature phase, the exponential tail persists. This is related to temperature being a dangerously irrelevant variable  near the zero temperature fixed point of a renormalisation group flow, between the lower and the upper critical dimension. The constant $`g_2`$ now falls to zero with system size but it does so more slowly than $`1/N`$ (see eqn. (100)). As a result of this slow decay, the ratio $`g_p/g_2^{p/2},p>2`$ in eqn. (28) is independent of $`N`$ and the distribution is non-Gaussian, despite $`g_p`$ and $`g_2`$ both being zero in the thermodynamic limit. At low temperature the magnetisation is finite, but the Goldstone mode influences the PDF sufficiently to produce an exponential tail. A physical consequence of this anomaly is that the longitudinal susceptibility

$`\chi {\displaystyle \frac{N}{T}}\left(m^2m^2\right)N^{(4D)/D}`$ (123)
is weakly divergent throughout the ordered phase . This is true for all magnetic systems with Heisenberg or XY symmetry. It could therefore be interesting to look for evidence of the departure from Gaussian behaviour experimentally in a non-critical three-dimensional system. Precision temperature control would not be required, however, as the ratio $`\sigma /m`$ falls off as $`1/N^{1/3}`$, the divergence in the susceptibility is very weak and this phenomenon may be out of experimental reach.

| $`D`$ | $`\gamma `$ | $`\kappa `$ |
| --- | --- | --- |
| $`\mathrm{\hspace{0.33em}\hspace{0.33em}\hspace{0.33em}1}`$ | $`1.807`$ | $`\mathrm{\hspace{0.33em}\hspace{0.33em}8.14}`$ |
| 2 | $`0.891`$ | $`4.41`$ |
| 3 | $`0.354`$ | $`3.31`$ |
| 4 | $`0`$ | $`3.0`$ |

TABLE 1 Variation of skewness $`\gamma `$ and kurtosis $`\kappa `$ with dimension $`D`$

Returning finally to critical systems; we have been able to exploit a system interacting via a quadratic Hamiltonian at exactly the lower critical dimension. In this particular case one has access to a critical point, with the fluctuation dominated behaviour that this implies, while retaining the benefit of Gaussian integration over phase space. As a result, all critical behaviour can be calculated microscopically, without the need for either the renormalisation group or the scaling hypothesis. The only price one pays for this simplicity is a critical system with a single independent exponent and the scaling relations satisfied through weak scaling only. In general, we believe that the analytic results that we have obtained are useful for the understanding of finite-size scaling and for the interpretation of experimental observations from more complex correlated systems. The examples we have discussed  point towards a behaviour analogous to criticality for an enclosed turbulent flow and for models showing self-organised criticality. However the detailed analysis presented here leaves many open questions and more experiment and simulation are clearly required if the generality and the limits of this proposition are to be tested further.

###### Acknowledgements.

This work was largely motivated by our collaboration with K. Christensen, H.J. Jensen, S. Lise, J. López and M. Nicodemi from Imperial College London and we are particularly grateful to H.J. Jensen and M. Nicodemi for bringing the theory of extremal statistics to our attention. In addition we have greatly benefitted from discussions with S. Fauve, N. Goldenfeld, J. Harte, A. Noullez, and Z. Rácz during the SIMU/CECAM planning meeting “Universal Statistics in Correlated Systems”, in Lyon, 29-31 March 2000 and from subsequent discussions with L. Berthier, E. Leveque, S. McNamara, P. Pujol and again Z. Rácz in connection with the $`1D`$ problem. It is our pleasure to thank all these people. MS acknowlegdes a Marie Curie fellowship of the European Commission (contract ERBFMBICT983561). The numerical work was supported by the Pôle Scientifique de Modélisation Numérique at the École Normale Supérieure de Lyon.

## A

### Some comments on the central limit theorem in critical systems

The central limit theorem is a powerful result of probability theory that provides the foundation for statistical thermodynamics . It states that the PDF of the sum $`Z=_{i=1}^Nz_i`$ of $`N`$ statistically independent variates $`z_i`$ tends, in the limit of large $`N`$ and for moderate values of the variate $`Z`$, to a Gaussian distribution. As well as the statistical independence of the $`z_i`$, another key criterion for the theorem to hold is that the $`z_i`$ are individually negligible . At a critical point, the first of these criteria is violated. The $`2D`$ XY model is of particular interest here as it is diagonalisable into statistically independent degrees of freedom and maps directly onto a problem where the second criterion is violated: the direct space variables, that is, the spins $`S_i`$, are certainly individually negligible for large system size $`N`$, but are strongly correlated. On the other hand, when diagonalised in reciprocal space the spin wave variables are statistically independent, but are no longer all individually negligible. In particular, the long wavelength modes make a significant impact on the fluctuations of the global measure; in this case the linearised order parameter (33). The PDF for the full and the linear order parameters are identical, even when the quantities themselves differ, which makes it an ideal system for the practical study of the breakdown of the central limit theorem. A conventional critical system cannot, in general be diagonalised in this way, as evidenced by the divergent specific heat.

Strictly speaking, the central limit theorem does not apply to the compound variate $`Z`$, but rather to the normalised quantity $`(ZZ)/N^{1/2}`$. This normalisation is essential for a reasonable PDF in the thermodynamic limit, as the standard deviation for fluctuations about the mean value $`Z`$ scales with system size in the same way. If a normalisation factor $`N^{1/2+\rho }`$, $`\rho 0`$, is chosen then one obtains a distribution that is concentrated either at zero or infinity . We illustrate this with an example from statistical thermodynamics. The total energy $`E`$ of an ideal gas of $`N`$ molecules has a PDF of the form $`P(E)E^{3N/21}\mathrm{exp}(\beta E)`$. It is straightforward to confirm that $`P(E)`$ tends to a delta function in the thermodynamic limit, while $`P(E/N^{1/2})`$ tends to a Gaussian function . One can see from this example that the function is never truly Gaussian - indeed it is always of the form $`\mathrm{ln}P(3N/21)\mathrm{ln}E\beta E`$, which can easily be made independent of $`N`$ by choosing appropriate units. The central limit theorem applies because the width of the distribution scales as $`N^{1/2}`$ which means that fluctuations with any physical significance are all concentrated near the turning point of the function $`\mathrm{ln}P`$. The theorem only has meaning because of the significance one attaches to values of the variate that differ by only a few standard deviations from the mean. In practical terms it is therefore essential to normalise fluctuations to the standard deviation in order to test the central limit theorem.

In the case of dependent variables, the limit distribution can be different from the Gaussian form. Two types of dependent random variable can be defined : (i) weakly dependent, in which the correlation function falls to a constant value in a finite range, and the standard deviation again varies as $`\sqrt{N}`$; (ii) strongly dependent, in which the fluctuations vary as a power of $`N`$ different from one half. Case (i) corresponds to a system with a finite correlation length. In case (ii), which includes systems with critical fluctuations, the central limit theorem does not hold, but a reasonable PDF can be obtained by normalising to the variance, hence to an appropriate power of $`N`$, with $`\rho 0`$. Defining the (scalar) order parameter to be the intensive quantity $`z=Z/N`$, and using the scaling relations for a finite system, one finds $`\rho =(1\eta /2)/D`$. The limit distribution is now expected to be non-Gaussian, as can be shown explicitly for the Ising model . Note however that, $`\rho `$ remains non-zero even at the upper critical dimension (taken as $`D=4`$ here), when $`\eta =0`$ and where one might legitimately expect a Gaussian PDF. The condition $`\rho 0`$ may therefore be a necessary but not a sufficient condition to ensure non-Gaussian order parameter fluctuations.

Case (ii) is not actually limited to critical fluctuations: the example of a dangerously irrelevant variable discussed in the text also falls into this category, with $`\rho =2/D1/2`$. Here, $`\rho `$ does go to zero as the upper critical dimension is reached and the danger of the irrelevant temperature variable disappears.

An ordinary critical point is more complicated than those of the $`2D`$-XY model. In this case the correlation length is only infinite precisely at the critical temperature. A non-Gaussian limit function can therefore only be found on a locus of points such that $`\xi /L`$ is a constant as the thermodynamic limit is taken. Thus, fixing the temperature $`TT_C`$ and varying $`N`$ will always cause a transition from non-Gaussian to Gaussian statistics. Conversely, fixing $`T=T_C`$ one will only arrive at the stable limit function in the thermodynamic limit. One can therefore imagine a set of loci of constant PDF in $`[T,L^1]`$ space that converge on $`[T_C,0]`$. We have suggested  that there is one such locus, $`[T^{}(L),L^1]`$, where the PDF has approximately the same form as that of the $`2D`$-XY model. Thus, to sit at the critical temperature and change $`L`$ is not the same as traveling along the locus $`[T^{}(L),L^1]`$. From scaling argument  one can check that the tails of the PDF at $`T_C`$ should have the form $`P(m)\mathrm{exp}(m^{\delta +1})`$ in order to yield the correct scaling relation in the presence of a weak magnetic field: $`mh^{1/\delta }`$. We do not find this, despite the same scaling relation holding for the $`2D`$-XY model with $`\delta =8\pi J/k_\mathrm{B}T1`$. This difference may come from the difference in trajectories in the space of variable $`T`$ and $`L`$.

A final point concerns the central limit theorem as applied to a vector order parameter, $`𝐦`$, such as the XY model. In the high temperature limit, the fluctuations in the vector $`𝐦`$ follow a two-dimension Gaussian centered on $`𝐦=0`$ and the PDF for the scalar $`m=|𝐦|`$ follows a “Maxwell speed distribution” for a two-dimensional gas. In an ordered regime and even in the critical regime for $`D=2`$ , $`\sigma m`$ which means $`m`$ behaves, to an excellent approximation, as a one-dimensional quantity. The symmetry breaking therefore induces a change in topology for the fluctuations in $`𝐦`$. This is generalisable to order parameters of higher dimension.

## B

The graphs $`g_k`$ can be written, in the large $`N`$ limit in terms of power series. For example:

$`g_2=\underset{N\mathrm{}}{lim}{\displaystyle \frac{4}{N^2}}{\displaystyle \underset{m=1}{\overset{Q}{}}}{\displaystyle \underset{n=1}{\overset{Q}{}}}{\displaystyle \frac{1}{(42\mathrm{cos}2\pi m/\sqrt{N}2\mathrm{cos}2\pi n/\sqrt{N})^2}}+{\displaystyle \frac{4}{N^2}}{\displaystyle \underset{m=1}{\overset{Q}{}}}{\displaystyle \frac{1}{(42\mathrm{cos}2\pi m/\sqrt{N})^2}},`$ (B1)
where $`Q=(\sqrt{N}1)/2`$. The sum is dominated by the contributions for small $`m`$ and $`n`$, but as the pole $`m=0,n=0`$ is explicitly excluded from the sum, it remains finite even in the limit $`N\mathrm{}`$. Taking only the first terms in a development of the cosines, which is exact in the thermodynamic limit one finds

$`g_2`$ $`=`$ $`{\displaystyle \frac{1}{4\pi ^2}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{m^4}}+{\displaystyle \frac{1}{4\pi ^4}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(m^2+n^2)^2}}`$ (B2)
$`=`$ $`{\displaystyle \frac{1}{360}}+{\displaystyle \frac{1}{4\pi ^4}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(m^2+n^2)^2}},`$ (B3)
and in general, for $`g_k`$

$`g_k={\displaystyle \frac{1}{4\pi ^2}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{m^{2k}}}+{\displaystyle \frac{1}{4\pi ^4}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{(m^2+n^2)^k}}.`$ (B4)
## C

For large and positive $`y`$ the functional form of $`\phi `$ is $`\phi (y)\frac{1}{8\pi }\mathrm{log}y+constant`$. To evaluate it in detail we use the results of the App. B to write:

$`\phi (y)=\underset{Q\mathrm{}}{lim}{\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{Q}{}}}\left({\displaystyle \frac{1}{m^2}}{\displaystyle \frac{1}{m^2+\widehat{y}}}\right)+{\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{Q}{}}}{\displaystyle \underset{n=1}{\overset{Q}{}}}\left({\displaystyle \frac{1}{m^2+n^2}}{\displaystyle \frac{1}{m^2+n^2+\widehat{y}}}\right),`$ (C1)
where $`\widehat{y}=y/4\pi ^2`$ The first two summations give, in the limit of large $`Q`$, a constant and a function of $`\widehat{y}`$ which tends to zero for large argument:

$`\underset{y\mathrm{}}{lim}{\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1}{m^2}}{\displaystyle \frac{1}{m^2+\widehat{y}}}\right)={\displaystyle \frac{1}{12}}`$ (C2)
The double sum can be rewritten as

$`{\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{Q}{}}}{\displaystyle \underset{n=1}{\overset{Q}{}}}\left({\displaystyle \frac{1}{m^2+n^2}}{\displaystyle \frac{1}{m^2+n^2+\widehat{y}}}\right)={\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{Q}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1}{m^2+n^2}}{\displaystyle \frac{1}{m^2+n^2+\widehat{y}}}\right)R(Q,y)`$ (C3)
where $`R`$ is a correction term which vanishes in the limit of large $`Q`$:

$`R(Q,y)={\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{Q}{}}}{\displaystyle \underset{n=Q+1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1}{m^2+n^2}}{\displaystyle \frac{1}{m^2+n^2+\widehat{y}}}\right)`$ (C4)
The sum can be evaluated in the limit $`Q\mathrm{}`$ using the Abel-Plana formula :

$`{\displaystyle \underset{i=p}{\overset{q}{}}}f(i)={\displaystyle _p^q}f(x)𝑑x+{\displaystyle \frac{1}{2}}f(p)+{\displaystyle \frac{1}{2}}f(q)+2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\mathrm{}[f(q+ix)f(p+ix)]}{\mathrm{exp}(2\pi x)1}}𝑑x,`$ (C5)
where $`f`$ is any real function that satisfied the assumptions in . Applying this to $`R(Q,y)`$, we have

$`R(Q1,y)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{Q1}{}}}{\displaystyle \frac{1}{m}}\left[\pi /2\mathrm{arctan}(Q/m)\right]+{\displaystyle \frac{1}{2(m^2+Q^2)}}+4Q{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{xdx}{(x^2m^2Q^2)^2+4Q^2x^2}}{\displaystyle \frac{1}{\mathrm{exp}(2\pi x)1}}`$ (C6)
$``$ $`{\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{Q1}{}}}{\displaystyle \frac{1}{\sqrt{m^2+\widehat{y}}}}\left[\pi /2\mathrm{arctan}\left(Q/\sqrt{m^2+\widehat{y}}\right)\right]+{\displaystyle \frac{1}{2(m^2+\widehat{y}+Q^2)}}`$ (C7)
$`+`$ $`4Q{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{xdx}{(x^2m^2\widehat{y}Q^2)^2+4Q^2x^2}}{\displaystyle \frac{1}{\mathrm{exp}(2\pi x)1}}`$ (C8)
The first term tends, in the large $`Q`$ limit, to the integral

$`{\displaystyle \underset{m=1}{\overset{Q1}{}}}{\displaystyle \frac{1}{m}}\left[\pi /2\mathrm{arctan}(Q/m)\right]{\displaystyle _0^1}{\displaystyle \frac{dx}{x}}\left[\pi /2\mathrm{arctan}(1/x)\right]={\displaystyle _0^1}{\displaystyle \frac{\mathrm{log}xdx}{1+x^2}}=\mathrm{Catalan}.`$ (C9)
A similar behaviour is found for the fourth term, since in this limit the dependence on $`\widehat{y}`$ of this term vanishes as $`\widehat{y}/Q^2`$. The other terms are corrections proportional to the inverse of some power of $`Q`$, so that $`R`$ vanishes in the large $`Q`$ limit. The double sum (C3) can thus be reduced to a simple sum, since

$`{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{n^2+z^2}}={\displaystyle \frac{1}{2z^2}}+{\displaystyle \frac{\pi }{2z}}\mathrm{coth}\pi z.`$ (C10)
We therefore have, for large $`y`$

$`{\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{Q}{}}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1}{m^2+n^2}}{\displaystyle \frac{1}{m^2+n^2+\widehat{y}}}\right)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi ^2}}{\displaystyle \underset{m=1}{\overset{Q}{}}}{\displaystyle \frac{1}{2m^2}}+{\displaystyle \frac{\pi }{2m}}\mathrm{coth}\pi m+{\displaystyle \frac{1}{2(m^2+\widehat{y})}}`$ (C11)
$``$ $`{\displaystyle \frac{\pi }{2\sqrt{m^2+\widehat{y}}}}\mathrm{coth}\pi \sqrt{m^2+\widehat{y}}`$ (C12)
$``$ $`{\displaystyle \frac{1}{24}}+{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{4\pi \sqrt{m^2+\widehat{y}}}}\left(1\mathrm{coth}\pi \sqrt{m^2+\widehat{y}}\right)+{\displaystyle \frac{1}{4\pi m}}(\mathrm{coth}\pi m1)`$ (C13)
$`+`$ $`{\displaystyle \frac{1}{4\pi }}\left({\displaystyle \frac{1}{m}}{\displaystyle \frac{1}{\sqrt{m^2+\widehat{y}}}}\right).`$ (C14)
The series containing the hyperbolic function of $`y`$ vanishes in the limit of large $`y`$ and the asymptotic behaviour of the last term can be evaluated with the Abel-Plana formula (C5)

$`{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}\left({\displaystyle \frac{1}{m}}{\displaystyle \frac{1}{\sqrt{m^2+\widehat{y}}}}\right)`$ $``$ $`{\displaystyle _1^{\mathrm{}}}𝑑x\left({\displaystyle \frac{1}{x}}{\displaystyle \frac{1}{\sqrt{x^2+\widehat{y}}}}\right)+{\displaystyle \frac{1}{2}}+2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{xdx}{(1+x^2)(\mathrm{exp}2\pi x1)}}`$ (C15)
$`=`$ $`\mathrm{log}\left(1+\sqrt{1+\widehat{y}}\right)\mathrm{log}2+\gamma `$ (C16)
where the constant $`\gamma `$ is equal to

$`\gamma ={\displaystyle \frac{1}{2}}+2{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{xdx}{(1+x^2)(\mathrm{exp}2\pi x1)}}.`$ (C17)
This can be proved by again applying the Abel-Plana formula to the function $`1/m`$, since we know that $`_{m=1}^n1/m\mathrm{log}n+\gamma `$. The constant $`_m(1\mathrm{coth}\pi m)/4\pi m`$ in (C14) can be rewritten as

$`{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{4\pi m}}(\mathrm{coth}\pi m1)`$ $`=`$ $`{\displaystyle \underset{m=1}{\overset{\mathrm{}}{}}}{\displaystyle \frac{1}{2\pi m}}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\mathrm{exp}(2\pi mn)`$ (C18)
$`=`$ $`{\displaystyle \frac{1}{2\pi }}\mathrm{log}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left[1\mathrm{exp}(2\pi n)\right]`$ (C19)
and finally, the results (C14,C16,C19) give the asymptotic behaviour of $`\phi `$ for large $`y`$:

$`\phi (y)={\displaystyle \frac{1}{8\pi }}\mathrm{log}y+{\displaystyle \frac{1}{24}}{\displaystyle \frac{1}{4\pi }}\mathrm{log}4\pi +{\displaystyle \frac{\gamma }{4\pi }}{\displaystyle \frac{1}{2\pi }}\mathrm{log}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left[1\mathrm{exp}(2\pi n)\right]+{\displaystyle \frac{1}{2y}}+\mathrm{}`$ (C20)
The last term comes from a further study of the Abel-Plana formula which gives the other correction terms in the inverse power of $`y`$. An identical analysis gives the finite size magnetisation

$`m=\mathrm{exp}\left({\displaystyle \frac{\tau }{2}}\mathrm{Tr}G/N\right)`$ (C21)
where $`\mathrm{Tr}G/N`$ can be expand as

$`{\displaystyle \frac{1}{N}}\mathrm{Tr}G`$ $`=`$ $`{\displaystyle \frac{1}{4\pi }}\mathrm{log}CN`$ (C22)
$`C`$ $`=`$ $`\mathrm{exp}\left\{{\displaystyle \frac{\pi }{3}}+2\mathrm{log}{\displaystyle \frac{\sqrt{2}}{\pi }}+2\gamma 4\mathrm{log}{\displaystyle \underset{n=1}{\overset{\mathrm{}}{}}}\left[1\mathrm{exp}(2\pi n)\right]\right\}=1.8456`$ (C23)