Sandpile Models in the Large

This contribution is a review of the deep and powerful connection between the large-scale properties of critical systems and their description in terms of a field theory. Although largely applicable to many other models, the details of this connection are illustrated in the class of two-dimensional Abelian sandpile models. Bulk and boundary height variables, spanning tree–related observables, boundary conditions, and dissipation are all discussed in this context and found to have a proper match in the field theoretic description.

1 Introduction

In statistical mechanics, critical points are very special points in the space of external parameters which control the state of a system. At such a point, the system is scale-invariant, and its thermodynamic functions and correlations are characterized by critical exponents and power laws. In many cases, physical systems have a finite number of critical points, most often only one. Typical examples include the endpoint of the liquid–gas coexistence line or the Curie point for ferromagnetic materials. In these cases, a system is brought to its critical point by tuning very precisely a few external parameters to their critical values.

In nature, however, power laws are commonplace and can be found in a large variety of different phenomena, like avalanches, earthquakes, solar flares, droplet formation. In all these cases, it is certainly not clear what parameters should be tuned, and even if they are perfectly tuned, it is unlikely that they would stay so over large periods of time. To solve this apparent paradox, Bak, Tang, and Wiesenfeld suggested in the 80s that the external parameters would tune themselves dynamically: even if the system is not initially in a critical state, its own dynamics will ineluctably drive it to criticality and maintain it in that state [1]. This attractive idea has led to the concept of self-organized criticality (SOC).

To support this idea, these authors proposed the sandpile model as a prototypical example of a system which shows a form of self-organized criticality. Since then, many other models showing SOC have been proposed, as abundantly illustrated in this volume and in introductory books and reviews [2–5].

The present review will be exclusively concerned with specific versions of two-dimensional sandpile models, formulated by Dhar [6], called Abelian sandpile models. Even though there are among the simplest and easiest sandpile models to handle, they show a large spectrum of interesting and difficult problems which have attracted considerable attention, in both the physical and mathematical communities. From the point of view taken here (like their scaling limit and the emerging conformal field theory), they are, to our knowledge, the only ones to have been studied. Yet, compared to many other equilibrium statistical models, a fair statement is that our present understanding of them is still very poor.

Our primary purpose is two-fold, namely, to give the unfamiliar reader an introduction of why and how the neighborhood of a critical point can be described by a Euclidean field theory, which, at first sight, appears to be a rather obscure statement, and also to show how this description can be worked out in practical terms. The second part will be illustrated in sandpile models, which lend themselves very well to this kind of analysis: they are simple enough that one can follow the steps in a clear and transparent way, yet they are rich enough to show the difficulties one sometimes has to face but also the elegance and the power of the approach. Understanding how of a field theory emerges from a stochastic lattice model enables to gain a probabilistic and intuitive view of what a field theory is in this context.

It turns out that the field theories which appear when analyzing critical systems are conformal field theories. The simple reason for this is that their large conformal symmetry integrates the fact that critical systems have a local scale invariance. Conformal theories in two dimensions have been tremendously successful since the 80s and have led to a deep understanding of the two-dimensional critical phenomena. It is certainly not our purpose to give an introduction to conformal field theories, and we will not go very deep into its technicalities, referring to the vast literature. We restrict to their most basic features, in the hope that these will be sufficient and useful to understand how conformal theories are so well suited for our study.

Section 2 starts with a brief review of the Abelian sandpile models, where the most basic features of the models are recalled. Section 3 is a general description, valid beyond the sandpile models, of what is called the scaling limit, which allows establishing the connection between the large-distance regime of a critical system and the associated field theory. A brief tour of conformal theories, and specifically logarithmic conformal theories, is presented in Section 4. The application of the conceptual ingredients is illustrated in the next three sections. Section 5 focuses on the bulk observables in the sandpile models, computes the first correlators, and explains how these should be understood in terms field theoretic quantities. Boundary conditions and boundary observables are examined in Section 6 as well as the way they should be thought of in conformal theories. Section 7 discusses a dissipative variant of the sandpile models and their description by a massive field theory, and also some universality aspects of the sandpile models. The last section summarizes the present status of the conformal theory at work in sandpile models.

The present text has some overlap with [7]. The latter was more concerned with the sandpile models as being described specifically by a logarithmic conformal field theory. Intended to a potentially wider readership, the present review is more devoted to the general connection between critical systems and field theories, illustrated in a specific class of models. The two are somehow complementary and, if combined, may provide a more complete overview.

2 Abelian Sandpile Models

The models we discuss are discrete stochastic dynamical systems. Their microscopic variables are attached to the vertices of a finite connected graph Γ = (V,E) (with V the set of vertices, or sites, and E the set of simple, unoriented edges) and evolve in discrete time as a random process. We label the vertices of Γ by Latin indices i, j,… and denote the microscopic variables by h_i. These are called height variables and simply give the height of the sandpile at vertex i (i.e., count the number of sand grains at i); they are integer-valued, with $h_i\mathrm{⩾}1$. A height configuration C is a set of height values ${\left\{h_i\right\}}_{i\in V}$.

We are not quite ready to define the dynamics. For reasons that will become clear in a moment, we need to extend Γ by adding one special vertex, noted s and called the sink, as well as a number of edges connecting s to the vertices of a non-empty subset $D\subset V$. Vertices in D are called dissipative or open, while those in $V\mathrm{\setminus }D$ are conservative or closed. If ${\mathrm{\Gamma }}^{*}=\left(V^{*},E^{*}\right)$ denotes the extended graph in an obvious notation, we define z_i to be the coordination number of i in Γ (the number of edges in E incident to i, or the number of its nearest neighbors in Γ) and, similarly, $z_i^{*}$ its coordination number in ${\mathrm{\Gamma }}^{*}$. Thus, $z_i^{*}=z_i$ if i is closed and $z_i^{*}>z_i$ if i is open. Finally, we say that a site i of V is stable¹ if its height satisfies $1\mathrm{⩽}{h}_i\mathrm{⩽}{z}_i^{*}$. A height configuration is stable if all sites are stable. Clearly, the number of stable configurations is equal to ${\displaystyle {\prod }_{i\in V}{z}_i^{*}}$.

The discrete, stochastic dynamics of the sandpile model is defined as follows. Assume that C_t = {h_i} is a stable configuration at time t. The stable configuration C_t+1 is obtained from the following two steps.

1) Deposition: One grain of sand is dropped on a random site j of V, selected with probability p_j, producing therefore a new configuration C^new with heights h_i^new = h_i + δ_i,j. If $h_j^{\text{new}}\mathrm{⩽}{z}_j^{*}$, then C^new is stable and defines C_t+1; if not, we proceed to step (ii).

2) Relaxation: If $h_j^{\text{new}}>z_j^{*}$ (it is in fact equal to $z_j^{*}+1$), we let the site j topple: its height is decreased by $z_j^{*}$, each of its neighbors in Γ receives one grain, and the remaining $z_j^{*}-z_j$ grains go to the sink. After this, one or more neighbors of j in Γ may become unstable, in which case they topple in the way explained above for the site j. The toppling process is pursued for all unstable sites until a stable configuration is obtained. That configuration defines C_t+1.

It is useful to introduce the toppling matrix Δ as it will play an important role in what follows:

${\mathrm{\Delta }}_{i,j}=\{\begin{array}{ll}{z}_i^{*}\hfill & \text{for}\hspace{0.17em}i=j,\hfill \\ -1\hfill & \text{if}\hspace{0.17em}i\hspace{0.17em}\text{and}\hspace{0.17em}j\hspace{0.17em}\text{are}\hspace{0.17em}\text{neighbours}\hspace{0.17em}\left(\text{i}.\text{e}.\text{connected}\right),\hfill \\ 0\hfill & \text{otherwise},\hfill \end{array}(1)$

for $i,j\in V$. The sand redistribution occurring when a site j topples can then be written as the update h_i→h_i−Δ_j,i for all $i\in V$. The matrix Δ is like a Laplacian on Γ, with mixed boundary conditions dictated by the open and closed sites, which induce, respectively, Dirichlet and Neumann boundary conditions (see Section 6).

The above dynamics is well-defined. We see that the total number of sand grains is conserved under the toppling of a closed site, whereas a nonzero number of grains are transferred to the sink under the toppling of an open site. The existence of at least one open site guarantees that the relaxation process terminates after a finite number of topplings and motivates the necessity of the extension of the graph Γ by the sink site. Moreover, if several sites are unstable during the relaxation process, the order in which they are toppled does not matter. More generally, one may define the operator a_i for each $i\in V$, whose action on a stable configuration returns the stable configuration resulting from the relaxation process after the deposition of a sand grain at i. One can then prove that the operators a_i commute [6], explaining the qualifier “Abelian” used to designate the models satisfying this property.

The dynamics described above is a discrete Markov chain on a finite configuration space: at each time step, one applies the operator a_i with probability p_i (it is the only stochastic element of the dynamics), going from C_t to C_t+1 = a_iC_t. An important question concerns the invariant measures, since they control the behavior of the model in the long run.

If there is no strong reason to favor certain sites, one takes all probabilities p_i equal (uniform distribution). In this case,² Dhar [6] has shown that there is unique invariant measure ${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}$, which is uniform on its support. In the Markov chain terminology, the configurations in the support of ${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}$ are called recurrent; the others are called transient. Being in the support of the unique invariant measure means that the recurrent configurations are those which are in the repeated image of the operators a_i. The transient ones either never appear (depending on the initial configuration) or cease to appear after some finite time.

If indeed the unique invariant measure is uniform, the situation appears to be deceptively simple. Not so. What makes the sandpile models nontrivial, fascinating, and rich is the support of the invariant measure. A generic recurrent configuration is really complicated because the height values are delicately correlated over the entire graph. In the general case, there is no simpler criterion characterizing the recurrent configurations than the following. Let C be a stable configuration, and let C_F be its restriction to a subgraph $F\subset \mathrm{\Gamma }$ (F can be assumed to be connected). We say that C_F is a forbidden subconfiguration if each vertex of F has a height smaller or equal to the number of its neighbors in F. It can be shown [8] that a forbidden subconfiguration cannot be in the repeated image of the dynamics (of the operators a_i). It follows that a configuration is recurrent if and only if it contains no forbidden subconfiguration. The simplest example of a forbidden subconfiguration is when F contains just two neighboring vertices with height values equal 1. The criterion also implies that the maximal configuration with heights $h_i=z_i^{*}$ is recurrent since a vertex i with height h_i > z_i cannot be in a forbidden subconfiguration. It is also clearly in the image of the iterated dynamics since it can be reached from any other stable configuration by an appropriate sequence of a_i’s.

The characterizing condition for recurrence shows that the heights of a recurrent configuration are not at all independent. They are not only correlated locally (think of two neighboring 1’s) but also globally because asserting that a given configuration is recurrent generally requires scanning the entire graph. For instance, the configuration having h_i = z_i for all i is not recurrent and possesses no other forbidden subconfiguration than the whole configuration itself. Moreover, the recurrent status is very sensitive to local changes and can be lost or gained by the change of a single height (for the configuration just discussed, the increase by one unit of the height at a single open site makes it recurrent). However, the increase in any height in a recurrent configuration preserves the recurrence.

The burning algorithm [8] (see also the review [5]) provides a convenient way to test whether a given stable configuration is recurrent. In addition to providing a completely automatic procedure, more importantly, it establishes a bijection between the set of recurrent configurations on Γ and the set of rooted spanning trees on ${\mathrm{\Gamma }}^{*}$, rooted at the sink site s. Let us recall that a spanning tree is a loopless connected subgraph $\left(V^{*},F\right)\subset {\mathrm{\Gamma }}^{*}=\left(V^{*},E^{*}\right)$ with $F\subset E^{*}$. This bijection is important and useful as most of the actual calculations use the spanning tree formulation. Interestingly, there is no canonical bijection between the two sets in the sense that there are in fact many burning algorithms (the detailed definition requires a certain prescription that is largely arbitrary), each giving rise to a different bijection. This freedom in the choice of a definite algorithm, a sort of huge gauge symmetry, has remained unexploited so far.

If the notion of recurrence remains somewhat elusive in the generic case, simple arguments lead to a remarkably simple and general formula for the number of recurrent configurations [6], naturally identified as the partition function Z since the invariant measure is uniform:

$Z=\#\left\{\text{recurrent}\hspace{0.17em}\text{configs}\right\}=\text{det}\text{Δ},(2)$

for Δ the toppling matrix introduced in (1). It is a standard result in combinatorics (Kirchhoff’s matrix-tree theorem) that det Δ also counts the number of spanning trees on Γ (see Section 5.7 for a proof). The previous formula usually implies that the recurrent configurations form an exponentially small fraction of the set of stable configurations (whose number is equal to ${\displaystyle \prod }_i{\mathrm{\Delta }}_{i,i}$). On a large grid in ${\mathrm{\mathbb{Z}}}^2$, for instance, for which the density of dissipative sites goes to 0 in the infinite volume limit, the effective number of degrees of freedom per site in a recurrent configuration is roughly 3.21 (as compared to 4 in a stable configuration), meaning that $\text{det}\Delta \simeq {\text{e}}^{\frac{4G}{\pi }N}={\left(3.21\dots \right)}^N$, with N the total number of sites and G the Catalan constant.

The definition of recurrence implies that all the operators a_i map recurrent configurations to recurrent configurations, implying that once the dynamics has brought the sandpile into a recurrent configuration, all subsequent configurations are recurrent. Therefore, the invariant measure is appropriate to study the long-term behavior of the sandpile.

The sandpile models summarized above have raised a large number of interesting and difficult questions. In the context of this review, most if not all of them focus on the stationary regime and study the statistical behavior of the sandpile when it runs over the recurrent configurations. In other words, all the probabilities we are interested in are induced by the invariant measure ${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}$. The use of ${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}$ is what makes most of the calculations fairly hard³ because as noted earlier, that measure is nonlocal in terms of the (local) height variables (equivalently the recurrence criterion is nonlocal).

We should remark that the measure ${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}$ fully depends on all the minute details which are necessary in order to completely specify the sandpile model under study. Not only the graph Γ itself but also the number and relative positions of closed and open vertices, and the values of the local thresholds $z_i^{*}$ affect the invariant measure. Many features which directly depend on these data will change if any of these parameters is modified, like the number of recurrent configurations, the structure of the sandpile group,⁴ the geometric structure of the identity configuration,⁵ or the average height at a given site for instance. All these features are mathematically interesting and challenging (hence interesting) but very sensitive to the underlying details.

One should however expect that more robust features would be shared by sandpile models that are “close enough.” The same situation prevails for other statistical models which, although having different microscopic descriptions, are considered to be essentially equivalent and grouped together to form a single universality class. Models belonging to the same universality class have identical behaviors “in the large,” a point of view made technically more precise by the renormalization group analysis.

In order to identify these common behaviors, one should not look at small scales as these are more likely to be determined by the local details. The probability that two vertices next to each other have a height 2, for instance, is not really interesting; in addition, it is a pure number, different for each different model. Robust behaviors are expected to be found at large scales as they are much less affected by the microscopic details. One convenient method to access the large-distance behaviors is by taking the scaling limit. (Readers familiar with the scaling limit and the ideas of the renormalization group can safely go straight to the next sections.)

3 The Scaling Limit and Continuum Field Theories

The simple idea underlying the scaling limit is this: if we want to concentrate on the large-scale behavior of a system, let us look at it from far away! The further away we look at the system, the larger our horizon is and the larger the distances we keep in sight. At the same time, when looking from a distance, the details get blurred and disappear: one can no longer recognize the type of graph, and its connectivities are no longer visible. What we see seems to become independent of the microscopic details of the model.

Rather than stepping back, an equivalent but more convenient way to proceed is to shrink the discrete structure (graph or grid or lattice) on which the microscopic variables live. This will involve a (real) small parameter ε such that the graph can be embedded in $\epsilon {\mathrm{\mathbb{Z}}}^d$ (or another shrinked regular lattice). For smaller and smaller ε, fixing a macroscopic distance $\overrightarrow{x}=\epsilon \overrightarrow{m}\in \epsilon {\mathrm{\mathbb{Z}}}^d$ amounts to probe larger and larger scales $\overrightarrow{m}$ in terms of lattice units and at the same time, allows to keep a macroscopic distance $r=\left|\overrightarrow{x}\right|$ under control. The scaling limit corresponds⁶ to take ε → 0.

We note that since the scaling limit is a way to focus on asymptotically large distances, we have to make sure that the system does have such asymptotic distances! Indeed the scaling limit requires that we also take the infinite volume limit, by allowing the system to remain finite but of increasing size, the growth being at least of order 1/ε.

The scaling limit has interesting consequences. The first most apparent one is that the substrate of the rescaled model goes to a continuum, either ${\mathrm{ℝ}}^d$ or a part of ${\mathrm{ℝ}}^d$, which may be bounded.⁷ This is the first sign that a continuum description ought to emerge in the scaling limit. This is confirmed by a second observation: the microscopic variables—the heights in the sandpile models—, which were attached to the vertices of a graph, or a grid, should, in some sense, converge to variables defined on a continuum. If indeed this is expected to happen, exactly what happens is quite subtle. To realize this, one may note that all the microscopic variables attached to sites contained in a ball of radius o(1/ε) will actually collapse to the same point in the scaling limit. Thus, every point in the continuum is the convergence point of an infinite number of vertices in the original discrete setting. The infinity of microscopic variables carried by these vertices will supposedly mix and fuse to generate some kind of degree of freedom located at a single point in the continuum. What is then the nature of the emerging continuum degree of freedom at that point, and how is it related to the lattice variables supposed to collectively generate it? The conceptual answer is provided by the renormalization group. It roughly goes as follows.⁸

The scaling limit, as explained at the beginning of this section, was carried out in one stroke: all distances are scaled by ε, which is then taken to 0. This limit was only designed to show how the large-distance behaviors can be assessed but is too rough to answer the question raised in the previous paragraph. The renormalization group is much better designed conceptually as it organizes the scaling limit scale by scale and keeps track, at each scale, of the degrees of freedom present in the system.

Let us suppose that we start with a statistical model defined on very large graph, or, to simplify and fix the ideas, on an infinite lattice. We fix a convenient scale Λ > 1, partition the lattice into boxes of size Λ, and shrink the lattice by a factor Λ. Each box is now of linear size 1 and contains of the order of Λ^d microscopic variables. Within each box, we associate an effective, coarse-grained degree of freedom which takes into account the overall behavior of the microscopic variables inside the box (it could be, f.i., their average value), and we then compute the sum over the microscopic variables conditioned by the values of the coarse-grained variables. The result is a statistical model for the coarse-grained variables, defined on a lattice similar to the original one. Once this is done (!), we iterate the process by defining a second generation of coarse-grained variables out of those of the first generation, and so on.

After the first iteration, each group of roughly Λ^d microscopic variables has collapsed to a single coarse-grained variable of the first generation; the statistical model obtained for these can be interpreted as the original model in which the fluctuations of scale smaller than Λ have been integrated out. The second iteration yields a statistical model for the coarse-grained variables of the second generation, each of which has integrated the fluctuations of Λ^2d microscopic variables over scales smaller than Λ², and so on, for the next iterations. In this way, each iteration, also called renormalization, yields a model where more small-scale fluctuations have been integrated out, and whose large-scale behavior should be identical to that of the original model, since the large-scale fluctuations have been preserved.

The continuum degrees of freedom we were asking about are what the coarse-grained variables of higher and higher generation should converge to when the number of iterations goes to infinity. Each of them is indeed what is left of the infinite collection of the microscopic variables that were located around it. Because the coarse-grained variables of one generation are representative of those of the previous generation, the continuum degrees of freedom should similarly carry the same characteristics as the original microscopic variables. In particular, the long-distance correlations should be identical, at dominant order.

The continuum degrees of freedom emerging in the scaling limit are called fields. Unlike their lattice ancestors, they usually take continuous values. Fields are all what remains when the short-ranged degrees of freedom have been integrated out: they form the complete set of variables which are relevant as far as the long-distance properties of the original model are concerned. It means that only the lattice degrees of freedom which have long range correlations, namely, with diverging correlation lengths, will survive the scaling limit and eventually give rise to a field; all the others progressively disappear in the renormalization process.

The microscopic variables in terms of which the discrete statistical model is defined usually give rise to fields, but they are not the only ones. Any lattice observable, that is, any function of the microscopic variables, can potentially give rise to a field in the scaling limit⁹ so that one is typically left with an infinite number of different fields. Each field has its own specific properties and should be interpreted as the scaling limit of one particular lattice observable (it may also happen that different lattice observables converge to fields with the same characteristics).

One last question must be addressed. The original statistical model was not only defined by its microscopic variables but also by a probability measure on the configuration space. That measure, which is a joint distribution for the (nonindependent) random microscopic variables, is usually given by a Gibbs measure and written, up to normalization, as $\mathrm{\mathbb{P}}\left(C\right)\sim \text{exp}\left(-H\left[C\right]\right)$, where H is the Hamiltonian of the system, that is, some given function of the microscopic variables which determines the relative probability of a configuration C. What is the equivalent of the Gibbs measure for the fields?

According to the discussion above, one starts from the original model and its Hamiltonian H₀ ≡ H. The first renormalization yields the coarse-grained variables of the first generation and a corresponding Hamiltonian H₁, computed (at least in principle) by summing exp(−H₀) over the microscopic variables inside the boxes. Similarly, the kth iteration will produce a Hamiltonian H_k, defining the statistical model for the coarse-grained variables of the kth generation. The appropriate measure for the fields should therefore be something like the formal limit lim_k→∞exp(−H_k). Physicists like to denote this formal object by exp(−S), where S, called the action, is a certain functional of the fields.

Thus, if the description of a statistical model is given, in the discrete lattice setting, in terms of a set of microscopic variables (h¹_i, h²_i,…) and a Hamiltonian H(h¹_i, h²_i,…), it is given in the scaling limit by a set of continuous fields $\left({\varphi }_1\left(\overrightarrow{x}\right),{\varphi }_2\left(\overrightarrow{x}\right),\dots \right)$ and an action $S\left[{\varphi }_1,{\varphi }_2,\dots \right]$. The pair $\left\{\left({\varphi }_1\left(\overrightarrow{x}\right),{\varphi }_2\left(\overrightarrow{x}\right),\dots \right),S\right\}$ is referred to as a continuum field theory.¹⁰ More precisely, specifying a set of fields and their action S is only one way to present a field theory; it is also the most comfortable one because it allows to compute the correlators of the various fields, at least in principle.

Needless to say, working out the successive renormalizations along with the Hamiltonians H₀, H₁,… is a formidable task that is, for all practical purposes, impossible to carry out explicitly, except on extremely rare occasions (and for tailored examples). As a consequence, the field theory describing the large distances of a statistical model cannot be obtained in a deductive way.

The situation however is not hopeless. Experience, heuristic arguments, or results obtained on the lattice can often give definite hints about the nature of the seeked field theory. More importantly, and even if one has no clue of what the correct field theory is, the relevance of a trial field theory, perhaps suggested by an educated guess, can be firmly tested by comparing correlation functions. If the lattice microscopic variable h_i=$\frac{x}{\epsilon }$ (at site i) converges in the scaling limit to the field ϕ(x), it must be true that the scaling limit of the lattice correlators is equal to field theoretic correlators, as follows:

$\underset{\epsilon \to 0}{\text{lim}}\hspace{0.17em}{\epsilon }^{-n\mathrm{\Delta }}\hspace{0.17em}{\langle h_{\frac{x_1}{\epsilon }}\hspace{0.17em}{h}_{\frac{x_2}{\epsilon }}\dots \hspace{0.17em}{h}_{\frac{x_n}{\epsilon }}\rangle }_{\text{lattice}}={\langle \varphi \left(x_1\right)\hspace{0.17em}\varphi \left(x_2\right)\dots \hspace{0.17em}\varphi \left(x_n\right)\rangle }_{\text{FT}},(3)$

where the exponent Δ is determined so that the limits on the left hand side exist: as shown below, it will eventually be related to the scale dimension of the field ϕ to which the lattice variable h_i converges. The previous identity must be satisfied for all n-point correlators, but also for any correlator of any number of lattice observables provided that for each observable O(i) around the site i inserted in the lattice correlator, the corresponding field Φ(x) to which it converges is inserted in the field theoretic correlator:

$\underset{\epsilon \to 0}{\text{lim}}{\epsilon }^{-{\displaystyle \sum }_i{\mathrm{\Delta }}_i}{\langle O_1\left(\frac{x_1}{\epsilon }\right)O_2\left(\frac{x_2}{\epsilon }\right)\dots O_n\left(\frac{x_n}{\epsilon }\right)\rangle }_{\text{lattice}}={\langle {\text{Φ}}_1\left(x_1\right){\text{Φ}}_2\left(x_2\right)\dots {\text{Φ}}_n\left(x_n\right)\rangle }_{\text{FT}}.(4)$

So, we can write the convergence of a lattice observable to a field as the formal identity:

$\underset{\epsilon \to 0}{\text{lim}}{\epsilon }^{-\mathrm{\Delta }}\hspace{0.17em}O\left(\frac{x}{\epsilon }\right)=\mathrm{\Phi }\left(x\right),(5)$

meant to be valid inside correlators.

If both types of correlators can be separately computed, the potential infinity of identities similar to the previous one will put very strong constraints on the field theory proposed and allow to validate it or, on the contrary, to discard it. The more identities we are able to test, the higher the level of confidence we gain for the conjectural field theory.

At this stage, we seem to be running in a vicious circle: we want to test the proposed field theory by comparing its correlators with the lattice quantities, but we cannot compute the field correlators if we do not know the field theory! If one thinks of a field theory as being given by a set of fields and an action S, this is indeed a serious problem because the action cannot be easily guessed, and even worse, there are many cases for which one has no clue as to what the action is. However, the action is just one convenient (and usually not simple) way to compute correlators. One could think of other ways to determine correlators, and one of them is the presence of symmetry: enough symmetry allows determining the correlators. It is precisely the principle underlying the conformal field theories, which therefore provides a field theoretic framework where no action is necessary. They are discussed in the next section.

Knowing the details of the field theory describing the long-distance properties of a statistical model is at the same time extremely powerful and immensely complicated. On the one hand, it is indeed powerful because it captures the very essential behavior of the statistical model without being cluttered with the many irrelevant lattice effects which make the lattice model so much more complex. On the other hand, it is also immensely complicated because every single element in the lattice model which affects the long distances must have a match in the field theory. Such elements include,

• of course, the bulk observables as discussed above;

• the boundary conditions, the changes of boundary conditions, and the boundary observables;

• the nonlocal observables (like disorder lines in the Ising model);

• the algebra of all the observables;

• the specific effects arising when the lattice is embedded in topologically nontrivial geometries (cylinder, torus, etc.);

• the symmetry, finite or other, that may be present in the model, and possibly many others. All this represents a huge amount of information that must be present and known in the field theory, and which can be only very rarely contemplated in full. A renown exception is when we consider critical statistical models, as we do here, which are in addition formulated on two-dimensional domains (d = 2).

4 Conformal Field Theories

Critical systems are primarily characterized by a scale invariance. The correlation lengths of the observables surviving the scaling limit diverge in the infinite volume limit so that there is no intrinsic length scale left: the fluctuation patterns appear to be the same at all scales. As a consequence, the correlation functions of those observables decay algebraically rather than exponentially. The large-distance 2-point correlator of a typical lattice observable O_i located around the site i takes the following form:

$\langle O_i\hspace{0.17em}{O}_j\rangle =\frac{A}{{|i-j|}^{2\mathrm{\Delta }}}+\dots ,(6)$

where A is a normalization, Δ is the exponent controlling the decay, and the dots indicate lower order terms.

The field theory emerging in the scaling limit inherits the scale invariance. Further assuming translation and rotation symmetries, the scale invariance is enhanced to the invariance under a larger group, namely, the group of conformal transformations, that is, the coordinate transformations which preserve angles.¹¹ In d dimensions, the conformal transformations include the transformations mentioned above, namely, the translations (d real parameters), dilations (1 parameter), and rotations (d(d − 1)/2 parameters), and the so-called special conformal transformations (or conformal inversions) which depend on an arbitrary vector $\overrightarrow{b}$ (d additional parameters) take the following general form:

$\frac{\overrightarrow{x}\prime }{{|\overrightarrow{x}\prime |}^2}=\frac{\overrightarrow{x}}{{|\overrightarrow{x}|}^2}+\overrightarrow{b}\hspace{1em}\iff \hspace{1em}\overrightarrow{x}\prime =\frac{\overrightarrow{x}+|x{|}^2\hspace{0.17em}\overrightarrow{b}}{1+2\overrightarrow{b}\cdot \overrightarrow{x}+|\overrightarrow{b}{|}^2\hspace{0.17em}|\overrightarrow{x}{|}^2}.(7)$

Together these transformations form a finite Lie group isomorphic to SO(d+1,1). They are all global conformal transformations because they are defined everywhere on ${\mathrm{ℝ}}^d{\displaystyle \cup }\left\{\infty \right\}$ and bijective. In dimension d > 2, a conformal transformation defined locally can be extended to a global transformation.

Typical spinless (i.e., rotationally invariant) fields transform tensorially under conformal transformations as follows:

$\mathrm{\Phi }\left(\overrightarrow{x}\right)\hspace{0.17em}\to \hspace{0.17em}{\left|\frac{\partial {x\prime }^i}{\partial x^j}\right|}^{\mathrm{\Delta }/d}\hspace{0.17em}\mathrm{\Phi }\left(\overrightarrow{x}\prime \right),(8)$

for some number Δ. Fields transforming that way under global conformal transformations are called quasi-primary. Global conformal invariance then fixes the average value of a quasi-primary field,

$\langle \mathrm{\Phi }(\overrightarrow{x})\rangle =0,\hspace{1em}\hspace{1em}\text{if}\hspace{0.17em}\mathrm{\Delta }\ne 0,(9)$

(a constant for Δ = 0 by translation invariance) and the 2-point correlator of two quasi-primary fields is given as follows:

$\langle {\mathrm{\Phi }}_1({\overrightarrow{x}}_1)\hspace{0.17em}{\mathrm{\Phi }}_2({\overrightarrow{x}}_2)\rangle =\{\begin{array}{ll}\frac{A_{12}}{{|{\overrightarrow{x}}_1-{\overrightarrow{x}}_2|}^{2{\mathrm{\Delta }}_1}}\hfill & \text{if }\hspace{0.17em}{\mathrm{\Delta }}_1={\mathrm{\Delta }}_2,\hfill \\ 0\hfill & \text{if }\hspace{0.17em}{\mathrm{\Delta }}_1\ne {\mathrm{\Delta }}_2.\hfill \end{array}(10)$

Specializing (8) to a dilation, $\overrightarrow{x}\prime =\alpha \overrightarrow{x}$, we have $\mathrm{\Phi }\left(\overrightarrow{x}\right)\to {\alpha }^{\mathrm{\Delta }}\hspace{0.17em}\mathrm{\Phi }\left(\alpha \overrightarrow{x}\right)$ so that Δ can be identified with the dimension of the field Φ (in units of inverse length).

Global conformal invariance also completely determines the correlator $\langle {\mathrm{\Phi }}_1\left({\overrightarrow{x}}_1\right)\hspace{0.17em}{\mathrm{\Phi }}_2\left({\overrightarrow{x}}_2\right)\hspace{0.17em}{\mathrm{\Phi }}_3\left({\overrightarrow{x}}_3\right)\rangle$ of three (and not more) quasi-primary fields as follows:

$\langle {\mathrm{\Phi }}_1({\overrightarrow{x}}_1)\hspace{0.17em}{\mathrm{\Phi }}_2({\overrightarrow{x}}_2)\hspace{0.17em}{\mathrm{\Phi }}_3({\overrightarrow{x}}_3)\rangle =\frac{A_{123}}{|{\overrightarrow{x}}_1-{\overrightarrow{x}}_2{|}^{{\mathrm{\Delta }}_1+{\mathrm{\Delta }}_2-{\mathrm{\Delta }}_3}\hspace{0.17em}|{\overrightarrow{x}}_1-{\overrightarrow{x}}_3{|}^{{\mathrm{\Delta }}_1+{\mathrm{\Delta }}_3-{\mathrm{\Delta }}_2}\hspace{0.17em}|{\overrightarrow{x}}_2-{\overrightarrow{x}}_3{|}^{{\mathrm{\Delta }}_2+{\mathrm{\Delta }}_3-{\mathrm{\Delta }}_1}}.(11)$

We see that the lattice 2-correlator (6) is consistent with the convergence of the observable O_ito a quasi-primary field$\mathrm{\Phi }\left(\overrightarrow{x}\right)$of dimension Δ upon setting$i=\overrightarrow{x}/\epsilon$since the matching identity (3) is satisfied as follows:

$\underset{\epsilon \to 0}{\text{lim}}\hspace{0.17em}{\epsilon }^{-2\mathrm{\Delta }}\hspace{0.17em}{\langle O_{{\overrightarrow{x}}_1/\epsilon }\hspace{0.17em}{O}_{{\overrightarrow{x}}_2/\epsilon }\rangle }_{\text{lattice}}=\frac{A}{|{\overrightarrow{x}}_1-{\overrightarrow{x}}_2{|}^{2\mathrm{\Delta }}}={\langle \mathrm{\Phi }({\overrightarrow{x}}_1)\hspace{0.17em}\mathrm{\Phi }({\overrightarrow{x}}_2)\rangle }_{\text{FT}}.(12)$

We note that all subdominant terms in the lattice correlator (6) drop out when taking the limit ε → 0, confirming once more that a field theory captures the large-distance behavior of a critical lattice model.

What has been just recalled is valid in any dimension $d\mathrm{⩾}2$ but is only the beginning of the story for d = 2. The global conformal group discussed above remains but is more conveniently presented in complex coordinates as the SL$\left(2,\mathrm{ℂ}\right)$ group of Möbius transformations $w=\frac{az+b}{cz+d}$, for $a,b,c,d\in \mathrm{ℂ}$ satisfying ad − bc = 1.

The two-dimensional world has however many more conformal transformations in store. Indeed, it is a well-known fact that any analytic map w(z) of the complex plane is conformal. Surely, an analytic function requires an infinite number of parameters to fix it (f.i., the coefficients of its Laurent expansion in some neighborhood) so that the conformal “group” is certainly infinite-dimensional. The term group is not really appropriate because the composition of analytic maps is generally not defined everywhere on the complex plane: unless it is a Möbius transformation, an analytic map is either not defined everywhere or its image is not the whole complex plane. For instance, the map $w=\frac{L}{2\pi \text{i}}\text{log}z$ maps the complex plane to a cylinder of circumference L. The discussion of the two-dimensional conformal group is thus usually carried out at the level of its algebra, for which infinitesimal transformations of the form $w=z+\mathit{ϵ}\hspace{0.17em}{z}^{n+1}$ are considered. The corresponding generators satisfy the famous infinite-dimensional Virasoro algebra,

$\left[L_m,L_n\right]=\left(m-n\right)L_{m+n}+\frac{c}{12}m\left(m^2-1\right){\delta }_{m+n,0},\hspace{1em}\hspace{1em}m,n\in \mathrm{\mathbb{Z}},(13)$

a central extension of the Witt algebra. The real number c is the central charge and is one of the most important data of a two-dimensional conformal field theory (CFT). The modes L₀ and L_±1, whose algebra is unaffected by the central charge, are the infinitesimal generators of the Möbius group, with L₋₁ and L₀ corresponding to translations and dilations, respectively. As it turns out, a second commuting copy of the Virasoro algebra, with modes ${\overline{L}}_n$, can formally be considered for the conformal transformations of the antiholomorphic variable $\overline{z}$.

It is not our purpose to give an introduction to the CFT, but one can easily conceive the huge difference between a finite symmetry algebra and an infinite one. A field theory that is to be invariant under an infinite algebra is immensely more constrained and therefore much more rigid, leaving the hope that one should be able to say a lot more about it. It is indeed the case.

For one thing, the field content of a CFT must be organized into representations of the Virasoro algebra, which are all infinite dimensional, and this opens up the possibility that an infinite number of fields be in fact accommodated in a finite number of representations (such CFTs are called rational). In this respect, the primary fields are particularly important. They are the strengthened version of quasi-primary fields in the sense that they transform tensorially under any conformal transformation. A primary field is an eigenfield of L₀ and ${\overline{L}}_0$ with real eigenvalues h and $\overline{h}$ and, more importantly, is annihilated by all positive modes $L_{n>0},{\overline{L}}_{n>0}$. It is in particular characterized by a total weight $\mathrm{\Delta }=h+\overline{h}$ (its eigenvalue under $L_0+{\overline{L}}_0$, the real dilation generator) and is, of course, quasi-primary. The action of any string of negative Virasoro modes $L_{n<0},{\overline{L}}_{n<0}$ on a primary field produces infinitely many new fields, called descendant fields, which include all derivatives of the primary field, since $L_{-1}={\partial }_z$ and ${\overline{L}}_{-1}={\partial }_{\overline{z}}$ act as derivatives on any field. All of them are eigenfields of L₀ and ${\overline{L}}_0$. Together, they form a highest weight representation of the Virasoro algebra whose structure is similar to highest weight representations of simple Lie algebras, the primary field playing the role of the highest weight state.

Like in higher dimensions, the forms of the 1-, 2-, and 3-point of quasi-primary fields are completely fixed by their invariance under Möbius transformations. They are more easily written in complex coordinates (z_ij = z_i − z_j) as follows:

$\langle \mathrm{\Phi }\left(z,\overline{z}\right)\rangle =A\hspace{0.17em}{\delta }_{h,0}{\delta }_{\overline{h},0},(14)$

$\langle {\mathrm{\Phi }}_1\left(z_1,{\overline{z}}_1\right)\hspace{0.17em}{\mathrm{\Phi }}_2\left(z_2,{\overline{z}}_2\right)\rangle =\frac{A_{12}}{z_{12}^{h_1+h_2}\hspace{0.17em}{\overline{z}}_{12}^{\hspace{0.17em}{\overline{h}}_1+{\overline{h}}_2}}\hspace{0.17em}{\delta }_{h_1,h_2}{\delta }_{{\overline{h}}_1,{\overline{h}}_2},(15)$

$\langle {\mathrm{\Phi }}_1\left(z_1,{\overline{z}}_1\right)\hspace{0.17em}{\mathrm{\Phi }}_2\left(z_2,{\overline{z}}_2\right)\hspace{0.17em}{\mathrm{\Phi }}_3\left(z_3,{\overline{z}}_3\right)\rangle =\frac{A_{123}}{z_{12}^{h_1+h_2-h_3}\hspace{0.17em}{z}_{13}^{h_1+h_3-h_2}\hspace{0.17em}{z}_{23}^{h_2+h_3-h_1}\hspace{0.17em}{\overline{z}}_{12}^{\hspace{0.17em}{\overline{h}}_1+{\overline{h}}_2-{\overline{h}}_3}\hspace{0.17em}{\overline{z}}_{13}^{\hspace{0.17em}{\overline{h}}_1+{\overline{h}}_3-{\overline{h}}_2}\hspace{0.17em}{\overline{z}}_{23}^{\hspace{0.17em}{\overline{h}}_2+{\overline{h}}_3-{\overline{h}}_1}}.(16)$

These forms suggest that the conformal weights $h_i,{\overline{h}}_i$ are positive so that the correlators decrease with the separation distances, as seems natural from a physical point of view. We will nonetheless encounter physical fields with negative weights, for which the correlators have a different meaning.

Occasionally, we will consider chiral correlators for which we only retain the dependence in the z_i variables of the full correlators (equivalently the action of the holomorphic modes L_n). Chiral correlators are appropriate for observables living on a boundary, like the real line bordering the upper-half plane, since a boundary is one-dimensional. In this case, only one copy of the Virasoro algebra remains so that the fields are characterized by a single conformal weight. Chiral correlators are also useful to compute the correlators of bulk variables on surfaces with boundaries (see Section 6).

The precise structure of a Virasoro highest weight representation (c, Δ) based on a primary field of weight Δ is crucial. In the good cases, it determines the properties of the primary field (and of its descendants) by fixing its correlators with itself or with other fields. The 2-point correlator of a primary field has the form (15) since it is quasi-primary, and the same is true for the 3-point correlator. To go beyond, the global conformal invariance is not enough.¹² It turns out to be often the case that the structure of a Virasoro highest weight representation implies that the correlators $\langle \mathrm{\Phi }\left(z,\overline{z}\right)\dots \hspace{0.17em}\rangle$ involving the primary field Φ obey differential equations. Four-point functions can be routinely computed in this way. All correlators can then be determined, at least in principle, without knowing anything of a possible Lagrangian realization of the underlying field theory (through its action).

The miracle of 2-dimensional CFTs can be paraphrased in the following way: to completely solve a CFT, that is, compute all its correlation functions, and thereby to know everything there is to know of the large-distance limit of a critical model; it is sufficient to know enough of the Virasoro representations making up that CFT. This methodology has been immensely successful since the mid-80’s and has led to a profound understanding of the many aspects of critical models listed at the end of the previous section. The Ising model is the prominent example of a model that can be treated that way, but the same is true of more general statistical models involving local interactions between the microscopic variables.

More recently, models showing some form of nonlocality have been examined at the conformal light. Sandpile models are in this class, since, as we have seen earlier, the height variables are in strong interaction over the entire domain to form global recurrent configurations. Other models with nonlocal interactions and/or nonlocal degrees of freedom include percolation, critical polymers, and more general loop models. It may sound surprising, but the conclusion seems to be that the conformal approach is still relevant. However, the CFTs underlying these models are more complex essentially because the representations of the Virasoro algebra that appear have a far more complicated structure. These special CFTs are called logarithmic conformal field theories (LCFTs). What follows is a very basic introduction to the salient features of the LCFT; various reviews and applications may be found in the special issue [18]. Let us also mention [19] which reviews the extension to the LCFT of the calculational tools used in CFT.

For the highest weight representations discussed above, the operators $L_0,{\overline{L}}_0$ are diagonalizable. LCFTs have the distinct feature to include Virasoro representations for which L₀ and ${\overline{L}}_0$ are no longer diagonalizable, but instead contain (infinitely many) Jordan blocks of finite rank. To have a rough idea of what these representations look like, one can think of a highest weight representation for which the highest weight is not a single primary field, but a pair of fields (Φ, Ψ), of which only Φ is primary. The action of L₀ on them would be typical of a rank 2 Jordan cell as follows:

$L_0\mathrm{\Phi }=h\mathrm{\Phi },\hspace{1em}\hspace{1em}{L}_0\text{Ψ}=h\text{Ψ}+\lambda \mathrm{\Phi },(17)$

where Ψ is called the logarithmic partner of the primary field Φ and a similar action of ${\overline{L}}_0$ (with $\overline{h}$). Under the action of the negative Virasoro modes, the Jordan block structure will propagate among the descendant fields. The presence of Jordan blocks is a sort of minimal ingredient to make a representation logarithmic; many mathematical complications can and do arise (see, f.i., [20]). Higher rank Jordan blocks can also appear.

A immediate consequence of the presence of Jordan blocks explains the use of the word “logarithmic”: the correlators of fields in an LCFT contain logarithmic terms in addition to the power laws encountered before. For instance, the 2-point correlators of the logarithmic pair {Φ, Ψ}, both of weights $\left(h,\overline{h}\right)$, read

$\langle \Phi \left(z_1,{\overline{z}}_1\right)\hspace{0.17em}\mathrm{\Phi }\left(z_2,{\overline{z}}_2\right)\rangle =0,\hspace{1em}\hspace{1em}\langle \mathrm{\Phi }\left(z_1,{\overline{z}}_1\right)\hspace{0.17em}\text{Ψ}\left(z_2,{\overline{z}}_2\right)\rangle =\frac{B}{{\left(z_1-z_2\right)}^{2h}{\left({\overline{z}}_1-{\overline{z}}_2\right)}^{2\overline{h}}},(\mathrm{18a})$

$\langle \text{Ψ}\left(z_1,{\overline{z}}_1\right)\hspace{0.17em}\text{Ψ}\left(z_2,{\overline{z}}_2\right)\rangle =\frac{C-2\lambda B\text{log}|z_1-z_2{|}^2}{{\left(z_1-z_2\right)}^{2h}{\left({\overline{z}}_1-{\overline{z}}_2\right)}^{2\overline{h}}}.(\mathrm{18b})$

For rank r Jordan blocks, the 2-point correlators would involve up to (r − 1)th powers of logarithms. The parameter λ is not intrinsic as it can be observed in the normalization of Φ or of Ψ; likewise, the logarithmic partner Ψ is defined up to a multiple of Φ without affecting the defining relations (17). The chiral version of the above 2-point functions reads

$\langle \mathrm{\Phi }\left(z_1\right)\hspace{0.17em}\mathrm{\Phi }\left(z_2\right)\rangle =0,\hspace{1em}\hspace{1em}\langle \mathrm{\Phi }\left(z_1\right)\hspace{0.17em}\text{Ψ}\left(z_2\right)\rangle =\frac{B}{{\left(z_1-z_2\right)}^{2h}},(\mathrm{19a})$

$\langle \text{Ψ}\left(z_1\right)\hspace{0.17em}\text{Ψ}\left(z_2\right)\rangle =\frac{C-2\lambda B\log \left(z_1-z_2\right)}{{\left(z_1-z_2\right)}^{2h}},(\mathrm{19b})$

It should not be too surprising that Jordan blocks and logarithms go hand in hand. Under dilation by a factor α, a logarithmic term transforms inhomogeneously log z → log z + log α, reflecting the inhomogeneous action of the dilation generator L₀ on Ψ. Under a finite dilation w = αz, the transformation laws of Φ and Ψ read

$\text{Φ′}\left(w,\overline{w}\right)=|\alpha {|}^{-\mathrm{\Delta }}\hspace{0.17em}\mathrm{\Phi }\left(z,\overline{z}\right),\Psi \prime \left(w,\overline{w}\right)=|\alpha {|}^{-\mathrm{\Delta }}\hspace{0.17em}\left\{\text{Ψ}\left(z,\overline{z}\right)-\lambda \log |\alpha {|}^2\hspace{0.17em}\mathrm{\Phi }\left(z,\overline{z}\right)\right\}.(20)$

One may check that the form of the correlators (18) is indeed invariant under the replacement $\mathrm{\Phi }\left(z,\overline{z}\right)\to \text{Φ′}\left(w,\overline{w}\right)$ and $\text{Ψ}\left(z,\overline{z}\right)\to \text{Ψ′}\left(w,\overline{w}\right)$. Let us also note that the scaling (5) must be redefined for the lattice observables described by logarithmic fields since it involves a dilation by a factor 1/ε, to which the field responds by an inhomogeneous term.

Despite all the efforts spent, LCFTs are generally much less understood than their non-logarithmic cousins, although a number of general features are known. On the statistical side, few models have been thoroughly studied as their nonlocal features make it hard to carry out exact calculations on the lattice. On the field theoretic side, it is not known what a generic LCFT looks like. The simplest of all (but nontrivial) and probably the only LCFT to be fully under control is the symplectic fermion theory with central charge c = −2, also called the triplet theory. It has been introduced in [21] and then investigated in greater detail in [22, 23]. It has the following Lagrangian realization in terms of a pair of free, massless, Grassmannian scalar fields $\theta ,\tilde{\theta }$,

$S=\frac{1}{\pi }{{\displaystyle \int }}^{\text{}}\text{d}z\text{d}\overline{z}\hspace{0.17em}\partial \theta \overline{\partial }\tilde{\theta }\hspace{0.17em},\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\partial ={\partial }_z,\hspace{0.17em}\overline{\partial }={\partial }_{\overline{z}}.(21)$

Several fields in this theory form logarithmic pairs, like the identity $\mathbb{I}$ and the composite field $\theta \tilde{\theta }$. We note that (18) then implies the somewhat unusual relation $\langle \mathbb{I}\rangle =0$, which indeed follows, using the rules of integration over Grassmannian variables, from the fact that the above action does not depend on the constant modes of θ and $\tilde{\theta }$. Since this is a free scalar theory, all correlators of fields that are local (i.e., product of derivatives) in $\theta ,\tilde{\theta }$ are polynomials in the derivatives of the Green function (the kernel of the inverse Laplacian $-4\partial \overline{\partial }$) given in complex coordinates by G(z, w) = −log ǀz−wǀ.

To finish, let us note that the statistical models which have a non-diagonalizable transfer matrix (when there is a proper one) are the natural candidates for being described by LCFTs in their scaling regime. Indeed, such a transfer matrix gives rise to a non-diagonalizable Hamiltonian, which itself is the lattice version of the field theoretic operator $L_0+{\overline{L}}_0$. As said above, the non-diagonalizability of $L_0,{\overline{L}}_0$ is the hallmark of LCFTs. The logarithmic minimal models form an infinite series of such lattice models [24].

The rest of this review is devoted to discussing the variables of the 2-dimensional sandpile models which have been successfully (i.e., with enough confidence) identified in the corresponding LCFT. These elements reveal some facets of the field theory at work in sandpile models: the big and complete picture is well out of reach for the moment.

5 Bulk Variables

The height variables are certainly the first and most natural variables to look at as they are the microscopic variables in terms of which the models are defined. The introduction we gave in Section 2 was for the Abelian sandpile on an arbitrary graph. If large-distance properties should be rather robust against local modifications of a graph, they are not expected to be the same on a graph with a high degree of connectivity (the extreme example being the complete graphs), a regular graph with a moderate degree of connectivity or a graph with a strong hierarchical structure (like Cayley trees). Most of the results reviewed here are obtained when the graph is a rectangular portion of the square lattice ${\mathrm{\mathbb{Z}}}^2$; varying the size of the grid is an easy way to approach the infinite volume limit, and this choice ensures that conservative sites away from the boundary have height variables taking the same number of values (namely, 4). The triangular and honeycomb lattices, for which the number of height values is, respectively, 6 and 3, will be briefly discussed as well in order to address universality issues (see Section 7.2).

In most cases, the only dissipative sites will be located on the boundary,¹³ except when we discuss the insertion of isolated dissipation. With one exception, we will exclusively consider open and closed boundary conditions, by which we mean that whole stretches of boundary sites are either dissipative or conservative, respectively. The choice of boundary conditions not only has an effect at finite volume but also in the infinite volume limit if some of the boundaries are kept at finite distance (e.g., on the upper-half plane or on a strip of finite width, see Section 6).

On a finite grid Γ, the heights assigned to the vertices form stable configurations, but only the recurrent ones have a nonzero (and uniform) weight with respect to the invariant measure ${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}$. So far, we have no clear idea of what a generic recurrent configuration looks like. Answers to questions like “What is the proportion of sites having height 1, height 2, …?” can certainly help figure out. Also, the heights must be correlated within a recurrent configuration. Can one characterize these correlations? Are they exponential or power-lawed? The computation of multisite height probabilities answers these questions and helps understand the statistics of recurrent configurations.

To be definite, let us consider Γ as an L × M rectangular grid in ${\mathrm{\mathbb{Z}}}^2$, with open boundary conditions: the non-boundary sites are conservative and have a maximal height value equal to $z_i^{*}=z_i=4$, whereas the boundary sites have a maximal height value chosen to be $z_i^{*}=4>z_i$ (boundary and corner sites dissipate 1 resp. 2 grains of sand under toppling; both types are connected to the sink). Thus, the toppling matrix is four times the identity minus the adjacency matrix of the grid, and the height at every site takes values in {1, 2, 3, 4}. In this section, all boundaries are sent off to infinity in the scaling limit so that the domain converges to ${\mathrm{ℝ}}^2$; in that limit, all multisite probabilities are fully invariant under translations.

5.1 One-Site Height Probabilities

As a warm up for what has to come, we ask the following: what is the probability ${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h_i=a\right)$ that in a recurrent configuration, a given site i has height equal to a, between 1 and 4? Because we are interested in the infinite volume limit of these numbers, we take i to be deep in the middle of the grid, well away from the boundaries.

If we pause for a while and ponder over that simple question, we feel a bit at a loss on how to handle it because the only means we have is the general criterion of recurrence, namely, the nonexistence of forbidden subconfigurations. Let us start with the height 1.

Since the total number of recurrent configuration is equal to detΔ_Γ (see Section 2), we can write

${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h_i=1\right)=\frac{\#\left\{\text{recurrent}\hspace{0.17em}\text{configs}\hspace{0.17em}\text{with}\hspace{0.17em}{h}_i=1\right\}}{{\text{det}\Delta }_{\mathrm{\Gamma }}}.(22)$

For h_i = 1 to be in a recurrent configuration C, the height of none of its neighbors N, E, S, or W can be equal to 1 (as they would form a forbidden subconfiguration). Following the clever trick proposed in [25], we consider a new grid ${\tilde{\mathrm{\Gamma }}}_i$ by deleting from Γ the vertex i and the four edges incident to it. We also define from C a new configuration $\tilde{C}$ on ${\tilde{\mathrm{\Gamma }}}_i$ by setting

${\tilde{h}}_j=\{\begin{array}{ll}{h}_j\hfill & \text{for}\hspace{0.17em}j\notin \left\{i,\text{N},\text{E},\text{S},\text{W}\right\},\hfill \\ h_j-1\mathrm{⩾}1\hfill & \text{for}\hspace{0.17em}j\in \left\{\text{N},\text{E},\text{S},\text{W}\right\}.\hfill \end{array}(23)$

Looking back at the criterion of recurrence for an arbitrary graph, it is not difficult to see that a configuration C with h_i = 1 is recurrent on Γ if and only if $\tilde{C}$ is recurrent on ${\tilde{\mathrm{\Gamma }}}_i$. We thus obtain

${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h_i=1\right)=\frac{{\text{det}\Delta }_{{\tilde{\mathrm{\Gamma }}}_i}}{{\text{det}\Delta }_{\mathrm{\Gamma }}}=\frac{\text{det}[{\mathrm{\Delta }}_{{\tilde{\mathrm{\Gamma }}}_i}\oplus 1_{ii}]}{{\text{det}\Delta }_{\mathrm{\Gamma }}},(24)$

where the matrix in the numerator has been extended by a one-dimensional diagonal block labeled by the vertex i, without changing the value of the determinant. One then can write

${\mathrm{\Delta }}_{{\tilde{\mathrm{\Gamma }}}_i}\oplus 1_{ii}={\mathrm{\Delta }}_{\mathrm{\Gamma }}+B\left(i\right),\hspace{1em}\hspace{1em}(25)$

with B(i), the defect matrix is given by

$B{\left(i\right)}_{k,k^{\prime }}=\left(\begin{array}{ccccc}-3& 1& 1& 1& 1\\ 1& -1& 0& 0& 0\\ 1& 0& -1& 0& 0\\ 1& 0& 0& -1& 0\\ 1& 0& 0& 0& -1\end{array}\right),\hspace{1em}\hspace{1em}k,k^{\prime }\in \left\{i,\text{N},\text{E},\text{S},\text{W}\right\},(26)$

and B(i) is zero everywhere else. We obtain

${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h_i=1\right)=\frac{\text{det}\left[{\mathrm{\Delta }}_{\mathrm{\Gamma }}+B\left(i\right)\right]}{{\text{det}\Delta }_{\mathrm{\Gamma }}}=\text{det}\left[\mathbb{I}+{\mathrm{\Delta }}_{\mathrm{\Gamma }}^{-1}\hspace{0.17em}B\left(i\right)\right].(27)$

It reduces to the computation of a finite determinant since B(i) has finite rank. In the infinite volume limit (both L, M → ∞), this probability converges to a constant ${\mathrm{\mathbb{P}}}_1$ (by translation invariance). As the matrix Δ_Γ becomes the discrete Laplacian on ${\mathrm{\mathbb{Z}}}^2$ in that limit,¹⁴ standard results yield [25]

${\mathrm{\mathbb{P}}}_1\equiv \underset{\left|\mathrm{\Gamma }\right|\to \infty }{\text{lim}}\hspace{0.17em}{\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h_i=1\right)=\frac{2\left(\pi -2\right)}{{\pi }^3}\simeq 0.073\hspace{0.17em}63.(28)$

It also means that a recurrent configuration has an average of about 7% of sites with a height equal to 1.

What about higher heights? We know for sure that the inequalities ${\mathrm{\mathbb{P}}}_4>{\mathrm{\mathbb{P}}}_3>{\mathrm{\mathbb{P}}}_2>{\mathrm{\mathbb{P}}}_1$ hold because adding one grain of sand to a recurrent configuration, at a site where h_i = a, yields a recurrent configuration if a < 4. However, to actually compute these numbers, can one use the same trick as for the height 1? The answer is definitely negative: no local modification of Γ like what we did above will allow computing the corresponding probabilities. To understand this, we turn to the description in terms of spanning trees.

As was briefly mentioned in Section 2, the burning algorithm yields a one-to-one correspondence between a recurrent configuration and a spanning tree rooted at the sink site s and growing into the interior of Γ. In a given spanning tree $\mathcal{T}$, a site j is called a predecessor of i if the unique path in $\mathcal{T}$ from j to the root passes through i. Let us also denote by X_k(i) the fraction of all spanning trees for which the site i has k predecessors among its nearest neighbors, for $0\mathrm{⩽}k\mathrm{⩽}3$. A careful analysis of the burning algorithm shows the following [26]:

${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h_i=a\right)={\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h_i=a-1\right)+\frac{X_{a-1}\left(i\right)}{5-a},\hspace{1em}\hspace{1em}1\mathrm{⩽}a\mathrm{⩽}4.(29)$

For a = 1, we see that ${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h_i=1\right)$ is related to X₀(i), namely, the fraction of spanning trees on Γ for which the site i is a leaf. All such trees can be obtained from arbitrary spanning trees on ${\tilde{\mathrm{\Gamma }}}_i$ by adding one edge between one neighbor of i and i itself (four different possibilities). Thus, both points of view coincide and lead to the same local modification $\mathrm{\Gamma }\to {\tilde{\mathrm{\Gamma }}}_i$.

The next case is ${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h_i=2\right)$, related to X₁(i). Here, the situation is dramatically different because the condition that i has only one predecessor among its nearest neighbors is highly nonlocal. The reason for this is that there are two manners for a neighbor of i to be a predecessor of i in a given tree. The first one is that the tree includes the edge between the two sites so that the neighbor of i is directly connected to i. In the second manner, the tree contains a potentially long chain of edges that forms a path between the two sites. The first one is a local connection and is easy to check, and the second one is nonlocal and more difficult. The same remark applies to the fractions X₂(i) and X₃(i) and makes the calculation of the corresponding probabilities much more complicated.

In fact, this first natural and simple-looking question we have raised, namely, the value of $\mathrm{\mathbb{P}}\left(h_i=a\right)$, turned into a fairly long warming up exercise as it took about twenty years before the completely explicit probabilities could be found. By using a rather heavy graph theoretical technology, Priezzhev [26] obtained the first expressions for ${\mathrm{\mathbb{P}}}_2,{\mathrm{\mathbb{P}}}_3$, and ${\mathrm{\mathbb{P}}}_4$, but these were given in the form of multivariate integrals. The problem was reconsidered in [27], where the following explicit values were conjectured:

${\mathrm{\mathbb{P}}}_2=\frac{1}{4}-\frac{1}{2\pi }-\frac{3}{{\pi }^2}+\frac{12}{{\pi }^3}\simeq 0.173\hspace{0.17em}90,(\mathrm{30a})$

${\mathrm{\mathbb{P}}}_3=\frac{3}{8}+\frac{1}{\pi }-\frac{12}{{\pi }^3}\simeq 0.306\hspace{0.17em}29,(\mathrm{30b})$

${\mathrm{\mathbb{P}}}_4=\frac{3}{8}-\frac{1}{2\pi }+\frac{1}{{\pi }^2}+\frac{4}{{\pi }^3}\simeq 0.446\hspace{0.17em}17.(\mathrm{30c})$

A few years later, three independent proofs were given. The first one was based on a relation with the probability of a loop-erased random walk (LERW) to visit a fixed nearest neighbor of its starting point, which was then computed in terms of dimer arrangements [28]. The second proof also used the relation with LERW passage probabilities but within a much more general approach [29]. Finally the third one [30] carried out the direct computation of the multiple integrals left open in [26]. Let us mention that the technique developed in [29] to enumerate the so-called cycle-rooted groves (which generalize spanning trees to spanning forests with marked points) currently provides by far the most efficient way to compute height probabilities, reducing the calculation of ${\mathrm{\mathbb{P}}}_2,{\mathrm{\mathbb{P}}}_3$ to just a few lines (see Ref [31]). Most of the height correlators presented below have been computed using this technique. Also noteworthy in this context is the work [32] which presents a direct and elementary derivation of the average height $\langle h\rangle ={\displaystyle \sum }_{a}a\hspace{0.17em}{\mathrm{\mathbb{P}}}_a$ on planar lattices (from the formulas above, it is equal to $\frac{25}{8}$ on ${\mathrm{\mathbb{Z}}}^2$) without computing the individual height probabilities.

Ironically, the four numbers ${\mathrm{\mathbb{P}}}_a$ are not very useful for a comparison with a field theory because they will have to be subtracted in correlators (see below). And indeed, some of the correlators have been determined exactly before the 1-site probabilities ${\mathrm{\mathbb{P}}}_a$ were found.

Even though the explicit expressions for ${\mathrm{\mathbb{P}}}_a$’s have the same level of simplicity, the far larger complexity of the combinatorial problem posed by the calculation of ${\mathrm{\mathbb{P}}}_{a\mathrm{⩾}2}$ hints at a striking difference of nature between the height 1 and the higher heights: the height 1 is essentially local, and the others are nonlocal. This will soon be confirmed.

5.2 Height Cluster Probabilities

Cluster height probabilities are a rather obvious generalization of one-site probabilities, by which we ask for the probability that a specific connected subconfiguration occurs in recurrent configurations, away from the boundaries and in the infinite volume limit. Examples of height clusters are shown below.

The three clusters on the left belong to the family of weakly allowed subconfigurations, or minimal height clusters, first introduced in [25], which contains the cluster made of a single height equal to 1. They are minimal subconfigurations in the sense that if one decreases any of its heights by 1, the clusters become (or contain) forbidden subconfigurations. As was done in the previous subsection for single height 1, their occurrence probabilities around position i can be computed by cutting off appropriate lattice sites and edges. They take the form of finite determinants ${\mathrm{\mathbb{P}}}_S\left(i\right)=\text{det}\left[\mathrm{I}+{\mathrm{\Delta }}_{\mathrm{\Gamma }}^{-1}{B}_S\left(i\right)\right]$, where the defect matrix B_s(i) depends on the cluster S considered [25, 33].

The three clusters on the right of (31) are not minimal and generalize the simple cluster made of a single height larger or equal to 2. Their level of complexity is comparable to the latter and is best computed using the methods of [29]. Explicit calculations become fairly tedious as the size of the cluster increases.

5.3 Height Correlations

In terms of the subtracted height variables,¹⁵

$h_a\left(i\right)\equiv {\delta }_{h\left(i\right),a}-{\mathrm{\mathbb{P}}}_a,(32)$

the n-point correlation functions are given by

${\sigma }_{a_1,a_2,\dots ,a_n}\left(i_1,i_2,\dots ,i_n\right)=\mathbb{E}\left[h_{a_1}\left(i_1\right)\hspace{0.17em}{h}_{a_2}\left(i_2\right)\dots h_{a_n}\left(i_n\right)\right].(33)$

These are the functions we are primarily interested in for a future comparison with a conformal field theory. To make the comparison sensible, we have to take the infinite volume limit and the limit of large separations ǀi_k − i_lǀ → +∞. In addition, to avoid the boundary effects—they will be studied later on, all insertion points i_k are to stay (infinitely) far from the boundaries. In practice, one first replaces Δ_Γ by the Laplacian Δ on ${\mathrm{\mathbb{Z}}}^2$ and then expand the Green matrix (i.e., the inverse Laplacian) for large separations.

The computation of correlations of heights 1 (or indeed any weakly allowed subconfigurations, see below) poses no particular problem. The argument used in Section 5.1 leading to consider new configurations $\tilde{C}$ on a locally modified lattice $\tilde{\mathrm{\Gamma }}$ is simply repeated for the neighborhood of each cluster. Thus, the probability to find a height h(i₁) = 1 at site i₁, a height h(i₂) = 1 at site i₂, and so on, is equal to

${\mathrm{\mathbb{P}}}_{\mathrm{\Gamma }}\left(h\left(i_1\right)=1,\hspace{0.17em}h\left(i_2\right)=1,\dots \right)=\frac{\text{det}\left[{\mathrm{\Delta }}_{\mathrm{\Gamma }}+B\left(i_1\right)+B\left(i_2\right)+\dots \right]}{{\text{det}\Delta }_{\mathrm{\Gamma }}}=\text{det}\left[\mathbb{I}+{\mathrm{\Delta }}_{\mathrm{\Gamma }}^{-1}\left\{B\left(i_1\right)+B\left(i_2\right)+\dots \right\}\right].(34)$

The correlators σ_1,1,…,1 are obtained by taking appropriate subtractions and the limits discussed above.

The first few n-point correlators can be easily computed for arbitrary configurations of insertion points [31, 33]. By construction, the 1-point function vanishes, σ₁(i₁) = 0 (the relation (9) is indeed the main motivation for the subtraction). The 2-point function is found to be ($i_1-i_2=\overrightarrow{r}=r{\text{e}}^{i\phi }$)

${\sigma }_{\mathrm{1,1}}\left(i_1,i_2\right)=-\frac{{\mathrm{\mathbb{P}}}_1^2}{2r^4}-\frac{4\left(\pi -2\right)\left[1+\left(\pi -2\right)\cos 4\phi \right]}{{\pi }^6\hspace{0.17em}{r}^6}+\dots (35)$

where the dots stand for lower order terms.

This first result is instructive for several reasons. First, for large separation distances, the dominant term indicates that the correlation decay is algebraic, which shows that the model is critical and makes room for a conformal field theoretic description. Second, choosing the scale dimension Δ = 2, the scaling limit (12) indeed retains the first term only, the form of which is the expected one (note that the second term, like all other subdominant ones, has only the lattice rotation invariance and is therefore not expected to survive the scaling limit). And third, the dominant term is negative, indicating an anticorrelation between heights 1. This is consistent with the fact that the presence of many heights 1 in a configuration makes it more likely to be nonrecurrent. Interestingly, the calculation can be carried out in ${\mathrm{\mathbb{Z}}}^d$, with the result that the correlation decays like r^−2d, giving a dimension Δ = d [25].

The mixed correlator of a height 1 and a height 2, 3, or 4 is harder. They have first been obtained in [34, 35] by using classical graph theoretic techniques, and then reconsidered and extended in Ref [31] using the results of Ref [29]. Whatever the method used, one has to evaluate the fractions ${\tilde{X}}_k\left(i_1\right)$ of spanning trees, as defined in Section 5.2, but on a lattice modified around the site i₂ where height 1 is located. This modification affects the toppling (Laplacian) matrix and its inverse and, consequently, the whole computation, heavily based on these two matrices. The result for a height 1 and a height 2 reads, at dominant order,

${\sigma }_{\mathrm{2,1}}\left(i_1,i_2\right)=-\frac{{\mathrm{\mathbb{P}}}_1^2}{2r^4}\left\{\text{log}r+\left(\gamma +\frac{3}{2}\text{log}2+\frac{16-5\pi }{2\left(\pi -2\right)}\right)\right\}+\dots (36)$

where γ = 0.577216… is the Euler constant. The first subdominant correction is of order r⁻⁶ and contains a nontrivial angular dependence, like in (35), but also a logr term [31]. The expressions of σ_3,1 and σ_4,1 are similar, with different coefficients.

The expressions σ_a,1 for a > 1 definitely establish the logarithmic character of the CFT underlying the sandpile model. The expressions (35) and (36) are strongly reminiscent of those in (18) but do not quite match. If in the scaling limit, heights 1 and 2 were to converge to a logarithmic pair {h₁(z), h₂(z)}, one would think that σ_1,1 and σ_2,1 ought to go over to the 2-point functions $\langle h_1\left(z_1\right)h_1\left(z_2\right)\rangle$ and $\langle h_2\left(z_1\right)h_1\left(z_2\right)\rangle$, respectively. However, conformal invariance implies that the former vanishes identically, whereas the latter is not logarithmic. We could think of computing σ_2,2 to see what comes out, but large-distance correlators with several heights strictly larger than 1 are far beyond our present computational capabilities. Let us add that the calculation of σ_a,b(i₁, i₂) for a, b > 1 does not merely reduce to the evaluation of numbers like X_a−1,_b−1(i₁, i₂) which would generalize the numbers X_a−1(i) defined earlier and enumerate the spanning trees with fixed numbers of predecessors among the nearest neighbors of i₁ resp. i₂. Indeed, the possibility that neighbors of i₁ are predecessors of i₂, or vice versa, substantially complicates the matter. Details on how to perform the correct counting have been given in [31].

To reconcile the previous lattice results and the LCFT predictions, we pause for a while to examine the effects of a seemingly unrelated observable.

5.4 Isolated Dissipation

In the previous section, the calculation of height probabilities started on a finite grid Γ, where the only dissipative sites are boundary sites. We did not pay too much attention to exactly which boundary sites are dissipative; in fact, since the infinite volume limit sends the boundaries off to infinity, there is no need to know precisely which boundary conditions are used (this is what we meant when we said that Δ_Γ becomes the Laplacian on ${\mathrm{\mathbb{Z}}}^2$). We show now that the situation changes if we make some of the bulk sites dissipative [36]. It is not difficult to understand why this is so in terms of spanning trees. We remember that dissipative sites are sites that are connected to the sink, the root of the trees, from which the branches of the tree are growing. Therefore, the existence of dissipative sites in the bulk makes it possible that branches grow from the middle of the grid, thereby affecting the macroscopic structure of the spanning trees.

To make a bulk site i₁ dissipative, one simply has to connect it to the sink. In the notations of Section 2, this amounts to increase the value $z_{i_1}^{*}$, for instance, from z_i1 = 4 (on ${\mathrm{\mathbb{Z}}}^2$) to 5 (a higher value would not make much difference). In turn, this changes by 1 the diagonal entry (ΔΓ)_i1,i1 of the toppling matrix, that is, Δ_Γ → Δ_Γ + D_i1, with ${\left(D_{i_1}\right)}_{ij}={\delta }_{i,i_1}{\delta }_{j,i_1}$. More generally, the new toppling matrix ${\tilde{\mathrm{\Delta }}}_n\equiv {\mathrm{\Delta }}_{\mathrm{\Gamma }}+D_{i_1}+D_{i_2}+\dots +D_{i_n}$ defines a new model in which several bulk sites i_k are dissipative. As a consequence, the height variables at these sites take values in the set {1, 2, 3, 4, 5}.

A simple and natural way to evaluate the effect of inserting isolated dissipation is to consider the change in the number of recurrent configurations, by computing the ratio ${\text{det}\tilde{\Delta }}_n/{\text{det}\Delta }_{\mathrm{\Gamma }}$, first at finite volume and then in the infinite volume limit.

We start by inserting dissipation at the single site i, far from the boundaries. The ratio is easy to compute since the defect matrix D_i has rank 1,

$\frac{{\text{det}\tilde{\Delta }}_1}{{\text{det}\Delta }_{\mathrm{\Gamma }}}=\text{det}\left(\mathbb{I}+{\mathrm{\Delta }}_{\mathrm{\Gamma }}^{-1}{D}_i\right)=1+{\left({\mathrm{\Delta }}_{\mathrm{\Gamma }}^{-1}\right)}_{i,i}.(37)$

It is a finite number at finite volume but diverges in the infinite volume limit, no matter where the site i is located. The divergence reflects the fact that the extra value h_i = 5 allows enormously more recurrent configurations in the modified model.¹⁶

The same divergence is present in the ratio ${\text{det}\tilde{\Delta }}_n/{\text{det}\Delta }_{\text{Γ}}$, which suggests to change the normalization and compare the effect of inserting n dissipative sites with respect to the situation where there is only one dissipative site, that is, to consider instead the ratio ${\text{det}\tilde{\Delta }}_n/{\text{det}\tilde{\Delta }}_1$, which is perfectly well-defined. Let us also remark that in the infinite volume limit, the denominator does not depend on the location of the (only) dissipative site so that the ratios are fully symmetric in the insertion points i_k and translation invariant.

The first two ratios read, with r = ǀi₁ − i₂ǀ for n = 2,

$\frac{{\text{det}\tilde{\Delta }}_1}{{\text{det}\tilde{\Delta }}_1}=1,\hspace{1em}\hspace{1em}\frac{{\text{det}\tilde{\Delta }}_2}{{\text{det}\tilde{\Delta }}_1}=\frac{1}{\pi }\text{log}r+2{\gamma }_0+\mathcal{O}\left(r^{-2}\right),(38)$

where ${\gamma }_0=\frac{1}{2\pi }\left(\gamma +\frac{3}{2}\text{log}2\right)+1$. If we denote by $\omega \left(z,\overline{z}\right)$ the field that describes, in the scaling limit, the insertion of dissipation at a bulk site, the previous two equations would imply

$\langle \omega \left(z,\overline{z}\right)\rangle =1,\hspace{1em}\hspace{1em}\langle \omega \left(z_1,{\overline{z}}_1\right)\hspace{0.17em}\omega \left(z_2,{\overline{z}}_2\right)\rangle =\frac{1}{\pi }\text{log}\left|z_1-z_2\right|+2{\gamma }_0.(39)$

Interestingly, they exactly match the last two equations of (18), with the logarithmic pair {Φ, Ψ} identified with $\left\{\mathbb{I},\omega \right\}$, both fields having the weights $h=\overline{h}=0$ (the identity field is primary). Moreover, the logarithmic term in the 2-point correlation fixes the coefficient λ of the logarithmic pair $\left(\mathbb{I},\omega \right)$ equal to $\lambda =-\frac{1}{4\pi }$ so that $L_0\hspace{0.17em}\omega ={\overline{L}}_0\hspace{0.17em}\omega =-\frac{1}{4\pi }\mathbb{I}$. The relation $\langle \mathbb{I}\rangle =0$, as noted for the free symplectic fermion theory, is here understood as being given by the inverse of the divergent quantity in (37).

The lattice calculation of ${\text{det}\tilde{\Delta }}_3/{\text{det}\tilde{\Delta }}_1$, corresponding to the insertion of three dissipative sites, is not difficult and yields the following 3-point correlation, with z_ij ≡ z_i − z_j,

$\langle \omega \left(1\right)\hspace{0.17em}\omega \left(2\right)\hspace{0.17em}\omega \left(3\right)\rangle =3{\gamma }_0^2+\frac{{\gamma }_0}{2\pi }\text{log}{\left|z_{12}{z}_{13}{z}_{23}\right|}^2+\frac{1}{16{\pi }^2}\left[\text{log}{\left|z_{12}\right|}^2\hspace{0.17em}\text{log}{\left|\frac{z_{13}{z}_{23}}{z_{12}}\right|}^2+\text{cyclic}\right].(40)$

It is fully consistent with the general 3-point correlators of fields in a logarithmic pair [19]. Many additional checks have been carried out [36] which all confirm the consistency of the above field assignment. It has been shown [27] that the bulk dissipation field can be realized in terms of symplectic free fermions as

$\omega \left(z,\overline{z}\right)=\frac{1}{2\pi }\theta \tilde{\theta }+{\gamma }_0,(41)$

in the sense that the correlators of this composite field, computed in the symplectic fermion theory, reproduce the above expressions.

5.5 Height Correlations Continued

The multisite height probabilities computed in Section 5.3 were obtained by taking the limit over a sequence of grids of increasing size. Because of the dissipation along the boundaries, the probabilities are well-defined for each finite grid and properly converge. On the field theoretic side, the CFT supposedly describing the scaling limit is defined right away on the infinite continuum and does not know about the dissipation of the finite systems. To make the CFT connect with the lattice description, we have to insert by hand the required dissipation in the correlators. Since on the lattice side, the boundary dissipation is pushed off to infinity when we take the infinite volume limit, the previous section suggests that we insert the additional field ω(∞) in the correlators. Thus, the proposal, first made in Ref [27], is that a lattice n-point height correlator is described in the scaling limit by an (n+1)-point field correlator as follows:

${\sigma }_{a_1,a_2,\dots ,a_n}\left(i_1,i_2,\dots ,i_n\right)\stackrel{\text{scalim}}{\to }\langle h_{a_1}\left(z_1\right)\hspace{0.17em}{h}_{a_2}\left(z_2\right)\dots h_{a_n}\left(z_n\right)\hspace{0.17em}\omega \left(\infty \right)\rangle .(42)$

It turns out that the proposed field correlations exactly reproduce the form of the lattice results obtained in Section 5.3. If {Φ, Ψ} are fields of weights $h=\overline{h}$ forming a logarithmic pair such that $\left(L_0-h\right)\text{Ψ}=\left({\overline{L}}_0-h\right)\text{Ψ}=\lambda \mathrm{\Phi }$, one finds with Δ = 2h [27] the following equations:

$\langle \mathrm{\Phi }\left(z_1,{\overline{z}}_1\right)\hspace{0.17em}\mathrm{\Phi }\left(z_2,{\overline{z}}_2\right)\hspace{0.17em}\omega \left(\infty \right)\rangle =\frac{A}{{\left|z_{12}\right|}^{2\mathrm{\Delta }}},\hspace{1em}\langle \mathrm{\Phi }\left(z_1,{\overline{z}}_1\right)\hspace{0.17em}\text{Ψ}\left(z_2,{\overline{z}}_2\right)\hspace{0.17em}\omega \left(\infty \right)\rangle =\frac{B-\lambda A\log {\left|z_{12}\right|}^2}{{\left|z_{12}\right|}^{2\mathrm{\Delta }}},(\mathrm{43a})$

$\langle \text{Ψ}\left(z_1,{\overline{z}}_1\right)\hspace{0.17em}\text{Ψ}\left(z_2,{\overline{z}}_2\right)\hspace{0.17em}\omega \left(\infty \right)\rangle =\frac{C-2\lambda B\log {\left|z_{12}\right|}^2+{\lambda }^{2}A{\log }^2{\left|z_{12}\right|}^2}{{\left|z_{12}\right|}^{2\mathrm{\Delta }}}.(\mathrm{43b})$

Comparing with equation (18), we see that the insertion of dissipation at infinity through ω(∞) allows a nonzero value of A and solves the problem encountered in Section 5.3.

From the dominant terms in equations (35) and (36) for the lattice correlations σ_1,1(i₁, i₂) and σ_2,1(i₁, i₂), we infer that the (subtracted) bulk lattice height 1 and height 2 variables converge to fields h₁(z) and h₂(z) that form a logarithmic pair of weight Δ = 2. Moreover, if we assign them the same normalization as their lattice companions ($A=-\frac{{\mathrm{\mathbb{P}}}_1^2}{2}$), we find the parameter of the logarithmic pair (h₁, h₂) equal to $\lambda =-\frac{1}{2}$. The explicit results for σ_3,1(i₁,i₂) and σ_4,1(i₁,i₂) [35] show that the fields h₃(z) and h₄(z) are also logarithmic partners of h₁(z), albeit of different normalizations and for different values of λ. As noted in Section 4, it means that they can be written as linear combinations¹⁷ of h₁(z) and h₂(z) with known coefficients,

$h_a\left(z\right)={\alpha }_a{h}_2\left(z\right)+{\beta }_a{h}_1\left(z\right),\hspace{1em}\hspace{1em}{\alpha }_1=0,\hspace{0.17em}{\alpha }_2=1,\hspace{0.17em}{\alpha }_3=\frac{8-\pi }{2\left(\pi -2\right)},\hspace{0.17em}{\alpha }_4=-\frac{\pi +4}{2\left(\pi -2\right)},(44)$

and other values for β_a [27]. These field assignments predict that the lattice correlation of heights larger or equal to 2 behaves asymptotically as follows:

${\sigma }_{a,b}\left(i_1,i_2\right)\simeq -{\alpha }_a{\alpha }_b\frac{{\mathrm{\mathbb{P}}}_1^2}{2}\hspace{0.17em}\frac{{\text{log}}^{2}r}{r^4}+\mathcal{O}\left(\frac{\log r}{r^4}\right),\hspace{1em}\hspace{1em}a,b\mathrm{⩾}2.(45)$

Because α₄ is the only negative coefficient among α_a, the height variables are all anticorrelated, except the height 4 which has a positive correlation with the other three heights. Numerical simulations have successfully confirmed the behavior (45) [27]. A lattice proof however remains one of the greatest challenges in the sandpile models.

The lattice 2-point correlators discussed above correspond to 3-point functions in the CFT. They are therefore completely generic, depending only on the weights of the fields involved and a few assumptions about their global conformal transformations. Higher correlators are not generic and depend on finer details of the nature of the fields and of the specific CFT at work, in particular its central charge. In this regard, the first hint for the value of the central charge was given in Ref [8] by looking at the finite-size corrections of the partition function (i.e., the number of recurrent configurations); the analysis yields the value c = −2.

The simplest higher correlators to consider on the lattice are the 3- and 4-point height 1 correlators. They have been computed in Ref [33] with the following results. Since the height 1 variable has weight Δ = 2 in the scaling limit, one would expect the dominant contribution to the 3-point correlator to be homogeneous of degree −6 in the separation distances. Surprisingly, the first nonzero term has degree −8,

${\sigma }_{\mathrm{1,1,1}}\left(i_1,i_2,i_3\right)=0+\dots (46)$

implying that its scaling limit, corresponding to the CFT 4-point function $\langle h_1\left(z_1\right)\hspace{0.17em}{h}_1\left(z_2\right)\hspace{0.17em}{h}_1\left(z_3\right)\hspace{0.17em}\omega \left(\infty \right)\rangle$, vanishes identically.

The lattice 4-point correlator has the expected dominant degree −8,

${\sigma }_{\mathrm{1,1,1,1}}\left(i_1,i_2,i_3,i_4\right)=\frac{{\mathrm{\mathbb{P}}}_1^4}{8}\hspace{0.17em}\left\{\frac{1}{{\left|z_{12}\hspace{0.17em}{z}_{34}\right|}^4}+\frac{1}{{\left|z_{13}\hspace{0.17em}{z}_{24}\right|}^4}+\frac{1}{{\left|z_{14}\hspace{0.17em}{z}_{23}\right|}^4}-\frac{1}{{\left(z_{12}\hspace{0.17em}{z}_{34}\hspace{0.17em}{\overline{z}}_{13}\hspace{0.17em}{\overline{z}}_{24}\right)}^2}-\frac{1}{{\left(z_{13}\hspace{0.17em}{z}_{24}\hspace{0.17em}{\overline{z}}_{14}\hspace{0.17em}{\overline{z}}_{23}\right)}^2}-\frac{1}{{\left(z_{14}\hspace{0.17em}{z}_{23}\hspace{0.17em}{\overline{z}}_{12}\hspace{0.17em}{\overline{z}}_{34}\right)}^2}+\text{c}.\text{c}.\right\}+\dots (47)$

and is much more instructive: it is precisely the expression we obtain for the 5-point CFT correlation function $\langle h_1\left(z_1\right)\hspace{0.17em}{h}_1\left(z_2\right)\hspace{0.17em}{h}_1\left(z_3\right)\hspace{0.17em}{h}_1\left(z_4\right)\hspace{0.17em}\omega \left(\infty \right)\rangle$ if we assume that the height 1 field h₁(z) is a primary field of dimensions $\left(h,\overline{h}\right)=\left(\mathrm{1,1}\right)$ in a CFT with central charge c = −2 and that it satisfies a certain degeneracy condition at level 2. This last condition furnishes a differential equation [16], from which the correlator can be fully determined, the result being exactly the function (47)! From this result, one can actually infer that the previous correlator $\langle h_1\left(z_1\right)\hspace{0.17em}{h}_1\left(z_2\right)\hspace{0.17em}{h}_1\left(z_3\right)\hspace{0.17em}\omega \left(\infty \right)\rangle$ must vanish if it is to be symmetrical in the three insertion points [31].

The lattice 3-point correlators σ_a,1,1(i₁,i₂,i₃), $a\mathrm{⩾}2$, have been computed more recently in [31]. For simplicity, the three points were assumed to be aligned horizontally in the plane, with real separations x_ij. The following result was obtained to dominant order:

${\sigma }_{a,\mathrm{1,1}}(i_1,i_2,i_3)={\alpha }_a\hspace{0.17em}\frac{{\mathrm{\mathbb{P}}}_1^3}{8}\hspace{0.17em}\frac{1}{x_{21}^3{x}_{31}^3}+\dots ,\hspace{1em}\hspace{1em}\left(a\mathrm{⩾}2\right),(48)$

where the coefficients α_a are those given in (44). In addition to being very simple, this expression is surprisingly non-logarithmic. Because its scaling limit should be given by $\langle h_a\left(z_1\right)\hspace{0.17em}{h}_1\left(z_2\right)\hspace{0.17em}{h}_1\left(z_3\right)\hspace{0.17em}\omega \left(\infty \right)\rangle$, it is particularly important to understand it from the CFT point of view. Indeed, the computation of such a correlator requires additional information on the height fields $h_{a\mathrm{⩾}2}$ as logarithmic partners of h₁. Since h₃ and h₄ can be regarded as linear combinations of h₁ and h₂, it is sufficient to consider h₂.

Inspired by the conformal representations appearing in the bosonic sector of the symplectic theory [22], the following proposal has been made in Ref [27] regarding the conformal nature of $h_2\left(z,\overline{z}\right)$. A more complete account will be presented in Section 8.

The field $h_2\left(z,\overline{z}\right)$ is not primary since it transforms into $h_1\left(z,\overline{z}\right)$ under dilations. It is also not quasi-primary because its L₁ and ${\overline{L}}_1$ transforms generate two new fields $\rho \left(z,\overline{z}\right)$ and $\overline{\rho }\left(z,\overline{z}\right)$, respectively, with weights (0,1) and (1,0). Moreover, the field $\rho \left(z,\overline{z}\right)$ is left primary, and its ${\overline{L}}_1$ transform is equal to $\kappa \mathbb{I}$; likewise, $\overline{\rho }\left(z,\overline{z}\right)$ is right primary, and its L₁ transform is also equal to $\kappa \mathbb{I}$. All this results in the following transformation law of $h_2\left(z,\overline{z}\right)$ under a general conformal transformation $z\to w\left(z\right)$ and $\overline{z}\to \overline{w}\left(\overline{z}\right)$ as follows:

$\begin{array}{l}{h}_2\left(z,\overline{z}\right)={\left|\frac{\text{d}w}{\text{d}z}\right|}^2\left[h_2\left(w,\overline{w}\right)+\text{log}{\left|\frac{\text{d}w}{\text{d}z}\right|}^2{h}_1\left(w,\overline{w}\right)\right]+\frac{1}{2}\left(\frac{{\text{d}}^{2}w}{\text{d}{z}^2}/\frac{\text{d}w}{\text{d}z}\right)\frac{\text{d}\overline{w}}{\text{d}\overline{z}}\rho \left(w,\overline{w}\right)\\ +\frac{1}{2}\frac{\text{d}w}{\text{d}z}\left(\frac{{\text{d}}^2\overline{w}}{\text{d}{\overline{z}}^2}/\frac{\text{d}\overline{w}}{\text{d}\overline{z}}\right)\overline{\rho }\left(w,\overline{w}\right)+\frac{\kappa }{4}\left(\frac{{\text{d}}^{2}w}{\text{d}{z}^2}/\frac{\text{d}w}{\text{d}z}\right)\left(\frac{{\text{d}}^2\overline{w}}{\text{d}{\overline{z}}^2}/\frac{\text{d}\overline{w}}{\text{d}\overline{z}}\right),\kappa =-\frac{{\mathrm{\mathbb{P}}}_1}{4}.\end{array}(49)$

This conformal transformation law of h₂ is sufficient to compute correlators involving h₂, but substantially complicates the calculations. Using this transformation and the left and right level 2 degeneracies of $h_1\left(z,\overline{z}\right)$, the required correlator can be nonetheless determined. The result reads [31]

$\langle h_2\left(z_1\right)\hspace{0.17em}{h}_1\left(z_2\right)\hspace{0.17em}{h}_1\left(z_3\right)\hspace{0.17em}\omega \left(\infty \right)\rangle =\frac{{\mathrm{\mathbb{P}}}_1^3}{16}\hspace{0.17em}\frac{1}{{\left|z_{12}{z}_{13}\right|}^2}\hspace{0.17em}\left[\frac{1}{z_{13}{\overline{z}}_{12}}+\frac{1}{z_{12}{\overline{z}}_{13}}\right].(50)$

When the three points z_i are aligned horizontally, it exactly reproduces the lattice result (48) for a = 2. The cases a = 3, 4 follow by multiplying by the proper coefficient α_a since $\langle h_1\left(z_1\right)\hspace{0.17em}{h}_1\left(z_2\right)\hspace{0.17em}{h}_1\left(z_3\right)\hspace{0.17em}\omega \left(\infty \right)\rangle =0$.

5.6 Minimal Height Cluster Correlations

The calculation of occurrence probabilities of minimal subconfigurations has been briefly discussed in Section 5.2. Their correlations can be computed very much like those of heights 1 by using a defect matrix. The calculation of mixed 2-point correlators for about a dozen different minimal subconfigurations has been reported in Ref [33]. It turns out that each such cluster S can be specified by a triplet (a, b₁, b₂) of real numbers.

We define as before subtracted variables

$h_S\left(i\right)={\delta }_{S\left(i\right)}-{\mathrm{\mathbb{P}}}_S,(51)$

where δ_S(i) denotes the event “a minimal subconfiguration S is found around site i” and ${\mathrm{\mathbb{P}}}_S\left(i\right)=\mathbb{E}[{\delta }_{S\left(i\right)}]$ is the probability of such an event. The mixed correlator of two such variables takes the form

${\sigma }_{S,S^{\prime }}\left(i_1,i_2\right)=\mathbb{E}\left[h_S\left(i_1\right)\hspace{0.17em}{h}_{S^{\prime }}\left(i_2\right)\right]=-\frac{1}{2r^4}\left\{a{a}^{\prime }+\left(b_1{b}_1^{\prime }-b_2{b}_2^{\prime }\right)\cos 4\phi \right\}+\dots (52)$

We see that the dominant contribution retains an angular dependence, which in this case is not surprising since the minimal clusters are generally not rotationally invariant. As a matter of illustration, the cluster reduced to a single height 1 is characterized by the triplet $\left(a,b_1,b_2\right)=\left({\mathrm{\mathbb{P}}}_1,\mathrm{0,0}\right)$, while for the second one in equation (31), one has

In the scaling limit, the subtracted cluster variables give rise to the fields h_S(z), whose mixed correlators are given by the terms displayed in equation (52). Interestingly, it has been observed [33] that these fields have a realization in terms of the symplectic free fermions, which is discussed at the end of Section 4. Indeed, one may check that the explicit fields given by

$h_S\left(z,\overline{z}\right)=-\left\{a\left(\partial \theta \overline{\partial }\tilde{\theta }+\overline{\partial }\theta \partial \tilde{\theta }\right)+\left(b_1+\text{i}{b}_2\right)\partial \theta \partial \tilde{\theta }+\left(b_1-\text{i}{b}_2\right)\overline{\partial }\theta \overline{\partial }\tilde{\theta }\right\},(54)$

reproduce the above 2-point correlators, as well as the higher order correlators computed in Ref [33], provided the dissipation field ω, proportional to $\theta \tilde{\theta }$, is inserted in the correlators, as explained earlier. In particular, the height 1 field $h_1\left(z,\overline{z}\right)$ is recovered upon setting $a={\mathrm{\mathbb{P}}}_1$ and b₁ = b₂ = 0. Let us note a generic field $h_S\left(z,\overline{z}\right)$ is a linear combination of three fields with different conformal weights, namely, (1,1), (2,0), and (0,2) and therefore different conformal transformations; the last two are responsible for inducing an angular dependence in the correlators. The field realization (54) has been proved in a much greater generality in Ref [37]: any lattice observable based on a conservative local bond modification¹⁸ converges in the scaling limit to a field of the form (54).

On general grounds, this should not be surprising. On the one hand, the multisite probabilities for minimal clusters can be computed by using defect matrices which implement the local bond modifications. On the other hand, defect matrices always yield contributions that are given by finite determinants of discrete Green matrix entries. In the limit of large separations, the determinants converge to polynomial expressions in the Green function and its derivatives. It is therefore not a complete surprise that the associated fields can be constructed out from the symplectic free fermions $\theta ,\tilde{\theta }$. Indeed, because the 2-point correlators of $\theta ,\tilde{\theta }$ are given by the Green function, the correlators of any fields that are local in $\theta ,\tilde{\theta }$ and their derivatives are necessarily polynomials in the Green function and its derivatives (Wick’s theorem). We expect this observation to extend to all the observables that correspond to local perturbations of the toppling matrix. Isolated dissipation and minimal cluster variables are among them; the arrow variables discussed in the next section are in this class too. These general remarks apply to the massive extension of the sandpile model (see Section 7.1).

What about the height 2, 3, and 4 variables? Can they also be accommodated in the free symplectic fermion theory? As explained earlier, these three variables cannot be handled with finite rank perturbations of the toppling matrix because they involve nonlocal constraints on the nearest neighbors (some of them should not be predecessors). Using the technique developed in Ref [29], the 1-site probabilities ${\mathrm{\mathbb{P}}}_{a\mathrm{⩾}2}$ can be efficiently computed. Surprisingly, the details show that the explicit values are given in terms of a few entries of the lattice Green matrix (at short distances), which explains why the values of ${\mathrm{\mathbb{P}}}_{a\mathrm{⩾}2}$ quoted in (30) are not much more complicated than for ${\mathrm{\mathbb{P}}}_1$. This is no longer the case for large-distance correlations ${\sigma }_{a\mathrm{⩾}\mathrm{2,1}}\left(i_1,i_2\right)$. The analysis of Ref [31] shows that those correlators are expressed in terms of sums of the product of Green matrix entries over a path connecting i₁ to i₂ and thus in terms of quantities that are not local in the Green matrix. It supports the view that the height fields $h_{a\mathrm{⩾}2}$ do not belong to the free symplectic fermion theory. A detailed analysis of this question has been carried out in Ref [27] and has reached the same conclusion. Section 8 below summarizes this somewhat strange situation.

5.7 Spanning Tree–Related Variables

A recurrent configuration of the sandpile model can be specified as a set of height values or as a spanning tree; the former has the local heights as natural variables, the latter has local connectivities as natural variables, namely, the existence or absence of specific bonds in the spanning tree.

We recall that a spanning tree is a connected subgraph with no loop which contains all vertices, including the sink. The latter is chosen to be the root of the tree, implying that there is a unique path connecting any vertex to the root and therefore any vertex to any other vertex. A rooted spanning tree can then naturally be oriented by deciding that the edges of the tree all point toward the root. As a consequence, in any rooted spanning tree, there is exactly one outgoing edge at each vertex but the root; there may however be more (or less) ingoing edges (a vertex with no ingoing edge is a leaf). A site j is then a predecessor of i if the unique path from j to i is consistently oriented (equivalently if the unique path form j to i does not pass through the root). As we have seen, the question of being predecessor is a nonlocal problem, even if i and j are close to each other, even nearest neighbors [38].

Connectivities between neighboring sites can be handled in much the same way as height 1 or minimal height cluster variables. To see this, we must first understand why the determinant of the toppling matrix on a graph counts the number of spanning trees on that graph. In the perspective of this section, we generalize the matrix by assigning arbitrary weights to the oriented edges of the graph ${\text{Γ}}^{*}=\left(V^{*},E^{*}\right)$. We define x_ij as the weight carried by the edge from the vertex i to the vertex j (x_ij = 0 means that there is no edge from i to j), and we set, for $i,j\in V$,

${\mathrm{\Delta }}_{i,j}=\{\begin{array}{ll}{y}_i=x_{i*}+{\displaystyle \sum }_{j\ne i}\hspace{0.17em}{x}_{ij}\hfill & \text{for}\hspace{0.17em}i=j,\hfill \\ -x_{ij}\hfill & \text{for}\hspace{0.17em}i\ne j,\hfill \end{array}(55)$

In the context of the sandpile model, the difference $x_{i*}=y_i-{\displaystyle \sum }_{j\ne i}\hspace{0.17em}{x}_{ij}$ can be viewed as the weight of the oriented edge from i to the sink $*$ so that the conservative vertices have this difference equal to 0 (no connection to the root).

If N = ǀVǀ is the number of vertices in the graph, let us write the determinant of Δ as a sum over the permutations σ of the symmetric group S_N, which we partition according to the number k of proper cycles they contain, that is, the cycles of length strictly larger than 1 as follows:

$\text{det}\Delta ={\displaystyle \sum }_{\sigma \in S_N}\hspace{0.17em}{\epsilon }_{\sigma }\hspace{0.17em}{\mathrm{\Delta }}_{1,\sigma \left(1\right)}{\mathrm{\Delta }}_{2,\sigma \left(2\right)}\dots {\mathrm{\Delta }}_{N,\sigma \left(N\right)}={\displaystyle \sum }_{k=0}^{\left[N/2\right]}{\left(-1\right)}^k{\displaystyle \sum }_{\begin{array}{c}\sigma \hspace{0.17em}\text{has}\hspace{0.17em}k\\ \text{proper}\hspace{0.17em}\text{cycles}\end{array}}\hspace{0.17em}{\mathrm{\Delta }}_{1,\sigma \left(1\right)}^{+}\dots {\mathrm{\Delta }}_{N,\sigma \left(N\right)}^{+}(56)$

where the matrix Δ⁺ is Δ without the minus signs in the non-diagonal part. The second equality follows by combining the signs in the non-diagonal entries of Δ with the parity of σ: every cycle of length $\ell \mathrm{⩾}2$ in a permutation σ brings a sign ${\left(-1\right)}^{\ell -1}$ coming from the parity ε_σ and another sign ${\left(-1\right)}^{\ell }$ from the product of non-diagonal entries of Δ, resulting in an overall sign −1 per proper cycle. A cycle of length $\ell =1$, corresponding to a vertex left invariant by σ, brings no sign.

The term k = 0 is simply equal to ${\displaystyle \prod }_i{y}_i$ as the only permutation with no proper cycle is the identity. In combinatorial terms, the product ${\displaystyle \prod }_i{y}_i$ is the weighted sum over all configurations of N arrows, where each vertex has exactly one arrow pointing to one of the other vertices or to the root, each configuration being weighted by the product of the weights carried by the arrows. The generic term k ≠ 0, apart for the sign (−1)^k, is a weighted sum of arrow configurations which contain at least k oriented loops. Indeed, the k cycles contained in a fixed σ give rise to k loops, and the arrows attached to the vertices left invariant by σ are unconstrained and possibly form more loops.

By using the inclusion–exclusion principle, one can see that the above alternating sum has the effect to subtract from the term k = 0 the weights of all the arrow configurations which contain at least one loop [26, 39]. Thus, detΔ is the sum over the oriented spanning trees on ${\text{Γ}}^{*}$, each tree being weighted by the product of the weights of the oriented edges present in the tree (Kirchhoff’s theorem). These oriented trees are also rooted spanning trees because one vertex at least must have its arrow oriented to the sink (a configuration of N arrows on N vertices necessarily contains a loop). When all weights x_ij are equal to 0 or 1, detΔ is simply the number of spanning trees.

Let us come back to the question of local connectivities on a rectangular grid in ${\mathrm{\mathbb{Z}}}^2$ and compute the probability that the outgoing arrow from the site i is oriented to its right neighbor. This amounts to compute the fraction of spanning trees with such an arrow given as follows:

${\mathrm{\mathbb{P}}}_{\to }\left(i\right)=\mathrm{\mathbb{P}}\left(\text{right}\hspace{0.17em}\text{arrow}\hspace{0.17em}\text{at}\hspace{0.17em}i\right)=\frac{\left\{\#\hspace{0.17em}\text{trees}\hspace{0.17em}\text{with}\hspace{0.17em}\text{right}\hspace{0.17em}\text{arrow}\hspace{0.17em}\text{at}\hspace{0.17em}i\right\}}{\text{det}\Delta },(57)$

where Δ is the discrete Laplacian. We take i to be a conservative, non-boundary site.

According to the general discussion above, in order to force an arrow from i to its right neighbor E, one could simply set to 0 the weights between i and its three neighbors S, W, and N. It is however computationally more efficient to set the weight from i to its right neighbor to x and take the limit x → + ∞ so as to give the edges to the other three neighbors a relative weight equal to 0. This implies that we define a new matrix $\tilde{\mathrm{\Delta }}$ which coincides with Δ, except on two entries, namely, ${\tilde{\mathrm{\Delta }}}_{\mathit{i,i}+{\widehat{e}}_1}=-x$ and ${\tilde{\mathrm{\Delta }}}_{ii}=x+3$. As before, we write $\tilde{\mathrm{\Delta }}=\mathrm{\Delta }+B_{\to }\left(i\right)$ for the defect matrix B_→(i) which is zero everywhere, except on the two sites $i,i+{\widehat{e}}_1$, where it reduces to $\left(\begin{array}{c}x-1\\ 0\end{array}\hspace{0.17em}\begin{array}{c}1-x\\ 0\end{array}\right)$. Since x is in any case large, we can simply set that part of B_→(i) to $\left(\begin{array}{c}x\\ 0\end{array}\hspace{0.17em}\begin{array}{c}-x\\ 0\end{array}\right)$. From Kirchhoff’s theorem, we obtain

${\mathrm{\mathbb{P}}}_{\to }\left(i\right)=\underset{x\to +\infty }{\text{lim}}\hspace{0.17em}\frac{1}{x}\hspace{0.17em}\text{det}\left[\mathbb{I}+{\mathrm{\Delta }}^{-1}{B}_{\to }\left(i\right)\right],(58)$

which reduces to a 2-by-2 determinant. In the infinite volume limit, one finds ${\mathrm{\mathbb{P}}}_{\to }\left(i\right)=\frac{1}{4}$, as expected. If we want to have the arrow at i oriented to its neighbor j, other than E, we use similar defect matrices B_↑, B_↓, B_←, which look the same as B_→ but with the nonzero 2-by-2 block $\left(\begin{array}{c}x\\ 0\end{array}\hspace{0.17em}\begin{array}{c}-x\\ 0\end{array}\right)$ placed on the sites i and j. The four orientations of the arrow at i yield the same result $\frac{1}{4}$.

Multipoint arrow probabilities can be computed in the now usual way, placing appropriate defect matrices at the different sites. For instance, the probability to find a right arrow at two sites i₁ and i₂ is given as follows:

${\sigma }_{\to ,\to }\left(i_1,i_2\right)=\underset{x\to +\infty }{\text{lim}}\hspace{0.17em}\frac{1}{x^2}\hspace{0.17em}\text{det}\left[\mathbb{I}+{\mathrm{\Delta }}^{-1}\left\{B_{\to }\left(i_1\right)+B_{\to }\left(i_2\right)\right\}\right].(59)$

In the infinite volume limit and for a large distance, the following two-point probabilities are found at dominant order:

${\sigma }_{\to ,\to }\left(i_1,i_2\right)-\frac{1}{16}=\frac{1}{16{\pi }^2}\hspace{0.17em}\frac{{\left(z+\overline{z}\right)}^2}{{\left|z\right|}^4}+\dots (\mathrm{60a})$

${\sigma }_{\uparrow ,\uparrow }\left(i_1,i_2\right)-\frac{1}{16}=-\frac{1}{16{\pi }^2}\hspace{0.17em}\frac{{\left(z-\overline{z}\right)}^2}{{\left|z\right|}^4}+\dots (\mathrm{60b})$

${\sigma }_{\to ,\uparrow }\left(i_1,i_2\right)-\frac{1}{16}=-\frac{\text{i}}{16{\pi }^2}\hspace{0.17em}\frac{z^2-{\overline{z}}^2}{{\left|z\right|}^4}+\dots (\mathrm{60c})$

Looking for symplectic fermion realizations of fields ρ_→ and ρ_↑ which reproduce these two-point correlators, one quickly sees that these have to include two parts with respective weights (1,0) and (0,1), leading in a natural way to the following forms:

${\rho }_{\to }\left(z,\overline{z}\right)=\frac{1}{2\pi }\left(\theta \hspace{0.17em}\partial \tilde{\theta }+\theta \hspace{0.17em}\overline{\partial }\tilde{\theta }\right),\hspace{1em}\hspace{1em}{\rho }_{\uparrow }\left(z,\overline{z}\right)=\frac{\text{i}}{2\pi }\left(\theta \hspace{0.17em}\partial \tilde{\theta }-\theta \hspace{0.17em}\overline{\partial }\tilde{\theta }\right),(61)$

in agreement with the fact that ${\rho }_{\uparrow }\left(z,\overline{z}\right)$ is formally the rotated form of ${\rho }_{\to }\left(z,\overline{z}\right)$ (under z → −iz). The first two correlators are indeed related by rotation. This observation also suggests that the other two orientations are described by fields which are the opposite of the previous two, ${\rho }_{\leftarrow }\left(z,\overline{z}\right)=-{\rho }_{\to }\left(z,\overline{z}\right)$ and ${\rho }_{\downarrow }\left(z,\overline{z}\right)=-{\rho }_{\uparrow }\left(z,\overline{z}\right)$. Explicit calculations confirm it.

In a similar way, probabilities that edges belong to a random spanning tree, irrespective of their orientation, can be computed. The probability that a single, fixed edge belonging to a tree is the sum of the probabilities to find it in either of the two possible orientations, and is thus equal to $\frac{1}{2}$.

Likewise, the probability to find m edges in a tree is the sum of the probabilities to find them in all possible orientations and so is the sum of 2^m probabilities of m oriented edges. That sum can however be obtained in one go by replacing the block $\left(\begin{array}{c}x\\ 0\end{array}\hspace{0.17em}\begin{array}{c}-x\\ 0\end{array}\right)$ used for the oriented edges by $\left(\begin{array}{c}x\\ -x\end{array}\hspace{0.17em}\begin{array}{c}-x\\ x\end{array}\right)$ and dividing as above the determinant by x^m. Because two arrows with opposite orientations cannot occupy the same edge (as they would form a loop), the summation over the 2^m terms is correctly realized.

Correlations of unoriented edges should decay faster than those of oriented edges because the sum of the two orientations is zero, in view of the relations ρ_← = −ρ_→ and ρ_↓= −ρ_↑, at least at the order that was dominant for the oriented edges (r⁻²). Indeed, explicit calculations yield a r⁻⁴ decay:

${\sigma }_{\leftrightarrow ,\leftrightarrow }\left(i_1,i_2\right)-\frac{1}{4}={\sigma }_{\updownarrow ,\updownarrow }\left(i_1,i_2\right)-\frac{1}{4}=-\frac{1}{16{\pi }^2}\hspace{0.17em}\frac{{\left(z^2+{\overline{z}}^2\right)}^2}{|z{|}^8}+\dots (\mathrm{62a})$

${\sigma }_{\leftrightarrow ,\updownarrow }\left(i_1,i_2\right)-\frac{1}{4}=\frac{1}{16{\pi }^2}\hspace{0.17em}\frac{{\left(z^2-{\overline{z}}^2\right)}^2}{|z{|}^8}+\dots (\mathrm{62b})$

From these correlators, the associated fields ϕ_↔ and ϕ_↕ must have components with conformal weights (2,0), (1,1), and (0,2). One finds the same form as the fields associated to the minimal clusters, in agreement with a previous remark since the defect matrix $\left(\begin{array}{c}x\\ -x\end{array}\hspace{0.17em}\begin{array}{c}-x\\ x\end{array}\right)$ is conservative. More precisely, the following fermionic expressions reproduce the above correlations:

${\varphi }_{\leftrightarrow }\left(z,\overline{z}\right)=\frac{1}{2\pi }\left(\partial \theta \overline{\partial }\tilde{\theta }+\overline{\partial }\theta \partial \tilde{\theta }+\partial \theta \partial \tilde{\theta }+\overline{\partial }\theta \overline{\partial }\tilde{\theta }\right)=\left(\partial +\overline{\partial }\right){\rho }_{\to }-\frac{1}{2\pi }\left(\theta {\partial }^2\tilde{\theta }+\theta {\overline{\partial }}^2\tilde{\theta }\right),(\mathrm{63a})$

${\varphi }_{\updownarrow }\left(z,\overline{z}\right)=\frac{1}{2\pi }\left(\partial \theta \overline{\partial }\tilde{\theta }+\overline{\partial }\theta \partial \tilde{\theta }-\partial \theta \partial \tilde{\theta }-\overline{\partial }\theta \overline{\partial }\tilde{\theta }\right)=\text{i}\left(\partial -\overline{\partial }\right){\rho }_{\uparrow }-\frac{1}{2\pi }\left(\theta {\partial }^2\tilde{\theta }+\theta {\overline{\partial }}^2\tilde{\theta }\right).(\mathrm{63b})$

In fact, given that an unoriented edge is a sum of two oriented edges with opposite orientation, or, from what we said above, a difference of two oriented edges with the same orientation, one would expect that the fields describing a horizontal resp. vertical unoriented edge are proportional to the horizontal resp. vertical derivative of the fields describing the oriented edges, namely, ${\varphi }_{\leftrightarrow }\sim \left(\partial +\overline{\partial }\right){\rho }_{\to }$ and ${\varphi }_{\updownarrow }\sim \text{i}\left(\partial -\overline{\partial }\right){\rho }_{\uparrow }$. It turns out not to be quite the case.

6 Boundaries, Boundary Conditions, and Boundary Variables

Formulating the sandpile model on a surface with boundaries is important to see how they affect the statistics of the model. The multisite correlations discussed in the previous section are likely to be modified by the presence of a boundary, and by the associated boundary conditions. Moreover, the microscopic variables on a boundary or very close to it will surely have a different behavior from their bulk versions. In the field theoretic description, the boundary fields have to be properly identified, and the way changes of boundary conditions are implemented must be clarified. All this adds to the known set of bulk fields a number of boundary-related fields and offers the opportunity to further test the consistency of their identification by computing mixed correlations combining both types of variables.

Surfaces with boundaries arise in the thermodynamic limit when some of the boundaries of the finite system are not sent off to infinity, unlike the situation considered in the previous section. The simplest case is when only one boundary of the rectangular grid is kept at finite distance, leading to a domain converging to the upper-half plane $\mathrm{ℍ}=\left\{\left(x,y\right)\in {\mathrm{ℝ}}^2:y\mathrm{⩾}0\right\}$. There is only one boundary to care about, and the invariance under horizontal translations is preserved.

From our earlier discussion of conservative vs. dissipative sites, we have already defined two possible boundary conditions: open and closed. Let us recall that the boundary condition is open resp. closed if the boundary sites are dissipative resp. conservative. As before, a boundary open site has $z_i^{*}=4$, while a boundary closed site has $z_i^{*}=z_i=3$. In the scaling limit, it endows $\mathrm{ℍ}$ with a homogeneous boundary condition, open or closed. We also can (and will) consider inhomogeneous boundary conditions, by alternating open and closed stretches on a single boundary.

In terms of height variables, the open boundary condition is equivalent¹⁹ to fix all the boundary heights to 4, whereas the closed condition amounts to constrain the boundary heights not to take the value 4. The fixed boundary condition with boundary heights equal to 1 is not possible (two neighboring 1’s form a forbidden subconfiguration); the fixed boundary conditions with boundary heights equal to 2 and/or 3 should be possible but seem to be difficult to handle in practice.

Two more boundary conditions, defined in the spanning tree description and previously called windy boundary conditions will be discussed at the end of this section (as we will see, they are not so far from the possibility just mentioned, namely, that of having height variables being equal to 2 or 3). No other boundary condition has been considered so far, though it would be very surprising that no other exist.²⁰

6.1 Bulk Variables With Homogeneous Open or Closed Boundary

In this section, we would like to reconsider the multisite height probabilities but on a domain with a boundary, the upper-half plane (UHP) being the simplest case. The principle underlying the calculations on the UHP stays the same as on the full plane. The most essential difference is that the toppling matrix becomes in the thermodynamic limit the Laplacian matrix on the discrete UHP with the appropriate boundary condition, open or closed. In this section, we consider homogeneous boundary conditions only.

To be specific, we choose the boundary row of sites to be located on the horizontal line y = 1 so that the discrete UHP we consider is $\left\{\left(x,y\right)\in {\mathrm{\mathbb{Z}}}^2\hspace{0.17em}|\hspace{0.17em}y\mathrm{⩾}1\right\}$. For either boundary condition, the Laplacian matrix Δ^op or Δ^cl is minus the adjacency matrix of the discrete UHP plus a diagonal matrix, everywhere equal to 4 for the open condition, equal to 4 and 3, respectively, for the bulk and boundary sites for the closed condition. By the method of images, the Green matrices G^op/cl = (Δ^op/cl)⁻¹ can be easily computed in terms of that on the full plane as follows:

$G_{\left(x_1,y_1\right),\left(x_2,y_2\right)}^{\text{op}}=G_{\left(x_1,y_1\right),\left(x_2,y_2\right)}-G_{\left(x_1,y_1\right),\left(x_2,-y_2\right)},(\mathrm{64a})$

$G_{\left(x_1,y_1\right),\left(x_2,y_2\right)}^{\text{cl}}=G_{\left(x_1,y_1\right),\left(x_2,y_2\right)}+G_{\left(x_1,y_1\right),\left(x_2,1-y_2\right)},(\mathrm{64b})$

for y₁ and y₂ > 0. As anticipated, we verify that G^op satisfies the Dirichlet condition, namely, it is odd under the reflection through the line y = 0 and therefore vanishes on it, and that G^cl satisfies the Neumann condition, namely, it is even under the reflection through the line $y=\frac{1}{2}$, inducing a vanishing normal derivative in the scaling limit. The calculations of the previous section can then be generalized to the UHP geometry by using these Green matrices. For height 1 and for the cluster variables, one merely has to use the appropriate Green matrix. For higher heights, the presence of a boundary makes the calculations more complicated because the combinatorics involved is heavier.

The simplest case is the 1-site height probability ${\mathrm{\mathbb{P}}}_1^{\text{op}/\text{cl}}\left(y\right)$ to find a height equal to 1 at a distance y from the boundary. It can be computed by using eq. (27), where Δ⁻¹_Γ is replaced, in the infinite volume limit, by one of the two Green matrices given above. This was historically the first calculations of boundary effects in sandpile models [40], with the result

${\sigma }_1^{\text{op}}\left(y\right)={\mathrm{\mathbb{P}}}_1^{\text{op}}\left(y\right)-{\mathrm{\mathbb{P}}}_1=\frac{{\mathrm{\mathbb{P}}}_1}{4y^2}+\dots ,{\sigma }_1^{\text{cl}}\left(y\right)={\mathrm{\mathbb{P}}}_1^{\text{cl}}\left(y\right)-{\mathrm{\mathbb{P}}}_1=-\frac{{\mathrm{\mathbb{P}}}_1}{4y^2}+\dots (65)$

The analogous results for higher heights were obtained somewhat later [27, 41] and were the first to firmly establish their logarithmic nature. They take the following form, valid for $a\mathrm{⩾}1$,

${\sigma }_a^{\text{op}}\left(y\right)=\frac{1}{y^2}\left(c_a+\frac{d_a}{2}+d_a\text{log}y\right)+\dots ,{\sigma }_a^{\text{cl}}\left(y\right)=-\frac{1}{y^2}\left(c_a+d_a\text{log}y\right)+\dots (66)$

up to terms of order $\mathcal{O}\left(y^{-4}\text{log}y\right)$, which have since been explicitly computed [31], as they enter the calculations of σ^op/cl_a,1 given below. The coefficients are explicitly known and are shown in Table 1. One may check that the relations (44) expressing h₃ and h₄ linearly in terms of h₁ and h₂ are confirmed. The distinctive change of sign between the two boundary conditions and the fact that for fixed a, both are controlled by the same constants c_a and d_a and the equality d₂ = c₁ are striking. As will be explained below and in one of the next sections, all three features will follow from the CFT picture.

TABLE 1

TABLE 1. Numerical coefficients for one-site height probabilities on the UHP, with $\tilde{\gamma }=\gamma +\frac{5}{2}\text{log}2$. They satisfy the relation ${\displaystyle \sum }_a{c}_a={\displaystyle \sum }_a{d}_a=0$.

Let us mention that these lattice calculations have been carried out on another lattice realization of the UHP, namely, on the diagonal upper-half plane $\left\{\left(x,y\right)\in {\mathrm{\mathbb{Z}}}^2\hspace{0.17em}|\hspace{0.17em}y>x\right\}$, for which the method of images allows to explicitly compute the Green matrices for the two boundary conditions. As expected, the dominant terms are exactly the same as above in terms of the Euclidean distance between height 1 and the diagonal boundary, while the subdominant terms are different [31].

The 2-site height correlators in the bulk of the UHP, at sites i₁ = (x₁, y₁) and i₂ = (x₂, y₂), and which involve the same subtractions as before,

${\sigma }_{a,1}^{\text{op}/\text{cl}}\left(x;y_1,y_2\right)={\mathrm{\mathbb{P}}}_{a,1}^{\text{op}/\text{cl}}\left(i_1,i_2\right)-{\mathrm{\mathbb{P}}}_a^{\text{op}/\text{cl}}\left(i_1\right)\hspace{0.17em}{\mathrm{\mathbb{P}}}_1-{\mathrm{\mathbb{P}}}_a\hspace{0.17em}{\mathrm{\mathbb{P}}}_1^{\text{op}/\text{cl}}\left(i_2\right)+{\mathrm{\mathbb{P}}}_a\hspace{0.17em}{\mathrm{\mathbb{P}}}_1,(67)$

have been computed in [31] when the two sites are far from the boundary and far from each other, again using the technique developed in [29]. They depend on three real variables, the horizontal distance x = x₁−x₂ between the two sites and their vertical positions y₁ and y₂. For simplicity however, the lattice calculations have been carried out for two vertically aligned sites, that is, for x = 0.

Defining the two bivariate functions,

$P\left(u,v\right)=\frac{1}{8u^2{v}^2}-\frac{1}{{\left(u-v\right)}^4}-\frac{1}{{\left(u+v\right)}^4}\hspace{0.17em},\hspace{1em}\hspace{1em}Q\left(u,v\right)=\frac{1}{{\left(u-v\right)}^4}-\frac{1}{{\left(u+v\right)}^4},(68)$

the results for a = 1, 2 take the following form, at dominant order:

${\sigma }_{\mathrm{1,1}}^{\text{op}}\left(0;y_1,y_2\right)={\sigma }_{\mathrm{1,1}}^{\text{cl}}\left(0;y_1,y_2\right)=\frac{{\mathrm{\mathbb{P}}}_1^2}{2}\hspace{0.17em}P\left(y_1,y_2\right)+\dots ,(\mathrm{69a})$

$\begin{array}{c}{\sigma }_{\mathrm{2,1}}^{\text{op}/\text{cl}}\left(0;y_1,y_2\right)=\frac{{\mathrm{\mathbb{P}}}_1^2}{2}\left[P\left(y_1,y_2\right)\left(\text{log}{y}_1+\gamma +\frac{5}{2}\text{log}2\right)+Q\left(y_1,y_2\right)\text{log}\left|\frac{y_2+y_1}{y_2-y_1}\right|\right]\\ +\frac{H^{\text{op}/\text{cl}}\left(y_1^2,y_2^2\right)}{y_1^2{y}_2^2{\left(y_1^2-y_2^2\right)}^4}+\dots ,\end{array}(\mathrm{69b})$

where H^op(u,v) and H^cl(u,v) are homogeneous polynomials of degree 4 in u, v, with explicitly known coefficients. The results for a = 3, 4 take the same form with different coefficients and confirm once more the linear relations (44).

Let us discuss these results in the CFT picture, using what we already know about the height fields $h_a\left(z,\overline{z}\right)$. For a homogeneous boundary condition like here, boundary CFT prescribes to compute bulk correlators on the UHP by viewing a field $\varphi \left(z,\overline{z}\right)$ with conformal weights $\left(h,\overline{h}\right)$ as the product ${\varphi }_h\left(z\right){\varphi }_{\overline{h}}\left(\overline{z}\right)$ of two chiral fields of weights h and $\overline{h}$, respectively, located at the points z and $\overline{z}$ (the latter being in the lower-half plane) [42]. A correlation function of n bulk (non-chiral) fields on the UHP can then be computed as a correlation of 2n chiral fields on the full plane; the correlation appropriate for the boundary condition under consideration is accordingly selected in the solution space of these 2n-correlators.

The above prescription must however be adapted in the case of logarithmic fields because the chiral factorization is not consistent with the non-diagonal action of L₀. Indeed, let us consider a logarithmic pair $\left(\mathrm{\Phi }\left(z,\overline{z}\right),\text{Ψ}\left(z,\overline{z}\right)\right)$. If we factorize the logarithmic partner as $\text{Ψ}\left(z,\overline{z}\right)={\psi }_h\left(z\right){\psi }_{\overline{h}}\left(\overline{z}\right)$, we find from the action of L₀,

$L_0\text{Ψ}=\left(L_0{\psi }_h\right){\psi }_{\overline{h}}=\left(h{\psi }_h+\lambda \hspace{0.17em}{\varphi }_h\right){\psi }_{\overline{h}}=h\text{Ψ}+\lambda \hspace{0.17em}{\varphi }_h{\psi }_{\overline{h}},(70)$

That the chiral factorization of the primary partner is $\mathrm{\Phi }={\varphi }_h{\psi }_{\overline{h}}$. The same argument with ${\overline{L}}_0$ shows that an equally good factorization is $\mathrm{\Phi }={\psi }_h{\varphi }_{\overline{h}}$.

Let us first see how this works for σ^op_a(y). Their dominant terms, made explicit in relations (65) and (66), should correspond to ${\langle h_a\left(z,\overline{z}\right)\rangle }_{\text{op}}$. Note that unlike the correlations on the plane discussed in Section 5, we do not insert dissipation at infinity since the whole boundary is open and therefore dissipative. From the prescription recalled above, these 1-point functions should have the general form of chiral 2-point functions. If ψ and ϕ denote chiral versions of the height 2 and height 1 fields, respectively, with $h=\overline{h}=1$, the chiral factorizations of the height fields h₁ and h₂ read $h_1\left(z,\overline{z}\right)=\psi \left(z\right)\varphi \left(\overline{z}\right)$ and $h_2\left(z,\overline{z}\right)=\psi \left(z\right)\psi \left(\overline{z}\right)$. The CFT formalism gives the general forms (19) as follows:

${\langle h_1\left(z,\overline{z}\right)\rangle }_{\text{op}}=\langle \varphi \left(z\right)\psi \left(\overline{z}\right)\rangle =\frac{B}{{\left(z-\overline{z}\right)}^2},{\langle h_2\left(z,\overline{z}\right)\rangle }_{\text{op}}=\langle \psi \left(z\right)\psi \left(\overline{z}\right)\rangle =\frac{C-2\lambda B\log \left(z-\overline{z}\right)}{{\left(z-\overline{z}\right)}^2}.(71)$

With the value $\lambda =-\frac{1}{2}$ noted in Section 5, these forms exactly reproduce the lattice results (66), including the relation d₂ = c₁ = B.

For σ^cl_a(y) and since the closed boundary is not dissipative, we insert by hand dissipation at infinity so that σ^op_a(y) should be given by ${\langle h_a\left(z,\overline{z}\right)\omega \left(\infty \right)\rangle }_{\text{cl}}$. Using the same chiral factorization as above leads to a 3-point function. The selection of a physically sensible solution leads to the same general form as for the open boundary condition [27].

The conformal calculations required to account for σ^op_a,1 are only technically more involved. The needed chiral correlators are $\langle \varphi \left(z_1\right)\psi \left({\overline{z}}_1\right)\varphi \left(z_2\right)\psi \left({\overline{z}}_2\right)\rangle$ for a = 1 and $\langle \psi \left(z_1\right)\psi \left({\overline{z}}_1\right)\varphi \left(z_2\right)\psi \left({\overline{z}}_2\right)\rangle$ for a = 2. Both can be computed from the primary nature of the chiral field ϕ, as established in Section 5, by solving a second-order differential equation and selecting the appropriate solution. As the details are given in [31], we merely give the results, valid for any relative positions of the two heights, z₁ = (x₁, y₁) and z₂ = (x₂, y₂)

${\langle h_1\left(z_1,{\overline{z}}_1\right)h_1\left(z_2,{\overline{z}}_2\right)\rangle }_{\text{op}}=\frac{{\mathrm{\mathbb{P}}}_1^2}{2}\left\{\frac{2}{{\left(z_1-{\overline{z}}_1\right)}^2{\left(z_2-{\overline{z}}_2\right)}^2}-\frac{1}{|z_1-z_2{|}^4}-\frac{1}{|z_1-{\overline{z}}_2{|}^4}\right\},(72)$

and

$\begin{array}{l}{\langle h_2\left(z_1,{\overline{z}}_1\right)h_1\left(z_2,{\overline{z}}_2\right)\rangle }_{\text{op}}=\frac{{\mathrm{\mathbb{P}}}_1}{32y_1^2{y}_2^2}\frac{t^4-2t^3+4t-2}{{\left(1-t\right)}^2}\left[\frac{3\left(3\pi -10\right)}{2{\pi }^3}-{\mathrm{\mathbb{P}}}_1\left(\log y_1+\gamma +\frac{5}{2}\log 2\right)\right]\\ +\frac{{\mathrm{\mathbb{P}}}_1^2}{64y_1^2{y}_2^2}\left[\frac{t^3\left(t-2\right)}{{\left(1-t\right)}^2}\left(\text{log}\left(1-t\right)+\frac{y_1}{2y_2}\right)-\frac{t^2}{1-t}\right],t=-\frac{4y_1{y}_2}{|z_1-z_2{|}^2}.\end{array}(73)$

One can check that setting x₁ = x₂ exactly reproduces the lattice results ${\sigma }_{\mathrm{1,1}}^{\text{op}}\left(0;y_1,y_2\right)$ and ${\sigma }_{\mathrm{2,1}}^{\text{op}}\left(0;y_1,y_2\right)$ reported above.

The analogous calculation for the closed boundary has been carried out in the case a = 1, yielding the same expression as for the open boundary. No calculation however has been successful for a = 2 as it involves a nontrivial 5-point chiral correlator (in this case, the dissipation field ω must be added).

Similar calculations with isolated bulk dissipation, instead of height variables, have been considered in [36]; it was found in all cases that the CFT predictions compare successfully with the lattice results.

6.2 Changing the Boundary Condition

We have considered so far two different boundary conditions, the open and closed conditions. This allows addressing a fundamentally new issue, namely, how to think of a change of boundary condition, both on the lattice and in the emerging field theory. Like in the previous section, we consider the UHP.

We have seen that the calculation of correlations on the UHP, of height or dissipation variables, involves the use of the appropriate Laplacian (toppling) matrix and its inverse. On the lattice, the way we can change the boundary condition at a boundary site i is thus fairly clear: since an open boundary site has ${\mathrm{\Delta }}_{ii}=z_i^{*}=4$ and a closed one has ${\mathrm{\Delta }}_{ii}=z_i^{*}=3$, we simply lower by 1 the diagonal entry Δ_i,i to close an open site, and we increase it by 1 to open a closed site (as we did in Section 5.4 to introduce dissipation at bulk sites). We do it either way for n consecutive boundary sites to change the boundary condition on an interval I of length n, that is, we do the following change on the toppling matrix Δ → Δ ± D_I, where D_I implements the diagonal shifts described above.

Let us examine the effect of closing n consecutive sites in an otherwise open boundary. We decide to measure this effect as in Section 5.4, namely, by comparing the number of recurrent configurations before and after the closing of n sites. So, we want to compute the ratio Z_op(n)/Z_op. At finite volume, the two partition functions can be computed as determinants of the corresponding toppling matrices on rectangular grids, with say four open boundaries in the case of Z_op, and with n closed sites inserted on the lower boundary for Z_op(n). As usual, we can readily write the infinite volume limit of the ratio as follows:

$\frac{Z_{\text{op}}\left(n\right)}{Z_{\text{op}}}=\frac{{\text{det}\Delta }^{\text{op}}\left(n\right)}{{\text{det}\Delta }^{\text{op}}}=\frac{\text{det}\left[{\mathrm{\Delta }}^{\text{op}}-D_{I_n}\right]}{{\text{det}\Delta }^{\text{op}}}=\text{det}\left[\mathbb{I}-G^{\text{op}}{D}_{I_n}\right]=\text{det}{\left[\mathbb{I}-G^{\text{op}}\right]}_{i,j\in I_n},(74)$

where ${\left(D_{I_n}\right)}_{ij}={\delta }_{ij}$ for $i,j\in I_n$ is zero elsewhere, and $I_n=\left\{\left(\ell ,1\right):1\mathrm{⩽}\ell \mathrm{⩽}n\right\}$ is the set of sites being closed. Using the relation (64a) expressing G^op in terms of the Green matrix on the full plane ${\mathrm{\mathbb{Z}}}^2$, the matrix in the determinant reads

${\left(\mathbb{I}-G^{\text{op}}\right)}_{i,j\in I}={\left({\delta }_{\ell ,{\ell }^{\text{′}}}-G_{\left(\ell ,1\right),\left({\ell }^{\text{′}},1\right)}+G_{\left(\ell ,1\right),\left({\ell }^{\text{′}},-1\right)}\right)}_{1\mathrm{⩽}\ell ,{\ell }^{\text{′}}\mathrm{⩽}n}.(75)$

By the horizontal translation invariance of G^op, this is a Toeplitz matrix of the form $a_{\ell -\ell \text{′}}$. Using standard results on the Green matrix on the plane, one finds that the entries a_m are the Fourier coefficients of the following symbol, which has a so-called Fisher–Hartwig singularity,

${\sigma }^{\text{op}}\left(k\right)=\sqrt{1-\cos k}\cdot \left\{\sqrt{3-\cos k}-\sqrt{1-\cos k}\right\}.(76)$

For large n, the asymptotics of such Toeplitz determinants is well-known (see f.i., [43]) and leads to the following result [44]:

$\frac{Z_{\text{op}}\left(n\right)}{Z_{\text{op}}}\simeq A\hspace{0.17em}{n}^{1/4}\hspace{0.17em}{\text{e}}^{-\frac{2\text{G}}{\pi }n},\hspace{1em}\hspace{1em}n\gg 1,(77)$

with G = 0.915965, the Catalan constant. The proportionality constant A is explicitly known but is unimportant here.

What if we consider the opposite situation in which we open n consecutive sites of a closed boundary? Reasoning as above, we quickly get the corresponding ratio,

$\frac{Z_{\text{cl}}\left(n\right)}{Z_{\text{cl}}}=\frac{{\text{det}\Delta }^{\text{cl}}\left(n\right)}{{\text{det}\Delta }^{\text{cl}}}=\frac{\text{det}\left[{\mathrm{\Delta }}^{\text{cl}}+D_{I_n}\right]}{{\text{det}\Delta }^{\text{cl}}}=\text{det}{\left[\mathbb{I}+G^{\text{cl}}\right]}_{i,j\in I_n},(78)$

which is also a Toeplitz determinant. However, this one is infinite—each entry is infinite—for the same reason we have pointed out in Section 5.4. Adopting the same point of view, we similarly evaluate the effect of opening n sites with respect to the situation where only one site is open. One therefore considers instead the ratio $\frac{Z_{\text{cl}}\left(n\right)}{Z_{\text{cl}}\left(1\right)}$, which one can write as

$\frac{Z_{\text{cl}}\left(n\right)}{Z_{\text{cl}}\left(1\right)}=\frac{1}{b_0}\text{det}{\left(b_{\ell -{\ell }^{\text{′}}}\right)}_{1\mathrm{⩽}\ell ,\ell \text{′}\mathrm{⩽}n},(79)$

where the entries b_m are the Fourier coefficients of a symbol σ^cl(k) given by

${\sigma }^{\text{cl}}\left(k\right)=\frac{1}{2}{\left(1-\cos k\right)}^{\alpha }\cdot \left\{\sqrt{3-\cos k}+\sqrt{1-\cos k}\right\},\hspace{1em}\hspace{1em}\alpha =-\frac{1}{2}.(80)$

Its Fourier coefficients are well-defined for $\alpha >-\frac{1}{2}$, diverge in the limit $\alpha \to -\frac{1}{2}$, but nonetheless keep the ratio (79) finite. Remarkably, for large n, it takes the form [44].

$\frac{Z_{\text{cl}}\left(n\right)}{Z_{\text{cl}}\left(1\right)}\simeq A\hspace{0.17em}{n}^{1/4}\hspace{0.17em}{\text{e}}^{\frac{2\text{G}}{\pi }\left(n-1\right)},\hspace{1em}\hspace{1em}n\gg 1,(81)$

for the same constant A as above.

Before discussing the CFT side, let us remark that the exponential factors in equations (77) and (81) are expected. On a finite N × N grid, all four partition functions (numerators and denominators) are asymptotically dominated by the bulk free energy, given by ${\text{e}}^{\frac{4\text{G}}{\pi }{N}^2}$, as mentioned in Section 2. These terms drop out in the ratios. The next correction is related to the boundary free energy f_b (per site) and takes the form ${\text{e}}^{4N{f}_{\text{b}}}$ in case the boundary condition b is the same at all boundary sites. For the partition functions considered above, the boundary conditions only differ on the lower edge of the grid so that for large $N\gg n\gg 1$, the ratios are asymptotic to

$\frac{Z_{\text{op}}\left(n\right)}{Z_{\text{op}}}\simeq \frac{{\text{e}}^{\left(N-n\right)f_{\text{op}}+n{f}_{\text{cl}}}}{{\text{e}}^{N{f}_{\text{op}}}}={\text{e}}^{-n\left(f_{\text{op}}-f_{\text{cl}}\right)},\hspace{1em}\hspace{1em}\frac{Z_{\text{cl}}\left(n\right)}{Z_{\text{cl}}\left(1\right)}\simeq \frac{{\text{e}}^{\left(N-n\right)f_{\text{cl}}+n{f}_{\text{op}}}}{{\text{e}}^{\left(N-1\right)f_{\text{cl}}+f_{\text{op}}}}={\text{e}}^{\left(n-1\right)\left(f_{\text{op}}-f_{\text{cl}}\right)}.(82)$

The free energies f_op and f_cl represent (the logarithm of) the effective number of values taken by the boundary heights in the set of recurrent configurations. The number of possible values taken by the height at an open boundary site is 4, and is 3 at a closed site. If these numbers of values get effectively reduced in the set of recurrent configurations, one should expect that the number of values at an open boundary site remains larger than that at a closed site, implying f_op−f_cl > 0. An explicit calculation [44] confirms this and yields $f_{\text{op}}-f_{\text{cl}}=\frac{2\text{G}}{\pi }$, in agreement with the above results. To fix the ideas, the effective number of values taken by a boundary height is ${\text{e}}^{f_{\text{op}}}=3.70$ at an open site, and ${\text{e}}^{f_{\text{cl}}}=2.07$ at a closed site.

In the CFT approach, a change of the boundary condition at x, from condition a to condition b, is implemented by the insertion in the correlators of a specific field ϕ^a,b(x). Such boundary condition changing fields²¹ are usually expected to be chiral primary fields and satisfy ϕ^a,b(x) = ϕ^b,a(x) when the boundary conditions a and b do not carry an intrinsic orientation (see Section 6.4 for counterexamples). The insertion of the product ϕ^a,b(x₁)ϕ^b,a(x₂) accounts for the change at x₁ from condition a to condition b and then back from b to a at x₂ but does not account for the exponential terms related to the difference of boundary free energies of condition a vs. condition b, namely, the terms we have just discussed in the previous paragraph. These are clearly nonuniversal, that is, depend on the specific model under consideration, and cannot be accounted for by the underlying CFT, which itself applies to all the models in the universality class to which the sandpile model belongs.

It follows that the effect of changing the boundary condition given above, in which we omit the exponential terms, should correspond to the 2-point function $\langle {\varphi }^{\text{op},\text{cl}}\left(0\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(n\right)\rangle =\langle {\varphi }^{\text{cl},\text{op}}\left(0\right)\hspace{0.17em}{\varphi }^{\text{op},\text{cl}}\left(n\right)\rangle$. The two are indeed equal on the lattice and asymptotic to $n^{\frac{1}{4}}$, and from this, we infer that the boundary condition changing field ϕ^op,cl(x) = ϕ^cl,op(x) is a chiral conformal field of weight $h=-\frac{1}{8}$, with a correlator given by

${\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_2\right)\rangle }_{\text{op}}=A\hspace{0.17em}{\left|x_1-x_2\right|}^{1/4}.(83)$

For physical reasons, we might worry about having a correlator that actually increases with the distance, suggesting somehow the existence of a strange interaction that would get stronger at larger distances. There is nothing strange however as it does not really correspond to the physical correlation of two observables. As said above, the field ϕ^op,cl is expected to be primary. As usual, this conjecture can be put to test: the consequences of this statement must have a match in the lattice properties of the model.

One of the strongest consequence of the primary nature of ϕ^op,cl and the assumed structure of the conformal module that contains it is that any correlator where this field appears must satisfy a second-order partial differential equation,²² the precise form of which depends on the other fields involved. A first and simple test is to look at a 4-point function,²³ for instance $\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right){\varphi }^{\text{cl},\text{op}}\left(x_2\right)\hspace{0.17em}{\varphi }^{\text{op},\text{cl}}\left(x_3\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_4\right)\rangle$, which should describe the effect of closing the sites on two disjoint interval [x₁, x₂] and [x₃, x₄] in the otherwise open boundary of the UHP. Using the global conformal invariance, one can reduce the partial differential equation to a second-order ordinary differential equation. In the two-dimensional solution space, we select the only solution which reduces to the product $\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_2\right)\rangle \langle {\varphi }^{\text{op},\text{cl}}\left(x_3\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_4\right)\rangle$ when the two intervals are infinitely distant. This unique solution reads (with x_ij = x_i − x_j)

${\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_2\right)\hspace{0.17em}{\varphi }^{\text{op},\text{cl}}\left(x_3\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_4\right)\rangle }_{\text{op}}=\frac{2A^2}{\pi }{\left(x_{12}\hspace{0.17em}{x}_{34}\right)}^{1/4}\hspace{0.17em}{\left(1-t\right)}^{1/4}\hspace{0.17em}K\left(t\right),\hspace{1em}t\equiv \frac{x_{12}{x}_{34}}{x_{13}{x}_{24}},(84)$

where $K\left(t\right)={{{\displaystyle \int }}^{\text{}}}_0^{\frac{\pi }{2}}\frac{\text{d}\theta }{\sqrt{1-t{\sin }^2\theta }}$ is the complete elliptic integral.

To compare with a lattice calculation, we take x_i integers, with x₂₁, x₃₂, and x₄₃ all large, and try to compute the determinant in eq. (74) with $I=I_1{\displaystyle \cup }{I}_2$, the union of the two intervals [x₁, x₂] and [x₃, x₄]. This determinant is no longer Toeplitz, which makes it difficult to compute its asymptotics analytically. Dividing it by the prefactor ${\left(x_{12}\hspace{0.17em}{x}_{34}\right)}^{1/4}$, it can however be evaluated numerically as a function of t by varying the lengths of the intervals and their separation distance. The agreement with eq. (84) is more than satisfactory [44].

The opposite situation—two open intervals in a closed boundary—has also been considered. The appropriate 4-point function can be obtained from eq. (84) by making a simple cyclic permutation (x₁,x₂,x₃,x₄) → (x₄,x₁,x₂,x₃), with the result that K(t) gets replaced by K(1 − t). An equally successful agreement was observed [36]. Many other cross-checks have been done, confirming that the open/closed boundary condition changing field is indeed a primary field with conformal weight $h=-\frac{1}{8}$. One of them, particularly convincing, is presented in the next section.

6.3 Bulk Variables With Inhomogeneous Boundary

In Section 6.1, we have computed the lattice 1- and 2-site height probabilities on the UHP, with either the open or the closed boundary condition. Here, we would like to revisit these results in light of what we have learned of the boundary condition field, in terms of which one should be able to relate the probabilities for the two boundary conditions. In particular, we would like to understand the 1-site probabilities σ^op_a(y) and σ^cl_a(y),

${\sigma }_a^{\text{op}}\left(y\right)=\frac{1}{y^2}\left(c_a+\frac{d_a}{2}+d_a\log y\right)+\dots ,\hspace{1em}\hspace{1em}{\sigma }_a^{\text{cl}}\left(y\right)=-\frac{1}{y^2}\left(c_a+d_a\log y\right)+\dots (85)$

One can do this by computing, on the CFT side, a more general probability, namely, we look for the probability to find a height equal to a at a distance y from the boundary, when the boundary condition is mixed, namely, open everywhere, except on the interval [x₁, x₂] where the condition is closed. The two homogeneous open and closed conditions can be recovered in the limits x₁ → x₂ and x₁ → −∞, x₂ → ∞. Let us denote by ${\langle h_a\left(z,\overline{z}\right)\rangle }_{\text{mix}}$ the corresponding quantity in the CFT, given by

${\langle h_a\left(z,\overline{z}\right)\rangle }_{\text{mix}}=\frac{{\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_2\right)\hspace{0.17em}{h}_a\left(z,\overline{z}\right)\rangle }_{\text{op}}}{{\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_2\right)\rangle }_{\text{op}}},(86)$

where the division by ${\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_2\right)\rangle }_{\text{op}}$ comes from the fact that we want to evaluate the probability to have a height 1 in front of a mixed boundary condition and not the combined effects of having a height 1 and the closing the boundary between x₁ and x₂. The denominator is known from eq. (83).

To compute the numerator, we represent the height fields in terms of the chiral fields as $h_1\left(z,\overline{z}\right)=\varphi \left(z\right)\psi \left(\overline{z}\right)$ and $h_2\left(z,\overline{z}\right)=\psi \left(z\right)\psi \left(\overline{z}\right)$ (as usual, considering the heights a = 1, 2 is enough) and write the differential equation satisfied by the two ensuing 4-point correlators, as a consequence of the primary nature of ϕ^op,cl. Because ψ is the chiral logarithmic partner of ϕ, the general solution for ${\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_2\right)\hspace{0.17em}\psi \left(z\right)\hspace{0.17em}\psi \left(\overline{z}\right)\rangle }_{\text{op}}$ in fact depends on that of ${\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_2\right)\hspace{0.17em}\varphi \left(z\right)\hspace{0.17em}\psi \left(\overline{z}\right)\rangle }_{\text{op}}$.

