The scaling window of the 5D Ising model with free boundary conditions

The five-dimensional Ising model with free boundary conditions has recently received a renewed interest in a debate concerning the finite-size scaling of the susceptibility near the critical temperature. We provide evidence in favour of the conventional scaling picture, where the susceptibility scales as $O(L^2)$ inside a critical scaling window of width $O(1/L^2)$. Our results are based on Monte Carlo data gathered on system sizes up to $L=79$ (ca. three billion spins) for a wide range of temperatures near the critical point. We analyse the magnetisation distribution, the susceptibility and also the scaling and distribution of the size of the Fortuin-Kasteleyn cluster containing the origin. The probability of this cluster reaching the boundary determines the correlation length, and its behaviour agrees with the mean field critical exponent $\delta=3$, that the scaling window has width $O(1/L^2)$.


I. INTRODUCTION
The Ising model in dimension d = 5 is of particular interest since it is the first case where the model is strictly above its upper critical dimension d c = 4. Rigorous results [1,2] establish that the critical exponents of the model assume their mean field values. Here the specific heat exponent α = 0, and the results of [2] also imply that the specific heat is bounded at the critical point. Recent simulation results [3] also indicate that, just as for the mean field case, the specific heat is discontinuous at the critical point.
In contrast to these asymptotic results there has been a long running debate over the finite size scaling for the model with free boundary conditions. We'll refer the reader to [4] for a fuller overview of the history and stick to the presently most relevant parts here. For d = 5 and cyclic boundary conditions there is agreement that e.g. χ ∝ L 5/2 for a lattice of side L. The conventional picture for the free boundary case is that χ ∝ L 2 . However, it has also been suggested [4] that the free boundary case should scale in the same way as the cyclic boundary case near the finite size susceptibility maximum. A computational study [5] of the, then, largest lattices possible supported the conventional picture, but in [4] it was suggested this was due to an underestimate of the influence of the large boundaries of the used systems. For systems exactly at the critical coupling this issue was resolved in [6] where a study of systems up to L = 160 demonstrated an increasingly good agreement with the conventional picture as the system size was increased. But, this left the behaviour in the rest of the critical scaling window open.
The aim of the current paper is to extend the study of large systems from [6] to the full critical window, including the location of the maximum of χ, and give the best possible estimates for the scaling behaviour in the coupling region discussed by all the previous papers. Apart from the susceptibility we also study properties of the Fortuin-Kasteleyn cluster containing the origin and use those to estimate both the susceptibility and the correlation length of the model.
To concretize, the predictions from [4] are that the location of the maximum for χ will scale as L −2 and the maximum value as L 5/2 . The more recent [7] agrees with these predictions, and also the prediction from [8] that the location of the maximum of the susceptibility should scale as L 2 , as does [9], but both are based on smaller system sizes than those considered in the present work.
In short, our conclusion is that the data is well fitted by the conventional scaling picture, both for the location and value of the susceptibility, and location of the finite size critical point for the magnetization.

II. DEFINITIONS AND DETAILS
For a given graph G on N vertices 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. Here the graph G is a 5-dimensional grid graph of linear order L with free boundary conditions, i.e. a cartesian product of 5 paths on L vertices, so that the number of vertices is N = L 5 and the number of edges is 5L 5 (1 − 1/L). We use K = 1/k B T as the dimensionless inverse temperature (coupling) and denote the thermal equilibrium mean by · · · . The critical coupling K c was recently estimated by us [3] to K c = 0.11391498(2). We will define a rescaled coupling as κ = L 2 (K −K c ) which gives a scaling window of width O(1/L 2 ). The standard definitions apply; the magnetisation is M = i S i (summing over the vertices i) and the energy is E = ij S i S j (summing over the edges ij). We let m = M/N , U = E/N and U = U . The susceptibility is χ = M 2 /N while we define the modulus susceptibility asχ = var(|M |) /N . The standard deviation is denoted σ, as is customary. We will refer to the point where the distribution of M changes from unimodal to bimodal as K * c (L), or, in its rescaled form, κ * c (L). Recall also that thermodynamic derivatives of e.g. log χ can be obtained through correlation mea- Let S denote the average size of a flipped cluster and recall [10] that χ = S . We use the subscript o to denote the origin (i.e. centre) vertex of the grid G so that S o is the average size of a cluster containing the origin. The event that a cluster containing the origin also reaches the boundary ∂G of the graph is denoted o ↔ ∂G.
We have collected data using Wolff-cluster spin updating for L = 15, 19, 23, 31, 39, 47, 55, 63, 79 for a wide and dense set of κ-values in 0 ≤ κ ≤ 0.92. The data include e.g. magnetisation, energy (L ≤ 63), size and radius of the origin cluster. The total number of measurements range from about 300000 for the small systems (L ≤ 39) down to 100000 for L = 63 and 50000 for L = 79.

