Logarithmic scaling for height variables in the Abelian sandpile model

We report on the exact computation of the scaling form of the 1-point function, on the upper-half plane, of the height 2 variable in the two-dimensional Abelian sandpile model. By comparing the open versus the closed boundary condition, we find that the scaling field associated to the height 2 is a logarithmic scalar field of scaling dimension 2, belonging to a c=-2 logarithmic conformal field theory. This identification is confirmed by numerical simulations and extended to the height 3 and 4 variables, which exhibit the same scaling form. Using the conformal setting, we make precise proposals for the bulk 2-point functions of all height variables.


I. INTRODUCTION
Sandpile models are open dynamical systems which generically show a wide variety of critical behaviours [1,2]. In two dimensions, one may be tempted to analyse the stationary state of a sandpile in terms of a local conformal field theory. As most sandpile models have a number of unusual features, the underlying conformal theories are expected to be peculiar as well. It is now established that at least some of them fall in the class of logarithmic conformal theories (LCFT) [3,4].
The Abelian sandpile model is most naturally formulated in terms of height variables. These random variables, sitting at the sites of a lattice, take the values 1, 2, 3 and 4, and evolve under a discrete time dynamics. When the lattice is finite, there is a unique stationary measure, which is uniform on its support, formed by the so-called recurrent configurations [6]. If the continuum limit of this measure is a c = −2 LCFT measure, each of the four height variables should converge to a definite conformal field in the scaling limit, in such a way that the multisite probabilities are reproduced by conformal correlation functions.
Early on, fundamental differences among the four height variables emerged. The Bombay trick [7] allowed Majumdar and Dhar to handle the height 1 variable rather easily, and led to the clear identification of its scaling form with a primary field φ(z,z) of dimensions (1,1) [11]. Technically the other three height variables are harder. To date, the only known result, due to Priezzhev, is an exact integral formula for the 1-site probabilities P 2 , P 3 , P 4 on the discrete plane Z 2 (or infinitely far from boundaries) [18]. As these numbers must be subtracted from the random variables before making contact with the conformal fields (which have a zero vacuum expectation value on the plane), they give no information on the fields themselves. Even though different arguments have suggested a genuine difference between the height 1 and the heights 2, 3 and 4 in the bulk, the exact scaling form of the latters remains an open question.
In this Letter, we provide an answer to this problem. The results which follow show that the height 2 scales like a logarithmic conformal field ψ(z,z) of dimensions (1,1) -the logarithmic partner of the primary field φ-, and relying on numerical simulations of the lattice model, we argue that this is also the scaling form for the heights 3 and 4. Here we content ourselves with giving the main results of the lattice calculations and simulations, and their comparison with the LCFT predictions. The details of the analysis will be presented elsewhere.