All calculations done, one finds that they depend on two integration constants c₂ and d₂ in such a way that the ratios in eq. (86) take the following forms, where y = Rez [27]:

${\langle h_1\left(z,\overline{z}\right)\rangle }_{\text{mix}}=\frac{d_2}{2y^2}\hspace{0.17em}\frac{1+t}{\sqrt{t}},\hspace{1em}\hspace{1em}\hspace{1em}t=\frac{\left(x_1-\overline{z}\right)\left(x_2-z\right)}{\left(x_1-z\right)\left(x_2-\overline{z}\right)},(\mathrm{87a})$

${\langle h_2\left(z,\overline{z}\right)\rangle }_{\text{mix}}=\frac{1}{2y^2}\hspace{0.17em}\frac{1+t}{\sqrt{t}}\hspace{0.17em}\left\{c_2+\frac{d_2}{8}\frac{{\left(1+\sqrt{t}\right)}^2}{\sqrt{t}}+d_2\log y-d_2\frac{\text{i}y\left(1-t\right)}{\left(x_1-x_2\right)}\left[\frac{1}{1+t}-\frac{1}{2\sqrt{t}}\right]\right\}.(\mathrm{87b})$

Although t is complex, both expressions are real on account of $t^{\text{*}}=1/t$.

Let us now discuss the above two limits x₁ → x₂ and x₁ → −∞, x₂ → ∞. For convenience, we set x₁ = −x₂ and examine the limits x₂ → 0⁺ and x₂ → +∞. To compute the two limits, the important thing to notice is that the complex variable t, now equal to