III. THE WIDTH OF THE SCALING WINDOW AND THE SHIFT EXPONENTS
We begin by establishing the natural width of the scaling window, i.e. the form of a natural rescaled temperature. We have already defined the rescaled coupling as κ = L 2 (K − K c ) with the expectation that this will be the relevant scaling. In principle, the location of different effective finite size critical points could scale with different exponents and we will compare a number of such possibilities.
We used a 2nd order interpolation of the data points to estimate the location of the maximum for three different parameters:χ, ∂ log |m| /∂K and ∂ log χ/∂K. For the last two we only have these data for L ≤ 63. The locations are plotted versus 1/L in Figure 1. The values clearly point to three different limit κ-values. Based on L ≥ 19 we estimate the limit values and an error bar by removing each of the data points in turn and fitting a line to the remaining points. Thus we find κ c (L) = 0.8763(5) − 1.21(2)/L for theχ-maximum, κ c (L) = 0.8742(5) − 1.49(1)/L for the maximum ∂ |m| /∂K and 0.8725(3) − 1.50(1)/L for the maximum ∂ log χ/∂K. We conclude that our rescaled coupling κ is the correct scaling, and each of the effective critical points scale with the same exponent.

IV. THE LOW AND HIGH-TEMPERATURE REGIONS
Next we estimate the location of the point where the magnetisation distribution switches from unimodal to bimodal, denoted κ * c (L). This point will be seen as the effective finite size separator between the high and lowtemperature regions of the model.
To do this we note that near this point the standard- . The lower inset of Figure 2 shows the distribution of M/σ for L = 23 at κ = 0.76 (K = 0.11535165) together with the fitted formula. The binning of the data was done with Mathematica's built-in procedures but the results are not very sensitive to this choice. At this κ (i.e. for L = 23) the distribution is very close to shifting from unimodal to bimodal, as is indicated by the very small coefficient c 2 = −0.002735. The distribution has kurtosis 2.2217 and when c 2 = 0 the kurtosis becomes Γ(1/4) 4 /8π 2 ≈ 2.1884. In fact, it works equally well to simply measure the kurtosis and then estimate through interpolation at which point the kurtosis is 2.1884. All this suggests that the distribution's shape is well captured by f (x). See also [3] where we apply this formula to 5D grids with periodic boundary conditions.
The upper inset of Figure 2 shows the measured coefficient c 2 plotted against κ for a range of L. The value of c 4 (not shown) at κ * c (L) is somewhat noisy but it seems to converge to −0.114 (3). For a cyclic boundary this value is close to −0.602 [11]. Figure 2 shows κ * c (L) (where c 2 = 0) versus 1/L. The fitted 2nd degree polynomial suggests κ * c (L) = 0.8561(6) − 2.52(3)/L + 7.6(4)/L 2 , with error bars obtained as above, so that κ * c ≈ 0.856. This point then represents where the true low-temperature region begins.
The point κ * c (L) constitutes an effective critical temperature but its scaling rule appears to depend strongly on the boundary conditions. Note that for periodic boundary conditions the corresponding point K * c (L) converges to K c with rate 1/L 3 [12], i.e. faster than the width of the scaling window, which is 1/L 5/2 in that case. For free boundary conditions we instead see that K * c (L) → K c with the same rate as the width of the scaling window, i.e. 1/L 2 .
We note that in [7] an estimate based in the kurtosis was used, but instead of using the universal kurtosisvalue given at the point where c 2 = 0 they hand picked a value for the kurtosis and then used the fact that all fixed such values should give the same scaling, with larger or smaller corrections to the scaling.

V. SCALING OF SUSCEPTIBILITY
Let us now return to the matter of the finite-size scaling of χ andχ. With the onset of the low-temperature region at (asymptotically) κ * c = 0.856 we look at the scaling behaviour of χ at fixed κ below and above this point. Consider Figure 3 where we show χ/L 2 versus L for six different fixed κ-values. Lines are fitted on L ≥ 31 to demonstrate the presence of correction terms. The behaviour is clearly linear for larger L, with the possible exception of κ = 0.86, i.e. above κ * c = 0.856 where we have entered the low-temperature region and have a bimodal magnetization, but even here the correction term is still rather weak. However, higher κ-values render us stronger corrections. The normalised susceptibility was also used in [7], but there the estimates, which lead to scaling proportional to L 5/2 were based on much smaller lattices L ≤ 36, where finite size effects are much stronger.
Moving on to the modulus susceptibilityχ/L 2 the corrections become clear and present, even for κ ≤ 0.856, but never more than can be captured by a simple 2nd degree polynomial. Figure 4 shows this effect. Note also that once we have gone beyond the maximum at κ = 0.876 the behavior quickly becomes linear again as is clearly seen at κ = 0.92, which is where our data end. Combining all the measuredχ/L 2 with their asymptotic  values, based on 2nd degree polynomials fitted on L ≥ 31, we now obtain Figure 5.
For this range of lattice sizes and at least some values of κ one can make a reasonable fit to the data with a L 5/2 scaling as well, by including sufficient correction terms. But, as we can see in Figure 5, the L 2 scaling gives values that decrease to their limit value for a large range of κ, possibly all κ below κ * c , making an even larger power of L seem unlikely here.

