Finite size scaling of the 5D Ising model with free boundary conditions

There has been a long running debate on the ﬁnite size scaling for the Ising model with free boundary conditions above the upper critical dimension, where the standard picture gives an L 2 scaling for the susceptibility and an alternative theory has promoted an L 5 / 2 scaling, as would be the case for cyclic boundary. In this paper we present results from simulation of the far largest systems used so far, up to side L = 160 and ﬁnd that this data clearly supports the standard scaling. Further we present a discussion of why rigorous results for the random-cluster model provide both supports for the standard scaling picture and a clear explanation of why the scalings for free and cyclic boundary should be different. ©

Much less has been written on the subject of free boundary conditions above the upper critical dimension, but see e.g.[10,11].As can be seen from the references in those two papers there has been some debate on whether the standard scaling picture, saying that e.g. the susceptibility scales as L 2 for free boundary, holds or whether an alternative theory proposing that it scales as L 5/2 is correct.In [10] the current authors simulated the 5-dimensional model on larger systems than previous authors and found that the data supported the standard scaling picture.In a reply [11] it was again suggested that the alternative picture is correct and that the results of [10] were due to too small systems, dominated by finite size effects stemming from their large boundaries.
The purpose of this paper is two-fold.We have extended the 5-dimensional simulations with free boundary from [10] to much larger systems, up to L = 160, where the boundary vertices make up less than 6.1% of the system.First, using the new data, we give improved estimates of the critical energy, specific heat and several other quantities at the critical point.Second, we compare how well the standard scaling and the alternative theory fit our new large system data, and discuss why, based on mathematical results on the random cluster model, there are good reasons for expecting the standard picture to be the correct one, as the data also suggests.

DEFINITIONS AND DETAILS
For a given graph G the Hamiltonian with interactions of unit strength along the edges is H = − ij S i S j where the sum is taken over the edges ij.As usual K = 1/k B T is the dimensionless inverse temperature and we denote the temperature equilibrium mean by • • • .The susceptibility is defined as χ = N m 2 , where m = (1/N ) i S i is the magnetisation per spin, and the specific heat as C = N U 2 − U 2 , where U = (1/N ) ij S i S j is the energy per spin, and for short we write U = U .
The underlying graph in question is an L×L×L×L×L grid graph with free boundary conditions, or equivalently, the cartesian product of five paths on L vertices.
For 3 ≤ L ≤ 63 we kept 64 separate systems at nine couplings K = 0.1139130, 0.1139135, . . .0.1139170.For L = 72, 96 we used 48 separate systems, 116 for L = 128 and 108 for L = 160.For these larger systems the number of measurements were just short of 30000 for L = 72 down to about 4000 for L = 160.Means and standard errors were estimated by exploiting the separate systems.For the larger systems (L ≥ 72), due to the comparably few measurements, we also used bootstrapping on the entire data set for estimating standard errors.To doublecheck for equilibration problems we compared with subsets of the data after rejecting early measurements.
The different couplings for 3 ≤ L ≤ 63 showed no discernible difference in their scaling behaviour for L ≤ 63.Based on the behaviour for L ≤ 63 we have designated K c = 0.1139150, which is a little higher than what we used in Ref. [10] and marginally lower than that used in e.g.Ref. [11].The lion's share of sampling were then made at 0.1139150 and for L ≥ 72 we have measured only at K c .
For all sizes we measured magnetisation and energy, storing their moment sums.For 3 ≤ L ≤ 63 we also measured many properties regarding the clusters that were generated and used for consistency checks, and one of them will be shown in the Discussion section.

Geometry and boundary effects
A potentially important issue for systems with free boundary condition is the size of the boundary, and in particular the fraction of vertices on the boundary.Of the N = L 5 vertices in the graph (L − 2) 5 are inner vertices and thus L 5 − (L − 2) 5 vertices sit on the boundary.The fraction of boundary vertices is then 1 − (1 − 2/L) 5 .For L = 8 this means that the boundary constitutes no less than 76% of all vertices.To continue, for L = 16 the boundary's share is 49%, for L = 32 it is 28%, for L = 64 it is 15%, for L = 128 it is 7.6% and for L = 160, our largest system studied here, it makes up 6.1%.So our largest system have a clear minority of their vertices on the boundary.
Another important measure is the number of vertices with a given minimum distance to the boundary.If we consider the cube of side cL around the centre vertex of the cube, i.e. the set of vertices with distance at least cL 2 to the boundary we find that it contains a c 5 fraction of the N vertices.That means that at least 50% of the vertices are at a distance of at most 0.065L from the boundary, for every L. Similarly, the central cube with side L/2 contains just 3.1% of the vertices of the cube.
This means that even in the limit the effect of vertices close to the boundary will always be large, and that the vertices close to the central vertex will also remain atypical for any property which depend both on the distance to the boundary and the majority of the vertices in the cube.In particular we should expect such properties to be bounded from above the corresponding values in the infinite system thermodynamic limit, if they tend to decrease with the distance to the boundary.