II. THE HEIGHT 2 ON THE UPPER-HALF PLANE
Information about the scaling form of the height 2 variable requires the knowledge of some of its correlations. The 1-site probability P 2 (z) on the upper-half plane is a disguised correlator of the height 2 with itself (its mirror image), and gives the simplest way to assess its scaling behaviour. It can be computed in essentially the same way the 1-site probability P 2 = lim z→i∞ P 2 (z) has been computed on the full plane by Priezzhev [18].
The height variables are the most basic microscopic variables, but they are strongly correlated over the whole lattice because they must make up recurrent configurations. A useful alternative is to treat the recurrent configurations as a global random variable, uniformly distributed, and in turn, the recurrent configurations can be traded for the set of spanning trees on the lattice one considers [8]. While the spanning tree picture may be avoided for calculations with heights 1, it seems to be unavoidable for the height 2 (and 3 and 4).
Being interested in height observables, one should translate the height values at certain sites into well-defined characteristics of the spanning trees. This can be done generally, but for what follows, we only discuss the case of a single height 2. In a spanning tree, any site belonging to the branch(es) grown from a site i is called a predecessor of i. It can then be shown [18] that a given recurrent configuration has a height 2 at site i if and only if, in the corresponding spanning tree, the site i has exactly one predecessor among its nearest neighbours. The distribution on spanning trees being uniform, one has to count the spanning trees which satisfy this (non-local) characterization at i.
Priezzhev has devised a method to compute this number [18]. It involves counting the trees with a certain local property, and also graphs with some other non-local properties. The former can be done by standard methods of graph theory, like the introduction of a defect matrix, the calculation then boiling down to a small determinant. The latter calculation requires to enumerate the graphs with a loop (a non-local constraint but easily computed using a defect matrix) and also what Priezzhev called Θ-graphs, which are graphs with an imposed Θ-shaped circuit. Both non-local contributions are divergent, but the combination which is relevant to the site probabilities is finite.
We have used this technique to compute the probability to have a height 2 somewhere on the upper-half plane. With respect to the plane, this situation brings a few changes.
The most obvious change is the presence of the boundary, on which boundary conditions must be prescribed. Two natural -and the only ones known-boundary conditions may be imposed along the boundaries, open (Dirichlet) and closed (Neumann). Their precise definition is not important for what follows and can be found for instance in [16]. More important is the fact that the operator which converts an open boundary condition into a closed one, or vice-versa, has been identified as a c = −2 primary field µ(z) of weights − 1 8 [12]. It will be a crucial ingredient in our proof that the height 2 variable scales like the logarithmic field ψ. Indeed the probability that a site have height 2 will be computed for both an open and a closed boundary condition, but the two are not independent since they can be related to each other by the use of the field µ. This is a non-trivial consistency check that involves the 4-point function with two µ's and two ψ's.
The counting of Θ-graphs is also a little bit more complicated on the upper-half plane than on the full plane. The rotational invariance is lost and the translational symmetry is broken in one direction. Moreover the way the Θ-circuit can touch the boundary has to be controlled. These circumstances increase the number of different graph configurations that have to be computed separately, and complicate the expressions.
The expressions we have obtained for the probabilities using the Θ-graph technique are rather cumbersome, to the extent that reproducing them here is unlikely to be helpful. The local and the loop contributions are simple but those coming from the Θ-graphs are fairly long. They involve six-fold integrals and a summation over the sites of the upper-half plane. The horizontal translational invariance enables one to carry out the infinite sum over the horizontal coordinate, which in turn allows to do one of the six integrations (on the full plane, the two infinite summations can be done, and then two of the six integrations, leaving a four-fold integral [18]). The integrand/summand involves the determinants of 4-by-4 matrices, and contain a rather complicated dependence in m, the distance to the boundary of the reference site.
The asymptotic analysis of these expressions is long and technical, but yields the asymptotic form of P op 2 (m) and P cl 2 (m), the probability that a site at a distance m from the boundary has a height equal to 2, when the boundary is all open or all closed. We have obtained the following results, The height 2 probability P 2 on the full plane is known exactly as a multiple integral, whose numerical evaluation yields P 2 = 0.1739 [18]. The constant a is also given in terms of integrals which we could not to evaluate analytically, but a numerical integration gives a = 0.0403. As to the other constant, we found the exact value b = 1 4 P 1 = 1 2π 2 (1 − 2/π), with P 1 = 0.07363 the probability that a site in the plane has height 1 [7].
These results give a clear signal that the scaling field describing the height 2 variable has scale dimension 2 and has a logarithmic component. The natural candidate in the c = −2 LCFT is the field ψ, namely the logarithmic partner of the primary field φ, with conformal dimensions (1, 1).
If ψ is indeed the scaling field of the height 2, the above dominant terms for the subtracted probabilities should be equal to ψ(z, z * ) op and ψ(z, z * ) cl with z − z * = 2im. The non-chiral 1-point functions on the upper-half plane are also equal to the chiral correlators on the full plane, of two chiral ψ inserted at z and z * , provided the coefficients defining the most general chiral 2-point function are properly chosen [19].
In order to be able to relate the open and the closed boundary conditions, we consider the more general expectation value of ψ on the upper-half plane, for which the boundary condition on the real axis is all open except on the interval I = [z 1 , z 2 ] where it is closed. It can be expressed as a chiral 4-point function on the plane as where N µ = 1.18894 normalizes the lattice 2-point functions of the fields µ [12]. In the sandpile model, it yields the (scaling limit of the) probability that the site z has a height equal to 2 in presence of a boundary which is open on R \ I and closed on I. The two lattice results (1) and (2) are then easily recovered upon taking z 1 , z 2 → 0 and −z 1 , z 2 → +∞ respectively. In these limits, one may invoke the operator product expansion of two µ's, which closes on different fields in the two cases, since the remaining boundary condition is different [16], where I D and I N are the identity fields on an open resp. closed boundary and ω N is the logarithmic partner of I N .
Thus we obtain the two expectation values on the upper-half plane as Because the primary field µ is degenerate at level 2, the chiral correlator µ(z 1 )µ(z 2 )ψ(z 3 )ψ(z 4 ) satisfies a secondorder differential equation, but the logarithmic nature of ψ complicates the matter. It is the logarithmic partner of a primary field φ of dimensions (1,1), into which it transforms under a dilatation, and it is not quasi-primary, because it transforms into primary fields ρ(z,z) andρ(z,z) of dimensions (0,1) and (1,0) under special conformal transformations [20]. The operator product expansion of ψ with the chiral stress-energy tensor reads As a consequence, the differential equation satisfied by µ(z 1 )µ(z 2 )ψ(z 3 )ψ(z 4 ) is inhomogeneous [3], where the inhomogeneity depends on the four other 4-point functions µµρψ , µµφψ , µµψρ and µµψφ . These 4-point correlators themselves satisfy inhomogeneous differential equations, and depend on yet four other correlators, µµρρ , µµφρ , µµρφ and µµφφ , which involve primary fields only, and hence may be computed independently, as solutions of homogeneous differential equations. Altogether, there is a sequence of nine correlators to compute, the last one being µµψψ . The five correlators which contain a φ satisfy a first order differential equation, because the φ is also degenerate at level 2, so each depends on one arbitrary coefficient. The other four satisfy second order equations, and depend on two arbitrary coefficients. Thus the general 4-point function µ(z 1 )µ(z 2 )ψ(z 3 )ψ(z 4 ) depends on 13 arbitrary coefficients.
To fix their values, we note, from the operator product expansion (4) in the limit z 1 , z 2 → 0, that none of the nine correlators mentioned above can have a logarithmic singularity, and that the resulting function ψ(z, z * ) op must be a function of z − z * only. These two facts together fix 10 coefficients. One more coefficient is fixed if one imposes that the limit −z 1 , z 2 → +∞ remains regular, i.e. the logarithmic singularity is absent. As is clear from (6), this implies that ψ(z, z * ) cl is related, through the image method, to a chiral field ψ satisfying ψ(z)ψ(z * ) = 0. The same is true for the height 1 variable and its field φ [16].
All this leaves just two free parameters, themselves determined by demanding that the limit z 1 , z 2 → 0 of (3) reproduces the lattice result (1). It leads to our final result for the expectation value of ψ on the UHP with a boundary closed on [z 1 , z 2 ] and otherwise open, where x = z12z34 z13z24 is a cross-ratio of the four insertion points (z 3 = z * 4 = z). One readily checks that in the limit z 1 , z 2 → 0, x → 0 and √ 1 − x → +1, the previous formula reduces to (1). To make the connection with the probability of a height 2 in presence of a closed boundary, one takes −z 1 , z 2 → +∞. Again x goes to 0, but √ 1 − x now tends to −1 [16]. Indeed if we set z 1 = −R and z 2 = R, we find that has unit modulus and draws a complete circle around 0 when R varies from 0 to +∞. This non-trivial monodromy factor brings a global change of sign in (8) with respect to the previous case (z 1 , z 2 → 0), and the non-logarithmic b terms drop out inside the curly brackets (the last term does not contribute in either case). The resulting 1-point function coincides exactly with the lattice result (2), thereby confirming in a non-trivial way the identification of ψ as the scaling field for the bulk height 2 variable in the sandpile model. We will add further support in the next section.
We finish this section by giving the expectation value of ψ in the upper-half plane, in presence of a boundary closed on R − and open on R + , a formula we will use in the next section. In this case, z 1 → −∞, z 2 → 0, and a simple calculation from (8) yields