$t=\frac{\left(x_2+\overline{z}\right)\left(x_2-z\right)}{\left(x_2+z\right)\left(x_2-\overline{z}\right)},(88)$

has complex norm equal to 1 and loops anticlockwise around the origin as x₂ varies from 0⁺ to +∞, starting from 1+ 0i to 1 − 0i. It follows that t itself goes to 1 in both limits, but $\sqrt{t}$ goes to +1 when x₂ → 0⁺ and goes to −1 when x₂ → +∞. The actual limits yield

${\langle h_1\left(z,\overline{z}\right)\rangle }_{\text{op}}=\underset{x_2\to 0^{+}}{\text{lim}}{\langle h_1\left(z,\overline{z}\right)\rangle }_{\text{mix}}=\frac{d_2}{y^2},{\langle h_1\left(z,\overline{z}\right)\rangle }_{\text{cl}}=\underset{x_2\to +\infty }{\text{lim}}{\langle h_1\left(z,\overline{z}\right)\rangle }_{\text{mix}}=-\frac{d_2}{y^2},(\mathrm{89a})$

${\langle h_2\left(z,\overline{z}\right)\rangle }_{\text{op}}=\frac{1}{y^2}\left(c_2+\frac{d_2}{2}+d_2\log y\right),{\langle h_2\left(z,\overline{z}\right)\rangle }_{\text{cl}}=-\frac{1}{y^2}\left(c_2+d_2\log y\right),(\mathrm{89b})$