ENERGY AND SPECIFIC HEAT
It is known [2] that in the limit ,for d ≥ 5, the specific heat, i.e. the energy variance, is bounded for all temperatures, but the value is not known, and likewise for the critical energy.In Fig. 1 the mean energy U is shown versus 1/L.The leading scaling term is here set to the order 1/L and the correction term to order 1/L 3/2 .This gave by far the most stable coefficients of the fitted curve among the simple exponents.We find that the best fitted curve is 0.675647(3) − 1.013 (1) x + 0.395(1) x 3/2 , where x = 1/L.The fit is excellent down to L = 3.The coefficients and their error estimates are here based on the median and interquartile range of the coefficients when deleting one of the data points from the fitting process.
In the inset picture in Fig. 1 we zoom into the plot by showing L(U − U c ) versus 1/L 1/2 together with the line −1.013+ 0.395x.The fit is vey good and hence we conclude that the correction term is of the order 1/L 1/2 .Note that the error bars in both plots are included but they are far too small to be seen at this scale.The limit energy U c = lim L→∞ U(K c , L) = 0.675647(3) is only marginally larger than the value we gave in [10] which may be explained by the slightly smaller K c .Note that for free boundary conditions the limit is reached from below whereas for periodic boundary conditions U approaches its limit from above, roughly as U(K c , L)−U c ∼ 5.5/L 5/2 [10].
For the specific heat we note that the error bars are noticeable, but this is to be expected.It is not known at which rate C(K c , L) approaches its asymptotic value C c but judging from the excellent line-up of the points in Fig. 2 it appears C(K c , L)− C c ∝ 1/L 1/3 .This is the only time we see a 1/3 in a scaling exponent for free boundary conditions and we have no theoretical basis for it.In Fig. 2 we show C(K c , L) versus 1/L 1/3 and the line 14.69(1) − 14.93 (2)x, where x = 1/L 1/3 .As before the error estimates of the coefficients are based on the variability of a fitted curve after deleting one of the points.We estimate thus that C c = lim L→∞ C(K c , L) = 14.69 (1).The inset picture of Fig. 2 zoom into the correction term by plotting L 1/3 (C − C c ) versus 1/L 1/3 together with the constant line −14.93.The error bars now become quite noticeable, especially for the larger L.There is no clear trend upwards or downwards in the data points which suggests that any further corrections to scaling must be truly negligible.

MAGNETISATION AND SUSCEPTIBILITY
The scaling of the modulus of the magnetisation |m| for free boundary conditions is very different from that of periodic boundary conditions.In the first case we find |m| ∝ L −3/2 whereas in the second it is wellknown that |m| ∝ L −5/4 .In the free boundary case we note the need for correction to scaling.Indeed, if we want perfectly fitted curves down to L = 3 we need two correction terms.We have instead chosen to ignore L ≤ 7 and stay with just one correction term for the remaining 13 points.In Fig. 3  As we mentioned above the susceptibility χ = N m 2 scales to leading order as χ ∝ L 5/2 for periodic boundary conditions but it is not known what the corresponding order is for free boundary conditions.We find here, as in Ref. [10], that χ ∝ L 2 is by far the best scaling rule.In Fig. 4 we show χ/L 2 versus 1/L for 7 ≤ L ≤ 160 together with the line 0.08269(2) + 0.8174 (3) x.The coefficients were determined after excluding L ≤ 5 from the fitting process.Had we included the two smallest systems an extra correction to scaling term would have been required.Again we zoom in and the inset of Fig.   shows (χ/L 2 − 0.08269)L versus 1/L together with the constant line 0.8174.Though the error bars are quite big for the largest systems the fit is quite acceptable.
In short we find χ(K c ) ∼ 0.08269(2)L 2 .We might add that the corresponding expression for periodic boundary is not known exactly but we suggested recently [8] that χ ∼ 1.742LIt has been suggested [13] that L 2 is in fact not the correct scaling for the susceptibility.The authors of [13] claim that the correct scaling should be L 5/2 , as for the case with cyclic boundary conditions, and that the exponent previously found by us, and other authors, are based on either finite size effects due to too many boundary vertices in small systems, for simulation studies, or incomplete theory.In our previous work the boundary did indeed contain a large fraction of the system's vertices but in our current study this fraction has been reduced to a lower value than in any previous study, including the truncated systems used in [13].
To avoid implicit bias in our scaling of the susceptibility to we can also test the ratio χ/L 5/2 .If the claims of [13] are correct this quantity should converge to a finite non-zero limit, at least for large enough systems, and if the standard scaling is correct it should to leading order converge to 0 as L −0.5 .We make a scaling ansatz c 0 + c 1 x λ1 + c 2 x λ2 , where x = 1/L, and let Mathematica find the five free parameters using a least squares fit, after excluding L = 3, 5.As usual, we let each remaining point be deleted in turn from the fitting data to obtain error bars of the parameters.We find on average the curve 0.0000(4) + 0.085(9) x 0.51(3) + 0.820(7) x 1.51 (3) .Clearly the parameters we estimated above for χ/L 2 falls inside these estimates, though the error bars are a magnitude larger here.Using the middle point values we plot the curve together with the data points in Fig. 5.The data for large systems is clearly consistent with the standard scaling, even for an unrestricted data fitting like this.