III. SIMULATIONS ON A STRIP
We have so far verified the identification of the height 2 variable with the field ψ in the two extreme cases, namely in the upper-half plane with a homogeneous boundary condition on the real line, either all open or all closed. It may be tested more deeply by seeing how the profile of the height 2 probability varies in a region where the boundary condition changes. The formula (9) corresponds to such a situation, with half the real line open and the other half closed, and would describe the way the probability P 2 (z) is interpolated between its open boundary value and its closed boundary value.
We could however not carry out the lattice analysis for this case as we did for the homogeneous boundary condition, so that the exact lattice result is not available. Thus we turned to numerical simulations of the model to measure the lattice profile of P 2 (z), with the extra advantage that the two other probabilities P 3 (z) and P 4 (z) can be obtained with no more effort.
Clearly the simulations can only be done on a finite array, in our case, a rectangle of base M and height N . We have imposed an open boundary condition on the bottom boundary, a closed boundary condition on the top boundary, and we have evaluated the four probabilities P 1 (n), P 2 (n), P 3 (n) and P 4 (n) on the vertical line joining the bottom to the top boundary, and located at mid-length (we have included P 1 (n) for completeness and for comparison purposes). For these results to be significant (with respect to the CFT approach), both the height N and the ratio τ = M/N must be large enough so that the scaling regime is approached, and so that the rectangle can be considered as a good approximation of the infinite strip. We found that the values N = 50 and τ = 4 were sufficient, by verifying that larger values did not bring any noticeable changes in the curves. To speed up the calculations, we chose open boundary conditions on the left and right boundaries.
We have determined the 1-site probabilities P i (n) on a sample of nearly 17 · 10 9 recurrent configurations, generated from 1680 random recurrent configurations in the following way. Applied to each of these, the dynamics of the model produces a new recurrent configuration at each time step. The first 10 5 ones are not sampled, and from then on, every 50th configuration is sampled, for a total of 10 7 sampled configurations for each initial configuration. At each point n of the line joining the two boundaries, the probabilities are the ratios of the number of occurrences by the size of the sample. From these numbers, we should subtract in each case the probability P k in the bulk (far from the boundaries) before making the comparison with the field-theoretic results. The most distant sites from the boundaries are the central sites, but because the strip is finite, the lattice probabilities at the central sites have not yet attained their limit value P k . So we have subtracted an effective probability P eff k , computed by a fitting procedure. The results are shown as colour dots on the left plots of Figures 1 and 2.
Despite the fact that the size of the sample (∼ 10 10 ) is desperately small compared to the volume of the phase space (∼ 10 5 045 !), the fluctuations in the data are surprisingly small, and follow a normal distribution to a good approximation. Based on an analysis of these, we expect errors less than 2% for P 1 (n), less than 1% for P 2 (n) and P 3 (n), and less than 0.5% for P 4 (n). The precision naturally increases for more frequent events (the bulk values are P 1 = 0.0736, P 2 = 0.1739, P 3 = 0.3063, P 4 = 0.4461 [18]).
The results pertaining to the height 2 probability P 2 (n) may be directly compared with the probability profile on the infinite strip obtained from the LCFT picture, on the basis that the height 2 variable goes to the field ψ in the scaling limit. For that, the corresponding curve on the upper-half plane is simply transformed to the infinite strip, of width L, by the conformal transformation w = L π log z. By integrating the infinitesimal transformation of ψ given by (7), one finds its finite transformation under z → w, Under the above mapping z = exp (πw/L) of the upper-half plane onto the infinite strip, and taking expectation values, one obtains ψ(w,w) strip in terms of ψ(w,w) UHP , φ(w,w) UHP , ρ(w,w) UHP and ρ(w,w) UHP . The  (12) and (13) to which they can be directly compared. On the right, the same plots, but multiplied by the factor L π 2 sin 2 (πn/L) cos(πn/L) , reveal the presence or the absence of a logarithmic dependence. expectation of ψ is given in (9), to which that of φ is related through a dilatation, whereas the expectations of ρ and ρ can be shown to vanish, Putting all together, we obtained the following formula for the expectation value of ψ on an infinite strip of width L, and coordinates w = u + iv, ψ(w, w * ) op,cl strip = π L 2 cos(πv/L) sin 2 (πv/L) The open boundary is at v = 0 and the closed one at v = L. The height 1 variable scales exactly to the primary field φ, normalized as in (11) [16] (that is why we introduced the factor − 1 2 in (7)). So one quickly obtains φ(w, w * ) op,cl strip = b π L 2 cos(πv/L) sin 2 (πv/L) .
The two functions of v in (12) and (13) have been plotted in the left graph of Figure 1, as solid lines. In both cases, we have chosen L = N + 1 2 , for the following reason. In the continuum, v = L is the position of the closed boundary, on which the Neumann condition is enforced. On the lattice, this condition is imposed through a mirror symmetry about the line n = N + 1 2 , so that it is at that point that the normal derivative vanishes. One sees that the conformal curves are very close to the profiles obtained from numerical simulations, adding much support to the identification of ψ as the scaling field for the height 2 variable. Because N = 50 is relatively small, the plots with L = N would be slightly shifted near the closed side, resulting in a visible discrepancy at the right end of the plots. The gross features of the two curves are essentially the same and do not reveal any qualitative difference between the height 1 and the height 2, namely the logarithmic nature of the latter. In order to make the difference to appear more clearly, we have plotted, in the right graph of Figure 1, the same curves and the same data but divided by the damping prefactor π L 2 cos(πn/L) sin 2 (πn/L) , which tends to make the function vanish precisely for the values of v where the logarithm becomes important. One then sees that the plot for the height 1 is flat, while the other one shows a non-trivial dependence, and agrees well with the simulations (the discrepancy in the central region is due to the difference between the true bulk value of P 1 and P 2 and the effective values used in the subtraction, a difference which is amplified by the factor L 2 /(π 2 cos (πv/L)); the same factor also amplifies any disparity between the dots and the curve). We note that for the height 2, the curve is asymmetric around the middle point, a consequence of the extra factor b 2 in the probability near an open boundary, see (1) and (2). The graphs in Figure 2 concern the height 3 and height 4 probabilities. As before the dots in the left plot represent the subtracted probabilities obtained from numerical simulations. Suspecting that the height 3 and 4 variables scale the same way as the height 2, we have fitted these data with the expressions they would follow if that were true, namely the expression (12). In other words, we have assumed that the probability profiles for the heights 3 and 4 across the finite discrete strip of width N = L − 1 2 have the form The fit of the three parameters, performed on the 46 most central values of n, yields the following values in each case, By comparison, a fit of the data for the height 2 yields b 2 = 0.01896 for the parameter controlling the logarithmic dependence (the exact value is b = P1 4 = 0.01841), and b 1 = −0.000388 for the height 1 (the exact value is 0). Using the fitted values of the parameters, the curves (14) are plotted on the left graph of Figure 2, as solid lines. The rather large values of b 3 and b 4 , as well as the good agreement between the solid curves and the data, are a very strong indication that the height 3 and 4 variables possess the same logarithmic scaling as the height 2 variable. The plots in the right part of Figure 2 bear the same relation to the left plots as those in Figure 1, and exhibit a clear logarithmic signature.
If the height 3 and 4 variables scale like the height 2 variable, their associated fields must be ψ i = α i ψ + β i φ, for i = 3, 4. Using the above formulae for the expectation values and the fitted values in (15) and (16), we obtain the relations b i = bα i and a i = aα i + bβ i , from which we deduce the values α 3 = 2.07, β 3 = −5.21, α 4 = −3.08, β 4 = 4.22.
The numerical values are consistent with the sum being zero.

IV. HEIGHT CORRELATIONS IN THE BULK
The results we obtained in the previous sections can finally be used to make definite statements about the height correlations on the plane. As far as we know, and apart from those involving heights 1 exclusively, these correlations have never been computed nor measured from simulations.
The correlation of two distant heights h i , h j , located at z 1 , z 2 on the plane, is given by the 3-point function of the corresponding fields with the field ω, the logarithmic partner of the identity, inserted at infinity, For this one needs the 3-correlators φφω , φψω and ψψω , the last one anyway requiring the other two.
Interestingly, the constant A governs the long distance behaviour of all 2-site height correlations, and is known form the calculation of the 2-site probability for heights 1, A = −P 2 1 /2 [7]. Granting the field assignments made above for the heights 3 and 4, we conclude that the leading terms of the 2-site probabilities on the plane are given by