P H Lundow

A Rosengren

On the p, q-binomial distribution and the Ising model

Keywords: Ising models, magnetic phase transition, theoretical, numerical modelling Ising, Magnetisation, Distribution

We employ p, q-binomial coefficients, a generalisation of the binomial coefficients, to describe the magnetisation distributions of the Ising model. For the complete graph this distribution corresponds exactly to the limit case p = q. We take our investigation to the simple d-dimensional lattices for d = 1, 2, 3, 4, 5 and fit p, q-binomial distributions to our data, some of which are exact but most are sampled. For d = 1 and d = 5 the magnetisation distributions are remarkably well-fitted by p, q-binomial distributions. For d = 4 we are only slightly less successful, while for d = 2, 3 we see some deviations (with exceptions!) between the p, q-binomial and the Ising distribution. However, at certain temperatures near T c the statistical moments of the fitted distribution agree with the moments of the sampled data within the precision of sampling. We begin the paper by giving results on the behaviour of the p, q-distribution and its moment growth exponents given a certain parameterization of p, q. Since the moment exponents are known for the Ising model (or at least approximately for d = 3) we can predict how p, q should behave and compare this to our measured p, q. The results speak in favour of the p, q-binomial distribution's correctness regarding their general behaviour in comparison to the Ising model. The full extent to which they correctly model the Ising distribution is not settled though.

The Lenz-Ising model is arguably the most studied model in theoretical physics. Though solved completely in one dimension [START_REF] Ising | Beitrag zur Theorie des Ferromagnetismus[END_REF] and for the special case without an external field in two dimensions [START_REF] Onsager | Crystal statistics I. A two-dimensional model with an order-disorder transition[END_REF], the three dimensional case has so far defied all attempts at an exact solution. Most studies have been devoted to the free energy, or rather, its derivatives with respect to either the temperature or the field. Series expansions and Monte Carlo studies in combination with finite-size scaling techniques have been quite successful methods for pinpointing the critical temperature or the critical exponents. This then has provided us with good knowledge of how the magnetisation cumulants, such as the susceptibility, scale with system size. There are also good estimates of shape parameters of the magnetisation distribution. One of these, the Binder ratio [START_REF] Binder | Finite size scaling analysis of Ising model block distribution functions[END_REF], gives us a measure of the kurtosis, or peakedness, of the distribution. However, less attention has been paid to the underlying magnetisation distribution. In this paper we focus on this distribution with the hope that our investigation may provide a new lead on the exact solution of the Ising model. We will take a very general perspective at first and only later concentrate on the distribution for the standard lattice systems.

Choose a graph, e.g. a square lattice, on n vertices and compute its Ising partition function Z, keeping track of its terms according to their magnetisation, so that Z = Z 0 +Z 1 +• • •+Z n . The quotient Z k /Z is the probability of having k negative spins, or magnetisation M = n -2 k. The controlling parameter of the partition function is the temperature, that is, for any given temperature the partition function provides us with a distribution of magnetisations. At infinite temperature this is simply the binomial distribution. At zero temperature, on the other hand, we receive a distribution with two peaks, one at k = 0 and one at k = n, both with 50% of the probability mass. The distribution is always symmetrical. What happens between these two extreme temperatures? Though there are exceptions to the rule, for most graphs the distribution begins its life at high temperatures as a unimodal distribution with the peak at the middle k = n/2. As we lower the temperature the distribution gets increasingly wider until we reach a temperature where the distribution changes from unimodal to bimodal. Near this temperature, slightly above and slightly below, the distribution is particularly wide. Lowering the temperature even further the distribution develops two sharp peaks, both essentially gaussian [START_REF] Martin-Löf | Mixing properties, differentiability of the free energy and the central limit theorem for a pure phase in the Ising model at low temperature[END_REF], and these peaks move outwards.

Early studies of the magnetisation distribution include [START_REF] Bruce | Probability density functions for collective coordinates in Ising-like systems[END_REF] and [START_REF] Binder | Finite size scaling analysis of Ising model block distribution functions[END_REF]. Distribution moment ratios were studied in for example [START_REF] Salas | Universal amplitude ratios in the critical twodimensional Ising model on a torus[END_REF] for 2-dimensional lattices. In [START_REF] Martins | Probability distribution of the order parameter[END_REF] the distribution was used to estimate the critical temperature of the model and in [START_REF] Kaneda | Effects of shape and boundary conditions on finite-size scaling functions for anisotropic threedimensional Ising systems[END_REF] they studied the effects of anisotropy on the distribution for the 3D-lattice. In the 5-dimensional case renormalisation group techniques have been used to give expressions for the distribution, see [START_REF] Brezin | Finite size effects in phase transitions[END_REF] and [START_REF] Mon | Finite-size scaling of the 5D Ising model[END_REF]. We should also mention [START_REF] Häggkvist | The Ising partition function for 2D grids with cyclic boundary: computation and analysis[END_REF] and [START_REF] Häggkvist | On the Ising model for the simple cubic lattice[END_REF] which include studies of the magnetisation distribution for the 2-dimensional and 3-dimensional case respectively. The distribution has also been studied for higher spins in two dimensions, see [START_REF] Plascak | Universality and double critical end points[END_REF]. Even for 1-dimensional lattices, for which the scenario described above does not hold, the distribution is complex enough to merit interest, see [START_REF] Garcia-Pelayo | Distribution of magnetization in the finite ising chain[END_REF]. The magnetisation distribution was recently studied for a quantum model in [START_REF] Lamacraft | Order parameter statistics in the critical quantum ising chain[END_REF].

Since the magnetisation distribution is binomial at infinite temperature, is it then possible to describe this distribution at all temperatures with a generalisation of the binomial distribution? The study of q-series and q-analogs of the classical special functions has popped up in such diverse quarters as Lie algebras, statistical mechanics, trancendental number theory and computer algebra [START_REF] Andrews | q-series: their development and application in analysis, number theory, combinatorics, physics, and computer algebra[END_REF]. We, in this article, try to model the distributions using p, q-binomial coefficients. In one case, the complete graph, they model exactly the distributions and for d-dimensional lattices the dimension seems to determine how well they fit to the Ising distributions. Especially in the case of d = 1 and d = 5 the Ising distributions are particularly well-fitted by the p, q-binomial coefficients.

These coefficients, controlled by only two parameters once the size of the system is set, constitute a very flexible family of distributions with a clear relationship to the Ising distribution under very general settings. It would of course be very beneficial to our understanding of the Ising model if we knew how close this relationship is. As it turns out, the p, q-binomial distribution, interesting enough in itself to merit a study, can be analysed theoretically in a way the Ising distribution so far has resisted. We believe this study will provide ample evidence of their belonging to the same family of distributions.

The problem now turns into that of determining how the parameters p, q depend on the temperature. We estimate this dependence for very high temperatures for d = 1. For d > 1 we give theory-based predictions of how the parameters should depend on n near T c to receive the correct finitesize scaling of the moments. This prediction even includes the logarithmic correction of the susceptibility for d = 4, in agreement with renormalisation group predictions.

The paper begins in section 2 by providing the basic definitions, such as the q-Pochhammer symbol, q-binomial coefficients and finally the p, qbinomial coefficients. Some results on their properties are also stated, intended to give the reader a feel for how p, q-binomials behave. In section 3 we define the p, q-binomial distribution and give an algorithm for finding values of p and q when the distribution is given as input. The problem is to determine an optimal choice for p and q. As it turns out we only have to focus on the value of the probability and the location of the distribution's peaks, at least for a bimodal distribution. The unimodal distribution always has its peak at the middle so in this case we instead take the quotient be- tween the two middle probabilities as controlling parameter. This quotient is unfortunately rather sensitive to noise, making it difficult to determine the parameters p and q for sampled data.

In section 4 we provide some useful tools for working with p, q-binomial coefficients. After this we are, at long last, ready to give some general results on the distribution of p, q-binomial coefficients in section 5. Using the parameterization p = 1 + y/n and q = 1 + z/n we find the asymptotic value of y, given z, where the distribution is flat in the middle. We also allow for a small change in y, using a higher order parameter a, so that we can follow properly how the distribution changes from unimodal to bimodal. However, the computations that we rely on involve some rather complicated series expansions that were made using Mathematica. These are much too long to fit into this paper. We have prepared a simplified Mathematica notebook that performs all the necessary computations. The interested reader can obtain it by contacting the first author.

In section 6 we give exact scaling formulae for the moments of the distributions depending on the parameters a and z. For a given moment of these distributions we always obtain the same exponent on n, regardless of a and z. In section 7 we try to remedy this by letting the previously fixed parameter z depend ever so slightly (at most logarithmically) on n. We can now adjust the exponent of n though this comes at the cost of an extremely slow convergence.

Section 8 defines the Ising model, laying the ground for studying distributions of magnetisations, the intended application of our endeavour. In section 9 we apply our tools to the d-dimensional lattice graphs for d = 1, 2, 3, 4, 5 fitting p, q-distributions to simulated distributions and comparing them. With statistical certainty we find that at least for d = 4, 5 the magnetisation distribution is correctly modelled by a p, q-binomial distribution.

The Appendix A gives detailed results for the special case when p = q. This case corresponds exactly to the complete graph and is the only case where we can give asymptotically exact expressions for the sum of the coefficients.

A condensed reading, more suitable to the reader who is pressed for time, should include a look at (1), ( 5), [START_REF] Corcino | On p, q-binomial coefficients[END_REF], [START_REF] Garcia-Pelayo | Distribution of magnetization in the finite ising chain[END_REF], ( 12), ( 28), ( 29), [START_REF] Rosengren | On the combinatorial solution of the Ising model[END_REF] for the necessary definitions and results concerning the basics. After that, the most important results are stated in equations (39), (45), (68), (79) and (82). After looking up the basic definitions regarding the Ising model the reader can skip to (91). In Section 9 the reader can now pick and choose his favourite lattice and look at the pictures. The q-binomial coefficient

n k q = k i=1 1 -q n-i+1 1 -q i , q = 1, 0 ≤ k ≤ n (1)
is a natural extension of the standard binomial coefficient. Note that (1) can be viewed as a formal polynomial in q of degree k (n -k) where the coefficient of q j counts the number of k-subsets of {1, . . . , n} with element sum j + k (k + 1)/2. It is thus easy to see that lim

q→1 n k q = n k (2)
The Pochhammer symbol (a) n , or shifted factorial, also has a q-deformed relative, the q-Pochhammer symbol, defined as