FOURTH MOMENT AND KURTOSIS
Our final property of interest is the fourth moment of the magnetisation at K c .We find that the fourth moment of the magnetisation scales as m 4 ∝ L −6 .The general rule would then be |m k | ∝ L −3k/2 whereas the corresponding rule for periodic boundary conditions is |m k | ∝ L −5k/4 .Proceeding in the same manner as before, we plot the normalised fourth moment's behaviour as m 4 L 6 versus 1/L in Fig. 6 together with the estimated polynomial 0.02051(3) + 0.4045 (8) x + 1.989(4) x 2 .We excluded L = 3 from the coefficient estimates.To leading order we thus find m 4 ∼ 0.02051(3)L −6 .
The moment ratio Q = m 4 / m 2 2 , or kurtosis, indicates the shape of the underlying magnetisation distribution.In Fig. 7 we show the kurtosis versus 1/L and the line 3−0.14 x.The error bars are based on the formula for the error of a quotient, d(x/y 2 ).The line is based on the coefficient estimates for m 2 and m 4 above by taking the quotient of their respective series expansions in the standard fashion.Inserting the coefficients and their error estimates gives the limit m 4 / m 2 2 → 3.000 (6) which is the characteristic value of a gaussian distribution.Recall that for periodic boundary conditions the kurtosis at K c takes the asymptotic value Γ(1/4) 4 /2π 2 = 2.1884 . .., see Refs.[7,8]

KURTOSIS AT AN EFFECTIVE CRITICAL POINT
That the kurtosis takes the asymptotical value 3 at K c does of course not mean that the kurtosis converges to 3 for every sequence of temperatures K c (L) that has K c as limit.We exemplify this by reexamining some of our data used in Ref. [10] where we relied on extremely detailed data on a wide temperature range for 4 ≤ L ≤ 20.For L = 4, 6, 8, 10 we also have magnetisation distributions.Let us say that K c (L) is the point where the variance of the modulus magnetisation, i.e. χ = N m 2 − |m| 2 , takes its maximum value.The distribution is here at its widest and on the threshold of breaking up into two parts, see Fig. 8 where we show a scaled distribution at K c (L) for L = 4, 6, 8, 10, in stark contrast to the distribution at K c of Fig. 9. Measuring the kurtosis at this point produces Fig. 10 which shows Q(K c (L), L) versus 1/L and a fitted 2nd degree polynomial which suggests the limit 1.520.The absence of error bars is due to the method by which the original data were produced.However, we expect the error to be smaller than the plotted points.In fact, repeating this exercise for periodic boundary conditions suggests the limit 1.517.There is a distinct possibility that these two limits are in fact the same, but that would require high-resolution data for larger systems to resolve than is at our disposal.In any case this subject falls outside the scope of this paper.

