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.


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 end-point 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, dropplet 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 80's 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 [BTK87]. 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 others models showing SOC have been proposed, as abundantly illustrated in this volume, and in introductory books and reviews [Ba96,Je98,Pr12,Dh06].
The present review will be exclusively concerned with specific versions of two-dimensional sandpile models, formulated by Dhar [Dh90] and 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 (namely 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 neighbourhood 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 themsleves 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 80's 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, refering 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 to establish 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 [Ru13]. 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.

Abelian sandpile models
The models we will 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 1. A height configuration C is a set of heights values {h i } i∈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 some vertices in a non-empty subset D ⊂ V . Vertices in D are called dissipative or open, while those in V \ D are conservative or closed. If Γ ⋆ = (V ⋆ , E ⋆ ) 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 neighbours in Γ), and similarly z ⋆ i its coordination number in Γ ⋆ . Thus Finally we say that a site i of V is stable 1 if its height satisfies 1 h i z ⋆ i . A height configuration is stable if all sites are stable. Clearly the number of stable configurations is equal to i∈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.
(i) 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 new i = h i + δ i,j . If h j z ⋆ j , then C new is stable and defines C t+1 ; if not, we proceed to step (ii).
(ii) Relaxation: if h new j > 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 neighbours in Γ receives one grain, and the remaining z ⋆ j − z j grains go to the sink. After this, one or more neighbours 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, for i, j ∈ 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 ∈ 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 non-zero 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 motivated 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 ∈ 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 [Dh90], 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 i C t . An important question concerns the invariant measures, since they control the behaviour of the model in the long run.
If there is no strong reason to favour certains sites, one takes all probabilities p i equal (uniform distribution). In this case 2 , Dhar [Dh90] has shown that there is unique invariant measure P Γ , which is uniform on its support. In the Markov chain terminology, the configurations in the support of P Γ 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 non-trivial, 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 ⊂ Γ (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 neighbours in F . It can be shown [MD92] 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 neighbouring 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 neighbouring 1's) but also globally because asserting that a given configuration is recurrent generally requires to scan 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 of any height in a recurrent configuration preserves the recurrence.
The burning algorithm [MD92] (see also the review [Dh06]) provides a convenient way to test whether a given stable configuration is recurrent. In addition to provide 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 Γ ⋆ , rooted at the sink site s. Let us recall that a spanning 2 The result holds in the more general case where p i = 0 for every i. tree is a loopless connected subgraph (V ⋆ , F ) ⊂ Γ ⋆ = (V ⋆ , E ⋆ ) with F ⊂ 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 [Dh90], naturally identified as the partition function Z since the invariant measure is uniform, for ∆ the toppling matrix introduced in (2.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 i ∆ i,i ). On a large grid in 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 det ∆ ≃ e 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 behaviour 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 behaviour 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 P Γ . The use of P Γ is what makes most of the calculations fairly hard 3 because as noted earlier, that measure is non-local in terms of the (local) height variables (equivalently the recurrence criterion is non-local).
We should remark that the measure P Γ 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 4 , the geometric structure of the identity configuration 5 , or 3 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 for instance [RS92,AD95]). 4 We have mentionned 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. 5 Recurrent configurations form an Abelian group under the sitewise 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 [Cr91,DRSV95,LBR02,CPS08].
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 behaviours 'in the large', a point of view made technically more precise by the renormalization group analysis.
In order to identify these common behaviours, 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 behaviours 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 behaviours 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.

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 behaviour 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 εZ d (or another shrinked regular lattice). For smaller and smaller ε, fixing a macroscopic distance x = ε m ∈ εZ d amounts to probe larger and larger scales m in terms of lattice units, and at the same time, allows to keep a macroscopic distance r = | x| under control. The scaling limit corresponds 6 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 R d or a part of R d , which may be bounded 7 . 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 originial 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 8 .
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 behaviours 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 behaviour 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 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 second generation, each of which has integrated the fluctuations of Λ 2d microscopic variables over scales smaller than Λ 2 , 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 behaviour 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 9 , 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 (non-independent) random microscopic variables, is usually given by a Gibbs measure, and written, up to normalization, as P where H is the Hamiltonian of the system, i.e. 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 0 ≡ H. The first renormalization yields the coarse-grained variables of the first generation and a corresponding Hamiltonian H 1 , computed (at least in principle) by summing exp (−H 0 ) over the microscopic variables inside the boxes. Similarly the k-th iteration will produce a Hamiltonian H k defining the statistical model for the coarse-grained variables of the k-th 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 1 i , h 2 i , . . .) and a Hamiltonian H(h 1 i , h 2 i , . . .), it is given in the scaling limit by a set of continuous fields (φ 1 ( x), φ 2 ( x), . . . ) and an action S[φ 1 , φ 2 , . . .]. The pair {(φ 1 ( x), φ 2 ( x), . . . ), S} is refered to as a continuum field theory 10 . 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 0 , H 1 , . . . is a formidable task that is, for all practical purposes, impossible to carry out explicitely, 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 9 F.i. the energy density in the Ising model, namely the product of two neighbouring 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).
10 One should add 'Euclidean' field theory because it is formulated on a Euclidean space R d .
trial field theory, perhaps suggested by an educated guess, can be firmly tested by comparing correlations functions. If the lattice microscopic variable h i= x ε (at site i) converge in the scaling limit to the field φ(x), it must be true that the scaling limit of the lattice correlators are equal to field theoretic correlators, namely where the exponent ∆ is determined so that the limits on the l.h.s. 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 site i inserted in the lattice correlator, the corresponding field Φ(x) to which it converges is inserted in the field theoretic correlator, So we can write the convergence of a lattice observable to a field as the formal identity, 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 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 to determine 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 behaviour 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 non-local 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 non-trivial geometries (cylinder, torus, ...), • 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).

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 site i takes the following form, 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, i.e. the coordinate transformations which preserve angles 11 . 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 b (d additional parameters) and take the following general form 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 R d ∪ {∞} 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, 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, (a constant for ∆ = 0 by translation invariance) and the 2-point correlator of two quasi-primary fields, (4.5) so that ∆ can be identified with the dimension of the field Φ (in units of inverse length). Global conformal invariance also completely determines the correlator Φ 1 ( x 1 ) Φ 2 ( x 2 ) Φ 3 ( x 3 ) of three (and not more) quasi-primary fields, We see that the lattice 2-correlator (4.1) is consistent with the convergence of the observable O i to a quasi-primary field Φ( x) of dimension ∆ upon setting i = x/ε since the matching identity (3.1) is satisfied, (4.7) We note that all subdominant terms in the lattice correlator (4.1) drop out when taking the limit ε → 0, confirming once more that a field theory captures the large distance behaviour of a critical lattice model. What has been just recalled is valid in any dimension d 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(2, C) group of Möbius transformations w = az+b cz+d , for a, b, c, d ∈ C 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 parameter to fix it (f.i. the coefficients of its Laurent expansion in some neighbourhood), so that the conformal 'group' is certainly infinitedimensional. 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 else its image is not the whole complex plane.
For instance the map w = L 2πi log z maps the complex plane to a cylinder of circumference L. The discussion of two-dimensional conformal group is thus usually carried out at the level of its algebra, for which infinitesimal transformations of the form w = z + ǫ z n+1 are considered. The corresponding generators satisfy the famous infinite-dimensional Virasoro algebra, 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 0 , L ±1 , whose algebra is unaffected by the central charge, are the infinitesimal generators of the Möbius group, with L −1 and L 0 corresponding to translations and dilations respectively. As it turns out, a second commuting copy of the Virasoro algebra, with modes L n , can formally be considered for the conformal transformations of the antiholomorphic variable z. It is not our purpose to give an introduction to 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 accomodated in a finite number of representations (such CFT 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 0 and L 0 with real eigenvalues h and h, and, more importantly, is annihiliated by all positive modes L n>0 , L n>0 . It is in particular characterized by a total weight ∆ = h + h (its eigenvalue under L 0 + L 0 , the real dilation generator) and is of course quasi-primary. The action of any string of negative Virasoro modes L n<0 , 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 = ∂ z and L −1 = ∂ z act as derivatives on any field. All of them are eigenfields of L 0 and 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 dimension, 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 ( (4.11) These forms suggest that the conformal weights h i , 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 has the form (4.5) 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 12 . It turns out to be often the case that the structure of a Virasoro highest weight representation implies that the correlators Φ(z, z) . . . 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 2d CFTs can be paraphrased in the following way: to completely solve a CFT, i.e. 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 non-locality 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 non-local interactions and/or non-local 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 (LCFT). What follows is a very basic introduction to the salient features of LCFT; various reviews and applications may be found in the special issue [GRR13]. Let us also mention [Fl03] which reviews the extension to LCFT of the calculational tools used in CFT.
For the highest weight representations discussed above, the operators L 0 , L 0 are diagonalizable. LCFTs have the distinct feature to include Virasoro representations for which L 0 and 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 0 on them would be typical of a rank 2 Jordan cell, namely where Ψ is called the logarithmic partner of the primary field Φ, and a similar action of L 0 (with 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 for instance [KR09]. 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 a 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 (h, h), read 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 absorbed in the normalization of Φ or of Ψ; likewise, the logarithmic partner Ψ is defined up to a multiple of Φ without affecting the defining relations (4.12). The chiral version of the above 2-point functions reads 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 0 on Ψ. Under a finite dilation w = αz, the transformation laws of Φ and Ψ read One may check that the form of the correlators (4.13) is indeed invariant under the replacement Let us also note that the scaling (3.3) 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 nonlogarithmic cousins, altough a number of general features are known. On the statistical side, few models have been thoroughly studied, as their non-local 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 non-trivial) 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 [Gu93] and then investigated in greater detail in [GK99,GR06]. It has the following Lagrangian realization in terms of a pair of free, massless, Grassmanian scalar fields θ,θ, Several fields in this theory form logarithmic pairs, like the identity I and the composite field θθ. We note that (4.13) then implies the somewhat unusual relation I = 0, which indeed follows, using the rules of integration over Grassmanian variables, from the fact that the above action does not depend on the constant modes of θ andθ. Since this is a free scalar theory, all correlators of fields that are local (i.e. product of derivatives) in θ,θ are polynomials in the derivatives of the Green function (the kernel of the inverse Laplacian −4∂∂) given in complex coordinates by 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 + L 0 . As said above, the nondiagonalizability of L 0 , L 0 is the hallmark of LCFTs. The logarithmic minimal models form an infinite series of such lattice models [PRZ06].
The rest of this review is devoted to discussing the variables of the 2d 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.

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 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 13 , 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 has clearly an effect at finite volume, but also in the infinite volume limit if some of the boundaries are kept at finite distance (for example 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 non-zero (and uniform) weight with respect to the invariant measure P Γ . 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 Γ to be an L × M rectangular grid in Z 2 , with open boundary conditions: the non-boundary sites are conservative and have maximal height value equal to z ⋆ i = z i = 4, whereas the boundary sites have 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 R 2 ; in that limit, all multisite probabilities are fully invariant under translations.

One-site height probabilities
As a warm up for what has to come, we ask the following: what is the probability P Γ (h i = a) 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 mean we have is the general criterion of recurrence, namely the non-existence 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 For h i = 1 to be in a recurrent configuration C, the height of none of its neighbours N, E, S or W can be equal to 1 (as they would form a forbidden subconfiguration). Following the clever trick proposed in [MD91], we consider a new gridΓ i by deleting from Γ the vertex i and the four edges incident to it. We also define from C a new configurationC onΓ i by setting 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 ifC is recurrent onΓ i . We thus obtain where the matrix in the numerator has been extended by a one-dimensional diagonal block labelled by the vertex i, without changing the value of the determinant. One then can write with B(i) the defect matrix given by and B(i) is zero everywhere else. We obtain 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 P 1 (by translation invariance). As the matrix ∆ Γ becomes the discrete Laplacian on Z 2 in that limit 14 , standard results yield [MD91] 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 P 4 > P 3 > P 2 > 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 to compute the corresponding probabilities. To undertand 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 T , a site j is called a predecessor of i if the unique path in 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 neighbours, for 0 k 3. A careful analysis of the burning algorithm shows the following [Pr94], For a = 1, we see that P Γ (h i = 1) is related to X 0 (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Γ i 14 The reader will legitimately point out that the Laplacian on Z 2 has a zero mode and is therefore not invertible. A closer look at the determinants (5.6) however reveals that they only depend on differences of the inverse matrix entries, which are perfectly well-defined. by adding one edge between one neighbour of i and i itself (four different possibilities). Thus both points of view coincide and lead to the same local modification Γ →Γ i . The next case is P Γ (h i = 2), related to X 1 (i). Here the situation is dramatically different because the condition that i has only one predecessor among its nearest neighbours is highly nonlocal. The reason for this is that there are two manners for a neighbour 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 neighbour 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, the second one is non-local and more difficult. The same remark applies to the fractions X 2 (i) and X 3 (i), and make 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 P(h i = a), 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 [Pr94] obtained the first expressions for P 2 , P 3 and P 4 , but these were given in the form of multivariate integrals. The problem was reconsidered in [JPR06], where the following explicit values were conjectured, 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 neighbour of its starting point, which was then computed in terms of dimer arrangements [PPR11]. The second proof also used the relation with LERW passage probabilities but within a much more general approach [KW15]. Finally the third one [CS12] carried out the direct computation of the multiple integrals left open in [Pr94]. Let us mention that the technique developed in [KW15] to enumerate 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 P 2 , P 3 to just a few lines, see [PR17]. Most of the height correlators presented below have been computed using this technique. Also noteworthy in this context is the work [KW16] which presents a direct and elementary derivation of the average height h = a a P a on planar lattices (from the formulas above, it is equal to 25 8 on Z 2 ) without computing the individual height probabilities.
Ironically, the four numbers 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 P a were found.
Even though the explicit expressions for the P a 's have the same level of simplicity, the far larger complexity of the combinatorial problem posed by the calculation of P a 2 hints at a striking difference of nature between the height 1 and the higher heights: the height 1 is essentially local, the others are non-local. This will soon be confirmed.

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 [MD91], and 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 the 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 The three clusters on the right of (5.10) 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 are best computed using the methods of [KW15]. Explicit calculations become fairly tedious as the size of the cluster increases.

Height correlations
In terms of the subtracted height variables 15 , the n-point correlation functions are given by 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 ℓ | → +∞. 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 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 configurationsC on a locally modified latticeΓ is simply repeated for the neighbourhood of each 15 In order to make contact with fields, we slightly change the notation h i → h(i) for the height at site i.
cluster. Thus the probability to find a height h(i 1 ) = 1 at site i 1 , a height h(i 2 ) = 1 at site i 2 , and so on, is equal to 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 [MR01,PR17]. By construction, the 1-point function vanishes, σ 1 (i 1 ) = 0 (the relation (4.4) is indeed the main motivation for the subtraction). The 2-point function is found to be 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 (4.7) indeed retains the first term only, the form of which has the expected form (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 the heights 1. This is consistent with the fact that the presence of many heights 1 in a configuration makes it more likely to be non-recurrent. Interestingly, the calculation can be carried out in Z d , with the result that the correlation decays like r −2d , giving a dimension ∆ = d [MD91].
The mixed correlator of a height 1 and a height 2, 3 or 4 is harder. They have first been obtained in [PGPR08,PGPR10] by using classical graph-theoretic techniques, and then reconsidered and extended in [PR17] using the results of [KW15]. Whatever the method used, one has to evaluate the fractionsX k (i 1 ) of spanning trees, as defined in Section 5.2, but on a lattice modified around the i 2 where the 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, where γ = 0.577 216... is the Euler constant. The first subdominant correction is of order r −6 and contains a non-trivial angular dependence, like in (5.14), but also a log r term [PR17]. 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 forms (5.14) and (5.15) are strongly reminiscent of those in (4.13), but do not quite match. If in the scaling limit, the heights 1 and 2 were to converge to a logarithmic pair {h 1 (z), h 2 (z)}, one would think that σ 1,1 and σ 2,1 ought to go over to the 2-point functions h 1 (z 1 )h 1 (z 2 ) and h 2 (z 1 )h 1 (z 2 ) 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 1 , i 2 ) for a, b > 1 does not merely reduce to the evaluation of numbers like X a−1,b−1 (i 1 , i 2 ) which would generalize the numbers X a−1 (i) defined earlier and enumerate the spanning trees with fixed numbers of predecessors among the nearest neighbours of i 1 resp. i 2 . Indeed the possibility that neighbours of i 1 are predecessors of i 2 , or vice-versa, substantially complicates the matter. Details on how to perform the correct counting have been given in [PR17].
To reconcile the previous lattice results and the LCFT predictions, we pause for a while to examine the effects of a seemingly unrelated observable.

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 Z 2 ). We show now that the situation changes if we make some of the bulk sites dissipative [PR04]. 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 make 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 1 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 i 1 = 4 (on Z 2 ) to 5 (a higher value would not make much difference). In turn this changes by 1 the diagonal entry (∆ Γ ) i 1 ,i 1 of the toppling matrix, that is, ∆ Γ → ∆ Γ + D i 1 , with (D i 1 ) i,j = δ i,i 1 δ j,i 1 . More generally, the new toppling matrix∆ n ≡ ∆ Γ + D i 1 + D i 2 + . . . + D in 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 det∆ n / det ∆ Γ , first at finite volume, then in the infinite volume limit.
We start by inserting dissipation at single site i, far from the boundaries. The ratio is easy to compute since the defect matrix D i has rank 1, 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 16 .
The same divergence is present in the ratios det∆ n / det ∆ Γ , 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 ratios det∆ n / det∆ 1 , 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 1 − i 2 | for n = 2, where γ 0 = 1 2π (γ + 3 2 log 2) + 1. If we denote by ω(z, z) the field that describes, in the scaling limit, the insertion of dissipation at a bulk site, the previous two equations would imply Interestingly, they exactly match the last two equations of (4.13), with the logarithmic pair {Φ, Ψ} identified with {I, ω}, both fields having the weights h = h = 0 (the identity field is primary). Moreover the logarithmic term in the 2-point correlation fixes the coefficient λ of the logarithmic pair (I, ω) equal to λ = − 1 4π so that L 0 ω = L 0 ω = − 1 4π I. The relation I = 0, as noted for the free symplectic fermion theory, is here understood as being given by the inverse of the divergent quantity in (5.16).
The lattice calculation of det∆ 3 / det∆ 1 , corresponding to the insertions of three dissipative sites, is not difficult and yields the following 3-point correlation, with z ij ≡ z i − z j , It is fully consistent with the general 3-point correlators of fields in a logarithmic pair [Fl03]. Many additional checks have been carried out [PR04] which all confirm the consistency of the above field assignment. It has been shown [JPR06] that the bulk dissipation field can be realized in terms of symplectic free fermions as in the sense that the correlators of this composite field, computed in the symplectic fermion theory, reproduce the above expressions.

Height correlations cont'd
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, 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 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 [JPR06], is that a lattice n-point height correlator is described in the scaling limit by an (n + 1)-point field correlator, a 2 ,...,an (i 1 , i 2 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 = h forming a logarithmic pair such that Comparing with (4.13), we see that the insertion of dissipation at infinity through ω(∞) allows a non-zero value of A, and cures the problem encountered in Section 5.3. From the dominant terms in (5.14) and (5.15) for the lattice correlations σ 1,1 (i 1 , i 2 ) and σ 2,1 (i 1 , i 2 ), we infer that the (subtracted) bulk lattice height 1 and height 2 variables converge to fields, h 1 (z) and h 2 (z), that form a logarithmic pair of weight ∆ = 2. Moreover if we assign them the same normalization as their lattice companions (A = − P 2 1 2 ), we find the parameter of the logarithmic pair (h 1 , h 2 ) equal to λ = − 1 2 . The explicit results for σ 3,1 (i 1 , i 2 ) and σ 4,1 (i 1 , i 2 ) [PGPR10] show that the fields h 3 (z) and h 4 (z) are also logarithmic partners of h 1 (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 17 of h 1 (z) and h 2 (z) with known coefficients, and other values for the β a [JPR06]. These field assignments predict that the lattice correlation of heights larger or equal to 2 behave asymptotically as Because α 4 is the only negative coefficient among the α 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 behaviour (5.24) [JPR06]. 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 [MD92] 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 [MR01] 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 seperation distances. Surprisingly the first non-zero term has degree −8, (5.25) implying that its scaling limit, corresponding to the CFT 4-point function and is much more instructive: it is precisely the expression we obtain for the 5-point CFT cor- if we assume that the height 1 field h 1 (z) is a primary field of dimensions (h, h) = (1, 1) 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 [DFMS97], from which the correlator can be fully determined, the result being exactly the function (5.26) ! From this result, one can actually infer that the previous correlator h 1 (z 1 ) h 1 (z 2 ) h 1 (z 3 ) ω(∞) must vanish if it is to be symmetrical in the three insertion points [PR17].
The lattice 3-point correlators σ a,1,1 (i 1 , i 2 , i 3 ), a 2, have been computed more recently in [PR17]. 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, where the coefficients α a are those given in (5.23). In addition to being very simple, this expression is surprisingly non-logarithmic. Because its scaling limit should be given by h a (z 1 ) h 1 (z 2 ) h 1 (z 3 ) ω(∞) , 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 2 as logarithmic partners of h 1 . Since h 3 and h 4 can be regarded as linear combinations of h 1 and h 2 , it is sufficient to consider h 2 . Inspired by the conformal representations appearing in the bosonic sector of the symplectic theory [GK99], the following proposal has been made in [JPR06] regarding the conformal nature of h 2 (z, z). A more complete account will be presented in Section 8.
The field h 2 (z, z) is not primary since it transforms into h 1 (z, z) under dilations. It is also not quasi-primary because its L 1 and L 1 transforms generate two new fields, ρ(z, z) and ρ(z, z) respectively, with weights (0, 1) and (1, 0). Moreover the field ρ(z, z) is left primary and its L 1 transform is equal to κI; likewise ρ(z, z) is right primary and its L 1 transform is also equal to κI. All this results in the following transformation law of h 2 (z, z) under a general conformal transformation z → w(z) and z → w(z), This conformal transformation law of h 2 is sufficient to compute correlators involving h 2 , but substantially complicates the calculations. Using this transformation and the left and right level 2 degeneracy of h 1 (z, z), the required correlator can be nonetheless determined. The result reads [PR17] When the tree points z i are aligned horizontally, it exactly reproduces the lattice result (5.27) for a = 2. The cases a = 3, 4 follow by multiplying by the proper coefficient

Minimal height cluster correlations
The calculation of occurrence probabilities of minimal subconfigurations have 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 [MR01]. It turns out that each such cluster S can be specified by a triplet (a, b 1 , b 2 ) of real numbers. We define as before subtracted variables where δ S(i) denotes the event "a minimal subconfiguration S is found around site i", and P S (i) = E δ S(i) is the probability of such an event. The mixed correlator of two such variables takes the form 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 (a, b 1 , b 2 ) = (P 1 , 0, 0) while for the second one in (5.10), 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 (5.31). Interestingly it has been observed [MR01] that these fields have a realization in terms of the symplectic free fermions discussed at the end of Section 4. Indeed one may check that the explicit fields given by h S (z, z) = − a ∂θ ∂θ + ∂θ ∂θ + (b 1 + ib 2 ) ∂θ ∂θ + (b 1 − ib 2 ) ∂θ ∂θ , (5.33) reproduce the above 2-point correlators, as well as the higher order correlators computed in [MR01], provided the dissipation field ω, proportional to θθ, is inserted in the correlators, as explained earlier. In particular, the height 1 field h 1 (z, z) is recovered upon setting a = P 1 and b 1 = b 2 = 0. Let us note a generic field h S (z, z) is a linear combination of three fields with different conformal weights, namely (1, 1), (2, 0) and (0, 2), and with therefore different conformal transformations; the last two are responsible for inducing an angular dependence in the correlators. The field realization (5.33) has been proved in a much greater generality in [Je05a]: any lattice observable based on a conservative local bond modification 18 converges in the scaling limit to a field of the form (5.33).
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 θ,θ. Indeed, because the 2-point correlators of θ,θ are given by the Green function, the correlators of any fields that are local in θ,θ 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 perturbation 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 accomodated 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 non-local constraints on the nearest neighbours (some of them should not be predecessors). Using the technique developed in [KW15], the 1-site probabilities P a 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 P a 2 quoted in (5.9) are not much more complicated than for P 1 . This is no longer the case for large distance correlations σ a 2,1 (i 1 , i 2 ). The analysis [PR17] shows that those correlators are expressed in terms of sums of product of Green matrix entries over a path connecting i 1 to i 2 , and thus in terms of quantities that are not local in the Green matrix. It supports the view that the height fields h a 2 do not belong to the free symplectic fermion theory. A detailed analysis of this question has been carried out in [JPR06], and has reached the same conclusion. Section 8 below summarizes this somewhat strange situation. 18 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 [MR01]. An example of a non-conservative bond modification is given in Section 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 towards 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 a 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 non-local problem, even if i and j are close to each other, even nearest neighbours [PP11].
Connectivities between neighbouring sites can be handled in much the same way as heights 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 can generalize the matrix by assigning arbitrary weights to the oriented edges of the graph Γ ⋆ = (V ⋆ , E ⋆ ). We define x ij as the weight carried by the edge from vertex i to vertex j (x ij = 0 means that there is no edge from i to j), and we set, for i, j ∈ V , In the context of the sandpile model, the difference x i⋆ = y i − j =i 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, 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 ℓ 2 in a permutation σ brings a sign (−1) ℓ−1 coming from the parity ε σ and another sign (−1) ℓ from the product of non-diagonal entries of ∆, resulting in an overall sign −1 per proper cycle. A cycle of length ℓ = 1, corresponding to a vertex left invariant by σ, brings no sign. The term k = 0 is simply equal to i y i , as the only permutation with no proper cycle is the identity. In combinatorial terms, the product 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 [Pr94,IPR07]. Thus det ∆ is the sum over the oriented spanning trees on Γ ⋆ , 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 Z 2 and compute the probability that the outgoing arrow from the site i is oriented to its right neighbour. This amounts to compute the fraction of spanning trees with such an arrow, namely 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 neighbour E, one could simply set to 0 the weights between i and its three neighbours S, W and N. It is however computationally more efficient to set the weight from i to its right neighbour to x, and take the limit x → +∞ so as to give the edges to the other three neighbours a relative weight equal to 0. This implies that we define a new matrix∆ which coincides with ∆ except on two entries, namely∆ i,i+ê 1 = −x and∆ i,i = x + 3. As before we write∆ = ∆ + B → (i) for the defect matrix  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 1 and i 2 is (5.38) In the infinite volume limit and for a large distance, the following two-point probabilities are found at dominant order, 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, in agreement with the fact that ρ ↑ (z, z) is formally the rotated form of ρ → (z, z) (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, ρ ← (z, z) = −ρ → (z, z) and ρ ↓ (z, z) = −ρ ↑ (z, z). 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 belongs to a tree is the sum of the probabilities to find it in either of the two possible orientations, is thus equal to 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 x Correlations of unoriented edges should decay faster than that 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 −2 ). Indeed explicit calculations yield a r −4 decay, 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 x −x −x x is conservative. More precisely, the following fermionic expressions reproduce the above correlations, 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 φ ↔ ∼ (∂ + ∂)ρ → and φ ∼ i(∂ − ∂)ρ ↑ . It turns out not to be quite the case.

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 behaviour from their bulk versions. In the field theoretic description, the boundary fields have to be properly identifed, 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 H = {(x, y) ∈ R 2 : y 0}. There is only one boundary to care about, and the invariance under horizontal translations is preserved.
From our earlier discussion of conservative versus 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 H 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 19 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 neighbouring 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 20 . Table 1: Numerical coefficients for one-site height probabilities on the UHP, with γ = γ + 5 2 log 2. They satisfy the relations a c a = a d a = 0.

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 {(x, y) ∈ Z 2 | y 1}. 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 ) −1 can be easily computed in terms of that on the full plane, y 1 ),(x 2 ,y 2 ) = G (x 1 ,y 1 ),(x 2 ,y 2 ) − G (x 1 ,y 1 ),(x 2 ,−y 2 ) , (6.1a) G cl (x 1 ,y 1 ),(x 2 ,y 2 ) = G (x 1 ,y 1 ),(x 2 ,y 2 ) + G (x 1 ,y 1 ),(x 2 ,1−y 2 ) , (6.1b) for y 1 , y 2 > 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 = 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 the 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 P op/cl 1 (y) to find a height equal to 1 at a distance y from the boundary. It can be computed by using the formula (5.6) where ∆ −1 Γ 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 [BIP93], with the result σ op 1 (y) = P op 1 (y) − P 1 = P 1 4y 2 + . . . , σ cl 1 (y) = P cl 1 (y) − P 1 = − The analogous results for higher heights were obtained somewhat later [PR05a,JPR06] and were the first to firmly establish their logarithmic nature. They take the following form, valid for boundary condition in the scaling limit. given below. The coefficients are explicitly known and are collected in Table 1. One may check that the relations (5.23) expressing h 3 and h 4 linearly in terms of h 1 , h 2 are confirmed. The distinctive change of sign between the two boundary conditions, the fact that for fixed a, both are controlled by the same constants c a and d a and the equality d 2 = c 1 are striking. As will be explained below and in one of the next sections, all three features will follow from the CFT picture.
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 {(x, y) ∈ Z 2 | y > x}, 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 the height 1 and the diagonal boundary, while the subdominant terms are different [PR17].
The 2-site height correlators in the bulk of the UHP, at sites i 1 = (x 1 , y 1 ) and i 2 = (x 2 , y 2 ), and which involve the same subtractions as before, have been computed in [PR17] when the two sites are far from the boundary and far from each other, again using the technique developped in [KW15]. They depend on three real variables, the horizontal distance x = x 1 − x 2 between the two sites and their vertical positions y 1 , y 2 . 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, the results for a = 1, 2 take the following form, at dominant order, σ op 1,1 (0; y 1 , y 2 ) = σ cl 1,1 (0; y 1 , y 2 ) = P 2 1 2 P (y 1 , y 2 ) + . . . , (6.6a) σ op/cl 2,1 (0; y 1 , y 2 ) = P 2 1 2 P (y 1 , y 2 ) log y 1 + γ + 5 2 log 2 + Q(y 1 , y 2 ) log y 2 + y 1 y 2 − y 1 + H op/cl (y 2 1 , y 2 2 ) y 2 1 y 2 2 (y 2 1 − y 2 2 ) 4 + . . . , 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 (5.23). Let us discuss these results in the CFT picture, using what we already know about the height fields h a (z, z). For a homogeneous boundary condition like here, boundary CFT prescribes to compute bulk correlators on the UHP by viewing a field φ(z, z) with conformal weights (h, h) as the product φ h (z)φ h (z) of two chiral fields of weight h and h respectively, located at the points z and z (the latter being in the lower-half plane) [Ca84]. 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 0 . Indeed let us consider a logarithmic pair (Φ(z, z), Ψ(z, z)). If we factorize the logarithmic partner as Ψ(z, z) = ψ h (z)ψ h (z), we find from the action of L 0 , that the chiral factorization of the primary partner is Φ = φ h ψ h . The same argument with L 0 shows that an equally good factorization is Φ = ψ h φ h . Let us first see how this works for σ op a (y). Their dominant terms, made explicit in (6.2) and (6.3), should correspond to h a (z, z) 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 = h = 1, the chiral factorizations of the height fields h 1 and h 2 read h 1 (z, z) = ψ(z)φ(z) and h 2 (z, z) = ψ(z)ψ(z). The CFT formalism gives the general forms (4.14), With the value λ = − 1 2 noted in Section 5, these forms exactly reproduce the lattice results (6.3), including the relation d 2 = c 1 = 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 h a (z, z)ω(∞) 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 [JPR06].
The conformal calculations required to account for σ op a,1 are only technically more involved. The needed chiral correlators are φ(z 1 )ψ(z 1 )φ(z 2 )ψ(z 2 ) for a = 1 and ψ(z 1 )ψ(z 1 )φ(z 2 )ψ(z 2 ) 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 [PR17], we merely give the results, valid for any relative positions of the two heights, z 1 = (x 1 , y 1 ) and z 2 = (x 2 , y 2 ), and h 2 (z 1 , z 1 )h 1 (z 2 , z 2 ) op = P 1 32y 2 1 y 2 (6.10) One can check that setting x 1 = x 2 exactly reproduces the lattice results σ op 1,1 (0; y 1 , y 2 ) and σ op 2,1 (0; y 1 , y 2 ) reported above. The analogous calculation for the closed boundary has been carried 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 non-trivial 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 [PR04]; it was found in all cases that the CFT predictions compare successfully with the lattice results.

Changing the boundary condition
We have considered so far two different boundary conditions, the open and closed conditions. This allows to address 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 ∆ i,i = z ⋆ i = 4 and a closed one has ∆ i,i = 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 a 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 usually, we can readily write the infinite volume limit of the ratio as where (D In ) i,j = δ i,j for i, j ∈ I n , is zero elsewhere, and I n = {(ℓ, 1) : 1 ℓ n} is the set of sites being closed. Using the relation (6.1a) expressing G op in terms of the Green matrix on the full plane Z 2 , the matrix in the determinant reads (6.12) By the horizontal translation invariance of G op , this is a Toeplitz matrix of the form a ℓ−ℓ ′ . 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, For large n, the asymptotics of such Toeplitz determinants is well-known (see f.i. [DIK13]), and leads to the following result [Ru02], π n , n ≫ 1, (6.14) with G= 0.915965... the Calatan 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, 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 Z cl (n) Z cl (1) , which one can write as where the entries b m are the Fourier coefficients of a symbol σ cl (k) given by Its Fourier coefficients are well-defined for α > − 1 2 , diverge in the limit α → − 1 2 but nonetheless keep the ratio (6.16) finite. Remarkably, for large n, it takes the form [Ru02] Z cl (n) Z cl (1) ≃ A n 1/4 e 2G π (n−1) , n ≫ 1, (6.18) for the same constant A as above. Before discussing the CFT side, let us remark that the exponential factors in (6.14) and (6.18) 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 e 4G π 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 e 4N f 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 ≫ n ≫ 1 the ratios are asymptotic to The An explicit calculation [Ru02] confirms this and yields f op − f cl = 2G π , in agreement with the above results. To fix the ideas, the effective number of values taken by a boundary height is e fop = 3.70 at an open site, and e f cl = 2.07 at a closed site.
In the CFT approach, a change of 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 21 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 1 ) φ b,a (x 2 ) accounts for the change at x 1 from condition a to condition b, and then back from b to a at x 2 , but does not account for the exponential terms related to the difference of boundary free energies of condition a versus condition b, namely the terms we have just discussed in the previous paragraph. These are clearly non-universal, i.e. 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 φ op,cl (0) φ cl,op (n) = φ cl,op (0) φ op,cl (n) . The two are indeed equal on the lattice and asymptotic to n 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 = − 1 8 , with a correlator given by 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 usually, this conjecture can be put to the 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 22 , the precise form of which depends on the other fields involved. A first and simple test is to look at a 4-point function 23 , for instance φ op,cl (x 1 ) φ cl,op (x 2 ) φ op,cl (x 3 ) φ cl,op (x 4 ) , which should describe the effect of closing the sites on two disjoint interval [x 1 , x 2 ] and [x 3 , x 4 ] 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 21 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 boundary condition but are not associated to a lattice observable, unlike the height fields h a for instance. 22 Again the technical assumption is that the field φ op,cl is degenerate at level 2, similarly to the height 1 field h 1 , see Section 5.5. 23 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 boundary condition thus involves an even number of insertions of φ op,cl . For instance φ op,cl (0) φ cl,op (∞) changes the boundary condition from open to closed on the positive real axis. which reduces to the product φ op,cl (x 1 ) φ cl,op (x 2 ) φ op,cl (x 3 ) φ cl,op (x 4 ) when the two intervals are infinitely distant. This unique solution solution reads (with is the complete elliptic integral. To compare with a lattice calculation, we take x i integers, with x 21 , x 32 and x 43 all large, and try to compute the determinant in (6.11) with I = I 1 ∪I 2 the union of the two intervals [x 1 , x 2 ] and [x 3 , x 4 ]. This determinant is no longer Toeplitz, which makes it difficult to compute its asymptotics analytically. Dividing it by the prefactor (x 12 x 34 ) 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 (6.21) is more than satisfactory [Ru02].
The opposite situation -two open intervals in a closed boundary-has also been considered. The appropriate 4-point function can be obtained from (6.21) by making a simple cyclic permutation (x 1 , x 2 , x 3 , x 4 ) → (x 4 , x 1 , x 2 , x 3 ), with the result that K(x) gets replaced by K(1 − x). An equally successful agreement was observed [PR04]. Many other crosschecks have been done, confirming that the open/closed boundary condition changing field is indeed a primary field with conformal weight h = − 1 8 . One of them, particularly convincing, is presented in the next section.

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 the light of what we have learned of the boundary condition changing 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), 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 1 , x 2 ] where the condition is closed. The two homogeneous open and closed conditions can be recovered in the limits x 1 → x 2 and x 1 → −∞, x 2 → ∞. Let us denote by h a (z, z) mix the corresponding quantity in the CFT, given by where the division by φ op,cl (x 1 ) φ cl,op (x 2 ) 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 1 and x 2 . The denominator is known from (6.20).
To compute the numerator, we represent the height fields in terms of the chiral fields as h 1 (z, z) = φ(z)ψ(z) and h 2 (z, z) = ψ(z)ψ(z) (as usually, considering the heights a = 1, 2 is enough) and write the differential equation statisfied 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 φ op,cl (x 1 ) φ cl,op (x 2 ) ψ(z) ψ(z) op in fact depends on that of φ op,cl (x 1 ) φ cl,op (x 2 ) φ(z) ψ(z) op .
All calculations done, one finds that they depend on two integration constants c 2 and d 2 in such a way that the ratios (6.23) take the following forms where y = Re z [JPR06], Although t is complex, both expressions are real on account of t * = 1/t. Let us now discuss the above two limits x 1 → x 2 and x 1 → −∞, x 2 → ∞. For convenience, we set x 1 = −x 2 and examine the limits x 2 → 0 + and x 2 → +∞. To compute the two limits, the important thing to notice is that the complex variable t, now equal to has complex norm equal to 1 and loops anticlockwise around the origin as x 2 varies from 0 + to +∞, starting from 1 + 0i to 1 − 0i. It follows that t itself goes to 1 in both limits but √ t goes to +1 when x 2 → 0 + and goes to −1 when x 2 → +∞. The actual limits yield 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 2 and d 2 , and the equality c 1 = d 2 are all accounted for ! The conformal approach however cannot fix the two coefficients c 2 and d 2 ; these must be determined by lattice calculations. The expressions (6.24) 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 1 → −∞ and x 2 → 0. The conformal transformation w = L π log z can be used to map the UHP onto an infinite strip of width L, with open boundary condition on the left side, closed on the right side. The conformal transformation rules of the fields involved being known, the expressions (6.24) 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 [JPR06] that the conformal predictions and the numerical plots match remarkably well, thereby confirming once more all the field identifications made so far.

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 [Ru07], 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 24 . In terms of height values, either condition means that none of the boundary sites has a height 1 (because each boundary site has an ingoing arrow) or a height 4 (because the burning algorithm would imply that the arrow is pointing down, towards 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 to evaluate 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 [Ru07] for the 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 a, b ∈ {op, cl, ←, →}, 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 a ∈ {op, cl} is unoriented and b ∈ {←, →} 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 •→ • ←• 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 upwards into the bulk of the UHP. In the two cases, the macroscopic configurations of arrows are different. Thus we should distinguish two different fields, φ → op , ← and φ → cl , ← . The detailed analysis confirm 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.
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 φ a,b ⋆ φ b,c ≃ φ a,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 non-trivial, 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 to [Ru07] for more details on these specific points.

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 Φ(z, z), characterized by a pair of conformal weights (h, h); 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 wih 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 case, only open and closed boundaries have been considered.
The case of isolated dissipation is simpler and has been examined in details in [PR04], where isolated dissipation has been considered on a closed boundary only. The calculation proceeds much like those 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 for its bulk version. 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 ω cl ∼ θθ [PR04] and ω op ∼ ∂θ∂θ [Je05b] 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 (x) ∼ (L −2 I)(x).
Boundary height variables are more complicated than dissipation but simpler than the bulk height variables. The first results have been derived by Ivashkevich in [Iv94], where the oneand 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 neighbours. So it leads essentially to the same problems of computing non-local contributions. Both in the bulk and on a boundary, one can write linear identities expressing combinations of non-local contributions in terms of local ones, themselves calculable with a defect matrix. In turn, the non-local contributions can be used to calculate probabilities. In the bulk, the linear system is underdetermined and cannot be inverted to provide the required non-local contributions, and then the probabilities themselves. The main observation made in [Iv94] was that on a boundary, the linear system can be inverted, and therefore allows to compute 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 [Iv94]) and recall those in the bulk, as given in Section 5.1, 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 1 , x 2 ) and τ cl a,b (x 1 , x 2 ) on an open or closed boundary were also computed in [Iv94]; all of them were found to decay like |x 1 −x 2 | −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 θ,θ. This was done independently in [Je05b] and in [PR05b], 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 to distinguish more efficiently different fields which would otherwise have the same 2-point correlators in the non-massive (critical) limit. where the θ,θ 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 (6.1a) and (6.1b). 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.

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 flavour. It has been considered and discussed in [IPRH05,IPR07] in the dimer model, intrinsically related, like the sandpile model, to spanning trees.
Let us consider a rectangular portion of 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 Γ ⋆ , obtained from Γ by adding the sink vertex and the edges connecting the open boundary sites to the sink. The graph Γ ⋆ corresponding to a 3 × 3 grid with three open edges and one closed egde is shown in Figure 1.
Once the graph Γ ⋆ is embedded in the plane 25 , the faces of Γ ⋆ 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 (Γ ⋆ ) * is standard: the vertices of (Γ ⋆ ) * are associated to the faces of Γ ⋆ , and two such vertices are connected if their corresponding two faces are separated from each other by an edge of Γ ⋆ . The dual graph of the example above is also shown is 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 Γ ⋆ , it becomes homogeneously closed resp. open in (Γ ⋆ ) * . In addition, the dual graph (Γ ⋆ ) * is the extension (Γ * ) ⋆ by a sink of a dual rectangular grid Γ * , of size slightly different from the original grid Γ. A classical result states that the number of spanning trees on Γ ⋆ is equal to the number of spanning trees on its dual (Γ ⋆ ) * . In fact, for every spanning tree T on Γ ⋆ , there is a unique dual spanning tree T * on (Γ ⋆ ) * such that the two are perfectly interdigitating: the edges of T * are exactly those of (Γ ⋆ ) * which cross the edges of Γ ⋆ not used in 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 on Γ ⋆ can be isomorphically described by those on (Γ * ) ⋆ . As far as the counting goes, the equality of their partition functions can be explicitly written for rectangular grids. If Γ is an L 1 × L 2 rectangular grid with k of its four boundaries being open, the other 4 − k being closed, the dual Γ * is an L ′ 1 × L ′ 2 rectangular grid with swapped boundary conditions, and the following identity holds, The dimensions are related as follows: 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.

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.