in complete agreement with the lattice results: the change of the overall sign between the open and closed boundary conditions, the specific dependence on the two coefficients c₂ and d₂, and the equality c₁ = d₂ are all accounted for! The conformal approach however cannot fix the two coefficients c₂ and d₂; these must be determined by lattice calculations.

The expressions (87) can also be tested in situations where the boundary condition along the real axis is no longer homogeneous. A particularly instructive case is when the boundary condition is closed on the negative part of the real axis and open on the positive part, corresponding to the limits x₁ → −∞ and x₂ → 0. The conformal transformation $w=\frac{L}{\pi }\text{log}\hspace{0.17em}z$ can be used to map the UHP onto an infinite strip of width L, with open boundary condition on the left side and closed on the right side. The conformal transformation rules of the fields involved being known, the expressions (87) can be transformed to the strip and compared with numerical simulations on a truncated (and large) strip (exact calculations on the lattice are not available). It was found [27] that the conformal predictions and the numerical plots match remarkably well, thereby confirming once more all the field identifications made so far.

6.4 Wind on the Boundary

The open and closed boundary conditions are very natural as the very definition of the sandpile model uses dissipative and conservative sites. One may wonder what other type of boundary condition could be thought of. Perhaps, we could think of alternating open and closed boundary sites; we expect however that such a boundary condition would flow to the open condition in the scaling limit, as numerical experiments confirm. We have already commented on the possibility to uniformly fix the boundary heights. Fixing the boundary heights to 2 or to 3, or even to 2 or 3, seems difficult. The two boundary conditions, different from open and closed, which have been considered in [45], are in fact closely related, but not quite identical, to the third possibility. They are fixed boundary conditions but in the language of spanning trees.