(a; q) n = n-1 i=0 1 -a q i , n ≥ 0 (3)
The q-binomial coefficient can be expressed as n k q = (q; q) n (q; q) k (q; q) n-k = q n-k+1 ; q k (q; q) k (4)

The q-binomials and the q-Pochhammer function have many interesting properties and we point the interested reader to the book [START_REF] Andrews | Special functions, vol. 71 of Encyclopedia of mathematics and its applications[END_REF] and especially the charming little book [START_REF] Andrews | q-series: their development and application in analysis, number theory, combinatorics, physics, and computer algebra[END_REF]. The p, q-binomial coefficient was defined in [START_REF] Corcino | On p, q-binomial coefficients[END_REF] as

n k p,q = k i=1 p n-i+1 -q n-i+1 p i -q i , p = q, 0 ≤ k ≤ n (5)
Clearly, in the case p = 1 this reduces to a q-binomial coefficient. Also, note that p and q are interchangable and that the coefficient takes the same value at k and n -k. It is an easy exercise to show the following identity and we leave this to the reader.

n k p,q = p k (n-k) n k q/p = q k (n-k) n k p/q (6)
Combining ( 6) with (2) we let the p, q-binomial coefficient be defined for p = q as n k q,q = q k (n-k) n k An identity involving q-binomial coefficients can be extended to a p, qbinomial identity by first replacing q with q/p and then using the identity [START_REF] Brezin | Finite size effects in phase transitions[END_REF]. As an example we consider the q-Vandermonde identity, see e.g. [START_REF] Corcino | On p, q-binomial coefficients[END_REF], which can be stated as

m + n k q = k ℓ=0 m k -ℓ q n ℓ q q ℓ (m-k+ℓ) (8)
and is an extension of the Chu-Vandermonde identity for binomial coefficients. Using (6) above we can now obtain a p, q-analog of this.

Theorem 2.1.

m + n k p,q = k ℓ=0 m k -ℓ p,q n ℓ p,q p (n-ℓ) (k-ℓ) q ℓ (m-k+ℓ) (9)
Proof. In the q-Vandermonde identity (8), replace q with q/p and multiply both sides with p k (m+n-k) . Using [START_REF] Brezin | Finite size effects in phase transitions[END_REF] the left hand side is now a pure p, qbinomial coefficient. The ℓth term of the right hand side becomes

m k -ℓ q/p n ℓ q/p p k (m+n-k) (q/p) ℓ (m-k+ℓ) = p (k-ℓ) (m-(k-ℓ)) m k -ℓ q/p p ℓ (n-ℓ) n ℓ q/p p (n-ℓ) (k-ℓ) q ℓ (m-k+ℓ) = m k -ℓ p,q n ℓ p,q p (n-ℓ) (k-ℓ) q ℓ (m-k+ℓ)
and the theorem follows.

The q-binomial coefficients have been shown to form a log-concave (and thus unimodal) sequence for q ≥ 0, see e.g. [START_REF] Butler | The q-log-concavity of q-binomial coefficients[END_REF] and [START_REF] Krattenthaler | On the q-log-concavity of Gaussian binomial coefficients[END_REF]. However, for the p, q-binomial coefficients this does not always hold. Rewriting them as a product like in [START_REF] Brezin | Finite size effects in phase transitions[END_REF] we have in fact a product of two sequences; that of p k (n-k) and n k q/p for k = 0, . . . , n. The first sequence is log-concave for p ≥ 1 and log-convex for p ≤ 1. It is well-known that the element-wise product of two log-concave positive sequences is also log-concave. So, if p ≥ 1 and q ≥ 0 then the sequence of p, q-binomial coefficients is log-concave. We conjecture that for p, q > 0 the sequence can be either unimodal, with the maximum at k = ⌊n/2⌋, or bimodal, with the maxima at k and n -k for some 0 ≤ k ≤ n/2, but not trimodal etc. We will assume this to be true in this paper but a formal proof is still lacking. It is, however, fairly easy to show in the special case p = q. Note that if we allow negative values of p the sequence can have a local maximum at every alternate index k.

(n) ∼ g(n) to denote that f (n)/g(n) → 1 as n → ∞. Anal- ogously f (n) ∝ g(n) denotes that f (n)/g(n) → A,
for some non-zero real number A, as n → ∞. Also, henceforth we will assume that n is even to simplify some calculations.

3 The p, q-binomial distribution Now we are ready to introduce the notation

Ψ p,q (n) = n k=0 n k p,q (10)
and define the p, q-binomial probability function

P p,q (n, k) = n k p,q Ψ p,q (n) (11)
The reader should here observe that the sum of the coefficients has, to the best of our knowledge, no simpler expression in the general case. Neither do the sum of the q-binomial coefficients have a simpler expression that we are aware of. Compare this with the case of the standard binomial coefficients for which the sum is simply 2 n . Having made our assumption of unimodality/bimodality we can now set up a simple computational scheme to find values of p given q. First we need to define a highly useful quantity; the ratio between two coefficients

R p,q (n, k, ℓ) = n k -ℓ p,q n k p,q (12)
In the special case when ℓ = 1 we are looking at two consecutive coefficients. The ratio then becomes

R p,q (n, k, 1) = p k -q k p n-k+1 -q n-k+1 = p -(n-2 k+1) 1 -(q/p) k 1 -(q/p) n-k+1 = q -(n-2 k+1) 1 -(p/q) k 1 -(p/q) n-k+1 (13)
Due to symmetry we assume now that 1 ≤ k ≤ n/2. Given q and a ratio r = R (n, k, 1) we wish to find the p, assumed to be larger than q, such that p and q satisfies (13). This is done through the iteration scheme p ← q k -r q n-k+1 + r p n-k+1 which is obtained from setting (13) to r. Choosing e.g. p = 1, or a number slightly larger than q, as start value works well for all practical purposes. To prove that this iteration scheme actually converges one would have to show that the derivative with respect to p is at most 1 using the start value.

We give no such proof but will leave it at that and just claim that this is a practical method. Given a p, q-binomial distribution, or any distribution that we wish to approximate by a p, q-binomial distribution, can we find the pair p, q > 0 that generated it such that the distribution of p, q-binomial coefficients have the correct probability and ratio r at coefficient k? A simple bisection procedure solves this problem practically under the assumption p > q. Suppose the input distribution has the probabilitites P (0) , P (1) . . . , P (n). Let k, P (k), r = P (k -1) /P (k) and an ǫ be given as input parameters.

Algorithm p, q-Find 1. Assign q min ← 0 and

q max ← k r (n-k+1) 1 n-2 k+1 .
2. Assign q ← (q min + q max )/2 3. Compute p using (14). 4. Compute P p,q (n, k) using [START_REF] Garcia-Pelayo | Distribution of magnetization in the finite ising chain[END_REF].

5. If P p,q (n, k) < P (k) then q min ← q, otherwise q max ← q.

6. If q max -q min < ǫ then exit loop, otherwise jump to step 2.

The initial value of q max in step 1 is due to that the ratio r is bounded by

0 < r < q -(n-2 k+1) k n -k + 1 (15)
which is a consequence of

0 < 1 -x k 1 -x n-k+1 < k n -k + 1 (16)
for x > 1 and 1 ≤ k ≤ n/2. The p, q-Find procedure seems to work best when P (k) is one of the maximum probabilities. However, if the distribution is unimodal so that the maximum probability is at k = n/2, then the scheme will depend heavily on the quality of r. On the other hand, if the distribution is bimodal then this problem goes away and we may simply set r = 1, unless n is too small. It is implied here, though we do not have a proof, that increasing q while keeping k and r fixed also increases the probability P p,q (n, k). It actually increases until p = q which then constitutes an interesting limit case, which we will deal with in the appendix A. In this section we set up bounds for the q-Pochhammer function and use them for giving bounds of quotients between q-binomial coefficients. Basically we mimic the upper bound in [START_REF] Kirousis | Upper bounds adn asymptotics for the q-binomial coefficients[END_REF] and [START_REF] Kaporis | The unsatisfiability threshold revisited[END_REF] but we extend it to obtain a lower bound as well. They are very useful bounds so we will do this in some detail though everything is based on standard elementary methods. First we need the integral estimate of a sum. Let f (x) be a continuous, positive, decreasing function on the interval m ≤ x ≤ n + 1 where m and n are integers. Then

n+1 m f (x) dx ≤ n k=m f (k) ≤ f (m) + n m f (x) dx (17)
Recall that the dilogarithm is defined as

Li (x) = ∞ n=1 x n n 2 = - x 0 log(1 -t) t dt (18)
Now let 0 < a < 1 and 0 < q < 1 and note that

-log (a; q) n = n-1 k=0 -log 1 -a q k (19)
Note also thatlog (1 -a q x ) is a positive and decreasing function for x ≥ 0. Take the series expansion

-log(1 -x) = ∞ k=1 x k k (20) so that -log(1 -a q x ) = ∞ k=1 (a q x ) k k (21)
Integration gives

u 0 -log(1 -a q x ) dx = ∞ k=1 u 0 (a q x ) k k dx = ∞ k=1 a k k q k x k log q u 0 = (22) 1 log q ∞ k=1 (a q u ) k k 2 - ∞ k=1 a k k 2 = Li (a q u ) -Li (a) log q (23)
Together with the integral estimates above we have Li (a q n ) -Li (a) log q ≤log (a; q) n ≤log(1 -a) + Li (a q n-1 ) -Li (a) log q (24) Reversing the signs and taking exponentials we finally obtain (a; q) n ≥ (1 -a) exp Li (a) -Li (a q n-1 ) log q (25)

(a; q) n ≤ exp Li (a) -Li (a q n ) log q (26)
Now we turn to the q-binomial coefficients. Let 0 ≤ ℓ ≤ k ≤ n/2 and use (4) to note that

R q (n, k, ℓ) = q k-ℓ+1 ; q ℓ (q n-k+1 ; q) ℓ (27)
Using the bounds for the q-Pochhammer function we can now bound the ratio. For the upper bound of the ratio we take the quotient of the upper bound and the lower bound. The ratio R q (n, k, ℓ) then has the upper and lower bounds

R q (n, k, ℓ) ≤ exp Li(q k-ℓ+1 )+Li(q n-k+ℓ )-Li(q k+1 )-Li(q n-k+1 ) log q 1 -q n-k+1 (28)
R q (n, k, ℓ) ≥ 1 -q k-ℓ+1 exp Li(q k-ℓ+1 )+Li(q n-k+ℓ+1 )-Li(q k )-Li(q n-k+1 ) log q (29) Using [START_REF] Brezin | Finite size effects in phase transitions[END_REF] we can now obtain bounds for the quotients of p, q-binomial coefficients. Note simply that