VI. THE FK-CLUSTER AT THE ORIGIN
Our next object of study is the Fortuin-Kasteleyn cluster containing the origin. This is the cluster which the Wolff random cluster algorithm builds when started at the central vertex of our lattice, and is one of the clusters in the FK-random cluster representation of the Ising model. From the FK-model the usual susceptibility can be computed as the expected cluster size, and with free boundary the expected size of the origin cluster will be larger than the average cluster size, thereby giving an upper bound on the susceptibility. This observation was used in [6] to give a susceptibility estimate which was hoped to be less affected by the effects of the free boundary, in part as a response to the truncation method of [4].
As described we expect S o to place an upper bound on χ, but we do not expect them to be of different scaling orders. Plotting S o /L 2 versus 1/L demonstrates a similar behaviour to that of χ/L 2 . In Figure 6 we show the scaling at a number of different κ-values below κ * c . In fact, a strong correction enters the picture already at κ = 0.85 (not shown in the figure). This is to be expected since above κ * c there are additional terms from the FKmodel affecting the Ising susceptibility, and those are not included in our sampled data.
The distribution of the origin cluster size S o is also an interesting object. It is a Pareto-like distribution with density proportional to x −1−1/δ , or, by definition [13], where δ is the critical exponent reflecting how an external field affects the magnetisation at the critical temperature,  here having the mean field value δ = 3. In Figure 7 we show a log-log plot of the distribution density (for small x = S o ) together with a line having slope −4/3, consistent with δ = 3. The line is clearly an excellent fit. The plot shows the distribution for L = 63 at κ = 0.85 but for this range of cluster sizes the plot is indistinguishable from that of other κ-values, and indeed of other L (except the smallest L).

VII. CORRELATION LENGTH
The origin cluster can also be used to estimate the correlation length of the model. In the high temperature region the correlation length ξ is given, see [13], by the limit For a fixed κ we expect not just ξ → ∞, but rather that the correlation length is comparable to L, i.e. ξ/L → a(κ) for some function a. Thus we propose to estimate a(κ) by first defining ξ/L = −1/ log Pr(o ↔ ∂G) and then use a simple projection rule to estimate the limit function a(κ).
The inset of Figure 8 shows ξ/L versus 1/L for a range of L at different κ-values. Figure 8 also shows a rough estimate of a(κ) based upon fitted 2nd degree polynomials (on L ≥ 15) for each κ. We expect the errors in the plot to be significant although the general behaviour as such is correct. In the low-temperature region the identity in Equation 2 is no longer valid, and the probability that the origin cluster reaches the boundary tends to 1, and this agrees well with the steep growth beginning at κ ≈ 0.85, the start of the effective finite size low-temperature region.

VIII. CONCLUSIONS
As we have seen here our data fits excellently with the conventional scaling picture for the free boundary case. Now, apart from this fact we believe that there are good reasons for expecting distinct behaviour from the free and cyclic boundary cases. As is well known the Ising model is equivalent to the Fortuin-Kasteleyn random cluster model for q = 2. The case q = 1 corresponds to ordinary percolation and q → 0 gives the random spanning tree for the lattice.
For the random spanning tree Pemantle [14] proved that for large enough d the distance between two points in the tree scales as L 2 with free boundary and as L d/2 with cyclic boundary. For percolation, Aizenmann [15] conjectured that for large enough d the size of the largest cluster should scale as L 4 with free boundary and L 2d/3 with cyclic boundary. This conjecture was later proven in [16,17]. So, for both q = 0 and q = 1 we have rigorous results showing that free and cyclic boundary conditions lead to distinct scalings.
Aizenman's conjectured scaling L 2d/3 came from the fact that the cyclic boundary case was expected to behave like the critical Erdős-Renyi random graph, where the maximum cluster has size N 2/3 , where N is the number of vertices, which would be L d for the lattice. Again, there are detailed rigorous results [18] concerning the random cluster model on the complete graph and in [19] it was demonstrated that Monte Carlo data for the FKmodel with q = 2 and d = 5 with cyclic boundary gives a detailed agreement with the scaling for the complete graph. In particular, the size of the largest cluster has the same scaling for both graphs for several different rescaled coupling ranges. Now, the L 2 scaling observed in this paper for the origin cluster can be seen as an indicator of the correct scaling of the largest cluster as well. This would give a scaling which is different for the free and cyclic boundary cases, just as one would expect from the known results for lower q. In [5] the authors conjectured that for every q there is a d(q) such that above this dimension the FK-model with cyclic boundary behaves like the model on a complete graph, giving a partial generalisation of Aizenmann's conjecture for all q. In fact, we also expect the model with free boundary to follow a simpler meanfield scaling, similar to that seen on an infinite tree, or Bethe lattice. The fact that there are several distinct basic model systems, e.g. the complete graphs and the Bethe lattices, which on one hand have mean field critical exponents but on the other hand differ on more detailed properties is probably under-appreciated in the older literature.