We recall that in a rooted spanning tree, there is exactly one outgoing arrow at each vertex. The two new boundary conditions, noted ← and →, force the outgoing arrows at the boundary sites to be uniformly left or uniformly right.²⁴ In terms of height values, either condition means that none of the boundary sites has height 1 (because each boundary site has an ingoing arrow) or height 4 (because the burning algorithm would imply that the arrow is pointing down, toward the root). The converse is however not true: recurrent configurations with height values equal to 2 or 3 on the boundary do not necessarily have boundary arrows uniformly oriented.

The way the orientation of an edge can be forced has been briefly discussed in Section 5.7. This allows evaluating the effects of inserting a stretch of left or right arrows into an open or a closed boundary, similarly to what we did in Section 6.2. We refer the reader to Ref [45] for details of the analysis and restrict here to a summary of the results.

An obvious but unusual feature of the boundary conditions ← and → is that they are intrinsically oriented. It implies that the boundary condition changing field ϕ^a,→ turning the boundary condition from a to → may not be the same as the field ϕ^→,a implementing the opposite change. With $\text{a},\text{b}\in \left\{\text{op},\text{cl},\leftarrow ,\to \right\}$, this makes potentially twelve distinct fields ϕ^a,b (the fields ϕ^a,a are just the identity). We already know ϕ^op,cl = ϕ^cl,op, and likewise, if $\text{a}\in \left\{\text{op},\text{cl}\right\}$ is unoriented and $\text{b}\in \left\{\leftarrow ,\to \right\}$ is oriented, we expect the identifications ϕ^a,→ = ϕ^←,a and ϕ^a,← = ϕ^→,a on the basis of a left–right reflection symmetry. These identifications have been confirmed and reduce the number of distinct fields to seven.

There is an additional subtlety for the field that changes the orientation from → to ←. Indeed the right and left arrows $\circ \to •\leftarrow \circ$ point to the same boundary site (in black), and whether that site is open or closed may be relevant. Indeed, if it is open, the flow of arrows, which eventually terminates at the root, can go directly to the root; if it is closed, it must necessarily go upward into the bulk of the UHP. In the two cases, the macroscopic configurations of arrows are different. Thus, we should distinguish two different fields, ${\varphi }^{\to \stackrel{\text{op}}{,}\leftarrow }$ and ${\varphi }^{\to \stackrel{\text{cl}}{,}\leftarrow }$. The detailed analysis confirms that they are distinct fields as their conformal weights are different.

We therefore have eight distinct boundary condition changing fields. A mix of analytical calculations and numerical simulations has been used to determine the conformal weights of these eight fields. The results are given in Table 2.

TABLE 2

TABLE 2. Conformal weights of the fields ϕ^a,b which implement a change of boundary condition from a (row label) to b (column).

The more delicate question of the exact nature of all these fields has been addressed by considering the fusion of the representations to which they belong. Loosely speaking, the fusion rules implement the composition law ${\varphi }^{\text{a},\text{b}}*{\varphi }^{\text{b},\text{c}}\simeq {\varphi }^{\text{a},\text{c}}$ of boundary condition changing fields in the limit where the insertion points coincide. The ensuing consistency conditions suggest that all of them are primary fields, except two, which could belong to logarithmic representations (i.e., reducible indecomposable with Jordan cells). Also, the fields of weight 0 are nontrivial, that is, not equal to the identity (they are found to be degenerate at level 3). Relying on these proposals, various 4-point correlators have been computed and successfully compared with numerical simulations. We refer the readers to Ref [45] for more details on these specific points.

6.5 Boundary Height Variables

The boundary condition changing fields are not the only ones to live on a boundary. The lattice model includes observables in the bulk as well as on the boundaries. Those in the bulk have been discussed at length and give rise in the scaling limit to non-chiral fields $\mathrm{\Phi }\left(z,\overline{z}\right)$, characterized by a pair of conformal weights $\left(h,\overline{h}\right)$; those on the boundaries give rise to boundary, chiral fields Φ^a(x), characterized by a single conformal dimension h_a. In general, the nature of the boundary field associated with a boundary observable and its conformal weight depend on the boundary condition.

In the Abelian sandpile model, only the boundary fields arising from the height variables and from the insertion of isolated dissipation have been studied on the UHP. In both cases, only open and closed boundaries have been considered.

The case of isolated dissipation is simpler and has been examined in detail in Ref [36], where isolated dissipation has been considered on a closed boundary only. The calculation proceeds much like that for the bulk, reviewed in Section 5.4, for which the same regularization is used. The results are similar: the dissipation field ω^cl(x) turns out to be a chiral field with conformal weight h_cl = 0 and is a logarithmic partner of the identity. The multipoint correlators involve various combinations of logarithms like their bulk versions. On an open boundary, already dissipative, the dissipation field ω^op(x) is expected to be a descendant of the identity. Isolated dissipation is the simplest observable that can be associated and computed in terms of a local defect matrix. This, from what we have said in Section 5.6 of the minimal clusters, suggests that both ω^cl and ω^op can be realized as local fields in the symplectic fermions. It was indeed shown that ${\omega }^{\text{cl}}\sim \theta \tilde{\theta }$ [36] and ${\omega }^{\text{op}}\sim \partial \theta \partial \tilde{\theta }$ [46] reproduce all known correlations. We note that the latter is proportional to the boundary stress–energy tensor T(x) of the symplectic theory, a non-primary chiral field of weight h_op = 2 and a descendant of the identity since $T\left(x\right)\sim \left(L_{-2}\hspace{0.17em}\mathbb{I}\right)\left(x\right)$.

Boundary height variables are more complicated than dissipation but simpler than the bulk height variables. The first results have been derived by Ivashkevich in Ref [47], where the one- and two-site height probabilities on open and closed boundaries were obtained. The probabilities involving heights 1 only are no more complicated than in the bulk and can be easily obtained by using a defect matrix. As could be expected, probabilities for higher heights are more difficult.

On a boundary, heights larger or equal to 2 are characterized as in the bulk, namely, in terms of the number of predecessors among their nearest neighbors. So, it leads essentially to the same problems of computing nonlocal contributions. Both in the bulk and on a boundary, one can write linear identities expressing combinations of nonlocal contributions in terms of local ones, themselves calculable with a defect matrix. In turn, the nonlocal contributions can be used to calculate probabilities. In the bulk, the linear system is underdetermined and cannot be inverted to provide the required nonlocal contributions and then the probabilities themselves. The main observation made in Ref [47] was that on a boundary, the linear system can be inverted and therefore allows computing the height probabilities and correlations in terms of local contributions only. The following results were obtained.

The 1-site height probabilities on the boundary, open and closed, of the infinite UHP were computed exactly. For comparison purposes, we reproduce here their numerical values (the exact values can be found in Ref [47]) and recall those in the bulk, as given in Section 5.1:

${\mathrm{\mathbb{P}}}_1=0.073\hspace{0.17em}63,\hspace{1em}{\mathrm{\mathbb{P}}}_2=0.173\hspace{0.17em}90,\hspace{1em}{\mathrm{\mathbb{P}}}_3=0.306\hspace{0.17em}29,\hspace{1em}{\mathrm{\mathbb{P}}}_4=0.446\hspace{0.17em}17,(90)$

${\mathrm{\mathbb{P}}}_1^{\text{op}}=0.103\hspace{0.17em}82,\hspace{1em}{\mathrm{\mathbb{P}}}_2^{\text{op}}=0.216\hspace{0.17em}57,\hspace{1em}{\mathrm{\mathbb{P}}}_3^{\text{op}}=0.316\hspace{0.17em}23,\hspace{1em}{\mathrm{\mathbb{P}}}_4^{\text{op}}=0.363\hspace{0.17em}38,(91)$

${\mathrm{\mathbb{P}}}_1^{\text{cl}}=0.113\hspace{0.17em}38,\hspace{1em}{\mathrm{\mathbb{P}}}_2^{\text{cl}}=0.318\hspace{0.17em}31,\hspace{1em}{\mathrm{\mathbb{P}}}_3^{\text{cl}}=0.568\hspace{0.17em}31.(92)$

On the open boundary, for which the comparison makes more sense, lower heights are thus more likely.

Mixed 2-site correlators τ^op_a,b(x₁, x₂) and τ^cl_a,b(x₁, x₂) on an open or a closed boundary were also computed in Ref [47]; all of them were found to decay like ${\left|x_1-x_2\right|}^{-4}$. Although logarithmic conformal field theory was in its infancy at the time, it indicates in hindsight that unlike their bulk cousins, boundary height fields are not logarithmic. This is also in agreement with the fact explained above that boundary height correlations can be fully computed in terms of local contributions.

The decay of the 2-site correlators strongly suggest that all boundary height fields, whatever the boundary condition, have a conformal dimension equal to h_op = h_cl = 2. But like for the other observables discussed so far, we are interested to know the precise nature of the associated fields. Since the multisite boundary height probabilities appear to be calculable in terms of local contributions using defect matrices, it suggests again to look for field candidates constructed out from the symplectic free fermions $\theta ,\tilde{\theta }$. This was done independently in Ref [46] and in Ref [48], following however two different approaches: the former computed various 3-point correlators, whereas the latter considered 2-point correlators only but in the massive extension of the sandpile model (see Section 7.1). The massive extension indeed allows distinguishing more efficiently different fields which would otherwise have the same 2-point correlators in the non-massive (critical) limit.

The results are as follows. The four height fields on an open boundary are all proportional to a single field,

$h_a^{\text{op}}\left(x\right)=O_a\hspace{0.17em}\partial \theta \hspace{0.17em}\partial \tilde{\theta },\hspace{1em}\hspace{1em}1\mathrm{⩽}a\mathrm{⩽}4,(93)$

with explicit normalization constants O_a and where the $\theta ,\tilde{\theta }$ fields satisfy the Dirichlet boundary condition. Thus, on an open boundary, the four height fields and the dissipation field turn out to be all proportional to each other. On a closed boundary, the three height fields are distinct and given by