The massive sandpile model
In the standard sandpile we have studied so far, the sites in the bulk of grid, that is, the vast majority of sites, are conservative. This meant that when such a site topples, it loses a certain number of sand grains which are all redistributed to its nearest neighbours. 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 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 [GLJ97, TK00,MR01]. A mathematically rigorous proof that all correlations decay exponentially has been provided in [MRS04]. Presumably a non-zero density of dissipative sites could be a 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, ∆ → ∆(t) = ∆ + t 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 ∆ −1 (t) has a kernel given at large distances by G i 1 ,i 2 (t) ≃ 1 2π K 0 |i 1 − i 2 | √ t + . . ., and decays exponentially like e −r √ t at large distances (K 0 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 examle, let us look at the correlation of two heights 1 in the dissipative has a 2-point correlator 26 in the massive fermionic theory that is precisely given by (7.2). The 3and 4-point correlators of the same field have been checked to reproduce the corresponding lattice results. The field h 1 (z, z; M) 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 [MR01], with an unexpectedly simple result. The field describing the minimal cluster variable S in the dissipative model appears to be simply given by where 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 (z, z) is still given by (5.33) 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 [PR05b]. For them too, it is found that the boundary fields given in Section 6.5 get additional terms proportional to M 2 θθ.
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.

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 neighbours (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 [DR89], a dynamics with stochastic toppling rules [Ma91], the formulation of the model on a hierarchical geometric structure like the Bethe lattice [DM90], 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 [PP97,LH02] and numerical simulations [HL03]. More recently, exact calculations of height correlations have been carried out on the honeycomb and triangular lattices.
In [ADMR10], 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 [PR18]. 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 [KW15], already used on the square lattice. The 1-site probabilities on the infinite honeycomb lattice are all rational, and very similar expressions for P 1 a 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 (h, h) = (1, 1), 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 27 exhibit the same bulk and boundary behaviours as on the square lattice. Thus for what concerns the type of the underlying lattice, universality has been explicitly and successfully verified.

Conformal summary
This last section is more specifically oriented towards 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 and not so well understood issues. 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 weight − 1 8 and 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 accomodated 28 in the symplectic fermion theory [GK99]. 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 2 not equal to 1, logarithmic partners of the height 1 field h 1 , as well as the two fields ρ and ρ, to which they transform under L 1 and L 1 respectively.
Although we have not given any physical interpretation of ρ and ρ, they appear to be related to the derivatives of the dissipation field ω [PR17], where δ is a constant which may depend on the lattice considered, and equal to δ = πP 1 2 on the square lattice. In addition, the primary field h 1 may be consistently identified as being proportional to the derivatives of ρ and ρ, where λ is defined from L 0 h 2 = L 0 h 2 = h 2 + λh 1 and depends on the normalizations of h 1 , h 2 . Combining these relations with the previous ones yields the somewhat surprising result that the height 1 field is proportional to the Laplacian of the dissipation field, h 1 ∼ ∂∂ω. The correlator (5.19) confirms this: applying ∂ 1 ∂ 1 ∂ 2 ∂ 2 on it indeed yields a multiple of 1/|z 1 − z 2 | 4 , itself proportional to h 1 (z 1 , z 1 )h 1 (z 2 , z 2 )ω(∞) . From these observations, it follows that all bulk fields encountered so far, namely h a>1 , h 1 , ρ, ρ, ρ → , ρ ↑ , φ S , φ ↔ , φ , ω, I, (8.3) belong to the same conformal representation, as they are all related to each other by the action of Virasoro modes L n or L n . Indeed ρ → and ρ ↑ are not quasi-primary and transform to a multiple of I under L 1 or L 1 , while φ S , φ ↔ and φ are a linear combinations of h 1 and the chiral and antichiral stress-energy tensors T and T . In fact in terms of fermions, all these fields, except h a>1 , are proportional to or are linear combinations of I, θθ, θ∂θ, θ∂θ, ∂θ∂θ, ∂θ∂θ, ∂θ∂θ and ∂θ∂θ. 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 θ∂θ or ∂θ∂θ, 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 φ αβ = ∂θ α ∂θ β for the primary fields, and 28 Very much like the spin field of the Ising model belongs naturally to the free Majorana fermionc theory with c = 1 2 , despite being non-local in the fermions. ψ αβ = θθ ∂θ α ∂θ β for their logarithmic partners, where θ α and θ β are independently either θ orθ, see [GK99] for more details. Conformal representations of the above type are called staggered modules, and have been first studied in [Ro96] in their chiral version. As far as we know, it has been first noticed in [GK96] 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 [KR09] 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 [DF08] and [Ri12]. 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 0 − 1)ψ = λφ and L 1 ψ = ρ. 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 −1 ρ is proportional to λφ, L −1 ρ = β(λφ), (8.4) 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 1 φ = 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 1 , h 2 ), with the same normalization as the height variables on the square lattice, in which case λ = − 1 2 , then one finds β = 1 2 [JPR06]. The field h 1 has been already identified in terms of fermions, and yields a natural choice for the primary field φ on the symplectic side, φ θ = −P 1 ∂θ ∂θ + ∂θ ∂θ . (8.5) As mentioned earlier, the lattice results in the scaling limit are consistent with h 1 being degenerate at level 2, namely (L 2 −1 − 2L −2 )h 1 = 0, see Section 5.5. The same equation is satisfied by φ θ . The only candidate for the logarithmic partner of φ θ is proportional to θθ (∂θ ∂θ + ∂θ ∂θ) 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, ψ θ = −P 1 θθ (∂θ ∂θ + ∂θ ∂θ), (8.6) satisfies L 0 ψ θ = ψ θ − 1 2 φ θ , for the same value λ = − 1 2 . The same OPE reveal in addition that ρ θ = L 1 ψ θ = − P 1 2 θ ∂θ + ∂θθ , (8.7) from which, upon using ∂∂θ = ∂∂θ = 0, one obtains L −1 ρ θ = ∂ρ θ = − P 1 2 ∂θ ∂θ + ∂θ ∂θ = 1 2 φ θ . (8.8) Comparing with (8.4), the value of the logarithmic coupling is found to be β θ = −1 in the fermionic realization. As a consequence, the symplectic fermion theory cannot accomodate 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 β = 1 2 versus the fermionic realization β θ = −1 has been made in [JPR06]; it was shown in particular that the correlations with a trial field h 2 corresponding to a value β = −1 do not match the lattice results 29 . 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 β [KR09], −6βL 2 −3 − 6(β − 2)L −4 L 2 −1 + 8βL −4 L −2 − 2 3 (5β + 2)L −5 L −1 + 4βL −6 ρ. (8.9) Using the relations L 1 ψ = ρ, (L 0 −1)ψ = λφ, as well as the degeneracy condition (L 2 −1 −2L −2 )φ = 0 (and the value c = −2), one can verify that the field ξ satisfies L 1 ξ = L 2 ξ = 0 provided the identity L −1 ρ = βλφ holds. A rather convincing confirmation for the value of β = 1 2 in the sandpile model is therefore to check that the various correlators involving h 2 indeed satisfy the condition ξ = 0 for β = 1 2 . It has been done for the correlator (6.24b). The situation seems therefore 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 β = 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 have been successfully identified.