DISCUSSION AND CONCLUSIONS
As we have seen the sampled data for cubes up to side L = 160 agree well with the standard scaling picture for free boundary conditions, and e.g.recent long series expansions [14] also appears to favour this version, nontheless without a rigorous bound for the rate of convergence a simulation study is always open to the claim that it is dominated by finite size effects.
However, the last decade has seen a number of rigorous mathematical results on the behaviour of the randomcluster model which leads us to believe that the standard scaling is indeed the right one.The Fortuin-Kasteleyn   random-cluster model has a parameter q which governs the properties of the model.We'll refer the reader to [15] for more details and history.We recall that for q = 1 the model is the standard bond percolation model, for q = 2 it is equivalent to the Ising model, and in the limit q → 0 we get the uniform random spanning tree, or the uniform spanning forest, depending on the parameter p.
For the random spanning tree on the d-dimensional lattice with different boundary conditions Pemantle [16] begun a study which related it to the loop erased random walk and demonstrated a strong dependency on both the dimension and boundary condition.In later papers [17][18][19] these results were refined to show, among other things, that for large enough d that for two points in a grid with side L and free boundary the distance between the two points within the tree scale as L 2 , but that for the torus of side L the distance scales as L d/2 .
Coming to the case q = 1, Aizenman [20] studied the behaviour of the largest crossing clusters, i.e clusters which contain vertices on opposite sides of the box, in percolation on grids in different dimensions.Among other things he conjectured that at the critical point p c , for all large enough dimensions d the largest cluster in the case with free boundary should scale as L 4 , and for the torus, or cyclic boundary, it should scale as L 2d/3 .This conjecture was proven in [21,22], and so we know that for percolation the boundary condition have a large and non-vanishing effect on the size of the largest clusters at the critical point.It is also important to point out that the results from [21,22] are for the scaling exactly at the critical point p c .For other sequences of points converging to p c different scaling behaviours can appear.
The reason for the drastic difference between the torus and free boundary case here is that for large d the clusters in the model become much more spread out than in low dimensions.For d = 2 clusters at p c are always finite and have a boundary which in the scaling limit is very close to a brownian motion [23].For a finite box of side L, with free boundary the number of crossing clusters, has a finite mean, bounded as L grows, and the probability that there are more than k such clusters is less than exp(−a 1 k 2 ), for some positive constant a 1 [20].For a box with free boundary in d > 6, the critical dimension for q = 1, the number of crossing clusters is at least a 2 L d−6 , for some constant a 2 [20].Further, the largest cluster and a positive proportion of the crossing clusters have size proportional to L 4 , independently of d.The clusters are here of much lower dimension relative to that of the lattice and much more tree like in their structure than for d = 2.If we instead consider a torus with d > 6 the large number of these more expanding clusters lead them to connect up with other clusters, which would have been separate in the free boundary case, and this merging leads to maximum clusters that are vastly larger than in the free boundary case.
The current authors expect a picture similar to that for q = 1 to hold for the Ising case q = 2 as well.Since the susceptibility in the Ising model is proportional to the average cluster size in the random-cluster model this leads to a prediction of L 2 as the correct scaling for the free boundary case and L 5/2 for cyclic boundaries, for d = 5.This would also lead us to expect the largest clusters in the free boundary case to scale as L 2 .For our small to medium sized systems we collected the size of the cluster containing the central vertex of our cubes, a property which was also considered in [13] for truncated systems.
In Fig. 11 we show the normalised mean cluster size S 0 /L 2 for 3 ≤ L ≤ 55 as used in [13], together with a linear function estimated to be 0.3505(2) + 0.600 (1)x.The inset shows the zoomed-in version S 0 L −5/4 − 0.3505L 3/4 versus 1/L which then should essentially take the constant value 0.600, the red line.As we can see we have an excellent fit to the prediction that this cluster size should scale as L 2 .
Aizenman's prediction [20] of L 2d/3 as the correct scaling for percolation on the torus, for high d, came from a comparison with the Erdös-Renyi random graph on N vertices, for which the largest connected component at the critical probability, scales as N 2/3 .In [24] we follow this analogy further by comparing the largest cluster for the random-cluster model on 5-dimensional tori with the detailed rigorous results on the random-cluster model for complete graphs from [25], and a good agreement is found.
To conclude, we find that both the data from our simulations and the current mathematical result for the random-cluster model gives good support for the standard scaling picture for the Ising model with free boundary conditions, as well as a framework predicting further properties for the case with cyclic boundary.
we show |m| L 3/2 versus 1/L for 11 ≤ L ≤ 160 together with the curve 0.22958(6) + 1.101(3) x − 1.63(3) x 2 .We test the fit of this curve by zooming into the picture and instead show ( |m| L 3/2 − 0.22958)L which then should be well fitted by the line 1.101 − 1.63 x.As the inset of Fig 3 shows, it is and we conclude that to leading order |m| ∼ 0.22958(6)L −3/2 .However, the error bars for the largest systems are now quite pronounced.