$h_1^{\text{cl}}\left(x\right)=C_1\hspace{0.17em}\partial \theta \hspace{0.17em}\partial \tilde{\theta },\hspace{1em}{h}_2^{\text{cl}}\left(x\right)=C_2\hspace{0.17em}\partial \theta \hspace{0.17em}\partial \tilde{\theta }+\frac{1}{2\pi }\hspace{0.17em}\theta \hspace{0.17em}\partial \partial \tilde{\theta },\hspace{1em}{h}_3^{\text{cl}}\left(x\right)=C_3\hspace{0.17em}\partial \theta \hspace{0.17em}\partial \tilde{\theta }-\frac{1}{2\pi }\hspace{0.17em}\theta \hspace{0.17em}\partial \partial \tilde{\theta },(94)$

where the $\theta ,\tilde{\theta }$ fields now satisfy the Neumann boundary condition. In both cases, the boundary condition means that the correlators are computed using the Wick theorem with the Wick contractions given by the Green functions G^op or G^cl (see eqs. (64a) and (64b)). Let us point out that for both boundary conditions, the 3-site correlations of three heights 1 do not vanish, unlike their bulk version.

The question of the nature of the height fields on the windy boundary conditions discussed in the previous section is definitely interesting but has not been considered so far.

6.6 Duality

This long section on boundaries has been largely devoted to a discussion of the open and closed boundary conditions, the best known and most studied ones. To finish, it is worth pointing out that a duality exists between these two boundary conditions, which has not been fully investigated nor exploited. This duality follows from a duality relation for planar graphs, well-known in graph theory, and acquires in the framework of the sandpile model an interesting flavor. It has been considered and discussed in [39, 49] in the dimer model, intrinsically related, like the sandpile model, to spanning trees.

Let us consider a rectangular portion of ${\mathrm{\mathbb{Z}}}^2$, that is, the graph Γ made of a rectangular array of vertices, in which two adjacent vertices are linked by a single edge. The boundary conditions chosen for the boundary vertices determine the extended graph ${\text{Γ}}^{*}$, obtained from Γ by adding the sink vertex and the edges connecting the open boundary sites to the sink. The graph ${\text{Γ}}^{*}$ corresponding to a 3 × 3 grid with three open edges and one closed edge is shown in Figure 1.

FIGURE 1

FIGURE 1. Drawing codes for the three figures are as follows. The open circle stands for the sink vertex, while the solid circles stand for the non-sink vertices. The solid lines represent true edges of the extended graphs, unlike the dashed lines connecting the open circles which indicate that these should be identified as the unique sink vertex. Let us note that the corner vertices which lie at the intersection of an open and a closed boundary have a single edge to the sink; those at the intersection of two open boundaries have two such edges. (A) The left figure shows the extended graph of a 3 × 3 grid with open boundary conditions on the left, lower, and right boundaries, and closed boundary condition on the upper boundary. (B) The second panel shows how the dual of the blue graph, in red, is constructed. The red sink is the vertex associated with the outer face of the blue graph. (C) For a better readability, the dual graph alone is reproduced on the third panel. (D) Two dual spanning trees are drawn on the far right.

Once the graph ${\text{Γ}}^{*}$ is embedded in the plane,²⁵ the faces of ${\text{Γ}}^{*}$ are the connected components of its complementary in the plane (for a finite graph, there is thus a large outer face, encircling the graph). The definition of the dual graph ${\left({\text{Γ}}^{*}\right)}^{\text{∗}}$ is standard: the vertices of ${\left({\text{Γ}}^{*}\right)}^{\text{∗}}$ are associated with the faces of ${\text{Γ}}^{*}$, and two such vertices are connected if their corresponding two faces are separated from each other by an edge of ${\text{Γ}}^{*}$. The dual graph of the example above is also shown in Figure 1.

By comparing the two graphs, one immediately notices that the boundary conditions are exchanged: if a boundary is homogeneously open resp. closed in ${\text{Γ}}^{*}$, it becomes homogeneously closed resp. open in ${\left({\text{Γ}}^{*}\right)}^{\text{∗}}$. In addition, the dual graph ${\left({\text{Γ}}^{*}\right)}^{\text{∗}}$ is the extension ${\left({\text{Γ}}^{\text{∗}}\right)}^{*}$ by a sink of a dual rectangular grid ${\text{Γ}}^{\text{∗}}$, of size slightly different from the original grid Γ.

A classical result states that the number of spanning trees on ${\text{Γ}}^{*}$ is equal to the number of spanning trees on its dual ${\left({\text{Γ}}^{*}\right)}^{\text{∗}}$. In fact, for every spanning tree $\mathcal{T}$ on ${\text{Γ}}^{*}$, there is a unique dual spanning tree ${\mathcal{T}}^{\text{∗}}$ on ${\left({\text{Γ}}^{*}\right)}^{\text{∗}}$ such that the two are perfectly interdigitating: the edges of ${\mathcal{T}}^{\text{∗}}$ are exactly those of ${\left({\text{Γ}}^{*}\right)}^{\text{∗}}$ which cross the edges of ${\text{Γ}}^{*}$ not used in $\mathcal{T}$, and vice versa. An example of this is given in Figure 1.

This dual picture implies that the recurrent configurations for the sandpile model defined using ${\text{Γ}}^{*}$ can be isomorphically described by those on ${\left({\text{Γ}}^{\text{∗}}\right)}^{*}$. As far as the counting goes, the equality of their partition functions can be explicitly written for rectangular grids. If Γ is an L₁ × L₂ rectangular grid with k of its four boundaries being open, the other 4−k being closed, the dual ${\text{Γ}}^{\text{∗}}$ is an $L{\text{′}}_1\times L{\text{′}}_2$ rectangular grid with swapped boundary conditions, and the following identity holds,

$Z_{k\hspace{0.17em}\text{op}\hspace{0.17em}|\hspace{0.17em}\left(4-k\right)\hspace{0.17em}\text{cl}}\hspace{0.17em}\left(L_1,L_2\right)=Z_{\left(4-k\right)\hspace{0.17em}\text{op}\hspace{0.17em}|\hspace{0.17em}k\hspace{0.17em}\text{cl}}\hspace{0.17em}\left(L{\text{′}}_1,L{\text{′}}_2\right).(95)$

The dimensions are related as follows: $L_i\text{′}=L_i+1$ resp. L_i − 1 if the opposites sides of length L_i of Γ are both open resp. closed, and $L{\text{′}}_i=L_i$ otherwise. If k = 4, the dual rectangle has all its boundaries closed with a single boundary site open.

The isomorphism of the two descriptions may be hard to formulate in concrete terms for the height variables as it is defined for the associated trees. Its practical utility remains to be seen.

7 More Developments

We would like to add a few more considerations about two further features of the sandpile model, namely, the dissipative sandpile model and some aspects of universality.

7.1 The Massive Sandpile Model

In the standard sandpile we have studied so far, the sites in the bulk of the grid, that is, the vast majority of sites, are conservative. This means that when such a site topples, it loses a certain number of sand grains which are all redistributed to its nearest neighbors. Sand moves in the grid but remains conserved. Dissipative sites must be present for the dynamics of the model to be well-defined; however, the dissipative sites were located most of the time on the boundaries.

The mostly conservative nature of the model is what drives it dynamically to a critical state: when enough sand is stored in the system, large avalanches become likely and span macroscopic parts of it, inducing strong correlations between distant heights. In the long run, the system enters a critical state described by the invariant measure $\mathrm{\mathbb{P}}$, characterized by infinite correlation lengths in the infinite volume limit, and algebraic decays of the correlation functions. The field theory emerging in the scaling limit is conformal and consequently massless.

From the above point of view, a natural way to take the sandpile model out of criticality is to introduce a fair amount of dissipation so as to make the range of the avalanches shorter. It is not completely clear what a fair amount means as there are several ways to introduce dissipation. In the most common version, every site is made dissipative, with a dissipation rate controlled by an external parameter. In this case, it has been argued that indeed criticality is broken, resulting in an exponential decay of the correlations [33, 50, 51]. A mathematically rigorous proof that all correlations decay exponentially has been provided in Ref [52, 53]. Presumably, a nonzero density of dissipative sites could be sufficient to break criticality, but to our knowledge, this possibility has not been investigated. In any case, the field theory emerging from the dissipative sandpile model must be massive, with mass(es) inversely proportional to the lattice correlation length(s).

To make all sites dissipative, one can simply add to the toppling matrix of the standard model an integer multiple of the identity matrix, $\mathrm{\Delta }\to \mathrm{\Delta }\left(t\right)=\mathrm{\Delta }+t\hspace{0.17em}\mathbb{I}$ with t an integer, while leaving all non-diagonal entries unchanged. According to the update of the heights after the toppling of site j, namely, h_i → h_i − Δ_j,i(t), a toppled site loses t sand grains more than what it used to lose (whether or not the toppled site is on a boundary). That this change makes the correlation functions decay exponentially should be clear, for the following simple reason.

The new toppling matrix Δ(t) is a massive Laplacian matrix. It is well-known that the inverse Laplacian Δ⁻¹(t) has a kernel given at large distances by $G_{i_1,i_2}\left(t\right)\simeq \frac{1}{2\pi }{K}_0\left(\left|i_1-i_2\right|\sqrt{t}\right)+\dots$ and decays exponentially like ${\text{e}}^{-r\sqrt{t}}$ at large distances (K₀ is the modified Bessel function). Thus, all multisite probabilities examined in the earlier sections, for observables like minimal cluster variables, arrow variables, isolated dissipation, or boundary heights, will similarly decay exponentially. Though technically less clear for bulk heights equal to 2, 3, or 4, the same decay is expected for the reason explained above: there is a loss of sand each time a site topples, which makes the typical avalanches short-ranged, which in turn induces correlations of heights on local scales only.

To take a concrete example, let us look at the correlation of two heights 1 in the dissipative model. The technique explained in Section 5.3 in terms of defect matrices goes through. At dominant order, the result, which is the off-critical extension of (35), reads [33]

${\sigma }_{\mathrm{1,1}}\left(i_1,i_2;t\right)=-t^2\frac{{\mathrm{\mathbb{P}}}_1^2}{2}\left\{K_0^{″}{}^2-K_0\hspace{0.17em}{K}_0^{″}+\frac{1}{\pi }{K}_0{\text{′}}^2+\frac{1+{\pi }^2}{2{\pi }^2}{K}_0^2\right\}+\dots (96)$

where the argument of the Bessel functions is $r\sqrt{t}$ and ${\mathrm{\mathbb{P}}}_1$ on the right hand side is the critical probability; the dots stand for higher orders in t. We see that the correlation decays exponentially, with a correlation length proportional to ξ ∼ t^−1/2.

How do we compute the scaling limit in the massive model? The general discussion in Section 3 suggested that setting $i=\frac{\overrightarrow{x}}{\epsilon }$ in the lattice correlator and taking the limit over ε (after multiplying the correlator by a suitable power of ε) yields the field theoretic correlator. This cannot be the right way to proceed in the dissipative model. Because the correlators decay exponentially, the limit for ε going to zero of $\text{exp}\left(-\left|{\overrightarrow{x}}_1-{\overrightarrow{x}}_2\right|\sqrt{t}/\epsilon \right)$ vanishes whatever the power of ε it is multiplied by.

The only way to get a nontrivial limit is to take a double limit: as we take the large distance limit by setting $i=\frac{\overrightarrow{x}}{\epsilon }$, we simultaneously take the large correlation length limit by accordingly adjusting the dissipation rate. In the present case, we should take the latter proportional to ε²: we therefore set t = M²ε², with M playing the role of a mass (inverse correlation length in the continuum field theory).

Looking at the lattice correlator (96), we see that the factor t² carries the overall dimension of the fields involved: t² is proportional to M⁴ and thus inversely proportional to a distance to the fourth power. It replaces the explicit dependence in r⁻⁴ in the non-dissipative model. Eventually, we find that the scaling limit of the correlator (96) is

$\underset{\epsilon \to 0}{\text{lim}}\hspace{0.17em}{\epsilon }^{-4}\hspace{0.17em}{\sigma }_{\mathrm{1,1}}\left(\frac{z}{\epsilon },\frac{w}{\epsilon };{\epsilon }^2{M}^2\right)=-M^4\hspace{0.17em}\frac{{\mathrm{\mathbb{P}}}_1^2}{2}\left\{{K_0^{″}}^2-K_0\hspace{0.17em}{K}_0^{″}+\frac{1}{\pi }{K_0^{\prime }}^2+\frac{1+{\pi }^2}{2{\pi }^2}{K}_0^2\right\},(97)$

where the argument of the Bessel function is now Mǀz – wǀ. It is straightforward to check that the M → 0 limit of the previous expression is equal to $-{\mathrm{\mathbb{P}}}_1^2/2{\left|z-w\right|}^4$, obtained in Section 5.

The last question is: the expression above is: the correlator of what field and in what field theory? The most obvious guess turns out to be correct: let us look into the massive extension of the free symplectic fermion theory. It contains the same two fields as before, which simply acquire a mass through a mass term in the action,

$S=\frac{1}{\pi }{{\displaystyle \int }}^{\text{}}\text{d}z\text{d}\overline{z}\hspace{0.17em}\left(\partial \theta \overline{\partial }\tilde{\theta }+\frac{M^2}{4}\theta \tilde{\theta }\right)\hspace{0.17em}.(98)$

The 2-point correlators of the two fundamental fields are now given by

$\langle \theta \left(z,\overline{z}\right)\theta \left(w,\overline{w}\right)\rangle =\langle \tilde{\theta }\left(z,\overline{z}\right)\tilde{\theta }\left(w,\overline{w}\right)\rangle =0,\hspace{1em}\langle \theta \left(z,\overline{z}\right)\tilde{\theta }\left(w,\overline{w}\right)\rangle =K_0\left(M\left|z-w\right|\right).(99)$

Using Wick’s theorem, it is a simple matter to check that the following local field,

$h_1\left(z,\overline{z};M\right)=-{\mathrm{\mathbb{P}}}_1\left[\partial \theta \hspace{0.17em}\overline{\partial }\tilde{\theta }+\overline{\partial }\theta \hspace{0.17em}\partial \tilde{\theta }+\frac{M^2}{2\pi }\hspace{0.17em}\theta \tilde{\theta }\right],(100)$

has a 2-point correlator²⁶ in the massive fermionic theory that is precisely given by equation (97). The 3- and 4-point correlators of the same field have been checked to reproduce the corresponding lattice results. The field $h_1\left(z,\overline{z};M\right)$ is therefore what the height 1 variable in the dissipative sandpile model converges to in the scaling limit.

Similar correlators have been computed for many minimal clusters in Ref [33], with an unexpectedly simple result. The field describing the minimal cluster variable S in the dissipative model appears to be simply given by

$h_S\left(z,\overline{z};M\right)=h_S\left(z,\overline{z}\right)-{\mathrm{\mathbb{P}}}_S\hspace{0.17em}{N}_S\hspace{0.17em}\frac{M^2}{2\pi }\hspace{0.17em}\theta \tilde{\theta },(101)$

where ${\mathrm{\mathbb{P}}}_S$ is the probability of S in the non-dissipative model and N_s is the size of the cluster S. The field $h_S\left(z,\overline{z}\right)$ is still given by relation (54) in terms in the (now massive) fermions.

Likewise, the mixed 2-point correlators for all boundary heights on open and closed boundaries have been explicitly evaluated in the dissipative model [48]. For them too, it is found that the boundary fields given in Section 6.5 get additional terms proportional to $M^2\hspace{0.17em}\theta \tilde{\theta }$.

The nature of the higher height fields remains elusive but is definitely worth investigating as it would add a most valuable and crucial element of understanding of the sandpile model.

7.2 Aspects of Universality

Universality is the statement that the large distance properties of statistical models should only depend on some gross features of the way they are defined; microscopic details which become invisible from large distances should not matter. The statement is admittedly not very precise but, in concrete instances, leads to an expected robustness with respect to local modifications. In sandpile models, these would include the precise way sand is deterministically redistributed among neighbors (provided some form of isotropy is preserved), or, to a certain extent, the specific graph or lattice on which the model is defined. Features that do matter are a substantial introduction of dissipation, as we have seen in the previous section; a directed redistribution of sand after toppling [54]; a dynamics with stochastic toppling rules [55]; the formulation of the model on a hierarchical geometric structure like the Bethe lattice [56]; and, of course, a change of dimensionality of the underlying lattice.

Very early on, universality with respect to the planar lattice on which the sandpile is being formulated has been tested via a renormalization group approach [57, 58] and numerical simulations [59]. More recently, exact calculations of height correlations have been carried out on the honeycomb and triangular lattices.

In Ref [60], all calculations of height 1 correlations presented in the previous sections have been worked out on the hexagonal lattice (in the non-dissipative model). These include the 2-, 3-, and 4-site probabilities for heights 1 in general positions, in the bulk and on open and closed boundaries, as well as 1-site probabilities on the UHP, again for both types of boundary conditions. The results show that although the subdominant contributions differ from those on the square lattice, the dominant terms are exactly identical, up to normalizations. The same distinctive features are found, like the fact that the 3-site bulk correlation vanishes in the scaling limit (the dominant term in the lattice result has dimension −7, instead of −8), and the change of sign for the UHP 1-site probabilities when changing the boundary condition from open to closed (see Section 6.3). Up to normalization, the field identifications of the height 1 variable in the bulk and on open and closed boundaries have been confirmed.

The results have been extended to higher heights on the honeycomb lattice and to all heights on the triangular lattice [61]. Interestingly, these two regular lattices have coordination numbers different from the square lattice, with the consequence that the height variables take in each case a different number of values: four for the square lattice, three for the honeycomb lattice and six for the triangular lattice. This naturally raises the question of which height variables scale to logarithmic fields and which do not.

The calculations have been carried out by using the technique developed in Ref [29], already used on the square lattice. The 1-site probabilities on the infinite honeycomb lattice are all rational,

${\mathrm{\mathbb{P}}}_1=\frac{1}{12},\hspace{1em}\hspace{1em}{\mathrm{\mathbb{P}}}_2=\frac{7}{24},\hspace{1em}\hspace{1em}{\mathrm{\mathbb{P}}}_3=\frac{5}{8},(102)$

while those on the infinite triangular lattice are somewhat more complicated, like

${\mathrm{\mathbb{P}}}_6=\frac{1175}{864}-\frac{365}{144\sqrt{3}\pi }-\frac{289}{12{\pi }^2}+\frac{30\sqrt{3}}{{\pi }^3}+\frac{45}{{\pi }^4}-\frac{54\sqrt{3}}{{\pi }^5}\simeq 0.286,(103)$

and very similar expressions for ${\mathrm{\mathbb{P}}}_{1\mathrm{⩽}a\mathrm{⩽}5}$.

Concerning the nature of the height variables in the scaling limit, the results confirm what the reader has probably already suspected: far from boundaries, the height 1 variable becomes a primary field with conformal weights $\left(h,\overline{h}\right)=\left(\mathrm{1,1}\right)$, while each of the higher heights scales to a logarithmic partner of the height 1, exactly like on the square lattice. On boundaries, all height fields are non-logarithmic. Moreover, all computed correlations²⁷ exhibit the same bulk and boundary behaviors as on the square lattice. Thus, for what concerns the type of the underlying lattice, universality has been explicitly and successfully verified.

8 Conformal Summary

This last section is more specifically oriented toward conformal aspects of the sandpile model. We will summarize what we believe is currently known of the conformal picture, and discuss some of the most peculiar issues that are not so well-understood. We will almost exclusively discuss the non-chiral bulk fields, but before coming to those, we briefly comment on the chiral boundary fields encountered so far.

The boundary fields have been somewhat less investigated than the bulk fields. We have encountered two types of boundary fields, those arising from boundary observables and the boundary condition changing fields. In the first class, we have considered the height fields on open and closed boundaries and the dissipation field. Except for the dissipation on a closed boundary, none of them is logarithmic, and no evidence of a logarithmic partner has been found. All can be expressed as local fields in the symplectic fermions.

In the second class, we found primary fields of weights $-\frac{1}{8}$ and $\frac{3}{8}$, which are both standard fields in a c = −2 CFT. Due to the values of their conformal weight, they cannot be local in the symplectic fermions but are naturally accommodated²⁸ in the symplectic fermion theory [22]. The status of the other boundary condition changing fields related to the windy boundary conditions is uncertain and should be further investigated before their exact nature can be reliably stated.

Thus, overall, the boundary fields raise no particular questions. They are fairly simple fields which fit well within the symplectic theory. From this point of view, the bulk fields are somehow more intriguing.

Most of the bulk fields we have encountered seem to have a realization in terms of symplectic fermions, by which we mean that the fermionic expressions reproduce the known correlators. A few have not been realized in this way so far, namely, the height variables $h_{a\mathrm{⩾}2}$ not equal to 1, logarithmic partners of the height 1 field h₁, and the two fields ρ and $\overline{\rho }$, to which they transform under L₁ and ${\overline{L}}_1$, respectively.

Although we have not given any physical interpretation of ρ and $\overline{\rho }$, they appear to be related to the derivatives of the dissipation field ω [31],

$\rho =\delta \hspace{0.17em}{\overline{L}}_{-1}\omega ,\hspace{1em}\hspace{1em}\overline{\rho }=\delta \hspace{0.17em}{L}_{-1}\omega ,(104)$

where δ is a constant which may depend on the lattice considered and equal to $\delta =\frac{\pi {\mathrm{\mathbb{P}}}_1}{2}$ on the square lattice. In addition, the primary field h₁ may be consistently identified as being proportional to the derivatives of ρ and $\overline{\rho }$,

$L_{-1}\rho ={\overline{L}}_{-1}\overline{\rho }=\beta \lambda h_1,\hspace{1em}\hspace{1em}\beta =\frac{1}{2},(105)$

where λ is defined from $L_0\hspace{0.17em}{h}_2={\overline{L}}_0\hspace{0.17em}{h}_2=h_2+\lambda h_1$ and depends on the normalizations of h₁ and h₂. Combining these relations with the previous ones yields the somewhat surprising result is that the height 1 field is proportional to the Laplacian of the dissipation field, $h_1\sim \partial \overline{\partial }\omega$. The correlator (40) confirms this: applying ${\partial }_1{\overline{\partial }}_1{\partial }_2{\overline{\partial }}_2$ on it indeed yields a multiple of 1/ǀz₁ – z₂ǀ⁴, itself proportional to $\langle h_1\left(z_1,{\overline{z}}_1\right)h_1\left(z_2,{\overline{z}}_2\right)\omega \left(\infty \right)\rangle$.

From these observations, it follows that all bulk fields encountered so far, namely,

$h_{a>1},\hspace{0.17em}{h}_1,\hspace{0.17em}\rho ,\hspace{0.17em}\overline{\rho },\hspace{0.17em}{\rho }_{\to },\hspace{0.17em}{\rho }_{\uparrow },\hspace{0.17em}{\varphi }_S,\hspace{0.17em}{\varphi }_{\leftrightarrow },\hspace{0.17em}{\varphi }_{\updownarrow },\hspace{0.17em}\omega ,\hspace{0.17em}\mathbb{I},(106)$

belong to the same conformal representation as they are all related to each other by the action of Virasoro modes L_n or ${\overline{L}}_n$. Indeed, ρ_→ and ρ_↑ are not quasi-primary and transform to a multiple of $\mathbb{I}$ under L₁ or ${\overline{L}}_1$, while ϕ_S, ϕ_↔ and ϕ_↕ are linear combinations of h₁ and the chiral and antichiral stress–energy tensors T and $\overline{T}$. In fact, in terms of fermions, all these fields, except h_a>1, are proportional to or are linear combinations of $\mathbb{I}$, $\theta \tilde{\theta }$, $\theta \partial \tilde{\theta }$, $\theta \overline{\partial }\tilde{\theta }$, $\partial \theta \overline{\partial }\tilde{\theta }$, $\overline{\partial }\theta \partial \tilde{\theta }$, $\partial \theta \partial \tilde{\theta }$ and $\overline{\partial }\theta \overline{\partial }\tilde{\theta }$. Clearly, the main question is: do the fields h_a>1 also have a realization in terms of symplectic fermions?

In the symplectic fermion theory, the conformal representation which contains the fields quoted above is larger because it contains many other fields, like $\theta \partial \theta$ or $\partial \theta \overline{\partial }\theta$, which have not yet been found in the sandpile model. Among other peculiarities, the fermionic theory also contains four logarithmic pairs (ϕ^αβ, ψ^αβ) of weight (1,1), given by ${\varphi }^{\alpha \beta }=\partial {\theta }^{\alpha }\overline{\partial }{\theta }^{\beta }$ for the primary fields, and ${\psi }^{\alpha \beta }=\theta \tilde{\theta }\hspace{0.17em}\partial {\theta }^{\alpha }\overline{\partial }{\theta }^{\beta }$ for their logarithmic partners, where θ^α and θ^β are independently either θ or $\tilde{\theta }$ (see [22] for more details).

Conformal representations of the above type are called staggered modules and have been first studied in Ref [62] in their chiral version. As far as we know, it has been first noticed in Ref [63] for the case of modules containing rank 2 Jordan blocks that these representations are characterized by an intrinsic complex parameter β, known as a logarithmic coupling, an indecomposability parameter or a beta-invariant. The parameter β is crucial because it specifies the equivalence class of such representations, whose general structure was further studied in Ref [20] in the rank 2 case. The non-chiral staggered modules are far less understood and documented and reflect the difficulty to formulate a consistent and local logarithmic CFT (see however [64, 65]). It is nonetheless believed that the parameter β present in the chiral representations plays the same role of equivalence class label in the non-chiral ones, even if the latter may have more than one such label.

Concentrating on the action of the chiral Virasoro modes, the parameter β arises when we consider the triangular relations satisfied by a generic logarithmic pair (ϕ, ψ) of weights (1,1) and the associated ρ.

The arrows coming out of ψ indicate the actions (L₀ − 1)ψ = λϕ and L₁ψ = ρ. It is important to note that if the normalization of ψ is fixed, those of λϕ and ρ are fixed as well (the value of λ depends on the way ϕ is normalized). The vertical arrow indicates that L₋₁ρ is proportional to λϕ,

$L_{-1}\rho =\beta \left(\lambda \varphi \right),(107)$

the proportionality factor β being intrinsic to the representation as all normalizations have already been fixed. In addition, these relations are invariant under the change ψ → ψ + αϕ because the field ϕ is primary (L₁ϕ = 0), and so they do not depend on which logarithmic partner is considered.

