Finite-size corrections to scaling of the magnetization distribution in the 2d XY -model at zero temperature
aa r X i v : . [ c ond - m a t . s t a t - m ec h ] A p r Finite-size corrections to scaling of the magnetization distributionin the d XY -model at zero temperature G. Palma, ∗ F. Niedermayer, Z. R´acz, A. Riveros, and D. Zambrano Departamento de F´ısica, Universidad de Santiago de Chile, Casilla 307, Santiago 2, Chile Albert Einstein Center for Fundamental Physics, Institute for Theoretical Physics,University of Bern, Switzerland MTA-ELTE Theoretical Physics Research Group, Budapest, Hungary Departamento de F´ısica, Universidad T´ecnica Federico Santa Mar´ıa,Av. Espa˜na 1680, Casilla 110-V, Valpara´ıso, Chile
The zero-temperature, classical XY -model on an L × L square-lattice is studied by exploringthe distribution Φ L ( y ) of its centered and normalized magnetization y in the large L limit. Anintegral representation of the cumulant generating function, known from earlier works, is used forthe numerical evaluation of Φ L ( y ), and the limit distribution Φ L →∞ ( y ) = Φ ( y ) is obtained withhigh precision. The two leading finite-size corrections Φ L ( y ) − Φ ( y ) ≈ a ( L ) Φ ( y ) + a ( L ) Φ ( y )are also extracted both from numerics and from analytic calculations. We find that the amplitude a ( L ) scales as ln( L/L ) /L and the shape correction function Φ ( y ) can be expressed through thelow-order derivatives of the limit distribution, Φ ( y ) = [ y Φ ( y ) + Φ ′ ( y ) ] ′ . The second finite-sizecorrection has an amplitude a ( L ) ∝ /L and one finds that a Φ ( y ) ≪ a Φ ( y ) already for smallsystem size ( L > XY -model at low temperatures, including T = 0. PACS numbers: 05.50.+q, 02.70.-c, 68.35.Rh, 75.40.Cx
I. INTRODUCTION
Studies of the scaling properties of fluctuations haveplayed an important role in developing the theory ofequilibrium critical phenomena, and they also provedto be instrumental in exploring systems driven far fromequilibrium where fluctuations often diverge in the ther-modynamic limit. In particular, finding critical expo-nents through finite-size (FS) scaling has become a stan-dard method for establishing universality classes [1]-[3]and, furthermore, the shapes of the distribution func-tions of fluctuations have also been used as hallmarks ofuniversality classes for both equilibrium [4]-[7] and non-equilibrium systems [8]-[12].One of the most studied distribution of critical order-parameter fluctuations is that of the two dimensionalclassical XY model. The model is important in it-self since it serves as a prime example of an equilib-rium phase transition with topological order [13, 14].Recent interest comes also from a suggestion that itscritical magnetization distribution, P L ( m ), may describethe fluctuations in a remarkable number of diverse far-from-equilibrium steady states. Examples range from theenergy-dissipation in turbulence [9, 10] and interface fluc-tuations in surface evolution [8, 12] to electroconvectionin liquid crystals [15], as well as to river-height fluctua-tions [16]. In some of the examples, such as the surfacegrowth problems described by the Edwards-Wilkinson ∗ Electronic address: [email protected] model [17], one can establish a rigorous link to the low-temperature limit (spin-wave approximation) of the XY model. Furthermore, the non-gaussian features of P L ( m )such as its exponential and double exponential asymp-totes appear to have generic origins [18, 19], thus ex-plaining the observed collapse in a remarkable set of ex-perimental and simulation data. In the majority of theexamples, however, the experiments or simulations pro-vide data of insufficient accuracy to decide unambigu-ously about the universality class. Indeed, it is not easyeven to observe the changes occurring in P L ( m ) if the XY model is considered at finite temperatures where thespin-wave approximation breaks down [20–22].In order to be more confident about universalityclaims, one would like to be able to examine the dis-tribution functions in more details. A well known andmuch investigated method of extracting additional infor-mation in critical systems such as the zero-temperature XY model is the study of the FS behavior [1]-[3]. Theimplication of FS scaling is that the appropriately scaleddistribution function has a well defined limit when thesystem size tends to infinity ( L → ∞ ). A frequently usedscaling variable is the centered and normalized magneti-zation y = ( m − h m i ) /σ m where σ m = h m i − h m i . Thischoice of scaling variable eliminates possible divergencesin h m i and, in general, produces a non-degenerate limitdistribution Φ ( y ) [23]:lim L →∞ Φ L ( y ) ≡ lim L →∞ σ m P L ( h m i + σ m y ) = Φ ( y ) . (1)The function Φ L ( y ) can be expressed in a compact formfor the zero-temperature XY model [24]. Its large- L limit, Φ ( y ), has been numerically evaluated and com-ared with simulations and experiments on a variety ofsystems (see e.g. [9, 16, 18]).Once the limit distribution is known, one can inves-tigate the approach to Φ ( y ) as the system size in-creases. As we shall show below, keeping the lead-ing and the next to leading terms, the FS corrections, δ Φ L ( y ) = Φ L ( y ) − Φ ( y ), to the limit distribution can bewritten as δ Φ L ( y ) = a ( L ) Φ ( y )+ a ( L ) Φ ( y )+ O (cid:0) ln L/L (cid:1) . (2)We kept two terms since the asymptotic L -dependenceof the leading term differs from the next one only by aslowly varying logarithmic factor. Indeed, our calculationyields the amplitudes a i ( L ) in the following form a ( L ) = α ln L + γL = α ln( L/L ) L , a ( L ) = γ ′ L . (3)Here, the coefficient α of the logarithmic term can be ex-pressed through the Catalan constant G as α = 3 π/ G =2 . ... . The coefficients of the 1 /L terms are γ = − .
803 (giving L = 2 . γ ′ = 2 π/ ( y ) charac-terizing the leading shape correction to the limit distri-bution is that it can be expressed through Φ ( y ) asΦ ( y ) = [ y Φ ( y ) + Φ ′ ( y ) ] ′ , (4)where the prime denotes the derivative over y . The sec-ond scaling function Φ ( y ) is more complex in the sensethat it cannot be expressed through a finite sum of thederivatives of the limit distribution. However, it can becalculated efficiently using integral representations as ex-plained in Sec.III [see Eq.(24)] and in Appendix B [seeEq.(B18)].Evaluating the scaling functions Φ i ( y ) numerically,one observes that a ( L )Φ ( y ) ≫ a ( L )Φ ( y ) already forsmall systems ( L > δ Φ L ( y ) ≈ α ln( L/L ) L [ y Φ ( y ) + Φ ′ ( y ) ] ′ . (5)The above expression can be easily calculated and com-pared with experiments and simulations using both thescaling of the amplitude and checking the shape of thecorrection.It is important to note that Φ ( y ) is uniquely deter-mined by the limit distribution, thus the leading shapecorrection has the same universality attributes as thelimit distribution itself. In renormalization group lan-guage, the meaning of the above result is that the eigen-function corresponding to the direction of slowest ap-proach to the fixed point distribution can be expressed through the limit distribution. A simple and transpar-ent derivation of an analogous result for the case of thecentral limit theorem can be found in [25].In order to show the workings of the method of FSscaling, we carried out Monte Carlo (MC) simulations ofthe XY model in its low-temperature limit, including its T = 0-limit. We examined the FS corrections describedabove for the lattice of size L = 10, as an illustration,and found close agreement between the theoretical andMC results.The paper is organized as follows. The XY modeland the zero-temperature magnetization distribution isdescribed in Sec.II. Next, the FS corrections are calcu-lated in Sec.III with the technicalities relegated to twoAppendices (the direct calculations of the needed cumu-lants are found in Appendix A, while an integral repre-sentation that simplifies the evaluation of finite-volumesums is presented in Appendix B). Finally in Sec.IV, andas an illustration, we have compared the FS results withthe MC simulations for a lattice of size L = 10. II. PROBABILITY DENSITY OF THEMAGNETIZATION FOR THE CLASSICAL XY -MODEL IN TWO-DIMENSIONS The two-dimensional, classical XY model on a squarelattice is defined by the Hamiltonian H = − J X h i,j i cos( θ i − θ j ) (6)where the angle variables θ i describe the orientation ofunit vectors in the plane, J > La × La square lattice. From now on we willset J = 1 and a = 1, which defines the energy- and lengthscales in the problem.The order parameter m whose probability distributionis of our interest is defined as m = 1 L X i cos( θ i − ¯ θ ) (7)where ¯ θ is the instantaneous average orientation. Theprobability density function (PDF) of the magnetization, P L ( m ), has been much studied [9, 16, 18–22, 24]. Itszero temperature limit ( T →
0) has been first calculated[24] by using the spin-wave approximation cos( θ i − θ j ) ≈ − ( θ i − θ j ) / P L ( m )is based on evaluating the Fourier transform of the par-tition function of the system in a loop expansion [20].It turns out that the n -loops contribution correspondsto the T n − contribution, where T is the system tem-perature. Correspondingly, the expansion up to one-loopgives exactly the zero-temperature limit of P L ( m ).oth of the above methods yield an identical result toleading order. The second method, however, shows ex-plicitly how to compute the higher temperature correc-tions to P L ( m ) through an effective, T -dependent latticepropagator, Γ = ( I + ikTL G ) − G , where G is independentof T [26], and is defined by its Fourier representation as G ( k ) = 1 / ˆk . (8)Here, the components of the vector ˆk are ˆ k i = 2 sin( k i / k i = 2 πn i /L with i = 1 , n i ∈ Z such thatthe lattice momentum k lies in the first Brillouin zone: − π < k i ≤ π .Hence the PDF at T = 0 can be obtained by the wellknown ‘1 / L ( y ) = Z ∞−∞ dq π r g (cid:18) i r g qy (cid:19) × exp − X k =0 (cid:26) ln (cid:18) − iqL G ( k ) (cid:19) + iqL G ( k ) (cid:27) (9)where y is the centered and normalized magnetization y = ( m − h m i ) /σ m , and the coefficient g is defined asthe particular n = 2 case of the general expression: g n ≡ g n ( L ) = X k =0 G ( k ) n /L n . (10)The PDF in Eq.(9) agrees with the corresponding for-mula obtained in [18], (see Eq.(26) therein).In Fig.1 we plot Φ L ( y ) as a function of y using Eq.(9)for different L values. On the right panel, we display azoom into the region close to the peak of the distributionfunctions. It can be clearly seen that Φ L ( y ) tends to-wards an asymptotic distribution Φ ( y ) as L → ∞ and,furthermore, one observes that the convergence is fast.Nevertheless, in experiments, the value of L is not knownand deviations from Φ ( y ) are observed. They are oftenexplained away as finite-size effects and thus leaving theuniversality claims not entirely justified. We would liketo emphasize that there is information in the FS cor-rections (shown on Fig.1 around the peak of the PDF)and, by evaluating these corrections, one may refine thereasoning for or against finding a universality class. −3 −2 −1 0 1 2 300.050.10.150.20.250.30.350.40.45 y = (m−
The scaled probability density of the magnetization Φ L ( y ) , computed numerically using Eq. (9) , is displayed forvarious lattice sizes L . The remarkably fast convergence tothe limit distribution Φ ( y ) is demonstrated in the left panelwhile, to illustrate the corrections to the limit distribution, thepeak region of the Φ L ( y ) is magnified in the right panel. III. FINITE-SIZE CORRECTIONS TO THELIMIT DISTRIBUTION
As explained in Sec.II, the PDF of the magnetizationfor the 2 d XY -model at T = 0 is given by the 1-loop an-alytic expression of Eq.(9), and the numerical evaluationof the L → ∞ limit distribution, Φ ( y ), can be carriedout with an excellent precision. The aim of this section isto present the steps of the calculation of the leading andnext-to-leading FS corrections to the limit distribution.We start by expanding the logarithm in the expo-nential on the r.h.s. of Eq.(9) which allows rewrit-ing the equation in terms of the coefficients g n definedby Eq.(10). After rescaling the integration variable by p g /
2, we obtain Φ L ( y ) as the following Fourier integralΦ L ( y ) = ∞ Z −∞ dk π exp (cid:26) iky − k + F L ( k ) (cid:27) , (11)where we have defined: F L ( k ) = X n ≥ g n n (cid:18) ik r g (cid:19) n . (12)The L dependence of the above sum is in the coeffi-cients g n which have a finite non-zero thermodynamiclimit g n ( L → ∞ ) = g ∞ n for n ≥
2. Thus, in order tocompute the FS behavior of Φ L ( y ), we shall have to de-termine the FS corrections to g ∞ n g n = g ∞ n + δg n . (13)Assuming that δg n is known, we can write F L ( k ) as F L ( k ) = F ( k ) + δF ( k ) + δF ( k ) (14)where F ( k ) is the thermodynamic limit of F L ( k ): F ( k ) = X n ≥ g ∞ n n ik s g ∞ ! n (15)nd the FS corrections, due to δg and δg n for n ≥
3, arewritten separately δF ( k ) = − δg g ∞ X n ≥ g ∞ n ik s g ∞ ! n , (16) δF ( k ) = X n ≥ δg n n ik s g ∞ ! n . (17)The separation of the δg and the δg n ≥ contributions ispartly motivated by their L → ∞ asymptotic behaviors.It will be shown in the Appendices that δg ∼ ln L/L while for n ≥ δg n ∼ /L . Thus the leading correctioncomes from δF ( k ) and the sum in Eq.(16) determinesthe shape (the functional form) of the leading correction.As it turns out, the same shape correction can be easilyseparated from the δg n contributions for n ≥
3. Indeed,for large n , one has δg n ∼ ng ∞ n , and one can write (cf.Eqs.(A6), (B4) and (B15)) δg n = π L ng ∞ n + 2 π L c n , (18)where the second term is suppressed relative to the firstone by a factor 4 − n . Substituting the above split of δg n into Eq.(17), one can see the emergence of the same sumas in Eq.(16), and it allows us to write δF ( k ) + δF ( k ) = a ( L )Ψ ( k ) + a ( L )Ψ ( k ) (19)where the L -dependent amplitudes are given by a ( L ) = δg ( L )2 g ∞ − π L , a ( L ) = 2 π L , (20)while the corresponding L -independent functions byΨ ( k ) = − X n ≥ g ∞ n ik s g ∞ ! n = − k ddk F ( k ) , Ψ ( k ) = X n ≥ c n n ik s g ∞ ! n . (21)Integral representations for F ( k ) and Ψ ( k ) are given inEqs.(B5) and (B18), respectively, in Appendix B.Inserting Eq.(19) into Eq.(11), we obtain the PDF ofthe 2 d XY -model at zero-temperature, including its lead-ing FS correctionsΦ L ( y ) = Φ ( y ) + a ( L ) Φ ( y ) + a ( L ) Φ ( y ) . (22)Here Φ ( y ) = ∞ Z −∞ dk π exp (cid:26) iky − k F ( k ) (cid:27) , (23) andΦ , ( y ) = ∞ Z −∞ dk π exp (cid:26) iky − k F ( k ) (cid:27) Ψ , ( k ) . (24)Since Ψ ( k ) = − kF ′ ( k ), one can evaluate Φ ( y )through Φ ( y ) by carrying out the appropriate integra-tions by part on the r.h.s. of Eq.(24). The outcome isone of the main results of our work, namely a simple ex-pression is obtained for the leading shape correction interms of the limit distribution (as quoted in Eq.(4)):Φ ( y ) = [ y Φ ( y ) + Φ ′ ( y )] ′ . (25)The second shape correction Φ ( y ) can not be relatedto Φ ( y ) in such a simple manner but can be readilyevaluated using Eqs.(24) and (B18) of the Appendix B.The functions Φ ( y ) and Φ ( y ) are displayed in Fig.2.A general property of these functions is that their 0 th , 1 st and 2 nd moments are zero. This follows from their defi-nition as being corrections to a centered and normalizedprobability distribution, and can be explicitly verified us-ing their definitions, e.g. Eq.(24). −5 −4 −3 −2 −1 0 1 2 3 4−0.5−0.4−0.3−0.2−0.100.10.20.3 y = (m−
Behavior of the scale independent shape correctionfunctions Φ ( y ) and Φ ( y ) . In order to determine the amplitudes in front of theshape corrections, we have to calculate the FS behaviorof δg n for n ≥
2. This problem is addressed in details,using two distinct methods, in the Appendices A and B.It is found there that the asymptotic L -dependence of δg has the form ( C ln L + C ) /L (cf. Eq.(A3)) while δg n behaves as ∝ /L for n > ( y ) given by a ( L ) = α ln L + γL = αL ln LL (26)where the coefficient α is obtained analytically α =3 π/ (4 G ) = 2 . G being the Catalan con-stant, while γ = − . L = 2 . a ( L ) is the second main result of ourwork since the FS corrections are dominated by the a ( L )Φ ( y ) term. Indeed, a ( L ) /a ( L ) is small, 0 .
67, al-ready for L = 10 and it decreases with increasing L . Fur-thermore, evaluating Φ ( y ) numerically, one finds thatapart from the neighborhood of the zeros of Φ ( y ), theinequality a ( L )Φ ( y ) ≫ a ( L )Φ ( y ) holds for L >
10 asseen in Fig.3. −2 −1 0 1 2−15−10−505x 10 −3 L = 10 a (L) Φ (y)a (L) Φ (y) −2 −1 0 1 2−505x 10 −3 L = 16 a (L) Φ (y)a (L) Φ (y)−2 −1 0 1 2−3−2−1012x 10 −3 y = (m−
Comparison of the two leading correction terms a ( L )Φ ( y ) and a ( L )Φ ( y ) for L = 10 , , and 128. It should be mentioned that there is some freedomin separating the two contributions a ( L )Φ ( y ) and a ( L )Φ ( y ) in Eq.(2). One can replace Φ ( y ) by Φ ( y ) + c Φ ( y ) changing simultaneously γ to γ − cγ ′ in a ( L ),Eq.(20). Our choice of separating the leading large- n asymptote of δg n in Eqs.(18) and (21) leads naturallyto a function proportional to Φ ( y ) and also results in aremnant that is small already for small values of L . Thischoice is convenient, since it allows to write the FS cor-rection in a compact form [see Eq.(5)] with a very goodaccuracy as demonstrated in Fig.4.It is remarkable that the dominant correction term a ( L )Φ ( y ) also emerges from a simple assumption aboutthe PDF written in its original variable. Namely, if weassume that one can write P L ( m ) = P ( m ) + ǫ ( L ) P ′′ ( m )with ǫ → L → ∞ , we find that ǫ ( L ) = a ( L ) σ m . Using then the scaled variable y , the expres-sion ǫ ( L ) P ′′ ( m ) becomes a ( L )Φ ( y ).As discussed above, there are other choices for sepa-rating a contribution proportional to Φ ( y ). It should beclear, however, that the freedom is irrelevant when thesum of the two contributions a ( L )Φ ( y ) + a ( L )Φ ( y )is used. As expected, and as can be seen in Fig.5, theconvergence is fastest when the sum of both correctionsare used. −2 0 2−0.015−0.01−0.00500.0050.01 L = 10 Φ L (y) − Φ (y)a (L) Φ (y) −2 0 2−3−2−1012x 10 −3 L = 32 Φ L (y) − Φ (y)a (L) Φ (y)−2 0 2−10−505x 10 −4 y = (m−
Comparison of the exact FS correction Φ L ( y ) − Φ ( y ) with the leading term a ( L )Φ ( L ) given by (25) and (26) .Apart from the regions close to the maxima and minima, goodagreement can be observed already for small system sizes. −2 0 2−0.02−0.0100.01 L = 10 Φ L (y) − Φ (y)a (L) Φ (y) + a (L) Φ (y) −2 0 2−4−202x 10 −3 L = 32 Φ L (y) − Φ (y)a (L) Φ (y) + a (L) Φ (y) −2 0 2−15−10−505x 10 −4 y = (m−
Comparison of the exact FS correction Φ L ( y ) − Φ ( y ) with the leading terms a ( L )Φ ( L ) + a ( L )Φ ( y ) for latticesizes L = 10 , , and 128. IV. MONTE CARLO SIMULATIONS
Here, we briefly demonstrate that the calculated FScorrections can be observed in simulations. It is clearthat increasing the system size and improving the statis-tics by increasing the number of Monte Carlo samples,the leading FS corrections, as seen in Fig.5, will emergefrom the analysis. The question is whether the leadingcorrections we calculated could be seen already in smallsystems with reasonable simulation effort.We have thus performed MC simulations on a 2 d XY -model of size L = 10 and computed Φ L ( y ). Since our an-lytic results Eqs.(22)-(24) pertain to the T = 0 limit ofthe system, we made simulations deep in the low T region( T = 0 .
04 and 0 .
02) using the over-relaxation Metropolis(ORM) algorithm [27]. In addition, simulations of the T = 0 limit itself could also be carried out since therethe spin-wave approximation applies which yields inde-pendent modes with Gaussian action whose simulation isstraightforward [28]. −3 −2 −1 0 1 2 3−15−10−505x 10 −3 y = (m−
Comparison of the exact FS correction Φ L ( y ) − Φ ( y ) (red line) and the leading terms a ( L )Φ ( L ) + a ( L )Φ ( y ) (blue dashed line) with the MC simulations for lattice size L = 10 . The region close to the minimum is magnified inthe right panel. It shows that the ORM results approach to Φ L ( y ) − Φ ( y ) as T → . In the simulation using ORM, we typically measuredobservables using 100 × sweeps, and the errors wereestimated by using a binning method. On Fig.6, the redline shows the FS corrections to Φ ( y ) displayed as thedifference Φ L =10 ( y ) − Φ ( y ) computed from the integralrepresentation of Eq.(9). It is compared with data ob-tained by ORM at very low temperatures T = 0 .
04 (greencircles) and T = 0 .
02 (asterisk markers). The magentapoints represent the results from the simulation of theGaussian action at T = 0. The statistical errors on allthe data points displayed are smaller than the point size,which is not surprising for the quoted large number of MCsweeps. From the actual statistical error ∼ × − onecan find that a relatively small number, 3 × sweepsare enough to reach an accuracy 10 − , one tenth of themaximal FS correction.As one can see, the simulation temperatures used aresmall enough for the temperature corrections to be smallcompared to the FS corrections at L = 10. It can alsobe seen that the the sum of the two leading FS correc-tion terms, a ( L )Φ ( y ) + a ( L )Φ ( y ) (blue dashed line)is quite close to the exact result. Of course, L = 10 is arather small system to expect full agreement of the calcu-lated leading FS corrections with the full FS correction.Looking at Fig.5, however, one notes that increasing thesize of the system by only a factor three would yield acomplete domination of the leading FS corrections overthe higher order FS corrections. Thus, we conclude that observing the leading FS cor-rections is feasible in relatively small systems at relativelysmall computational coast. Based on this observation,we expect that a meaningful analysis of FS correctionsin experimental XY systems is also possible. V. CONCLUSIONS
We have computed the leading FS corrections to thePDF of the magnetization of the 2 d XY -model at zerotemperature. Two scale-independent functions Φ ( y )and Φ ( y ) were found with their amplitudes behavingwith system size as a ( L ) ∼ α ln( L/L ) /L and a ( L ) ∼ /L . The function Φ ( y ) can be expressed throughthe limit distribution Φ ( y ) and its low-order derivatives.This makes it a candidate for identifying universality fea-tures hidden in FS corrections.The leading and next to leading corrections were foundto describe the FS behavior very accurately already forsmall system size. Thus, as our MC simulations demon-strated, the observation of the calculated FS correctionsis possible in model systems. We expect that their ex-perimental observation may also be feasible. Acknowledgments
This work was partially supported by Dicyt-USACHGrant No. 041531PA and pai-conicyt 79140064, and bythe Hungarian Research Fund (OTKA NK100296). G. P.would like to thank the Institute for Theoretical Physicsat E¨otv¨os University, Budapest, for the invitation in thesummer of 2014, where this collaboration started and alsothanks to their members for the kind hospitality.
Appendix A: Finite-size corrections to δg n We begin by deriving the large- L asymptotic expansionfor g defined in Eq.(10). It is convenient to write thisequation in the form: g = 1(2 π ) X l,m ′ " l + m − (cid:18) πL (cid:19) ( l + m ) + · · · − = g ∞ + δg (A1)where the sum goes from l, m = − L/ L/ l = m = 0 term is left out. Theasymptotic value g ∞ is given by g ∞ = 1(2 π ) ∞ X l,m = −∞′ l + m ) = G/ (24 π ) ≈ . G is the Catalan’s constant.In the FS correction δg , the sum giving the coefficientof 1 /L diverges logarithmically with L . This leadingterm is obtained as δg ∼ π ) L X l + m Behavior of g ( L ) as a function of lattice size L upto L = 10 in logarithmic scale. The red dots correspond tothe direct numerical evaluation of Eq.(10) and the blue linecorresponds to the analytic expression of Eq.(A4) for γ =0 . . Similarly to Eq.(A2) the asymptotic value g ∞ n is given by g ∞ n = 1(2 π ) n ∞ X l,m = −∞′ l + m ) n = 4(2 π ) n ζ ( n ) β ( n ) (A5)where ζ ( n ) and β ( n ) are Riemann’s zeta function andDirichlet’s beta function, respectively [31].For n ≥ /L correctionterm δg n converges for L → ∞ hence one can extend thesummation to ±∞ (up to an error decreasing faster than1 /L ). One has then δg n ng ∞ n = (2 π ) L P l,m ′ ( l + m )( l + m ) − ( n +1) P l,m ′ ( l + m ) − n = π L (cid:18) U n V n (cid:19) (A6)where U n = X l,m ′ l − l + m − m ( l + m ) n +1 V n = X l,m ′ ( l + m ) − n = 4 ζ ( n ) β ( n ) . (A7)For large n the dominant terms in these sums come fromsmallest | l | , | m | with non-vanishing contributions. Thenumerator in U n vanishes for the two shells ( l = ± , m =0), ( l = 0 , m = ± 1) and ( l = ± , m = ± n is coming from ( l = ± , m = 0), ( l =0 , m = ± 2) and is given by 12 · − n . Since V n = 4(1 +2 − n + . . . ) one finds that U n /V n ∼ · − n . This is a smallcorrection – even for n = 3 it is just ≈ . δg n ng ∞ n = π L + 1 L O (cid:0) − n (cid:1) . n ≥ Appendix B: Integral representation for leadingfinite-size corrections For calculating finite-volume sums for a cubic box ofsize L in d dimensions in the continuum, like P ′ k / ( k ) n ,where k = P di =1 (2 πn i /L ) , it is useful to introduce thefunction S ( x ) (see e.g. [32]), (related to Jacobi’s thetafunction) defined as S ( x ) = ∞ X n = −∞ exp (cid:0) − πxn (cid:1) . (B1)It satisfies the relation S ( x ) = 1 √ x S (cid:18) x (cid:19) (B2)hich allows to calculate S ( x ) very precisely by takingonly a few terms in the sum, both for x < x > S ( x ) = x − / (1+2 exp( − π/x )+ . . . ) for x → S ( x ) = 1 + 2 exp( − πx ) + . . . for large x .As an illustration, it is easy to show that g ∞ n given byEq.(A5) has an integral representation g ∞ n = 1(4 π ) n Γ( n ) ∞ Z d x x n − (cid:0) S ( x ) − (cid:1) . (B3)The leading term for n → ∞ is given by the large- x behavior of the integrand. Separating it, one obtains anexpression g ∞ n = 4(2 π ) n + 1(4 π ) n Γ( n ) × ∞ Z d x x n − (cid:0) S ( x ) − − − πx (cid:1) (B4)which can be evaluated and shown to be in agreementwith Eq.(A5). With the help of this one can perform thesummation in Eq.(15) yielding F ( k ) = Z ∞ d u u (cid:2) S ( uw ) − (cid:3) × (cid:26) e iku − (cid:18) iku − k u (cid:19)(cid:27) (B5)where w = 4 π p g ∞ / 2. The Fourier transformation ap-pearing here can be performed efficiently by a fast Fouriertransform (FFT).This technique can be generalized to finite-volume lat-tice sums by introducing [33] Q L ( z ) = 1 L L − X l =0 exp (cid:16) − z ˆ k l (cid:17) = φ ( z ) + 2 ∞ X m =1 φ mL ( z ) (B6)where ˆ k l = 2(1 − cos(2 πl/L )), l = 0 , . . . , L − 1, and φ n ( z ) = e − z I n (2 z ) (B7)with I n ( z ) being the modified Bessel function. For large z one has φ ( z ) ∼ / √ πz .For fixed z with increasing L the function Q L ( z ) ap-proaches φ ( z ) exponentially fast. The approach be-comes slower with increasing z , but even when the ar-gument increases slower than L one still haslim L →∞ ( Q L ( cL α ) − φ ( cL α )) = 0 for α < L . This is not true for z ∝ L , and for this case one obtains another scaling function. Rescaling Q L ( z ),we introduce the lattice counterpart [33] of S ( x ) by S L ( x ) = LQ L (cid:18) xL π (cid:19) (B9)By expanding Eq.(B6) for large L one finds the asymp-totic expansion S L ( x ) = S ( x ) + π L xS ′′ ( x ) + O (cid:18) x L (cid:19) . (B10)As the error term indicates, the approach to L = ∞ isnot uniform in x .Using S L ( x ) one has for g n = g n ( L ) two integral rep-resentations g n = L − n +2 Γ( n ) Z ∞ d z z n − (cid:20) Q L ( z ) − L (cid:21) = 1(4 π ) n Γ( n ) Z ∞ d x x n − (cid:2) S L ( x ) − (cid:3) . (B11)We outline below the calculation of δg to O (cid:0) /L (cid:1) . Dueto the non-uniform convergence for L → ∞ it is usefulto split the integration region and write g ( L ) = 1 L Z z d z z (cid:20) Q L ( z ) − L (cid:21) +1(4 π ) Z ∞ x d x x (cid:2) S L ( x ) − (cid:3) (B12)where x = 4 πz /L . Choosing z = z ( L ) = cL − ǫ withsome fixed small ǫ > Q L ( z ) by φ ( z ) up to exponentially small corrections.Similarly, in the second integral one can use the expan-sion Eq.(B10). Note that z ( L ) → ∞ and x ( L ) → L → ∞ . Using Eq.(B10) and neglecting terms vanishingfaster than 1 /L one obtains δg = 1 L Z z d z zφ ( z ) + 2 π L (4 π ) Z ∞ x d x x S ( x ) S ′′ ( x )(B13)Separating the asymptotic behavior of the integrands forlarge z , and small x , respectively, one obtains the log-arithmic contribution 1 / (32 π ) ln( z /x ) = 1 / (16 π ) ln L ,and in the remaining terms one can make the substitu-tions z = ∞ and x = 0. Evaluating the correspondingintegrals one reproduces the fit result Eq.(A4) to all dig-its (cf. Fig.7).The leading correction of δg n for n > δg n = 2 π L π ) n Γ( n ) Z ∞ d x x n S ( x ) S ′′ ( x ) , n ≥ x term of the integrand one obtains δg n = π L n (2 π ) n + 2 π L π ) n Γ( n ) × Z ∞ d x x n (cid:2) S ( x ) S ′′ ( x ) − π e − πx (cid:3) (B15)here n ≥ 3. The leading term has the same form asfor ng ∞ n (cf. Eq.(B4)). Subtracting this way the leadingterm one can define c n by δg n = π L ng ∞ n + 2 π L c n (B16)where for n ≥ c n = 1(4 π ) n Γ( n ) ∞ Z d x x n S ( x ) ( S ′′ ( x ) + πS ′ ( x )) . (B17)For large n one has c n ∼ πn (4 π ) − n , i.e. it is sup- pressed by a factor of 4 − n compared to ng ∞ n .Inserting Eq.(B17) into Eq.(21) we obtain an integralrepresentation for Ψ ( k ),Ψ ( k ) = w ∞ Z d u S ( uw ) { S ′′ ( uw ) + πS ′ ( uw ) } × (cid:20) e iku − (cid:18) iku − k u (cid:19)(cid:21) . (B18)Using FFT one can evaluate this and finally the corre-sponding correction Φ ( y ) to the PDF. [1] M. E. Fisher in Critical Phenomena, Proceedings ofthe 1970 International School of Physics Enrico Fermi,Course 51, Edited by M. S. Green (Academic, New York,1971).[2] J. L. Cardy, Finite Size Scaling, (North-Holland, Ams-terdam, 1988).[3] Finite Size Scaling and Numerical Simulation of Statis-tical Physic s, edited by V. Privman (World Scientific,Singapore, 1990).[4] A. D. Bruce, J. Phys. C , 3667 (1981).[5] K. Binder, Z. Phys. B , 119 (1981).[6] D. Nicolaides, A. D. Bruce, J. Phys. A , 233 (1988).[7] Examples of wide-ranging applications are: A. D. Bruceand N. B. Wilding, Phys. Rev. Lett. , 193 (1992)(liquid-gas transition), D. Nicolaides and R. Ewans Phys.Rev. Lett. , 778 (1989) (wetting); N. B. Wilding and P.Niebala, Phys. Rev. E , 926 (1996) (tricritical point);M. M¨uller and N. B. Wilding, Phys. Rev. E , 2079(1995) (polymers); S. L. A. de Queiroz and R. B. Stinch-comb, Phys. Rev. E , 036117 (2001) (random-fieldIsing model); M. M. Tsypin, Phys. Rev. Lett. , 2015(1994) (field theory).[8] G. Foltin, K. Oerding, Z. R´acz, R. L. Workman and R.K. P. Zia, Phys. Rev. E , R639 (1994); Z. R´acz andM. Plischke, Phys. Rev. E , 3530 (1994); E. Marinari,A. Pagnani, G. Parisi, and Z. R´acz, Phys. Rev. E ,026136 (2002).[9] S. T. Bramwell, P. C. W. Holdsworth, J. -F. Pinton, Na-ture , 552 (1998); R. Labb´e, J. -F. Pinton, P. C. W.Holdsworth, Phys. Rev. E , R2452 (1999).[10] V. Aji and N. Goldenfeld, Phys. Rev. Lett. , 1007(2001)[11] S. T. Bramwell et al., Phys. Rev. Lett. , 3744 (2000).[12] G. Korniss, Z. Toroczkai, M. A. Novotny, and P. A.Rikvold, Phys. Rev. Lett. , 1351 (2000); S. Lubeckand P. C. Heger, Phys. Rev. Lett. , 230601 (2003);T. Halpin-Healy and G. Palasantzas, EPL , 50001(2014).[13] J. M. Kosterlitz and D. J. Thouless, J. Phys. C , 1181(1973).[14] P. M. Chaikin and T. C. Lubensky, Principles of Con-densed Matter Physics (Cambridge University Press,Cambridge, 1996).[15] T. T´oth-Katona and J. Gleeson, Phys. Rev. Lett. , 180601 (2008).[16] S. T. Bramwell, T. Fennel, P. C. W. Holdsworth, andB. Portelli, EPL , 310 (2002); K. Dahlstedt and H.Jensen, Physica A , 596 (2005).[17] S. F. Edwards and D. R. Wilkinson, Proc. R. Soc. Lon-don, Series A , 17 (1982).[18] S. T. Bramwell, J.-Y. Fortin, P. C. W. Holdsworth, S.Peysson, J.-F. Pinton, B. Portelli, and M. Sellito, Phys.Rev. E , 041106 (2001).[19] S. T. Bramwell, Nature Physics , 444 (2009).[20] G. Mack, G. Palma and L. Vergara, Phys. Rev. E ,026119 (2005).[21] G. Palma, Phys. Rev. E. , 046130 (2006).[22] S. T. Banks and S. T. Bramwell, J. Phys. A , 5603(2005).[23] T. Antal, M. Droz, G. Gy¨orgyi, and Z. R´acz, Phys. Rev.Lett. , 240601 (2001); Phys. Rev. E. , 046140 (2002).[24] P. Archambault, S. T. Bramwell, J.-Y. Fortin, P. C. W.Holdsworth, S. Peysson, J.-F. Pinton, J. Appl. Phys. ,7234 (1998).[25] G. Jona-Lasinio, Phys. Rep. , 439 (2001).[26] J. Cardy, Scaling and Renormalization in StatisticalPhysics, Cambridge University Press (1996).[27] M. Creutz, Phys. Rev. D , 515 (1987); K. Kanki, D.Loison and K.D. Schotte, Eur. Phys. J. B , 309?315(2005)[28] Integrating out the gaussian variables one obtains theanalytic expression in Eq.(9).[29] R. Kenna and A.C. Irving, Phys. Lett B , 273 (1995);W. Janke, Phys. Rev. B (1997); S.G. Chung, Phys.Rev. B , 16 (1999)[30] Note that in the double sum appearing in (10) the sum-mation over one of the integers can be done analytically,which makes possible to reach large values of L .[31] See, e.g., I. S. Gradstein and I. M. Ryshik, Table of In-tegral series and products, seventh edition, Elsevier Aca-demic Press (2007)[32] P. Hasenfratz and H. Leutwyler, Nucl. Phys. B , 241(1990).[33] F. Niedermayer and P. Weisz, Massless sunset diagramsin finite asymmetric volumes ,,