R p,q (n, k, ℓ) = R q/p (n, k, ℓ) p ℓ (n-2 k+ℓ) (30)
and use the bounds from ( 28) and [START_REF] Plascak | Universality and double critical end points[END_REF].

5 Controlling the p, q-binomial distribution

Choosing a k for a given n such that R p,q (n, k, 1) = 1 defines a set of pairs p, q; an isocurve. If we choose a value of q, then what value should p have to result in a distribution of coefficients which has a peak at k (and n -k), i.e. with R p,q (n, k, 1) = 1? The iterative method (14) above produces the correct p for any given q but reveals no information on p. To obtain this we need to parameterize q properly and one way to do this is to set q = 1 + z/n for some z ≤ 0. This parameterization was also used in [START_REF] Kirousis | Upper bounds adn asymptotics for the q-binomial coefficients[END_REF] and [START_REF] Kaporis | The unsatisfiability threshold revisited[END_REF] for computing upper bounds on q-binomial coefficients.

Wide and flat distributions

Let us begin with the particular distribution which has its peak in the middle, i.e. k = n/2. Due to symmetry we are interested only in the case p > q.

Recall that in the special case p = q the distribution has its peak at n/2

when q = n n + 2 = 1 - 2 n + 4 n 2 - 8 n 3 + • • • (31)
which would correspond to z = -2. Assuming that p has the expansion

p = 1 + y 1 n + y 2 n 2 + • • • (32)
we shall now determine y 1 , y 2 , . . .. According to (13), R p,q (n, n/2, 1) = 1 is equivalent to p n/2+1 -p n/2 = q n/2+1 -q n/2 (33)

Setting p = 1 + y 1 /n + y 2 /n 2 + • • • and q = 1 + z/n and performing a series expansion using (126), we find the right hand side to be

z e z/2 n - z 3 e z/2 4 n 2 + z 4 (3 z + 16) e z/2 96 n 3 + • • • (34)
and the first two terms of the left hand side are

y 1 e y 1 /2 n - y 3 1 -4 y 2 -2 y 1 y 2 e y 1 /2 4 n 2 + • • • (35)
We solve this term by term. The first equation

y 1 e y 1 /2 = z e z/2 (36)
has the solution y 1 = 2 w where

w = W z 2 e z/2 (37)
Here s = W (x) is the Lambert function solving s e s = x. Note that for z ≥ -2 we have y 1 = z but this is not the case for z < -2 where we have y 1 > z. It is easy to solve the corresponding equations for y 2 , y 3 , . . . but for the case in hand we actually only need y 1 . We will henceforth drop the subscript and refer to it as simply y.

What shape does the distribution have at this particular p and q? What we are seeking is an expression for R p,q (n, n/2, ℓ) and this is where we start using ( 28), ( 29) and [START_REF] Rosengren | On the combinatorial solution of the Ising model[END_REF]. We define ℓ = x n 3/4 , just as we did in the case of p = q. Let also where w is defined as in (37). We are now ready to compute the limits of the upper and lower bound of (30) by using ( 28) and [START_REF] Plascak | Universality and double critical end points[END_REF]. It turns out that the limits of these bounds coincide and we receive