To answer the above question thus amounts to check whether the sandpile representation and the symplectic representation have the same value of β. The value of β in the sandpile model has been given above: if the pair (ϕ, ψ) is chosen to be (h₁, h₂), with the same normalization as the height variables on the square lattice, in which case $\lambda =-\frac{1}{2}$, then one finds $\beta =\frac{1}{2}$ [27].

The field h₁ has been already identified in terms of fermions and yields a natural choice for the primary field ϕ on the symplectic side,

${\varphi }_{\theta }=-{\mathrm{\mathbb{P}}}_1\left(\partial \theta \hspace{0.17em}\overline{\partial }\tilde{\theta }+\overline{\partial }\theta \hspace{0.17em}\partial \tilde{\theta }\right).(108)$

As mentioned earlier, the lattice results in the scaling limit are consistent with h₁ being degenerate at level 2, namely, (L²₋₁ – 2L₋₂)h₁ = 0 (see Section 5.5). The same equation is satisfied by ϕ_θ.

The only candidate for the logarithmic partner of ϕ_θ is proportional to $\theta \tilde{\theta }\hspace{0.17em}\left(\partial \theta \hspace{0.17em}\overline{\partial }\tilde{\theta }+\overline{\partial }\theta \hspace{0.17em}\partial \tilde{\theta }\right)$ up to an irrelevant multiple of ϕ_θ. By computing its conformal transformations via its OPE with the chiral stress–energy tensor, one finds that the following normalization,

${\psi }_{\theta }=-{\mathrm{\mathbb{P}}}_1\hspace{0.17em}\theta \tilde{\theta }\hspace{0.17em}\left(\partial \theta \hspace{0.17em}\overline{\partial }\tilde{\theta }+\overline{\partial }\theta \hspace{0.17em}\partial \tilde{\theta }\right),(109)$

satisfies $L_0{\psi }_{\theta }={\psi }_{\theta }-\frac{1}{2}{\varphi }_{\theta }$, for the same value $\lambda =-\frac{1}{2}$. The same OPE reveals in addition that

${\rho }_{\theta }=L_1{\psi }_{\theta }=-\frac{{\mathrm{\mathbb{P}}}_1}{2}\hspace{0.17em}\left(\theta \hspace{0.17em}\overline{\partial }\tilde{\theta }+\overline{\partial }\theta \hspace{0.17em}\tilde{\theta }\right),(110)$

from which, upon using $\partial \overline{\partial }\theta =\partial \overline{\partial }\tilde{\theta }=0$, one obtains

$L_{-1}{\rho }_{\theta }=\partial {\rho }_{\theta }=-\frac{{\mathrm{\mathbb{P}}}_1}{2}\left(\partial \theta \hspace{0.17em}\overline{\partial }\tilde{\theta }+\overline{\partial }\theta \hspace{0.17em}\partial \tilde{\theta }\right)=\frac{1}{2}{\varphi }_{\theta }.(111)$

Comparing with equation (107), the value of the logarithmic coupling is found to be β_θ = −1 in the fermionic realization. As a consequence, the symplectic fermion theory cannot accommodate the height fields h_a>1 and therefore does not appear to be the correct CFT to describe the scaling limit of the sandpile model.

As one might suspect, the value of β has strong consequences on correlation functions involving ψ. A detailed comparison between $\beta =\frac{1}{2}$ and the fermionic realization ${\beta }_{\theta }=-1$ has been made in Ref [27]; it was shown in particular that the correlations with a trial field h₂ corresponding to a value $\beta =-1$ do not match the lattice results.²⁹ On general grounds, this can also be understood from the fact that the value of β determines the singular descendant of ψ, which, if set to zero, yields a β-dependent differential equation satisfied by any correlator containing ψ. In the present case, the singular logarithmic field is a combination of a descendant of ψ at level 5 and a descendant of ρ at level 6, with the following explicit dependence on β [20]:

$\begin{array}{c}\xi =\left(L_{-1}^3-8L_{-2}{L}_{-1}+12L_{-3}\right)\left(L_{-1}^2-2L_{-2}\right)\psi -\frac{1}{\beta }[-\frac{16}{3}\left(\beta +1\right)L_{-2}^2{L}_{-1}^2+\frac{4}{3}\left(14\beta +5\right)L_{-3}{L}_{-2}{L}_{-1}\\ -6\beta L_{-3}^2-6\left(\beta -2\right)L_{-4}{L}_{-1}^2+8\beta L_{-4}{L}_{-2}-\frac{2}{3}\left(5\beta +2\right)L_{-5}{L}_{-1}+4\beta L_{-6}]\rho .\end{array}(112)$

Using the relations $L_1\psi =\rho ,\hspace{0.17em}\left(L_0-1\right)\psi =\lambda \varphi$, as well as the degeneracy condition $\left(L_{-1}^2-2L_{-2}\right)\varphi =0$ (and the value c = −2), one can verify that the field ξ satisfies L₁ξ = L₂ξ = 0, provided the identity L₋₁ρ = βλϕ holds. A rather convincing confirmation for the value of $\beta =\frac{1}{2}$ in the sandpile model is therefore to check that the various correlators involving h₂ indeed satisfy the condition ξ = 0 for $\beta =\frac{1}{2}$. It has been done for the correlator (87b).

The situation therefore seems to be the following. The sandpile model contains a conformal logarithmic representation whose structure is very similar to the one appearing in the symplectic fermion theory, but which is nevertheless inequivalent to it. As far as the logarithmic partner ψ is not brought in, the two representations look the same; this explains why some of the fields can be realized in terms of symplectic fermions. However, the fermionic theory does not contain the $\beta =\frac{1}{2}$ representation found in the sandpile model, from which one concludes that it does not describe its scaling limit.

To characterize the CFT that does describe the sandpile model, even if a Lagrangian realization of it cannot be found, remains an enormous challenge. At the moment, this looks to be an extremely ambitious question in view of the (very) small number of fields which has been successfully identified.

Author Contributions

The author confirms being the sole contributor of this work and has approved it for publication.

Funding

The author is a senior research associate of FRS-FNRS (Belgian Fund for Scientific Research). This work was supported by the Fonds de la Recherche Scientifique (FNRS) and the Fonds Wetenschappelijk Onderzoek-Vlaanderen (FWO) under EOS (project no 30889451).

Conflict of Interest

The author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Footnotes

¹There is no need to keep track of the number of sand grains in the sink, and so we do not assign it a height variable.

²The result holds in the more general case where p_i≠0 for every i.

³A notable exception concerns the linear or almost linear graphs, for which the recurrence property usually takes a simpler form and allows for a larger number of explicit results (see f.i., [9, 10]).

⁴We have mentioned that the operators a_i generate an Abelian algebra. But when acting on recurrent configurations, they are invertible and therefore generate an Abelian group, called the sandpile group. The sandpile group, of order equal to detΔ, has been determined for a number of finite graphs.

⁵Recurrent configurations form an Abelian group under the site-wise addition of the heights, followed by relaxation. This group is isomorphic to the sandpile group. In particular, one of the recurrent configurations is the identity in the group and shows remarkable geometric patterns [11–14].

⁶For the scaling to be nontrivial, some external parameters may need to be appropriately scaled with ε. One example of this is discussed in Section 7.1.

⁷It is bounded if all the linear sizes of the finite systems in the sequence defining the infinite volume limit grow exactly like 1/ε.

⁸Among the many books and reviews on the renormalization group in statistical mechanics, see, for instance, the book by Cardy [15].

⁹For instance, the energy density in the Ising model, namely, the product of two neighboring spins, gives rise to a field that is different from the one obtained from the spin variable itself. Later, we will give examples of this in the sandpile models (cluster variables).

¹⁰One should add “Euclidean” field theory because it is formulated on a Euclidean space ${\mathrm{ℝ}}^d$.

¹¹The material recalled in this section is completely standard; useful references include Ref [16] (rather comprehensive) and [17] (more focused on critical statistical systems).

¹²A general 3-point correlator $\langle {\mathrm{\Phi }}_1\left(z_1\right)\hspace{0.17em}{\mathrm{\Phi }}_2\left(z_2\right)\hspace{0.17em}{\mathrm{\Phi }}_3\left(z_3\right)\rangle$ is a function of three complex numbers; if all three fields are quasi-primary, that function can be determined by trading z₁, z₂, and z₃ for the three complex parameters of a general Möbius transformation.

¹³For rectangular grids $\text{Γ}\subset {\mathrm{\mathbb{Z}}}^2$, the notion of boundary is clear: when Γ is embedded in ${\mathrm{\mathbb{Z}}}^2$, the boundary sites are those which are connected to sites of ${\mathrm{\mathbb{Z}}}^2$ not in Γ.

¹⁴The reader will legitimately point out that the Laplacian on ${\mathrm{\mathbb{Z}}}^2$ has a zero mode and is therefore not invertible. A closer look at the determinants (27) however reveals that they only depend on differences of the inverse matrix entries, which are perfectly well-defined.

¹⁵In order to make contact with fields, we slightly change the notation h_i → h(i) for the height at site i.

¹⁶We note that the inverse ratio ${\text{det}\Delta }_{\text{Γ}}/{\text{det}\tilde{\Delta }}_1$ is equal to $\text{Prob}\left[\text{all}\hspace{0.17em}{h}_j\mathrm{⩽}4\right]=\text{Prob}\left[h_i\mathrm{⩽}4\right]=1-\text{Prob}\left[h_i=5\right]$ where the probabilities are evaluated in the modified model. The divergence mentioned in the text therefore implies that h_i = 5 with probability 1.

¹⁷The four height fields h_a(z) also satisfy the trivial identity h₁(z) + h₂(z) + h₃(z) + h₄(z) = 0.

¹⁸The qualifier “conservative” means that the defect matrix that implements the bond modifications has zero row and column sums. The defect matrix used in Section 5.1 to compute the height 1 probability does not have this property. It can however be replaced by another one that does have it [33]. An example of a nonconservative bond modification is given in Section 5.7.

¹⁹Indeed, the burning algorithm implies that the boundary sites, all with a height equal to 4, will burn at the first step of the burning process. The sites on the next layer all have $z_i^{*}=4$, corresponding to the open condition.

²⁰Of course, we talk here of no other universality class of boundary conditions. Many boundary conditions may differ in the way they are microscopically defined on the lattice and nevertheless renormalize to the same continuum boundary condition in the scaling limit.

²¹In the correspondence between statistical system and field theory, the boundary condition changing fields are somehow special. They describe the effects of a change of the boundary condition but are not associated with a lattice observable, unlike the height fields h_a, for instance.

²²Again, the technical assumption is that the field ϕ^op,cl is degenerate at level 2, similarly to the height 1 field h₁ (see Section 5.5).

²³The CFT is really defined on the UHP plus the point at infinity. The boundary must therefore be thought of as the real line plus the two points ±∞ identified, and forming a loop closing at infinity. Any change of the boundary condition thus involves an even number of insertions of ϕ^op,cl. For instance, $\langle {\varphi }^{\text{op},\text{cl}}\left(0\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(\infty \right)\rangle$ changes the boundary condition from open to closed on the positive real axis.

²⁴In contrast, the outgoing arrow of a closed boundary site of the UHP can point left, up, or right, while that of an open boundary site can point in any of the four directions, a down arrow pointing to the root.

²⁵This requires ${\text{Γ}}^{*}$ to be planar and therefore excludes that some of the bulk vertices and some of the boundary vertices be open (dissipative) at the same time, except in a few very special cases.

²⁶The insertion by hand of the dissipation field ω(∞) at infinity in the field theoretic correlator is not required in the dissipative model as dissipation is present everywhere in the bulk.

²⁷Some of the calculations done on the square lattice could not be worked out. For instance, we could not find a proper method of images to compute the Green matrix on the triangular half-plane with the closed boundary condition and therefore could not investigate that boundary condition.

²⁸Very much like the spin field of the Ising model belongs naturally to the free Majorana fermionc theory with $c=\frac{1}{2}$, despite being nonlocal in the fermions.

²⁹As an example, the correlator ${\langle h_2\left(z,\overline{z}\right)\rangle }_{\text{mix}}$ displayed in (87b), and corresponding to the four-point function $\langle {\varphi }^{\text{op},\text{cl}}\left(x_1\right)\hspace{0.17em}{\varphi }^{\text{cl},\text{op}}\left(x_2\right)\hspace{0.17em}{h}_2\left(z,\overline{z}\right)\rangle$, can be computed upon assuming that h₂ is a logarithmic partner of h₁ carrying a generic value of β ≠ 0. Its general form is given in Ref [7].

References

1. Bak P, Tang C, Wiesenfeld K. Self-organized criticality: an explanation of the 1/fnoise. Phys Rev Lett (1987) 59:381. doi:10.1103/physrevlett.59.381

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Bak P. How nature works: the science of self-organised criticality. New York: Springer-Verlag (1996).

3. Jensen HJ. Self-organized criticality. Cambridge: Cambridge University Press (1998).

4. Pruessner G. Self-organised criticality: theory, models and characterisation. Cambridge: Cambridge University Press (2012).

CrossRef Full Text

5. Dhar D. Theoretical studies of self-organized criticality. Physica A: Stat Mech its Appl (2006) 369:29. doi:10.1016/j.physa.2006.04.004

CrossRef Full Text | Google Scholar

6. Dhar D. Self-organized critical state of sandpile automaton models. Phys Rev Lett (1990) 64:1613. doi:10.1103/physrevlett.64.1613

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Ruelle P. Logarithmic conformal invariance in the Abelian sandpile model. J Phys A: Math Theor (2013) 46:494014. doi:10.1088/1751-8113/46/49/494014

CrossRef Full Text | Google Scholar

8. Majumdar SN, Dhar D. Equivalence between the Abelian sandpile model and the q→0 limit of the Potts model. Physica A: Stat Mech its Appl (1992) 185:129. doi:10.1016/0378-4371(92)90447-x

CrossRef Full Text | Google Scholar

9. Ruelle P, Sen S. Toppling distributions in one-dimensional Abelian sandpiles. J Phys A: Math Gen (1992) 25:L1257. doi:10.1088/0305-4470/25/22/006

CrossRef Full Text | Google Scholar

10. Ali AA, Dhar D. Breakdown of simple scaling in Abelian sandpile models in one dimension. Phys Rev E (1995) 51:R2705. doi:10.1103/physreve.51.r2705

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Creutz M. Abelian sandpiles. Comput Phys (1991) 5:198. doi:10.1063/1.168408

CrossRef Full Text | Google Scholar

12. Dhar D, Ruelle P, Sen S, Verma D-N. Algebraic aspects of Abelian sandpile models. J Phys A: Math Gen (1995) 28:805. doi:10.1088/0305-4470/28/4/009

CrossRef Full Text | Google Scholar

13. Le Borgne Y, Rossin D. On the identity of the sandpile group. Discrete Math (2002) 256:775. doi:10.1016/s0012-365x(02)00347-3

CrossRef Full Text | Google Scholar

14. Caracciolo S, Paoletti G, Sportiello A. Explicit characterization of the identity configuration in an Abelian sandpile Model. J Phys A: Math Theor (2008) 41:495003. doi:10.1088/1751-8113/41/49/495003

CrossRef Full Text | Google Scholar

15. Cardy J. Scaling and Renormalization in statistical Physics, cambridge lecture notes in Physics. Cambridge: Cambridge University Press (1996).

16. Di Francesco P, Mathieu P, Sénéchal D. Conformal field theory. New York: Springer (1997).

17. Henkel M. Conformal invariance and critical phenomena. New York: Springer (1999).

18. Gainutdinov A, Ridout D, Runkel I. Logarithmic conformal field theory. J Phys A: Math Theor (2013) 46:490301.

CrossRef Full Text | Google Scholar

19. Flohr MAI. Bits and pieces in logarithmic conformal field theory. Int J Mod Phys A (2003) 18:4497. doi:10.1142/s0217751x03016859

CrossRef Full Text | Google Scholar

20. Kytölä K, Ridout D. On staggered indecomposable Virasoro modules. J Math Phys (2009) 50:123503. doi:10.1063/1.3191682

CrossRef Full Text | Google Scholar

21. Gurarie V. Logarithmic operators in conformal field theory. Nucl Phys B (1993) 410:535. doi:10.1016/0550-3213(93)90528-w

CrossRef Full Text | Google Scholar

22. Gaberdiel MR, Kausch HG. A local logarithmic conformal field theory. Nucl Phys B (1999) 538:631. doi:10.1016/s0550-3213(98)00701-9

CrossRef Full Text | Google Scholar

23. Gaberdiel MR, Runkel I. The logarithmic triplet theory with boundary. J Phys A: Math Gen (2006) 39:14745. doi:10.1088/0305-4470/39/47/016

CrossRef Full Text | Google Scholar

24. Pearce PA, Rasmussen J, Zuber J-B. Logarithmic minimal models. J Stat Mech (2006) 2006:P11017. doi:10.1088/1742-5468/2006/11/p11017

CrossRef Full Text | Google Scholar

25. Majumdar SN, Dhar D. Height correlations in the Abelian sandpile model. J Phys A: Math Gen (1991) 24:L357. doi:10.1088/0305-4470/24/7/008

CrossRef Full Text | Google Scholar

26. Priezzhev VB. Structure of two-dimensional sandpile. I. Height probabilities. J Stat Phys (1994) 74:955. doi:10.1007/bf02188212

CrossRef Full Text | Google Scholar

27. Jeng M, Piroux G, Ruelle P. Height variables in the Abelian sandpile model: scaling fields and correlations. J Stat Mech (2006) 2006:P10015. doi:10.1088/1742-5468/2006/10/p10015

CrossRef Full Text | Google Scholar

28. Poghosyan VS, Priezzhev VB, Ruelle P. Return probability for the loop-erased random walk and mean height in the Abelian sandpile model: a proof. J Stat Mech (2011) 2011:P10004. doi:10.1088/1742-5468/2011/10/p10004

CrossRef Full Text | Google Scholar

29. Kenyon RW, Wilson DB. Spanning trees of graphs on surfaces and the intensity of loop-erased random walk on planar graphs. J Amer Math Soc (2015) 28:985. doi:10.1090/S0894-0347-2014-00819-5

CrossRef Full Text | Google Scholar

30. Caracciolo S, Sportiello A. Exact integration of height probabilities in the Abelian sandpile model. J Stat Mech (2012) 2012:P09013. doi:10.1088/1742-5468/2012/09/p09013

CrossRef Full Text | Google Scholar

31. Poncelet A, Ruelle P. Multipoint correlators in the Abelian sandpile model. J Stat Mech (2017) 2017:123102. doi:10.1088/1742-5468/aa9a59

CrossRef Full Text | Google Scholar

32. Kassel A, Wilson DB. The looping rate and sandpile density of planar graphs. Amer Math Monthly (2016) 123:19. doi:10.4169/amer.math.monthly.123.1.19

CrossRef Full Text | Google Scholar

33. Mahieu S, Ruelle P. Scaling fields in the two-dimensional Abelian sandpile model. Phys Rev E Stat Nonlin Soft Matter Phys (2001) 64:66130. doi:10.1103/PhysRevE.64.066130

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Poghosyan VS, Grigorev SY, Priezzhev VB, Ruelle P. Pair correlations in sandpile model: a check of logarithmic conformal field theory. Phys Lett B (2008) 659:768. doi:10.1016/j.physletb.2007.12.002

CrossRef Full Text | Google Scholar

35. Poghosyan VS, Grigorev SY, Priezzhev VB, Ruelle P. Logarithmic two-point correlators in the Abelian sandpile model. J Stat Mech (2010) 2010:P07025. doi:10.1088/1742-5468/2010/07/p07025

CrossRef Full Text | Google Scholar

36. Piroux G, Ruelle P. Pre-logarithmic and logarithmic fields in a sandpile model. J Stat Mech Theor Exp (2004) 2004:P10005. doi:10.1088/1742-5468/2004/10/p10005

CrossRef Full Text | Google Scholar

37. Jeng M. Conformal field theory correlations in the Abelian sandpile model. Phys Rev E Stat Nonlin Soft Matter Phys (2005) 71:16140. doi:10.1103/PhysRevE.71.016140

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Poghosyan VS, Priezzhev VB. The problem of predecessors on spanning trees. Act Polytech (2011) 51:59. doi:10.14311/1364

CrossRef Full Text | Google Scholar

39. Izmailian NS, Priezzehv VB, Ruelle P. Non-local finite-size effects in the dimer model. SIGMA (2007) 3:1. doi:10.3842/SIGMA.2007.001

CrossRef Full Text | Google Scholar

40. Brankov JG, Ivashkevich EV, Priezzhev VB. Boundary effects in a two-dimensional Abelian sandpile. J Phys France (1993) 3:1729. doi:10.1051/jp1:1993212

CrossRef Full Text | Google Scholar

41. Piroux G, Ruelle P. Logarithmic scaling for height variables in the Abelian sandpile model. Phys Lett B (2005) 607:188. doi:10.1016/j.physletb.2004.12.045

CrossRef Full Text | Google Scholar

42. Cardy JL. Conformal invariance and surface critical behavior. Nucl Phys B (1984) 240:514. doi:10.1016/0550-3213(84)90241-4

CrossRef Full Text | Google Scholar

43. Deift P, Its A, Krasovsky I. Toeplitz matrices and Toeplitz determinants under the impetus of the ising model: some history and some recent results. Comm Pure Appl Math (2013) 66:1360. doi:10.1002/cpa.21467

CrossRef Full Text | Google Scholar

44. Ruelle P. A c=−2 boundary changing operator for the Abelian sandpile model. Phys Lett B (2002) 539:172. doi:10.1016/s0370-2693(02)02069-5

CrossRef Full Text | Google Scholar

45. Ruelle P. Wind on the boundary for the Abelian sandpile model. J Stat Mech (2007) 2007:P09013. doi:10.1088/1742-5468/2007/09/p09013

CrossRef Full Text | Google Scholar

46. Jeng M. Four height variables, boundary correlations, and dissipative defects in the Abelian sandpile model. Phys Rev E Stat Nonlin Soft Matter Phys (2005) 71:036153. doi:10.1103/PhysRevE.71.036153

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Ivashkevich EV. Boundary height correlations in a two-dimensional Abelian sandpile. J Phys A: Math Gen (1994) 27:3643. doi:10.1088/0305-4470/27/11/014

CrossRef Full Text | Google Scholar

48. Piroux G, Ruelle P. Boundary height fields in the Abelian sandpile model. J Phys A: Math Gen (2005) 38:1451. doi:10.1088/0305-4470/38/7/004

CrossRef Full Text | Google Scholar

49. Izmailian NS, Priezzehv VB, Ruelle P, Hu C-K. Logarithmic conformal field theory and boundary effects in the dimer model. Phys Rev Lett (2005) 95:260602. doi:10.1103/physrevlett.95.260602

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Ghaffari P, Lise S, Jensen HJ. Nonconservative sandpile models. Phys Rev E (1997) 56:6702. doi:10.1103/physreve.56.6702

CrossRef Full Text | Google Scholar

51. Tsuchiya T, Katori M. Proof of breaking of self-organized criticality in a nonconservative Abelian sandpile model. Phys Rev E (2000) 61:1183. doi:10.1103/physreve.61.1183

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Maes C, Redig F, Saada E. The infinite volume limit of dissipative Abelian sandpiles. Commun Math Phys (2004) 244:395. doi:10.1007/s00220-003-1000-8

CrossRef Full Text | Google Scholar

53. Járai A, Redig F. Approaching Criticality via the Zero Dissipation Limit in the Abelian Avalanche Model. J Stat Phys (2015) 159:1369–1407.

54. Dhar D, Ramaswamy R. Exactly solved model of self-organized critical phenomena. Phys Rev Lett (1989) 63:1659. doi:10.1103/physrevlett.63.1659

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Manna SS. Two-state model of self-organized criticality. J Phys A: Math Gen (1991) 24:L363. doi:10.1088/0305-4470/24/7/009

CrossRef Full Text | Google Scholar

56. Dhar D, Majumdar SN. Abelian sandpile model on the Bethe lattice. J Phys A: Math Gen (1990) 23:4333. doi:10.1088/0305-4470/23/19/018

CrossRef Full Text | Google Scholar

57. Papoyan VV, Povolotsky AM. Renormalization group study of sandpile on the triangular lattice. Physica A: Stat Mech its Appl (1997) 246:241. doi:10.1016/s0378-4371(97)00347-6

CrossRef Full Text | Google Scholar

58. Lin CY, Hu CK. Renormalization-group approach to an Abelian sandpile model on planar lattices. Phys Rev E Stat Nonlin Soft Matter Phys (2002) 66:21307. doi:10.1103/PhysRevE.66.021307

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Hu C-K, Lin C-Y. Universality in critical exponents for toppling waves of the BTW sandpile model on two-dimensional lattices. Physica A: Stat Mech its Appl (2003) 318:92. doi:10.1016/s0378-4371(02)01411-5

CrossRef Full Text | Google Scholar

60. Azimi-Tafreshi N, Dashti-Naserabadi H, Moghimi-Araghi S, Ruelle P. The Abelian sandpile model on the honeycomb lattice. J Stat Mech (2010) 2010:P02004. doi:10.1088/1742-5468/2010/02/p02004

CrossRef Full Text | Google Scholar

61. Poncelet A, Ruelle P. Sandpile probabilities on triangular and hexagonal lattices. J Phys A: Math Theor (2018) 51:15002. doi:10.1088/1751-8121/aa9255

CrossRef Full Text | Google Scholar

62. Rohsiepe F. On reducible but indecomposable representations of the Virasoro algebra. arXiv (1996).

Google Scholar

63. Gaberdie MR, Kausch HG. Indecomposable fusion products. Nucl Phys B (1996) 477:293. doi:10.1016/0550-3213(96)00364-1

CrossRef Full Text | Google Scholar

64. Do A-L, Flohr M. Towards the construction of local logarithmic conformal field theories. Nucl Phys B (2008) 802:475. doi:10.1016/j.nuclphysb.2008.05.001

CrossRef Full Text | Google Scholar

65. Ridout D. Non-chiral logarithmic couplings for the Virasoro algebra. J Phys A: Math Theor (2012) 45:255203. doi:10.1088/1751-8113/45/25/255203

CrossRef Full Text | Google Scholar