q p = 1 + z n 1 + 2 w n + • • • = 1 + z -2 w n + • • • (38
R p,q n, n/2, x n 3/4 ∼ exp w z (2 w + z) x 4 6 ( 39
)
Note that the special case z = -2 corresponds to p = q and gives the coefficient -4/3 of x 4 . The calculations were performed with Mathematica and are much to unwieldy to fit in this paper. We have prepared a Mathematica notebook that performs the calculations step by step, using some practical transformation rules. What we compute is actually the limit of the logarithm of the upper and lower bound. The steps are as follows; compute series expansions of the different powers of q/p, use them inside the dilogarithms and then compute their series expansions, add the dilogarithms and the series expansions of the logarithms of the remaining factors. Some transformations of this expression helps Mathematica to take the limit that gives the result. Figure 1 demonstrates how the asymptotic ratio is achieved with increasing n. It shows R p,q n, n/2, x n 3/4 for z = -9 at y = 2 w = -0.10539 . . . and the asymptotic ratio is given by (39), that is, e -0.719718 x 4 . The red curve is the asymptote and the blue curves are for finite n where the curves for larger n are closer to the asymptote.

Using (4) in combination with exactly the same technique as above we can also determine the growth rate of the middle coefficient. With r = q/p we have

n n/2 p,q = p n 2 /4 r n/2+1 ; r n/2 (r; r) n/2 (40)
and the q-Pochhammer bounds gives the growth base

θ(z) = lim n→∞ log n n/2 p,q n = π 2 + 6 w 2 -3 w z + 12 Li (-2 w/z) 12 w -6 z (41)
Note that when z → -2 -then θ(z) → -1/2 + log 2, corresponding to the case p = q in lemma A.2. With a slightly improved form of the integral estimate [START_REF] Kaporis | The unsatisfiability threshold revisited[END_REF] we can also show that

n n/2 p,q ∝ exp (n θ(z)) √ n (42)
Compare this with lemma A.2 where the exact order is given. However, in the manner of the proof of theorem A.4 we can determine the sum

Ψ p,q (n) ∝ n 1/4 exp (n θ(z)) (43

Wide and double-peaked distributions

Changing p only slightly, say on the order of 1/n 3/2 , allows us to move around in the region where the distribution is wide. With a proper choice of the coefficient of 1/n 3/2 we can get a nice ratio expression. Let q = 1 + z/n and

p = 1 + 2 w n + a w z 2 -4 w 2 3 (1 + w) n 3/2 (44)
and the end result is

R p,q n, n/2, x n 3/4 ∼ exp w z (2 w + z) 6 x 4 -2 a x 2 (45)
This puts the maximum at x = ± √ a for a > 0 and at x = 0 for a ≤ 0.

Figure 2 works like figure 1 but here with the parameter a set to 1. It shows R p,q n, n/2, x n 3/4 for z = -9 with p = 1 -0.105391/n -1.50171/n 3/2 and the asymptotic ratio is given by (45), that is, e -0.719718 (x 4 -2 x 2 ) . The red curve is the asymptote and the blue curves are for finite n where the larger n are closer to the asymptote. Had we set a < 0 the distributions would still be wide but with a single peak in the middle. (Colour online) R p,q n, n/2, x n 3/4 (blue curves) versus x at z = -9 and a = 1 so that q = 1-9/n and p = 1-0.105391/n-1.50171/n 3/2 for n = 2 6 , 2 8 , 2 10 , 2 12 and 2 14 . The red curve is the asymptote e -0.719718 (x 4 -2 x 2 ) .

Peakish distributions

If we instead choose a middle ratio of 1 + a/n then we receive a sharply peaked distribution. Say that we want

R p,q (n, n/2, 1) = p n/2 -q n/2 p n/2+1 -q n/2+1 = 1 + a n (46)
for p = 1 + y 1 /n + • • • and q = 1 + z/n as before. Note that a = -2 corresponds to a binomial distribution, so we are usually interested in the case -2 < a < 0. After expanding the equation we receive as before, from the first term of the left and right hand side, the equation

(a + y 1 ) e y 1 /2 = (a + z) e z/2 (47)
which has the solution y 1 = 2 w -a where

w = W a + z 2 e (a+z)/2 (48)
This gives a distribution with a width of the order √ n. Computing the limit ratio gives us thus giving us an essentially gaussian distribution. Note that it does not depend on z in its current form. Of course, we expect other terms to depend on z but these vanish when n → ∞.

R p,q n, n/2, x √ n ∼ e a x 2 (49
The growth base θ(z) of the middle coefficient can be determined to

θ(z) = (2 π 2 + 3 (a -2 w) (a + z -2 w) + 24 Li -2 w a+z 12 (2 w -a -z) (50)
as we did before in equation ( 41). More precisely the coefficient is

n n/2 p,q ∝ exp (n θ(z)) √ n (51)
and analogous to theorem A.10 this gives the sum

Ψ p,q (n) ∝ exp (n θ(z)) (52)

Two separate peaks

Suppose that we want the peaks located outside the middle. Defining k = n 2 (1 + µ) for 0 < |µ| < 1 means that we move the peaks out from the middle. We keep q = 1 + z/n and p = 1 + y 1 /n + y 2 /n 2 + • • • and solve R p,q (n, k, 1) = 1, i.e.

p n-k+1 -p k = q n-k+1 -q k (53) With k = n 2 (1 + µ)
we use (126) on both sides and find the equation

e y 1 /2 sinh µ y 1 2 = e z/2 sinh µ z 2 (54)
At this point it would be appropriate to define the function s = Ω µ (x), for 0 < µ < 1, as the maximum solution to the equation

x = e s sinh (s µ) (55)
Note here that the function e x sinh(µ x) has a minimum at (-atanh µ)/µ for 0 < µ < 1 and, due to symmetry, a maximum at the same point for -1 < µ < 0. The function Ω returns a value in the interval

- atanh µ µ < s < 0 (56)
The solution sought in our equation is thus y 1 = 2 w where We should mention that Ω µ (x) is a natural extension of W (x). In fact, if we let w = W (x e x ) then

w = Ω µ e z/
Ω µ (e x sinh µ x) = w + w (x 2 -w 2 ) 6 (1 + w) µ 2 + • • • (58)
giving a good approximation for small values of µ. Having computed p = 1 + y 1 /n we compute the upper and lower bounds of the ratio as before.

Unfortunately, the ratio has the rather ghastly expression

R p,q n, n 2 (1 + µ), x √ n ∼ exp x 2 2 2 w + z - e µw+ µz 2 e 2w -e z (z -2w) -e 2µw+w+ z 2 + e µw+2w+ µz 2 + e µw+ µz 2 +z -e w+µz+ z 2 ( 59
)
where w is defined by (57).

In figure 3 we show how the finite cases approach their asymptote for µ = 1/3 and z = -9. This gives p = 1 -0.153208/n and the asymptotic ratio is given by (59), that is, e -0.103551 x 2 . The red curve is the asymptote and the blue curves are for finite n where the larger n are closer to the asymptote.

The growth base θ(z) of the maximum coefficient is given by

θ(z) = π 2 + 3 w (z -2 w) (µ 2 -1) 6 (z -2 w) + Li (e z-2 w ) -Li (e (z-2 w) (1-µ)/2 ) -Li (e (z-2 w) (1+µ)/2 ) z -2 w ( 60
)
The coefficient is then proportional to

n n/2 p,q ∝ exp (n θ(z)) √ n (61)
and analogous to theorem A.14 this gives the sum

Ψ p,q (n) ∝ exp (n θ(z)) (62)
6 Moments

Once we have the ratios (39), (45), (49), (59) it is an easy task to compute moments of the distributions. Let us do this for the most interesting case of (45). First, to make the notation somewhat simpler, denote φ = Φ(z) = -w z (2 w + z)/6, i.e. φ > 0, with w as in (37), so that R p,q n, n/2, x n

1.0 R Figure 3: (Colour online)R p,q (n, (n/2) (1 + µ), x √ n) (blue curves) versus
x at µ = 1/3 and z = -9 so that q = 1 -9/n and p = 1 -0.153208/n for n = 2 6 , 2 8 , 2 10 , 2 12 , 2 14 and 2 16 . The red curve is the asymptote e -0.103551 x 2 .

where a is defined by (44). We use the notation

σ m = k - n 2 m (64)
for the mth moment of the probability distribution of k = 0, . . . , n. Define also

̺ m = +∞ -∞ |x| m exp -φ (x 4 -2 a x 2 ) dx, m ≥ 0 (65)
so that the mth moment becomes |x| m = ̺ m /̺ 0 . For m = 0 we have In general we have for m ≥ 0 that

̺ 0 =          π √ a 2 exp φ a 2 /2 I 1/4 (φ a 2 /2) + I -1/4 (φ a 2 /2) , for a > 0 √ -a √ 2 exp φ a 2 /2 K 1/4 (φ a 2 /2) , for a < 0 Γ(1/4) 2 φ 1/
σ m = k - n 2 m ∼ n 3 m/4 ̺ m ̺ 0 (68)
where ̺ m is given by

         φ -(m+1)/4 2 Γ( m+1 4 ) 1 F 1 ( m+1 4 , 1 2 , φ a 2 ) + 2 a √ φ Γ( m+3 4 ) 1 F 1 ( m+3 4 , 3 2 , φ a 2 ) , a > 0 φ -(m+1)/4 2 (m+1)/2 Γ( m+1 2 ) U( m+1 4 , 1 2 , φ a 2 ) , a < 0 φ -(m+1)/4 2 Γ( m+1 4
) , a = 0 (69) Here 1 F 1 (a, b, c) and U(a, b, c) denote the confluent hypergeometric functions of the first and second kind respectively. If we want to compute cumulant ratios we first need moment ratios which of course is easy now. For example, in the case of a = 0 we have

σ 4 σ 2 2 ∼ ̺ 0 ̺ 4 ̺ 2 2 = Γ(1/4) 4 8 π 2 = 2.1884 . . . (70)
This moment ratio is also obtained in the 5-dimensional Ising model, see e.g. [START_REF] Brezin | Finite size effects in phase transitions[END_REF].

7 Fine-tuning the exponents Note (45) and, as before, keep Φ(z) = -w z (2 w + z)/6 where w is defined by (37). Recall that the first and second absolute moments obtained from (69) for a = 0 are

σ 1 ∼ n 3/4 ̺ 1 ̺ 0 = n 3/4 √ π Γ(1/4) 1 Φ(z) 1/4 ∝ n 3/4 Φ(z) 1/4 (71) σ 2 ∼ n 3/2 ̺ 2 ̺ 0 = n 3/2 π √ 2 Γ(1/4) 2 1 Φ(z) ∝ n 3/2 Φ(z) (72)
The argument z is allowed to depend on n but probably not to a high order. At this point it is not clear how z may depend on n for the calculations leading to (45) and Φ(z) to work. We will assume, for the moment (see the end of this section), that the expressions for the moments of (69) are valid when z = O (log n).

The series expansion of W (x) is

w = W (x) = x -x 2 + 3 x 3 2 - 8 x 4 3 + • • • (73)
With x = e z/2 z/2, and note that z is negative, we have Set z = λ 0 + λ 1 log n + λ 2 log log n + λ 3 log log log n with λ 1 , λ 2 ≤ 0 and focus on the first term of (75).

w = W (x) = W z e z/2 2 = z e z/2
Φ(z) ∼ -1 12 (λ 0 + λ 1 log n + λ 2 log log n + λ 3 log log log n) 3 e λ 0 /2 n λ 1 /2 log λ 2 /2 n (log log n) λ 3 /2 (76)
We are interested in two special cases. First choose λ 1 < 0, λ 2 = -6 and λ 3 = 0. This gives

Φ(z) ∼ (-λ 1 ) 3 12 e λ 0 /2 n λ 1 /2 (77)
Combining this with (71) we receive

σ 1 ∼ 3 1/4 √ 2 π Γ(1/4) n 3/4-λ 1 /8 (-λ 1 ) 3/4 e λ 0 /8 (78)
and

σ 2 ∼ 2 √ 6 π Γ(1/4) 2 n 3/2-λ 1 /4 (-λ 1 ) 3/2 e λ 0 /4
(79)

Had we let λ 2 = 0, instead of λ 2 = -6, then we would have ended up with a factor log 3/4 n in the denominator of (78) and a factor log 3/2 n in the denominator of (79). For the second case we choose λ 1 = 0, λ 2 < 0 and λ 3 = -6. We get

Φ(z) ∼ (-λ 2 ) 3 12 e λ 0 /2 log λ 2 /2 n (80)
This together with (71) gives us

σ 1 ∼ 3 1/4 √ 2 π Γ(1/4) n 3/4 log -λ 2 /8 n (-λ 2 ) 3/4 e λ 0 /8 (81)
and

σ 2 ∼ 2 √ 6 π Γ(1/4) 2 n 3/2 log -λ 2 /4 n (-λ 2 ) 3/2 e λ 0 /4 (82)
These expressions obviously converge extremely slowly and are probably not of any use for n that might occur in practical situations.

We have managed to verify (78) and (79) by using the method described in section 5, again using Mathematica, only in the special case z =log n, i.e. λ 1 = -1, λ 0 = λ 2 = λ 3 = 0. We could then confirm that R p,q n, n/2, x n 7/8 Computing the moments of this distribution produces the same result as setting λ 1 = -1 and λ 0 = λ 2 = λ 3 = 0 in (76) and then computing the moments in the same the way we obtained (78) and (79). A more general computation seems not to be within reach with our current set of tools though. To conclude this section we note that the exponent of n in (83) is 3/4 -λ 1 /8. For this exponent to stay less than one we thus need λ 1 > -2, giving us a bound on z.

The Ising model

A state τ on a graph G is a function from the set of vertices to {±1}. There are thus 2 n states for a graph on n vertices. We define the energy of a state as E(τ ) = ij τ i τ j where the sum is taken over all edges ij of G. The magnetisation is defined as M (τ ) = i τ i with the sum taken over all vertices i of G. Note that -n ≤ M ≤ n and it only takes every alternate value, i.e. M ∈ {-n, -n + 2, -n + 4, . . . , n -4, n -2, n}. We will often need to refer to it in terms of how many negative spins the state has. If k spins are negative then M = n -2 k.

The partition function of the Ising model is defined for any graph G as

Z(G; x, y) = τ x E(τ ) y M (τ ) = E,M a(E, M ) x E y M (84)
The coefficients a(E, M ) then are defined as the number of states with energy E and magnetisation M . Denote the number of states at energy E by a(E) = M a(E, M ). Note that the number of states at magnetisation M is just n k , where k = (n -M )/2 is the number of negative spins. Let also Z k denote the terms of Z with magnetisation M = n -2 k for a graph on n vertices, so that Z = Z 0 + Z 2 + • • • + Z n , i.e. Z k are the terms corresponding to k negative spins.

If we evaluate the partition function in x = e K and y = e H with K the dimensionless coupling, or inverse temperature J/k B T , and H = h/k B T as the dimensionless external magnetic field, we obtain the physical partition function denoted Z = Z(G; K, H) = Z(G; e K , e H ), though we are usually interested only in the case when H = 0 (or y = 1). Analogously, we write

Z = Z 0 + • • • + Z n .
The dimensionless and normalised free energy is defined as F = (log Z) /n. From the derivatives of the free energy we can now obtain other physical quantities such as the internal energy ∂F/∂K and the specific heat ∂ 2 F/∂K 2 though we shall not be needing the latter for this investigation.

We assume the Boltzmann distribution on the states so that (with H = 0) the probability for state τ is We have then especially that the probability for energy E is

P (τ ) = e K E(τ ) Z (
P (E) = a(E) e K E Z (86)
and the probability for magnetisation M is

P (M ) = 1 Z E a(E, M ) e K E = Z k Z ( 87
)
where

M = n -2 k.
Denote by K * the coupling where

Z n/2-1 = Z n/2 = Z n/2+1
. This coupling will correspond to the relation R p,q (n, n/2, 1) = 1 for some choice of p, q. We define the (spontaneous) normalised magnetisation μ = |M | /n and the (spontaneous

) susceptibility χ = Var (|M |) /n = M 2 -|M | 2 /n.
The pure susceptibility is simply χ = Var (M ) /n = M 2 = 4 σ 2 /n. Since M = n-2 k we thus have μ = 2 σ 1 /n and χ = 4 (σ 2 -σ 2 1 )/n. Recall the traditional finite-size scaling laws which claim that in the critical region, i.e. near K c , μ ∝ L -β/ν and χ ∝ L γ/ν . Being near K c means that |K -K c | ∝ L -1/ν and we especially expect K * to belong to this region. Though the highand low-temperature exponents may or may not be equal for three dimensions, see [START_REF] Häggkvist | On the Ising model for the simple cubic lattice[END_REF] for an in-depth numerical investigation of this matter, the details of these exponents are not important for our present investigation. What matters is that there are exponents that guide the growth of e.g. the susceptibility near K c .

The complete graph

For a complete graph, denoted K n on n vertices and n 2 edges the partition function is easy to compute. Suppose k of the vertices are assigned spin -1 and the other n-k have spin +1. The magnetisation is obviously M = n-2 k and the energy is

E = k 2 + n -k 2 -k (n -k) = n 2 -2 k (n -k) (88)
The partition function is then

Z(K n ; x, y) = x ( n 2 ) y n n k=0 n k 1 x 2 k (n-k) 1 y 2 k ( 89
)
and with y = 1 we have where q = 1/x 2 . Thus we have

Z(K n ; x, 1) = x ( n 2 ) n k=0 n k q,q = x ( n 2 ) Ψ q,q (n) (90
Z(K n ; K, 0) = exp K n 2 Ψ q,q (n) ( 91
)
where q = exp(-2 K). Obviously we have

P (M = n -2 k) = P q,q (n, k) (92)
Since we have defined the critical temperature as the K = K * where the middle ratio is 1, i.e. P (M = -2) = P (M = 0) = P (M = +2) then this corresponds to the point where R q,q (n, n/2, 1) = 1, which takes place at q = n/(n + 2), see lemma A.1. Thus K * = 1 2 log 1 + 2 n for K n . In short, the partition function and the magnetisation distribution for K n can be expressed in terms of p, q-binomial coefficients. Does this hold for all graphs? No. In fact, it seems to only be true for K n . However, it does seem to hold asymptotically as the order of the graphs increase, for some interesting families of graphs. The precise formulation of such a statement remains and falls outside this paper.

The average graph

Let us compute the sum of all partition functions taken over all graphs on n vertices.

Zn (x, y) = G⊆Kn Z(G; x, y) = (93) n i=0 n i y n-2 i ( i 2 ) j=0 i 2 j x j ( n-i 2 ) k=0 n-i 2 k x k i (n-i) ℓ=0 i (n -i) ℓ x -ℓ = (94) n i=0 n i y n-2 i (1 + x) ( i 2 ) (1 + x) ( n-i 2 ) 1 + 1 x ( i (n-i) 2 ) = (95) (1 + x) ( n 2 ) y n n i=0 n i 1 y 2 i 1 x i (n-i) (96)
and for y = 1 we have

Zn (x, 1) = (1 + x) ( n 2 ) n i=0 n i 1 x i (n-i) = ( 97
) (1 + x) ( n 2 ) n i=0 n i q,q = (1 + x) ( n 2 ) Ψ q,q (n) ( 98
)
where q = 1/x so that K * = log(1+ 2/n). Again we have P (M = n -2 k) = P q,q (n, k). The mean magnetisation distribution can then be modelled by p, q-binomial coefficients, though with p = q, just as for the complete graph.

Z(K u,v ; x, y) = (99) u i=0 v j=0 u i v j y u+v-2 i-2 j x i j+(u-i) (v-j)-i (v-j)-j (u-i) = (100) x u v y u+v u i=0 v j=0 u i v j 1 y 2 i+j 1 x 2 i (v-j)+j (u-i) (101)
which for y = 1 gives us

Z(K u,v ; x, 1) = (102) x u v u i=0 v j=0 u i v j 1 x 2 i (v-j)+j (u-i) = ( 103
) u+v k=0 ℓ u ℓ v k -ℓ x (u-2 ℓ) (v-2 (k-ℓ)) (104)
which defines the partial sums for y = 1 as

Z k = ℓ u ℓ v k -ℓ x (u-2 ℓ) (v-2 (k-ℓ)) (105)
Data suggests that

P (M = n -2 k) = Z k Z ≈ P p,q (n, k) (106)
given an appropriate choice of p and q and for a rather wide range of temperatures. It does however not seem to hold if u differ from v. In the left panel of figure 4 we show a sample of magnetisation distributions together with fitted p, q-distributions for a K 32,32 . To find the appropriate p and q we used the method described in the p, q-find algorithm in section 3. The fit is excellent. The right panel of figure 4 shows y = n (p -1) versus z = n (q -1) for a range of temperatures and for complete bipartite graphs of different sizes. High temperatures are in the upper right corner and K = 0 gives p = q = 1, i.e. y = z = 0. As the temperature decreases, i.e. with increasing K, we move along the curves. The points are where the distribution is exactly flat, i.e. the inverse temperature K * where the middle probabilities are equal. We have no exact closed form expression for K * but one can show that the series expansion of this inverse temperature for K n/2,n/2 , i.e. a total of n vertices is The calculations behind this are rather long and were done with Mathematica. Compare this with the expansion in the previous subsection for the complete graph on n vertices, K n , which begins

K * = 2 n - 3
K * = 1/n -1/n 2 + 4/3n 3 + • • • .
In the lower left corner the distribution has P (M = n) = P (M = n -2), i.e. at K = log n/n. The points in the right panel of figure 4 should approach z = y = -2. One can show that the distribution of a balanced bipartite graph K n/2,n/2 has the shape R n, n/2, x n 3/4 ∼ exp -4 3

x 4 (108) at K * , just like the complete graph on n vertices. Since Φ(-2) = 4/3 we then assume that z, and thus also y, will approach -2.

The free energy

If the magnetisation distribution were indeed an exact p, q-binomial distribution then we could also express the free energy as

F(G; K) = K m n + log Ψ p,q (n) n ( 109
)
for a graph on n vertices and m edges and it is here implied that p and q depend on K. Why this expression? Note that Z 0 = a(m, n) e K m = e K m and thus we have 0, where the distributions are close to binomial, and for very high K where all the probability mass is concentrated on the extreme magnetisations. In the left plot of figure 5 we show the exactly computed free energy (red curve) for a complete bipartite graph on 16 + 16 vertices together with the p, q-approximation (109) (points). The fit is indeed very good for the whole temperature range. Taking a derivative of the points with respect to K produces a good approximation to the internal energy as the right plot shows.

P p,q (n, 0) = n 0 p,q Ψ p,q (n) = 1 Ψ p,q (n) = Z 0 Z = e K m Z (

Lattices

The plots in figure 4 are very representative for several graphs of interest. We intend to focus on the graphs that traditionally are studied in statistical physics; lattice graphs. We will take a look at the simple lattices in 1, 2, 3, 4 and 5 dimensions. Just to be clear, a 1D-lattice is a cycle C n on n vertices and it is 2-regular, i.e. 2 neighbours for each vertex. The 2-dimensional L × L-lattice is the cartesian product of two cycles on L vertices. The product thus has n = L 2 vertices and it is 4-regular. The d-dimensional L × L × • • • × L-lattice is a product of d cycles on L vertices, thus having a total of n = L d vertices. It is obviously 2 d-regular. Assuming finite-size scaling to hold then for a d-dimensional lattice we have

σ 1 = n μ/2 ∝ n L -β/ν = n n 1/d -β/ν = n 1-β/d ν (111)
and correspondingly for the second moment

σ 2 = n χ/4 ∝ n L γ/ν = n n 1/d γ/ν = n 1+γ/d ν (112)
Note also that for an r-regular triangle-free graphs we have

P (M = n) = P (M = n -2) when K = log n 2 r .
25 For 1D-lattices we can compute the coefficients a(E, M ) exactly. It is an exercise to show that the number of states with k negative spins and ℓ negative spin products (over the edges) is

a(E, M ) = k k -ℓ/2 n -k -1 n -k -ℓ/2 + k -1 k -ℓ/2 n -k n -k -ℓ/2 (113)
where M = n -2 k and E = n -2 ℓ. The distribution of magnetisations do not behave in a way representative for lattices of higher dimension. However, for extremely low temperatures the probabilities P (M = -n) = P (M = n) will dominate the other probabilities. The two outermost probabilities, P (M = n) and P (M = n -2), are equal when K = log n 4 . For the 1D-lattice the distribution is here sharply unimodal, while for higher dimensions the distribution is bimodal and has its peaks at the extreme magnetisations. For K larger than (log n)/4 the distribution actually has three peaks, i.e. a local maximum at M = 0. For 1D-lattices the p, q-approximation of the distribution thus breaks down beyond this K since it can not model a local maximum in the middle as well as peaks at the ends; they are at most bimodal. For K less than this point the p, q-distribution is a very good approximation. Figure 6 demonstrates this clearly; for the flattest distribution (low temperature) the fitted p, q-distribution starts to deviate from the actual distribution. In figure 7 we plot y = n (p -1) and z = n (q -1) versus K for a range of different n. Clearly there is some limit curve here, though we have not established what the limit function is. In figure 8 we see y = n (p-1) versus z = n (q -1) for different n. The right plot of figure 8 shows the value at K that gives the maximum value of y. The fitted straight line gives the limit 0.1333, very close to 2/15. What about the values of y and z? Indeed they converge beautifully as figure 9 indicates. The limit for y is about 1.010 and z approaches a value of -3.537. Though we can not exactly solve what y and z should be at K = 2/15 we can at least see how y and z relate at this point. For an infinite 1-dimensional lattice we have that χ = e 2 K , see e.g. [START_REF] Baxter | Exactly solved models in statistical mechanics[END_REF]. The second moment then should behave as

σ 2 ∼ n χ 4 = n e 2 K 4 ( 114
) Let ℓ = k -n/2 and σ = √ σ 2 .
For high temperatures we expect ℓ/σ to be normally distributed and thus

P (ℓ) ∼ exp -(ℓ/σ) 2 /2 σ √ 2 π (115)
The probability ratio is then and for ℓ = 1 this simplifies to

R (n, n/2, ℓ) = P (ℓ) P (0) = exp -ℓ 2 /
R (n, n/2, 1) = exp -1/2 σ 2 = exp -2 e -2 K /n ∼ 1 - 2 e -2 K n (117)
Compare this with (46). We thus have a = -2 e -2 K . Now y and z are related as y = 2 w -a where w is defined by (48). If we set K = 2/15 then a = -1.531857, and choosing z = -3.537 indeed gives us y = 1.01002. To actually solve z as a function of K seems harder though. -1) versus 1/n for C n . Right: value of z (q -1) versus 1/n for C n when y is at its maximum. In both cases n = 8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256. experimentation suggests that y and z for small K behave as

y(K) ≈ c 1 √ K + c 2 K (118) z(K) ≈ -c 1 √ K + c 2 K (119)
where c 1 ≈ 6.164 and c 2 ≈ -10.33.

If we compare the magnetisation distribution with the fitted p, q-binomial distribution we will of course see some small deviations. Say that we choose K = 1/4 and see how the two distributions differ with increasing n. For n = 128 we show how the difference in probabilities look in the left plot of Figure 10. We see a maximum which is located near ±1.5 σ and a minimum at ±3 σ, where σ = e K √ n/2 is the standard deviation. A simple scaling projection gives that the value at these peaks decrease at an order of 1/n 5/2 . Comparing the moments of the distributions we can detect no error at all, up to numerical precision.

However, even though the difference in probabilities is small, and decreases quickly with n, this is not enough to get the free energy as an ap- proximation through (109), except for very small K. Recall here that the asymptotic free energy for the 1D-lattice is given by F(K) = log (2 cosh K).

Let the p, q-approximation of the free energy, i.e. the right hand side of (109), be denoted F pq . For C n we thus have F pq = K + (log Ψ p,q (n)) /n. If the magnetisation distribution of C n were perfectly fitted by a p, q-binomial distribution then F = F pq . In the right plot of Figure 10 we show the asymptotic difference F -F pq versus K. This was obtained by fitting a second degree polynomial to the difference versus 1/n for n = 256, 384, 512, 768, 1024 and taking the polynomial at x = 0, i.e. n = ∞, as the asymptotic value.

2D-lattices

For the 2-dimensional lattices we can rely on exact data only for up to L = 16 and they were computed according to the method in [START_REF] Häggkvist | The Ising partition function for 2D grids with cyclic boundary: computation and analysis[END_REF]. We have sampled data for L = 32, 64, 128, 256, 512, collected with the methods described in [START_REF] Häggkvist | A Monte Carlo sampling scheme for the Ising model[END_REF] and [START_REF]Reconstruction of the finite size canonical ensemble from incomplete micro-canonical data[END_REF]. These methods gave us the energy distribution and then it is just a matter of combining this with the distribution of magnetisations for each given energy as described in [START_REF] Häggkvist | On the Ising model for the simple cubic lattice[END_REF]. Figure 11 shows an example of some distributions for the 128 × 128-lattice together with their fitted p, qbinomial distributions. The fit is fairly good, but hardly excellent near K * . However, as the figure shows, at K = 0.439 (i.e. for L = 128) the fit is practically spot on. For the 2D-lattices we have studied there is always one such temperature where the p, q-distribution fit particularly well. This point is located between K * and K c and is very close to, but not exactly equal to, the point where the susceptibility is at its maximum. Of course, for high temperatures (small K) and low temperatures (high K) the fit is typically very good but in the high-temperature region the measured y and z are unfortunately extremely sensitive to noise. As we get closer to the critical region where the distribution becomes bimodal The deviation of the fitted distribution from the sampled Ising distribution is of course reflected in the quotient of their variances. This quotient stays close to 1 for high temperatures (modulo noise) but jumps, first down, to about 0.6, and then up, to about 1.2 in the region where the distribution changes shape. It does, however, take the value 1 at some point and this is where the two distributions appear (almost) indistinguishable.

A chi-square test is more discerning than the human eye though. We tested the hypothesis that the fitted p, q-distribution is consistent with the distribution obtained from sampled data for a few temperatures between K * and K c for L = 128. The best result was actually found at K = 0.4384. Running a standard Metropolis sampling algorithm we collected about 2.6 million samples. Between each measurement of |M | we made about five sweeps to ensure at least an expected n successful spin flips.

The magnetisations at the tail which gave less than six samples were lumped into one bin. The test statistic then was about 6549 and we received 6039 bins. The reduced test statistic then is 1.08 which indicates a good fit. Unfortunately at the 5% significance level. However, the variance quotient between fitted and measured distribution is 1.004 and for higher moments the quotient is even closer to 1. Our fitted distribution thus have correct moments even though the probabilities are not correct.

Recall from section 7 how the exponents of the moment growth rates could be computed if we allow z to depend on n. For the 2D-lattices it is known that β = 1/8, γ = 7/4 and ν = 1, see [START_REF] Onsager | Crystal statistics I. A two-dimensional model with an order-disorder transition[END_REF], [START_REF] Yang | The spontaneous magnetization of a two-dimensional Ising model[END_REF] and [START_REF] Abraham | Susceptibility of the rectangular Ising ferromagnet[END_REF]. Thus the first moment σ 1 should scale as n 15/16 and the second moment σ 2 as n 15/8 . From equation ( 78) and (79) this would be achieved by choosing λ 1 = -3/2, λ 2 = -6 and λ 3 = 0. The left plot of Figure 12 shows z versus log n at K * together with the curve 3 -1.5 log n -6 log log n. The constant λ 0 is chosen only to make the curve look plausibly near the points. The point for L = 512 deviate slightly but we suspect that noise in the sampled data explains this. With λ 0 = 3 the coefficient of n 15/8 obtained from (79) would be 0.301 though the measured σ 2 divided by n 15/8 are closer 0.08. To get this we have to choose λ 0 ≈ 8.3. In that case the convergence is extremely slow. Note also that the fitted p, q-distribution is far from perfect which would contribute some amount of error as well. The right plot of figure 12 shows y vs z for a range of temperatures. The points representing K * may appear to lie on the z-axis but they are are slightly below it. In the 1D-case we suspected that there is a limit curve for the high-temperature region, but we suspect that the exact data that produced this part of the plot rely on far too small lattices to give any conclusive evidence. Also, the p, q-find algorithm is rather sensitive to noise in this region to be useful for sampled data. However, as we said before, this problem goes away once K ≥ K * . Figure 13 shows y and z versus K for all the lattices though for the sampled data we only show low-temperature data. The red line is located at

K c = atanh √ 2 -1 ≈ 0.44068.

3D-lattices

For these lattices we only have exact data for L = 4 and sampled data for L = 6, 8, 12, 16, 32, 64. The situation is actually somewhat better for 3D-lattices. Figure 14 shows some distributions in the vicinity of K * for L = 32 together with the fitted p, q-distributions. For K ≥ K * , just when the distributions become bimodal, the fit is certainly less than perfect, but near K * the p, q-approximation is actually rather good.

If we take the quotient between the variance of the fitted p, q-binomial distribution and the variance of the Ising distribution we find that it stays between 0.85 and 1.02 which is slightly better than for the 2D-lattices.

Running the chi-square test as we did for the 2D-lattice the best result with L = 32 was found at K = 0.2208 which gave the test statistic 5917 and, at 5550 degrees of freedom, the reduced statistic 1.066. Again, as for the 2D-lattice, the null hypothesis fails at the 5% significance level since χ 2 0.05 (5550) ≈ 5724. However, the first four moment quotients stays between 1.00 and 1.01 at this temperature.

In the left plot of figure 15 we show z versus log n at K * . The fitted line through the points corresponds to z = -5.3log n and is not too bad an approximation. However, in [START_REF] Häggkvist | On the Ising model for the simple cubic lattice[END_REF] it was estimated that the growth rate exponent at K c of the susceptibility is γ/ν = 1.978 ± 0.009 (assuming γ = γ ′ and ν = ν ′ ). For the magnetisation it was estimated β/ν = 0.5147 ± 0.0007. Translated into exponents of n this means 1.657 ≤ 1 + γ/3ν ≤ 1.663 and 0.8282 ≤ 1-β/3ν ≤ 0.8287. If we choose λ 1 = -5/8 in (78) and (79) the first moment exponent would be 53/64 = 0.828125 and 53/32 = 1.65625 for the second moment, slightly below the lower bound of the estimate intervals. Choosing λ 1 = -2/3 would give exponents 5/6 = 0.8333 . . . and 5/3 = 1.666 . . . respectively, slightly above the upper bound of the intervals. Let us suggest, as an example, that λ 0 = 6.8, λ 1 = -2/3, λ 2 = -6 and λ 3 = 0 in the expression (79). In figure 15 i.e. z = 6.8 -(2/3) log n -6 log log n. Will the points eventually converge to the curve? It would take considerably larger lattices to shed any light on this. We also have the problem what λ 0 should be. Using λ 0 = 6.8 means that the coefficient in (79) is about 0.393. Comparing the measured σ 2 with n 5/3 gives a factor of roughly 0.16 though the data are certainly far from conclusive. Since the distribution fit is not perfect a different constant is perhaps to be expected. Also, slow convergence is to be expected here.

The right plot of figure 15 shows y versus z for K > K * . Note the peculiar backwards movement of z getting more and more pronounced for larger L. The curves for 16, 32 and 64 show signs of approaching some limit curve. We don't have data for very low temperatures for the smaller lattices though, except for L = 4. The plots in figure 16 shows y and z versus K for K > K * . The red lines show location of K c ≈ 0.2216546, found in [START_REF] Häggkvist | On the Ising model for the simple cubic lattice[END_REF], but see also [START_REF] Rosengren | On the combinatorial solution of the Ising model[END_REF] for a theoretical estimate of K c .

4D-lattices

In the case of 4-dimensional lattices we have sampled data of magnetisation distributions for L = 4, 6, 8, 10, 12, 16. Figure 17 shows some of these magnetisation distributions for L = 12 near K * together with fitted p, q- binomial distributions. The fit is quite good, considerably better than for 2D and 3D, in the whole range of selected temperatures. Though it is hard to distinguish the fitted curves from the magnetisation curves, there is a small deviation near the middle. The quotient between the variances stays between 0.97 and 1.01 for all K and L where we have sampled data. The best result from a chi-square test for L = 12 was obtained at K = 0.1494. A Monte Carlo sampling here gave us 3.6 million measurements of |M |. Comparing the sampled data with the p, q-binomial distribution we received the test statistic 3205 with 3137 degrees of freedom. This gives the reduced test statistic 1.02 and since χ 2 0.05 (3137) ≈ 3268 the null hypothesis stands. The quotients between the first four moments stayed between 0.999 and 1.003 at this temperature. However, it should be mentioned that at other temperatures, say at K = 0.1498, the chi-square test fails the p, q-distribution. The best fit is thus found near K * . How should z at K * depend on n? Actually, taking the data at facevalue they are rather well-fitted to the simple formula z = -6.5 -0.45 log n. However, for the 4D-lattice we have γ = γ ′ = 1, β = 1/2 and ν = ν ′ = 1/2. This gives that 1 + γ/d ν = 3/2 and 1 -β/d ν = 3/4. Moreover, according to [START_REF] Lai | Finite-size scaling of the Ising model in four dimensions[END_REF] there should be a correction to this. They calculated, using renormalization group techniques, that the susceptibility should scale as L 2 √ log L near K c . This means that σ 2 should scale as n 3/2 √ log n. From (82) we see that we have to choose λ 2 = -2, with λ 1 = 0 and λ 3 = -6, to obtain this. In the left plot of figure 18 we have set z = -1.2 -2 log log n -6 log log log n and plotted it versus log log n. The curve would then behave as a limit curve rather than as a fitted curve. The choice of coefficient λ 0 = -1.2 is only supported by the human eye as a guide rather than any theory and herein lies a problem. With this choice the coefficient of (82) is about 0.558. However, dividing the measured σ 2 at the different K * with n 3/2 √ log n gives values close to 0.15. This discrepancy could be due to several sources; e.g. the expression in (82) could be incorrect or our data could be suffering from very slow convergence. In the right plot of figure 18 35

5D-lattices

For the 5-dimensional lattices we have sampled data of magnetisation distributions only for L = 4, 6, 8, 10, 12. The distributions in figure 20 are extremely well fitted by p, q-binomial distributions; it is almost impossible to tell them apart with the naked eye.

Looking at the quotient between the variances we note that it stays between 0.99 and 1.005 for all L and K. statistic 4001 at 3996 degrees of freedom. The reduced test statistic is thus 1.001, a very good fit. The test statistic is compared with χ 2 0.05 (3996) ≈ 4144, giving us a safe margin to let the null hypothesis stand. The quotients between the first four moments are all between 1.000 and 1.004 at this temperature.

In five dimensions the susceptibility near K c scales as L 5/2 , see [START_REF] Cardy | Scaling and renormalization in statistical physics[END_REF]. Thus σ 2 should scale as n 3/2 which is exactly what we receive when keeping z fixed. So, for z constant we obtain σ 1 ∝ n 3/4 and σ 2 ∝ n 3/2 . The left plot of figure 21 shows z at K * for L = 4, 6, 8, 10, 12. If z approaches a constant then what is the limit value? Extracting the limit z from this plot is futile of course. The right plot of figure 21 shows y vs z for the different lattices together with the points K * and the curve y = 2 w. In figure 22 we show y and z versus K for K ≥ K * with an estimated K c marked as a red line. Despite the noise in the plots it seems plausible that z stays essentially constant very close to K * (and K c ) and that only y moves. Let us assume this and see where this leads us. We employ the moment expressions in section 6 in terms of the parameter a to model the behaviour near K * . A normalised first cumulant of the absolute magnetisation |M | /2 n 3/4 = σ 1 /n 3/4 should approach ̺ 1 /̺ 0 when plotted as a function of a for a fixed z. Analogously, the second cumulant (normalised) should behave as

σ 2 -σ 2 1 n 3/2 → ̺ 2 ̺ 0 (120)
where the ̺ m were defined in section 6. Note that for a fixed z the ̺ m now depend only on a. The third and fourth cumulants of the absolute magnetisation, divided by respectively 8 n 9/4 and 16 n 3 , quite analogously approach their corresponding limits Through a simple scaling analysis based on our sampled data we have found that the normalised third cumulant has a limit maximum of about 0.0205 and a minimum of -0.0500. The fourth normalised cumulant has a limit maximum of 0.0229 and a minimum of -0.0278, based upon our sampled data. Choosing z = -9.87 puts the maximums and minimums of the limit curves at appropriate values. Now we identify the coupling K where the minimum of the fourth cumulant occurs with the point a where the minimum of the corresponding limit curve occurs and likewise for the maximum, thus providing us with a rescaling translating K into a. In figure 23 and 24 the first four cumulants are shown together with their estimated limit curves for z = -9.87. Indeed the red curve may provide us with a limit. Given a lattice size L we denote by K min (L) the location of the minimum fourth cumulant and by K max (L) the location of the maximum. Analogously for the limit curve, given a z we denote by a min (z) and a max (z) the location of the minimum and maximum fourth cumulant. For z = -9.87 we have a min ≈ 1.06965 and a max ≈ 2.51275. A simple scaling analysis gives that 39 roughly K max (L) ≈ K c + 0.22/L 5/2 and K max (L) -K min (L) ≈ 0.093/L 5/2 . Also K c ≈ 0.113915, see [START_REF] Luijten | Finite-size scaling above the critical dimension revisited: The case of the five-dimensional Ising model[END_REF]. Thus, in principle at least, the rescaling between a and K is

̺ 3 ̺ 0 -3 ̺ 1 σ 2 ̺ 2 0 + 2 ̺ 3 1 ̺ 3 0 (121) and ̺ 4 ̺ 0 -4 ̺ 1 σ 3 ̺ 2 0 -3 σ 2 2 ̺ 2 0 + 12 ̺ 2 1 ̺ 2 ̺ 3 0 -6 ̺ 4 1 ̺ 4
K(a) ∼ K max (L) -K min (L) a max (z) -a min (z) (a -a max (z)) + K max (L) (123)
However, this kind of expression is somewhat too simplistic to get figure 24.

It would take higher-order corrections to scaling to produce it but this would probably take a more involved numerical study of the 5D-model. Other investigations of the 5D-lattice includes e.g. [START_REF] Brezin | Finite size effects in phase transitions[END_REF], [START_REF] Mon | Finite-size scaling of the 5D Ising model[END_REF] and [START_REF] Luijten | Finite-size scaling above the critical dimension revisited: The case of the five-dimensional Ising model[END_REF].

Conclusions

The magnetisation distribution for the complete graph is exactly described by the p, q-binomial distribution, corresponding to the special (or limit) case of p = q. For balanced complete bipartite graphs this is most likely also true in some limit sense, yet to be made precise. Actually, it appears that for most graphs, at least those which are more or less regular, the magnetisations are well-fitted by a p, q-binomial distribution for some choice of p and q. The exact extent to which the p, q-binomial approximation is good we do not yet know (e.g. convergence in moment or probability) nor the exact class of graphs that would satisfy this. We have investigated the matter more closely for lattices of dimension one through five. In general they are always well-fitted by p, q-binomial distributions for high-and low-temperatures but the problems arise near K c , or rather K * where the distribution changes from unimodal to bimodal.

For the 1-dimensional lattices (having no such bounded K * ) the situation is basically always that of high temperatures. It seems possible to give expressions for p and q in terms of K in this case though we have not done so. The difference in probabilities between the magnetisation distribution and the fitted p, q-binomial distribution vanishes at a rapid rate and the moments of the distribution are equal up to numerical precision given large enough n. This is, however, not enough to give the correct free energy using only the fitted distribution.

For 2-dimensional lattices the distributions near K * are least well-fitted by the p, q-binomials. Though there are temperatures close to K * where the moments of the fitted p, q-binomial distribution are close to correct, the fitted distribution fails a chi-square test. In the 3-dimensional case we get a somewhat better fit but not enough to pass a chi-square test. For 4-dimensional lattices the distributions are clearly much better fitted by p, q-binomials, though some discrepancy still remains just above K * . For 5dimensional lattices even this small discrepancy is gone, leaving us perfectly fitted (that is, to the human eye) p, q-binomial distributions. In fact, for 40 both the 4-and 5-dimensional case the fitted distributions now also pass a chi-square test. In the 5-dimensional case the values of z at K * should approach a limit value. We estimated this limit to -9.87 and, using this limit value, compared the first four normalised cumulants for finite lattices with the (possible) limit curves.

Our theoretical investigation gives that the p, q-binomial distribution has correct moment growth exponents given the correct parameter scaling of z. We can in fact fine-tune this exponent to contain the logarithmic correction for d = 4, as predicted by renormalisation group theory. Also, both for d = 4 and d = 5 our numerical investigation gives that the sampled data are correctly modelled by a p, q-binomial distribution, not just the moments. We suggest thus that the magnetisation distribution is in fact a p, q-binomial distribution for d = 4, 5, at least near K * .

To obtain the fitted distributions we described and used a rather simple method to determine p and q given a distribution. Possibly this method is not optimal since it simply forces the distribution to be correct at a single point rather than providing a good overall-fit. It is also sensitive to noise when the distributions are unimodal, thus making it difficult to determine p and q. On the other hand it works extremely well for bimodal distributions where the noise sensitivity problem vanishes.

The p, q-binomial coefficients are just a tweaked form of q-binomials, i.e. they are multiplied by a power of p. It is possible that a different choice of factor would produce better results in the case of 2-and 3-dimensional lattices.

It would be interesting to see if the magnetisation distribution for quantum spin models or for spin-glass models can be modeled by p, q-binomial distributions. What is the sum of the coefficients at this point? To answer this we compare the middle coefficient with a coefficient situated at some carefully chosen distance from the middle. How big is the middle coefficient? Note first that n n/2 q,q = q n 2 4

n n/2 (125)

Lemma A.2. For q = n n+2 we have n n/2 q,q ∼ 2 e π n

2 √ e n meaning that the quotient between the left-and right-hand side goes to 1 as n → ∞. The proof follows from an easy application of the identity 1 +

x n n = e x 1 -

x 2 2 n + x 3 3 n 2 + x 4 8 n 2 + • • • (126)
and we leave it to the reader. A somewhat more involved application of (126) is the following lemma Lemma A.3. Let x be some real number. For q = n n+2 we have R q,q n, n/2, x n 3/4 ∼ exp -4 3

x 4

This allows us to give the exact order of the coefficient sum.

Theorem A.4. For q = n n+2 we have

Ψ q,q (n) ∼ Γ(1/4) 3 1/4 n 1/4 √ π 2 √ e n-1
Proof. The calculations goes as follows though we leave out some details.

Ψ q,q (n) = n n/2 q,q n/2 k=-n/2 R q,q (n, n/2, k) ∼ n 3/4 n n/2 q,q +∞ -∞ exp - This is the coefficient sum at the point where the distribution becomes flat in the middle region. We can do even better if we allow ourselves to move around in the vicinity of this point.

Lemma A.5. Let q = n n+2 + a n 3/2 for some real number a. Then

n n/2 q,q ∼ 2 e π n 2 √ e n exp a √ n 4
This can be verified using (126) as can the following lemma.

Lemma A.6. Let a and x be real numbers. For q = n n+2 + a n 3/2 we have R q,q n, n/2, x n 3/4 ∼ exp -a x 2 -4 3

x 4

We now have the resources to estimate the sum of the coefficients for a whole spectrum of values of q near n/(n + 2). The next theorem can be shown using the same technique as Theorem A.4, though the result gets slightly more complicated due to the integral on the right hand side in the previous lemma.

Theorem A.7. Let q = n n+2 + a n 3/2 . For a > 0 the asymptotic order of Ψ q,q (n) is For a < 0 the asymptotic order of Ψ q,q (n) is Here I α (x) and K α (x) denote the modified Bessel functions of the first and second kind respectively. The mth moment is simpler to express using an integral formulation. and the asymptotic behaviour of the middle coefficient is given by lemma A.5.

n 1/4
n 1/4
The same technique allows us to repeat this for points farther away from the critical point q = n/(n + 2). If we increase q by a/n then the coefficients get sharply concentrated in the middle like that of standard binomial coefficients. Again (126) to the rescue. Note that above, when the coefficients had a rather wide distribution, we examined their behaviour at x n 3/4 from the middle. Under our current assumption of q the distribution gets more sharply concentrated around the middle (basically they become gaussian), thus we study their behaviour at x √ n from the middle. Lemma A.9. Let q = n+a n+2 and a > 0. Then R q,q n, n/2, x √ n ∼ exp -a x 2

Combining lemma A.9 with lemma A.8 gives our next theorem on the coefficient sum.

Theorem A.10. Let q = n+a n+2 and a > 0. Then

Ψ q,q (n) ∼ 2 e a 2 √ e n exp a n 4 - a 2 8
If we instead decrease a below zero the coefficient sequence becomes sharply bimodal, with all its mass concentrated around two peaks. We can of course connect the position of the peaks with the parameter a. Suppose that we want one of the peaks to have its maximum located at k and k -1 where k = (n/2) (1 + µ) so that their ratio becomes 1. A simple calculation shows that lim n→∞ R q,q n, n 2 (1 + µ) , 1 = 1 + µ 1 -µ exp (µ (a -2)) (129)

Setting the limit to 1 and solving the equation gives the next lemma.

Lemma A.11. Let q = n+a n+2 and 0 < |µ| < 1. If a = 2 (1 -(1/µ) atanh µ) then the limit ratio in (129) is 1.

We continue as before and estimate the growth rate of the peak coefficient. The result (and the proof) is somewhat more complicated but follows from an application of (126).

Lemma A.12. Let q, µ and a be defined as in lemma A.11. Then Again we take x √ n steps away from the peak and find the shape of the distribution.

n n 2 (1 + µ) q,q ∼ √ 2 exp  n 2 " log 4
Lemma A.13. Let q, µ and a be defined as in lemma A.11. Then R q,q n, n 2 (1 + µ) , x √ n ∼ exp 2 x 2 1 µ 2 -1 + atanh µ µ

Finally we get the sum by multiplying the integral of the ratio with the peak coefficient and 2

√ n, where the factor 2 is due to that we have two peaks.

Theorem A.14. Let q, µ and a be defined as in lemma A.11. Then Ψ q,q (n) ∼ √ 2 π n n n 2 (1 + µ) q,q µ (1 -µ 2 ) µ + (µ 2 -1) atanh µ

2

 2

  notations, the very basics

) 5 Page

 5

8

 8

Figure 1 :

 1 Figure1: (Colour online) R p,q n, n/2, x n 3/4 (blue curves) versus x at z = -9 and y = 2 w = -0.105391 for n = 2 6 , 2 8 , 2 10 and 212 . The red curve is the asymptote e -0.719718 x 4 .

13

 13

Figure 2 :

 2 Figure2: (Colour online) R p,q n, n/2, x n 3/4 (blue curves) versus x at z = -9 and a = 1 so that q = 1-9/n and p = 1-0.105391/n-1.50171/n 3/2 for n = 2 6 , 2 8 , 2 10 , 2 12 and 214 . The red curve is the asymptote e -0.719718 (x 4 -2 x 2 ) .

22 http 8 . 3

 2283 The complete bipartite graphWhat about K u,v , i.e. the complete bipartite graph on n = u + v vertices? Now the partition function is

Figure 5 :

 5 Figure 5: (Colour online) Left: free energy (red curve) compared to the formula (109) (points) for K 16,16 . Right: internal energy (red curve) compared to the derivative of the points produced by (109) (points) for K 16,16 .

Figure 7 :

 7 Figure 7: (Colour online) Left: y = n (p -1) versus K for C n . Right: z = n (q -1) versus K for C n . Both plots are for n = 16, 32, 64, 128 (larger cycles stretch farther to the right).

Figure 9 :

 9 Figure 9: (Colour online) Left: Maximum value of y = n (p-1) versus 1/n for C n . Right: value of z (q -1) versus 1/n for C n when y is at its maximum. In both cases n = 8,12, 16, 24, 32, 48, 64, 96, 128, 192, 256.

Figure 10 :

 10 Figure 10: (Colour online) Left: The difference in probabilities vs k -n/2 between the fitted p, q-binomial distribution and the magnetisation distribution for C n at K = 1/4 for n = 128. Right: The asymptotic difference F pq -F versus K from scaling on n = 256, 384, 512, 768, 1024.

29

 29

Figure 11 :

 11 Figure 11: (Colour online) Magnetisation distributions for the 128 × 128lattice (red) together with fitted p, q-binomial distributions P p,q (n, k) (blue) vs k -n/2 at K = 0.432, K * = 0.43467, K = 0.437, K = 0.439 and K c = 0.44068 (downwards at y-axis).

30 Page

 30 χ 2 0.05 (6038) ≈ 6220 so we are obliged to reject the hypothesis

Figure 12 :

 12 Figure12: (Colour online) Left: z = n (q -1) vs log n at K * for the L × Llattice, L = 6, 8, 10, 12, 14, 16, 32, 64, 128, 256, 512. The curve is 3 -1.5 log n -6 log log n. Right: y = n (p -1) vs z = n (q -1) for the L × L-lattice, L = 6, 8, 10, 12, 14, 16, 32, 64, 128, 256, 512 (512 barely visible near the z-axis). The points represent K * .

31 httpFigure 13 :

 3113 Figure 13: (Colour online) y = n (p -1) vs K (left) and z = n (q -1) vs K (right) for the L × L-lattice, L = 6, 8, 10, 12, 14, 16, 32, 64, 128, 256, 512. The points represent K * and the red line is at K c . The larger lattices have their points farther to the right in the plots.

32 PageFigure 14 :

 3214 Figure 14: (Colour online) Magnetisation distributions for the 32 × 32 × 32lattice (red) and the fitted P p,q (n, k) (blue) vs k -n/2 for K = 0.2204, K * = 0.22066, K = 0.2210, K c = 0.2216546 and K = 0.2220 (downwards at the y-axis).

Figure 15 :

 15 Figure 15: (Colour online) Left: z = n (q-1) vs log n at K * for the L×L×Llattice, L = 4, 6, 8, 12, 16, 32, 64. The line though the points is -5.3log n and the curve is 6.8 -(2/3) log n -6 log log n. Right: y = n (p -1) vs z for the L × L × L-lattice, L = 4, 6, 8, 12, 16, 32, 64 (leftwards) for K > K * . Higher values of K when we move downwards left. The red curve is y = 2 w with w defined by (37).

Figure 16 :

 16 Figure 16: (Colour online) Left: z = n (q -1) vs K with K > K * for the L× L× L-lattice, L = 4, 6, 8, 12, 16, 32, 64 (downwards). Right: y = n (p -1) vs K with K > K * for the L × L × L-lattice, L = 4, 6, 8, 12, 16, 32, 64 (upwards). In both plots the red line indicates location of K c and the points are the locations of K * .

  34

Figure 17 :

 17 Figure 17: (Colour online) Magnetisation distributions for the 12× 12× 12× 12-lattice (red) and the fitted P p,q (n, k) (blue) vs k -n/2 for K = 0.1490, K * = 0.149255, K c = 0.149695, K = 0.1500 and K = 0.1505 (downwards at the y-axis).

Figure 18 :Figure 19 :

 1819 Figure 18: (Colour online) Left: z = n (q -1) vs log log n at K * for the L × L × L × L-lattice, L = 4, 6, 8, 10, 12, 16. The curve is -1.2 -2 log log n -6 log log log n. Right: y = n (p -1) vs z = n (q -1) for the L × L × L × Llattice, L = 4, 6, 8, 10, 12, 16 (leftwards). Higher temperatures (low K) begin at the upper right part of the plot and with lower temperatures we move down to the left. The red curve is y = 2 w with w defined by (37).

Figure 20 :

 20 Figure 20: (Colour online) Magnetisation distributions for 8 × 8 × 8 × 8 × 8lattice (red) and the fitted P p,q (n, k) (blue) vs k -n/2 for K = 0.1137, K * = 0.113786, K c = 0.113914, K = 0.1143 and K = 0.1147.

Figure 21 :

 21 Figure 21: (Colour online) Left: z = n (q -1) vs log n at K * for the L × L × L × L × L-lattice, L = 4, 6, 8, 10, 12. The straight line is constant at z = -9.87. Right: y = n (p -1) vs z = n (q -1) for the L × L × L × L × Llattice, L = 4, 6, 8, 10, 12 (leftwards). Higher temperatures (low K) begin at the upper right part of the plot and with lower temperatures we move down to the left. The red curve is y = 2 w with w defined by (37).

Figure 22 :

 22 Figure 22: (Colour online) Left: z = n (q -1) vs K with K > K * for the L×L×L×L×L-lattice, L = 4, 6, 8, 10, 12 (downwards). Right: y = n (p-1) vs K with K > K * for the L × L × L × L × L-lattice, L = 4, 6, 8, 10, 12 (upwards). In both plots the red line indicates location of K c and the points are the locations of K * .

Figure 23 :

 23 Figure 23: (Colour online) Normalised first (left) and second (right) cumulants for the L × L × L × L × L-lattice, L = 4, 6, 8, 10, 12, 16, 20, 24 (blue) versus a for z = -9.87 together with the limit curve (red).

Figure 24 :

 24 Figure 24: (Colour online) Normalised third (left) and fourth (right) cumulants for the L × L × L × L × L-lattice, L = 4, 6, 8, 10, 12, 16, 20, 24 (blue) versus a for z = -9.87 together with the limit curve (red).

  the factor n 3/4 in front of the integral comes from the change of variables k = x n 3/4 .

8 .

 8 Let q = n+a n+2 and a > 0.

1 -µ 2 -1+µ 2 µ

 12 atanh µ « + 1-µ 2 2 µ 2 (2 µ-atanh µ) atanh µ ff π n (1 -µ 2 )

10 Page 10 of 47 http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

  Figure 8: (Colour online) Left: y = n (p -1) versus z = n (q -1) for C n , with n = 8, 16, 32, 64, 128 with larger n extending farther to the left. Right: K giving the maximum y vs 1/n for C n , n =8, 12, 16, 24, 32, 48, 64, 96, 128, 192, 256.

								y		K
										0.135
								1.0
										0.130
								0.8
										0.125
								0.6		0.120
								0.4		0.115
		15 F o 10		5	0.2	z	0.105 0.110	256	64	32	24	16	12	8	1 n
				r
					P
	1.06 1.08	y				e r e
	1.04 1.02 1.00	256	64	32	24	16	12	R e e w i v 8 1 n 16 24 32 64 256 3.6	12	8	1 n
														O n
														l y
														However, numerical
									27

http://mc.manuscriptcentral.com/pm-pml Philosophical Magazine & Philosophical Magazine Letters

Acknowledgments

One of the authors (AR) wishes to thank the Swedish Research Council (VR) for financial support. This research was conducted using the resources of High Performance Computing Center North (HPC2N).

A The special case p = q Recalling (7) we observe that

It is then trivial to show the next lemma.

Lemma A.1. For q = n n+2 we have R q,q (n, n/2, 1) = 1. 41