aa r X i v : . [ c ond - m a t . s t a t - m ec h ] F e b Sandpile models in the large
Philippe Ruelle
Institut de Recherche en Mathématique et Physique
Université catholique de Louvain, Louvain-la-Neuve, B-1348, BelgiumInvited contribution to the Research Topic
Self-Organized Criticality, Three Decades Later (Frontiers in Physics)
Abstract
This contribution is a review of the deep and powerful connection between the large scaleproperties of critical systems and their description in terms of a field theory. Although largelyapplicable to many other models, the details of this connection are illustrated in the classof two-dimensional Abelian sandpile models. Bulk and boundary height variables, spanningtree related observables, boundary conditions and dissipation are all discussed in this contextand found to have a proper match in the field theoretic description.
Contents
Boundaries, boundary conditions and boundary variables 31
In statistical mechanics, critical points are very special points in the space of external parameterswhich control the state of a system. At such a point, the system is scale invariant and its ther-modynamic functions and correlations are characterized by critical exponents and power laws. Inmany cases, physical systems have a finite number of critical points, most often only one. Typicalexamples include the end-point of the liquid-gas coexistence line or the Curie point for ferromag-netic materials. In these cases, a system is brought to its critical point by tuning very precisely afew external parameters to their critical values.In Nature however, power laws are commonplace, and can be found in a large variety of differentphenomena, like avalanches, earthquakes, solar flares, dropplet formation, ... In all these cases, itis certainly not clear what parameters should be tuned, and even if they are perfectly tuned, it isunlikely 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 themselvesdynamically: even if the system is not initially in a critical state, its own dynamics will ineluctablydrive it to criticality and maintain it in that state [BTK87]. This attractive idea has led to theconcept of self-organized criticality (SOC).To support this idea, these authors proposed the sandpile model as a prototypical example ofa system which shows a form of self-organized criticality. Since then, many others models showingSOC have been proposed, as abundantly illustrated in this volume, and in introductory books andreviews [Ba96, Je98, Pr12, Dh06].The present review will be exclusively concerned with specific versions of two-dimensionalsandpile models, formulated by Dhar [Dh90] and called Abelian sandpile models. Even thoughthere are among the simplest and easiest sandpile models to handle, they show a large spectrum ofinteresting and difficult problems which have attracted considerable attention, in both the physicaland mathematical communities. From the point of view taken here (namely their scaling limit andthe emerging conformal field theory), they are to our knowledge the only ones to have been studied.2et, compared to many other equilibrium statistical models, a fair statement is that our presentunderstanding of them is still very poor.Our primary purpose is two-fold, namely to give the unfamiliar reader an introduction of whyand 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 canbe worked out in practical terms. The second part will be illustrated in sandpile models, whichlend themsleves very well to this kind of analysis: they are simple enough that one can follow thesteps in a clear and transparent way, yet they are rich enough to show the difficulties one sometimeshas to face, but also the elegance and the power of the approach. Understanding how of a fieldtheory emerges from a stochastic lattice model enables to gain a probabilistic and intuitive viewof what a field theory is in this context.It turns out that the field theories which appear when analyzing critical systems are conformalfield theories. The simple reason for this is that their large conformal symmetry integrates thefact that critical systems have a local scale invariance. Conformal theories in two dimensionshave been tremendously successful since the 80’s and have led to a deep understanding of thetwo-dimensional critical phenomena. It is certainly not our purpose to give an introduction toconformal field theories, and we will not go very deep into its technicalities, refering to the vastliterature. We restrict to their most basic features, in the hope that these will be sufficient anduseful 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 basicfeatures of the models are recalled. Section 3 is a general description, valid beyond the sandpilemodels, of what is called the scaling limit, which allows to establish the connection between thelarge distance regime of a critical system and the associated field theory. A brief tour of conformaltheories, and specifically logarithmic conformal theories, is presented in Section 4. The applicationof the conceptual ingredients is illustrated in the next three sections. Section 5 focuses on the bulkobservables in the sandpile models, computes the first correlators and explains how these should beunderstood in terms field theoretic quantities. Boundary conditions and boundary observables areexamined in Section 6 as well as the way they should be thought of in conformal theories. Section7 discusses a dissipative variant of the sandpile models and their description by a massive fieldtheory, and also some universality aspects of the sandpile models. The last section summarizesthe 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 thesandpile models as being described specifically by a logarithmic conformal field theory. Intendedto a potentially wider readership, the present review is more devoted to the general connectionbetween critical systems and field theories, illustrated in a specific class of models. The two aresomehow complementary, and, if combined, may provide a more complete overview.
The models we will discuss are discrete stochastic dynamical systems. Their microscopic variablesare attached to the vertices of a finite connected graph Γ = (
V, E ) (with V the set of vertices, orsites, 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 . These3re called height variables and simply give the height of the sandpile at vertex i (i.e. count thenumber of sand grains at i ); they are integer-valued, with h i >
1. A height configuration C is aset 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 numberof edges connecting s to some vertices in a non-empty subset D ⊂ V . Vertices in D are calleddissipative or open, while those in V \ D are conservative or closed. If Γ ⋆ = ( V ⋆ , E ⋆ ) denotes theextended graph in an obvious notation, we define z i to be the coordination number of i in Γ (thenumber of edges in E incident to i , or the number of its nearest neighbours in Γ), and similarly z ⋆i its coordination number in Γ ⋆ . Thus z ⋆i = z i if i is closed, z ⋆i > z i if i is open. Finally we saythat a site i of V is stable if its height satisfies 1 h i z ⋆i . A height configuration is stable ifall sites are stable. Clearly the number of stable configurations is equal to Q 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 thefollowing 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 isdecreased 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 processis pursued for all unstable sites until a stable configuration is obtained. That configurationdefines C t +1 .It is useful to introduce the toppling matrix ∆ as it will play an important role in what follows,∆ i,j = z ⋆i for i = j, − i and j are neighbours (i . e . connected) , , (2.1)for i, j ∈ V . The sand redistribution occurring when a site j topples can then be written as theupdate h i → h i − ∆ j,i for all i ∈ V . The matrix ∆ is like a Laplacian on Γ, with mixed boundaryconditions dictated by the open and closed sites, which induce respectively Dirichlet and Neumannboundary conditions (see Section 6).The above dynamics is well-defined. We see that the total number of sand grains is conservedunder the toppling of a closed site whereas a non-zero number of grains are transferred to thesink under the toppling of an open site. The existence of at least one open site guarantees thatthe relaxation process terminates after a finite number of topplings and motivated the necessityof the extension of the graph Γ by the sink site. Moreover if several sites are unstable during therelaxation process, the order in which they are toppled does not matter. More generally, one may There is no need to keep track of the number of sand grains in the sink, and so we do not assign it a heightvariable. a i for each i ∈ V , whose action on a stable configuration returns the stableconfiguration resulting from the relaxation process after the deposition of a sand grain at i . Onecan then prove that the operators a i commute [Dh90], explaining the qualifier ‘Abelian’ used todesignate the models satisfying this property.The dynamics described above is a discrete Markov chain on a finite configuration space: ateach time step, one applies the operator a i with probability p i (it is the only stochastic elementof the dynamics), going from C t to C t +1 = a i C t . An important question concerns the invariantmeasures, 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 (uniformdistribution). In this case , 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 supportof P Γ are called recurrent; the others are called transient. Being in the support of the uniqueinvariant measure means that the recurrent configurations are those which are in the repeated imageof 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 deceptivelysimple. Not so. What makes the sandpile models non-trivial, fascinating and rich is the support ofthe invariant measure. A generic recurrent configuration is really complicated because the heightvalues are delicately correlated over the entire graph. In the general case, there is no simplercriterion characterizing the recurrent configurations than the following. Let C be a stable config-uration and let C F be its restriction to a subgraph F ⊂ Γ ( F can be assumed to be connected). Wesay that C F is a forbidden subconfiguration if each vertex of F has a height smaller or equal to thenumber of its neighbours in F . It can be shown [MD92] that a forbidden subconfiguration cannotbe in the repeated image of the dynamics (of the operators a i ). It follows that a configurationis recurrent if and only if it contains no forbidden subconfiguration. The simplest example of aforbidden subconfiguration is when F contains just two neighbouring vertices with height valuesequal 1. The criterion also implies that the maximal configuration with heights h i = z ⋆i is recurrentsince a vertex i with height h i > z i cannot be in a forbidden subconfiguration. It is also clearly inthe image of the iterated dynamics since it can be reached from any other stable configuration byan appropriate sequence of a i ’s.The characterizing condition for recurrence shows that the heights of a recurrent configurationare not at all independent. They are not only correlated locally (think of two neighbouring 1’s) butalso globally because asserting that a given configuration is recurrent generally requires to scan theentire graph. For instance the configuration having h i = z i for all i is not recurrent and possessesno other forbidden subconfiguration than the whole configuration itself. Moreover the recurrentstatus 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 sitemakes it recurrent). However the increase of any height in a recurrent configuration preserves therecurrence.The burning algorithm [MD92] (see also the review [Dh06]) provides a convenient way to testwhether a given stable configuration is recurrent. In addition to provide a completely automaticprocedure, 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 The result holds in the more general case where p i = 0 for every i . V ⋆ , F ) ⊂ Γ ⋆ = ( V ⋆ , E ⋆ ) with F ⊂ E ⋆ . This bijectionis 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 arein fact many burning algorithms (the detailed definition requires a certain prescription that islargely arbitrary), each giving rise to a different bijection. This freedom in the choice of a definitealgorithm, 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 leadto 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, Z = { recurrent configs } = det ∆ , (2.2)for ∆ the toppling matrix introduced in (2.1). It is a standard result in combinatorics (Kirchhoff’smatrix-tree theorem) that det ∆ also counts the number of spanning trees on Γ (see Section 5.7for a proof). The previous formula usually implies that the recurrent configurations form anexponentially small fraction of the set of stable configurations (whose number is equal to Q i ∆ i,i ).On a large grid in Z for instance, for which the density of dissipative sites goes to 0 in the infinitevolume limit, the effective number of degrees of freedom per site in a recurrent configuration isroughly 3.21 (as compared to 4 in a stable configuration), meaning that det ∆ ≃ e Gπ N = (3 . ... ) N ,with N the total number of sites and G the Catalan constant.The definition of recurrence implies that all the operators a i map recurrent configurations torecurrent configurations, implying that once the dynamics has brought the sandpile into a recurrentconfiguration, all subsequent configurations are recurrent. Therefore the invariant measure isappropriate to study the long term behaviour of the sandpile.The sandpile models summarized above have raised a large number of interesting and difficultquestions. 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 because as noted earlier, thatmeasure is non-local in terms of the (local) height variables (equivalently the recurrence criterionis non-local).We should remark that the measure P Γ fully depends on all the minute details which arenecessary 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 thelocal thresholds z ⋆i affect the invariant measure. Many features which directly depend on thesedata will change if any of these parameters is modified, like the number of recurrent configurations,the structure of the sandpile group , the geometric structure of the identity configuration , or A notable exception concerns the linear or almost linear graphs, for which the recurrence property usually takesa simpler form and allows for a larger number of explicit results, see for instance [RS92, AD95]). We have mentionned that the operators a i generate an Abelian algebra. But when acting on recurrent con-figurations, they are invertible and therefore generate an Abelian group, called the sandpile group. The sandpilegroup, of order equal to det ∆, has been determined for a number of finite graphs. Recurrent configurations form an Abelian group under the sitewise addition of the heights, followed by relax-ation. This group is isomorphic to the sandpile group. In particular, one of the recurrent configurations is theidentity in the group, and shows remarkable geometric patterns [Cr91, DRSV95, LBR02, CPS08].
The simple idea underlying the scaling limit is this: if we want to concentrate on the large scalebehaviour of a system, let us look at it from far away ! The further away we look at the system, thelarger our horizon is and the larger the distances we keep in sight. At the same time, when lookingfrom a distance, the details get blurred and disappear: one can no longer recognize the type ofgraph and its connectivities are no longer visible. What we see seems to become independent ofthe microscopic details of the model.Rather than stepping back, an equivalent but more convenient way to proceed is to shrinkthe discrete structure (graph or grid or lattice) on which the microscopic variables live. This willinvolve a (real) small parameter ε such that the graph can be embedded in ε Z d (or another shrinkedregular lattice). For smaller and smaller ε , fixing a macroscopic distance ~x = ε ~m ∈ ε Z d amountsto probe larger and larger scales ~m in terms of lattice units, and at the same time, allows to keepa macroscopic distance r = | ~x | under control. The scaling limit corresponds to take ε → /ε .The scaling limit has interesting consequences. The first most apparent one is that the substrateof the rescaled model goes to a continuum, either R d or a part of R d , which may be bounded .This is the first sign that a continuum description ought to emerge in the scaling limit. This isconfirmed 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 variablesdefined on a continuum. If indeed this is expected to happen, exactly what happens is quite subtle. For the scaling to be non-trivial, some external parameters may need to be appropriately scaled with ε . Oneexample of this is discussed in Section 7.1. It is bounded if all the linear sizes of the finite systems in the sequence defining the infinite volume limit growexactly like 1 /ε .
7o realize this, one may note that all the microscopic variables attached to sites contained in aball of radius o (1 /ε ) will actually collapse to the same point in the scaling limit. Thus every pointin the continuum is the convergence point of an infinite number of vertices in the originial discretesetting. The infinity of microscopic variables carried by these vertices will supposedly mix and fuseto generate some kind of degree of freedom located at a single point in the continuum . What isthen the nature of the emerging continuum degree of freedom at that point, and how is it relatedto the lattice variables supposed to collectively generate it ? The conceptual answer is providedby the renormalization group. It roughly goes as follows .The scaling limit as explained at the beginning of this section was carried out in one stroke:all distances are scaled by ε , which is then taken to 0. This limit was only designed to show howthe large distance behaviours can be assessed, but is too rough to answer the question raised inthe previous paragraph. The renormalization group is much better designed conceptually as itorganizes the scaling limit scale by scale and keeps track, at each scale, of the degrees of freedompresent in the system.Let us suppose that we start with a statistical model defined on very large graph, or, to simplifyand fix the ideas, on an infinite lattice. We fix a convenient scale Λ >
1, partition the lattice intoboxes of size Λ and shrink the lattice by a factor Λ. Each box is now of linear size 1 and containsof the order of Λ d microscopic variables. Within each box, we associate an effective, coarse-graineddegree of freedom which takes into account the overall behaviour of the microscopic variables insidethe box (it could be f.i. their average value), and we then compute the sum over the microscopicvariables conditioned by the values of the coarse-grained variables. The result is a statistical modelfor 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 thoseof the first generation, and so on.After the first iteration, each group of roughly Λ d microscopic variables has collapsed to asingle coarse-grained variable of first generation; the statistical model obtained for these can beinterpreted as the original model in which the fluctuations of scale smaller than Λ have beenintegrated out. The second iteration yields a statistical model for the coarse-grained variables ofsecond generation, each of which has integrated the fluctuations of Λ d microscopic variables overscales smaller than Λ , and so on for the next iterations. In this way each iteration, also calledrenormalization, yields a model where more small scale fluctuations have been integrated out, andwhose large scale behaviour should be identical to that of the original model, since the large scalefluctuations have been preserved.The continuum degrees of freedom we were asking about are what the coarse-grained variablesof 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 werelocated around it. Because the coarse-grained variables of one generation are representative ofthose of the previous generation, the continuum degrees of freedom should similarly carry the samecharacteristics as the original microscopic variables. In particular the long distance correlationsshould be identical, at dominant order.The continuum degrees of freedom emerging in the scaling limit are called fields . Unlike theirlattice ancestors, they usually take continuous values. Fields are all what remains when the short- Among the many books and reviews on the renormalization group in statistical mechanics, see for instance thebook by Cardy [Ca96]. , so that one istypically left with an infinite number of different fields. Each field has its own specific propertiesand should be interpreted as the scaling limit of one particular lattice observable (it may alsohappen 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 byits microscopic variables, but also by a probability measure on the configuration space. Thatmeasure, which is a joint distribution for the (non-independent) random microscopic variables, isusually given by a Gibbs measure, and written, up to normalization, as P ( C ) ∼ exp (cid:16) − H [ C ] (cid:17) ,where H is the Hamiltonian of the system, i.e. some given function of the microscopic variableswhich determines the relative probability of a configuration C . What is the equivalent of the Gibbsmeasure for the fields ?According to the discussion above, one starts from the original model and its Hamiltonian H ≡ H . The first renormalization yields the coarse-grained variables of the first generation and acorresponding Hamiltonian H , computed (at least in principle) by summing exp ( − H ) over themicroscopic 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 appro-priate 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 certainfunctional of the fields.Thus if the description of a statistical model is given, in the discrete lattice setting, in termsof a set of microscopic variables ( h i , h i , . . . ) and a Hamiltonian H ( h i , h i , . . . ), it is given in thescaling limit by a set of continuous fields ( φ ( ~x ) , φ ( ~x ) , . . . ) and an action S [ φ , φ , . . . ]. The pair { ( φ ( ~x ) , φ ( ~x ) , . . . ) , S } is refered to as a continuum field theory . More precisely, specifying a setof fields and their action S is only one way to present a field theory; it is also the most comfortableone because it allows to compute the correlators of the various fields, at least in principle.Needless to say, working out the successive renormalizations along with the Hamiltonians H , H , . . . is a formidable task that is, for all practical purposes, impossible to carry out ex-plicitely, except on extremely rare occasions (and for tailored examples). As a consequence thefield theory describing the large distances of a statistical model cannot be obtained in a deductiveway.The situation however is not hopeless. Experience, heuristic arguments or results obtainedon the lattice can often give definite hints about the nature of the seeked field theory. Moreimportantly, and even if one has no clue of what the correct field theory is, the relevance of a F.i. the energy density in the Ising model, namely the product of two neighbouring spins, gives rise to a fieldthat is different from the one obtained from the spin variable itself. Later we will give examples of this in thesandpile models (cluster variables). One should add ‘Euclidean’ field theory because it is formulated on a Euclidean space R d . h i = xε (at site i ) converge in the scalinglimit to the field φ ( x ), it must be true that the scaling limit of the lattice correlators are equal tofield theoretic correlators, namelylim ε → ε − n ∆ h h x ε h x ε . . . h xnε i lattice = h φ ( x ) φ ( x ) . . . φ ( x n ) i FT , (3.1)where the exponent ∆ is determined so that the limits on the l.h.s. exist: as shown below, it willeventually 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 ofany number of lattice observables provided that for each observable O ( i ) around site i inserted inthe lattice correlator, the corresponding field Φ( x ) to which it converges is inserted in the fieldtheoretic correlator,lim ε → ε − P i ∆ i D O ( x ε ) O ( x ε ) . . . O n ( x n ε ) E lattice = h Φ ( x ) Φ ( x ) . . . Φ n ( x n ) i FT . (3.2)So we can write the convergence of a lattice observable to a field as the formal identity,lim ε → ε − ∆ O ( xε ) = Φ( x ) , (3.3)meant to be valid inside correlators.If both types of correlators can be separately computed, the potential infinity of identitiessimilar to the previous one put very strong constraints on the field theory proposed and allow tovalidate it or, on the contrary, to discard it. The more identities we are able to test, the higherthe 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 fieldtheory by comparing its correlators with the lattice quantities, but we cannot compute the fieldcorrelators if we do not know the field theory ! If one thinks of a field theory as being given by aset of fields and an action S , this is indeed a serious problem, because the action cannot be easilyguessed, 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. Onecould 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 theconformal field theories, which therefore provides a field theoretic framework where no action isnecessary. They are discussed in the next section.Knowing the details of the field theory describing the long distance properties of a statisticalmodel is at the same time extremely powerful and immensely complicated. On the one hand, it isindeed powerful because it captures the very essential behaviour of the statistical model withoutbeing cluttered with the many irrelevant lattice effects which make the lattice model so much morecomplex. On the other hand, it is also immensely complicated because every single element in thelattice model which affects the long distances must have a match in the field theory. Such elementsinclude• of course the bulk observables as discussed above,• the boundary conditions, the changes of boundary conditions, and the boundary observables,10 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 presentand known in the field theory, and which can be only very rarely contemplated in full. A renownexception is when we consider critical statistical models, as we do here, which are in additionformulated on two-dimensional domains ( d = 2). Critical systems are primarily characterized by a scale invariance. The correlation lengths of theobservables surviving the scaling limit diverge in the infinite volume limit, so that there is nointrinsic length scale left: the fluctuation patterns appear to be the same at all scales. As aconsequence, the correlation functions of those observables decay algebraically rather than expo-nentially. The large distance 2-point correlator of a typical lattice observable O i located aroundsite i takes the following form, h O i O j i = A | i − j | + . . . , (4.1)where A is a normalization, ∆ is the exponent controlling the decay and the dots indicate lowerorder terms.The field theory emerging in the scaling limit inherits the scale invariance. Further assumingtranslation and rotation symmetries, the scale invariance is enhanced to the invariance under alarger group, namely the group of conformal transformations, i.e. the coordinate transformationswhich preserve angles . In d dimensions, the conformal transformations include the transforma-tions mentioned above, namely the translations ( d real parameters), dilations (1 parameter) androtations ( d ( d − parameters), and the so-called special conformal transformations (or conformalinversions) which depend on an arbitrary vector ~b ( d additional parameters) and take the followinggeneral form ~x ′ | ~x ′ | = ~x | ~x | + ~b ⇐⇒ ~x ′ = ~x + | x | ~b ~b · ~x + | ~b | | ~x | . (4.2)Together these transformations form a finite Lie group isomorphic to SO( d + 1 , R d ∪ {∞} and bijective.In dimension d >
2, a conformal transformation defined locally can be extended to a globaltransformation. The material recalled in this section is completely standard; useful references include [DFMS97] (rather com-prehensive) and [He99] (more focused on critical statistical systems). ~x ) −→ (cid:12)(cid:12)(cid:12)(cid:12) ∂x ′ i ∂x j (cid:12)(cid:12)(cid:12)(cid:12) ∆ /d Φ( ~x ′ ) , (4.3)for some number ∆. Fields transforming that way under global conformal transformations arecalled quasi-primary . Global conformal invariance then fixes the average value of a quasi-primaryfield, h Φ( ~x ) i = 0 , if ∆ = 0 , (4.4)(a constant for ∆ = 0 by translation invariance) and the 2-point correlator of two quasi-primaryfields, h Φ ( ~x ) Φ ( ~x ) i = A | ~x − ~x | if ∆ = ∆ , = ∆ . (4.5)Specializing (4.3) to a dilation, ~x ′ = α~x , we have Φ( ~x ) → α ∆ Φ( α~x ) so that ∆ can be identifiedwith the dimension of the field Φ (in units of inverse length).Global conformal invariance also completely determines the correlator h Φ ( ~x ) Φ ( ~x ) Φ ( ~x ) i of three (and not more) quasi-primary fields, h Φ ( ~x ) Φ ( ~x ) Φ ( ~x ) i = A | ~x − ~x | ∆ +∆ − ∆ | ~x − ~x | ∆ +∆ − ∆ | ~x − ~x | ∆ +∆ − ∆ . (4.6)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,lim ε → ε − h O ~x /ε O ~x /ε i lattice = A | ~x − ~x | = h Φ( ~x ) Φ( ~x ) i FT . (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 criticallattice model.What has been just recalled is valid in any dimension d > d = 2. The global conformal group discussed above remains, but is more convenientlypresented in complex coordinates as the SL(2 , C ) group of Möbius transformations w = az + bcz + d , for a, b, c, d ∈ C satisfying ad − bc = 1.The two-dimensional world has however many more conformal transformations in store. Indeedit is a well-known fact that any analytic map w ( z ) of the complex plane is conformal. Surely ananalytic function requires an infinite number of parameter to fix it (f.i. the coefficients of itsLaurent expansion in some neighbourhood), so that the conformal ‘group’ is certainly infinite-dimensional. The term group is not really appropriate because the composition of analytic mapsis 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.12or instance the map w = L π 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 itsalgebra, for which infinitesimal transformations of the form w = z + ǫ z n +1 are considered. Thecorresponding generators satisfy the famous infinite-dimensional Virasoro algebra ,[ L m , L n ] = ( m − n ) L m + n + c m ( m − δ m + n, , m, n ∈ Z , (4.8)a central extension of the Witt algebra. The real number c is the central charge, and is one ofthe most important data of a two-dimensional conformal field theory (CFT). The modes L , L ± ,whose algebra is unaffected by the central charge, are the infinitesimal generators of the Möbiusgroup, with L − and L 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 forthe conformal transformations of the antiholomorphic variable z .It is not our purpose to give an introduction to CFT, but one can easily conceive the hugedifference between a finite symmetry algebra and an infinite one. A field theory that is to beinvariant 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 Vira-soro algebra, which are all infinite dimensional, and this opens up the possibility that an infinitenumber of fields be in fact accomodated in a finite number of representations (such CFT are calledrational). In this respect, the primary fields are particularly important. They are the strengthenedversion of quasi-primary fields in the sense that they transform tensorially under any conformaltransformation. A primary field is an eigenfield of L and L with real eigenvalues h and h , and,more importantly, is annihiliated by all positive modes L n> , L n> . It is in particular characterizedby a total weight ∆ = h + h (its eigenvalue under L + L , the real dilation generator) and is ofcourse quasi-primary. The action of any string of negative Virasoro modes L n< , L n< on a primaryfield produces infinitely many new fields, called descendant fields , which include all derivatives ofthe primary field, since L − = ∂ z and L − = ∂ z act as derivatives on any field. All of them are ei-genfields of L and L . Together they form a highest weight representation of the Virasoro algebrawhose structure is similar to highest weight representations of simple Lie algebras, the primaryfield 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 com-pletely fixed by their invariance under Möbius transformations. They are more easily written incomplex coordinates ( z ij = z i − z j ), h Φ( z, z ) i = A δ h, δ ¯ h, , (4.9) h Φ ( z , z ) Φ ( z , z ) i = A z h + h z ¯ h +¯ h δ h ,h δ ¯ h , ¯ h , (4.10) h Φ ( z , z ) Φ ( z , z ) Φ ( z , z ) i = A z h + h − h z h + h − h z h + h − h z ¯ h +¯ h − ¯ h z ¯ h +¯ h − ¯ h z ¯ h +¯ h − ¯ h . (4.11)These forms suggest that the conformal weights h i , h i are positive, so that the correlators decreasewith the separation distances, as seems natural from a physical point of view. We will nonethelessencounter physical fields with negative weights, for which the correlators have a different meaning.13ccasionally 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 ). Chiralcorrelators are appropriate for observables living on a boundary, like the real line bordering theupper-half plane, since a boundary is one-dimensional. In this case, only one copy of the Virasoroalgebra remains, so that the fields are characterized by a single conformal weight. Chiral correlatorsare also useful to compute the correlators of bulk variables on surfaces with boundaries, see Section6. The precise structure of a Virasoro highest weight representation ( c, ∆) based on a primaryfield 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-pointcorrelator of a primary has the form (4.5) since it is quasi-primary, and the same is true for the3-point correlator. To go beyond, the global conformal invariance is not enough . It turns outto be often the case that the structure of a Virasoro highest weight representation implies thatthe correlators h Φ( z, z ) . . . i involving the primary field Φ obey differential equations. Four-pointfunctions can be routinely computed in this way. All correlators can then be determined, at leastin principle, without knowing anything of a possible Lagrangian realization of the underlying fieldtheory (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 largedistance limit of a critical model, it is sufficient to know enough of the Virasoro representationsmaking up that CFT. This methodology has been immensely successful since the mid-80’s and hasled to a profound understanding of the many aspects of critical models listed at the end of theprevious section. The Ising model is the prominent example of a model that can be treated thatway, but the same is true of more general statistical models involving local interactions betweenthe microscopic variables.More recently, models showing some form of non-locality have been examined at the conformallight. Sandpile models are in this class, since, as we have seen earlier, the height variables are instrong interaction over the entire domain to form global recurrent configurations. Other modelswith non-local interactions and/or non-local degrees of freedom include percolation, critical poly-mers and more general loop models. It may sound surprising but the conclusion seems to be thatthe conformal approach is still relevant. However the CFTs underlying these models are morecomplex, essentially because the representations of the Virasoro algebra that appear have a farmore 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 reviewsand applications may be found in the special issue [GRR13]. Let us also mention [Fl03] whichreviews the extension to LCFT of the calculational tools used in CFT.For the highest weight representations discussed above, the operators L , L are diagonalizable.LCFTs have the distinct feature to include Virasoro representations for which L and L are nolonger diagonalizable, but instead contain (infinitely many) Jordan blocks of finite rank. To have arough idea of what these representations look like, one can think of a highest weight representationfor which the highest weight is not a single primary field, but a pair of fields (Φ , Ψ), of which only A general 3-point correlator h Φ ( z ) Φ ( z ) Φ ( z ) i is a function of three complex numbers; if all three fields arequasi-primary, that function can be determined by trading z , z , z for the three complex parameters of a generalMöbius transformation.
14 is primary. The action of L on them would be typical of a rank 2 Jordan cell, namely L Φ = h Φ , L Ψ = h Ψ + λ Φ , (4.12)where Ψ is called the logarithmic partner of the primary field Φ, and a similar action of L (with h ). Under the action of the negative Virasoro modes, the Jordan block structure will propagateamong the descendant fields. The presence of Jordan blocks is a sort of minimal ingredient tomake a representation logarithmic; many mathematical complications can and do arise, see forinstance [KR09]. Higher rank Jordan blocks can also appear.A immediate consequence of the presence of Jordan blocks explains the use of the word ‘log-arithmic’: the correlators of fields in a LCFT contain logarithmic terms in addition to the powerlaws encountered before. For instance the 2-point correlators of the logarithmic pair { Φ , Ψ } , bothof weights ( h, h ), read h Φ( z , z ) Φ( z , z ) i = 0 , h Φ( z , z ) Ψ( z , z ) i = B ( z − z ) h ( z − z ) h , (4.13a) h Ψ( z , z ) Ψ( z , z ) i = C − λB log | z − z | ( z − z ) h ( z − z ) h . (4.13b)For rank r Jordan blocks, the 2-point correlators would involve up to ( r − λ 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 h Φ( z ) Φ( z ) i = 0 , h Φ( z ) Ψ( z ) i = B ( z − z ) h , (4.14a) h Ψ( z ) Ψ( z ) i = C − λB log ( z − z )( z − z ) h , (4.14b)It should not be too surprising that Jordan blocks and logarithms go hand in hand. Underdilation by a factor α , a logarithmic term transforms inhomogeneously log z → log z + log α ,reflecting the inhomogeneous action of the dilation generator L on Ψ. Under a finite dilation w = αz , the transformation laws of Φ and Ψ readΦ ′ ( w, w ) = | α | − ∆ Φ( z, z ) , Ψ ′ ( w, w ) = | α | − ∆ { Ψ( z, z ) − λ log | α | Φ( z, z ) } . (4.15)One may check that the form of the correlators (4.13) is indeed invariant under the replacementΦ( z, z ) → Φ ′ ( w, w ) and Ψ( z, z ) → Ψ ′ ( w, w ). Let us also note that the scaling (3.3) must beredefined for the lattice observables described by logarithmic fields since it involves a dilation bya factor 1 /ε , to which the field responds by an inhomogeneous term.Despite all the efforts spent, LCFTs are generally much less understood than their non-logarithmic cousins, 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 outexact calculations on the lattice. On the field-theoretic side, it is not known what a generic LCFTlooks like. The simplest of all (but non-trivial) and probably the only LCFT to be fully undercontrol is the symplectic fermion theory with central charge c = −
2, also called the triplet theory.15t has been introduced in [Gu93] and then investigated in greater detail in [GK99, GR06]. It hasthe following Lagrangian realization in terms of a pair of free, massless, Grassmanian scalar fields θ, ˜ θ , S = 1 π ˆ d z d z ∂θ∂ ˜ θ , ∂ = ∂ z , ∂ = ∂ z . (4.16)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 h I i = 0, which indeed follows,using the rules of integration over Grassmanian variables, from the fact that the above action doesnot depend on the constant modes of θ and ˜ θ . Since this is a free scalar theory, all correlatorsof fields that are local (i.e. product of derivatives) in θ, ˜ θ are polynomials in the derivatives ofthe Green function (the kernel of the inverse Laplacian − ∂∂ ) given in complex coordinates by G ( z, w ) = − log | z − w | .To finish, let us note that the statistical models which have a non-diagonalizable transfer matrix(when there is a proper one) are the natural candidates for being described by LCFTs in theirscaling 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 + L . As said above, the non-diagonalizability of L , L is the hallmark of LCFTs. The logarithmic minimal models form aninfinite series of such lattice models [PRZ06].The rest of this review is devoted to discussing the variables of the 2d sandpile models whichhave been successfully (i.e. with enough confidence) identified in the corresponding LCFT. Theseelements reveal some facets of the field theory at work in sandpile models: the big and completepicture is well out of reach for the moment. The height variables are certainly the first and most natural variables to look at, as they are themicroscopic variables in terms of which the models are defined. The introduction we gave in Section2 was for the Abelian sandpile on an arbitrary graph. If large distance properties should be ratherrobust against local modifications of a graph, they are not expected to be the same on a graphwith a high degree of connectivity (the extreme example being the complete graphs), a regulargraph with a moderate degree of connectivity or a graph with a strong hierarchical structure (likeCayley trees). Most of the results reviewed here are obtained when the graph is a rectangularportion of the square lattice Z ; varying the size of the grid is an easy way to approach the infinitevolume limit and this choice ensures that conservative sites away from the boundary have heightvariables 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 inorder to address universality issues, see Section 7.2.In most cases, the only dissipative sites will be located on the boundary , except when wediscuss the insertion of isolated dissipation. With one exception, we will exclusively consider openand closed boundary conditions, by which we mean that whole stretches of boundary sites are For rectangular grids Γ ⊂ Z , the notion of boundary is clear: when Γ is embedded in Z , the boundary sitesare those which are connected to sites of Z not in Γ. P Γ . Sofar, we have no clear idea of what a generic recurrent configuration looks like. Answers to questionslike “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 thesecorrelations ? Are they exponential or power-lawed ? The computation of multisite height prob-abilities 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 , with open boundaryconditions: 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 (boundaryand corner sites dissipate 1 resp. 2 grains of sand under toppling; both types are connected tothe sink). Thus the toppling matrix is four times the identity minus the adjacency matrix of thegrid, and the height at every site takes values in { , , , } . In this section, all boundaries are sentoff to infinity in the scaling limit, so that the domain converges to R ; in that limit, all multisiteprobabilities are fully invariant under translations. 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 weare interested in the infinite volume limit of these numbers, we take i to be deep in the middle ofthe 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 howto handle it because the only mean we have is the general criterion of recurrence, namely thenon-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 canwrite P Γ ( h i = 1) = { recurrent configs with h i = 1 } det ∆ Γ , (5.1)For h i = 1 to be in a recurrent configuration C , the height of none of its neighbours N, E, S or Wcan be equal to 1 (as they would form a forbidden subconfiguration). Following the clever trickproposed in [MD91], we consider a new grid ˜Γ i by deleting from Γ the vertex i and the four edgesincident to it. We also define from C a new configuration ˜ C on ˜Γ i by setting˜ h j = h j for j
6∈ { i, N , E , S , W } ,h j − > j ∈ { N , E , S , W } . (5.2)Looking back at the criterion of recurrence for an arbitrary graph, it is not difficult to see that aconfiguration C with h i = 1 is recurrent on Γ if and only if ˜ C is recurrent on ˜Γ i . We thus obtain P Γ ( h i = 1) = det ∆ ˜Γ i det ∆ Γ = det[∆ ˜Γ i ⊕ ii ]det ∆ Γ , (5.3)17here the matrix in the numerator has been extended by a one-dimensional diagonal block labelledby the vertex i , without changing the value of the determinant. One then can write∆ ˜Γ i ⊕ ii = ∆ Γ + B ( i ) , (5.4)with B ( i ) the defect matrix given by B ( i ) k,k ′ = − − − − − , k, k ′ ∈ { i, N , E , S , W } , (5.5)and B ( i ) is zero everywhere else. We obtain P Γ ( h i = 1) = det[∆ Γ + B ( i )]det ∆ Γ = det[ I + ∆ − B ( i )] . (5.6)It reduces to the computation of a finite determinant since B ( i ) has finite rank. In the infinitevolume limit (both L, M → ∞ ), this probability converges to a constant P (by translation invari-ance). As the matrix ∆ Γ becomes the discrete Laplacian on Z in that limit , standard resultsyield [MD91] P ≡ lim | Γ |→∞ P Γ ( h i = 1) = 2( π − π ≃ .
073 63 . (5.7)It also means that a recurrent configuration has an average of about 7% of sites with a heightequal to 1.What about higher heights ? We know for sure that the inequalities P > P > P > P holdbecause adding one grain of sand to a recurrent configuration, at a site where h i = a , yields arecurrent configuration if a <
4. However to actually compute these numbers, can one use thesame trick as for the height 1 ? The answer is definitely negative: no local modification of Γ likewhat we did above will allow to compute the corresponding probabilities. To undertand this, weturn to the description in terms of spanning trees.As was briefly mentioned in Section 2, the burning algorithm yields a one-to-one correspondencebetween a recurrent configuration and a spanning tree rooted at the sink site s and growing intothe interior of Γ. In a given spanning tree T , a site j is called a predecessor of i if the unique pathin T from j to the root passes through i . Let us also denote by X k ( i ) the fraction of all spanningtrees for which the site i has k predecessors among its nearest neighbours , for 0 k
3. A carefulanalysis of the burning algorithm shows the following [Pr94], P Γ ( h i = a ) = P Γ ( h i = a −
1) + X a − ( i )5 − a , a . (5.8)For a = 1, we see that P Γ ( h i = 1) is related to X ( i ), namely the fraction of spanning trees onΓ for which the site i is a leaf. All such trees can be obtained from arbitrary spanning trees on ˜Γ i The reader will legitimately point out that the Laplacian on Z 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 matrixentries, which are perfectly well-defined.
18y adding one edge between one neighbour of i and i itself (four different possibilities). Thus bothpoints of view coincide and lead to the same local modification Γ → ˜Γ i .The next case is P Γ ( h i = 2), related to X ( i ). Here the situation is dramatically differentbecause the condition that i has only one predecessor among its nearest neighbours is highly non-local. The reason for this is that there are two manners for a neighbour of i to be a predecessorof i in a given tree. The first one is that the tree includes the edge between the two sites so thatthe neighbour of i is directly connected to i . In the second manner, the tree contains a potentiallylong chain of edges that forms a path between the two sites. The first one is a local connectionand is easy to check, the second one is non-local and more difficult. The same remark applies tothe fractions X ( i ) and X ( i ), and make the calculation of the corresponding probabilities muchmore 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 beforethe completely explicit probabilities could be found. By using a rather heavy graph-theoreticaltechnology, Priezzhev [Pr94] obtained the first expressions for P , P and P , but these were given inthe form of multivariate integrals. The problem was reconsidered in [JPR06], where the followingexplicit values were conjectured, P = 14 − π − π + 12 π ≃ .
173 90 , (5.9a) P = 38 + 1 π − π ≃ .
306 29 , (5.9b) P = 38 − π + 1 π + 4 π ≃ .
446 17 . (5.9c)A few years later, three independent proofs were given. The first one was based on a relationwith the probability of a loop-erased random walk (LERW) to visit a fixed nearest neighbour ofits starting point, which was then computed in terms of dimer arrangements [PPR11]. The secondproof also used the relation with LERW passage probabilities but within a much more generalapproach [KW15]. Finally the third one [CS12] carried out the direct computation of the multipleintegrals left open in [Pr94]. Let us mention that the technique developed in [KW15] to enumerateso-called cycle-rooted groves (which generalize spanning trees to spanning forests with markedpoints) currently provides by far the most efficient way to compute height probabilities, reducingthe calculation of P , P to just a few lines, see [PR17]. Most of the height correlators presentedbelow 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 h i = P a a P a onplanar lattices (from the formulas above, it is equal to on Z ) without computing the individualheight probabilities.Ironically, the four numbers P a are not very useful for a comparison with a field theory, becausethey will have to be subtracted in correlators (see below). And indeed some of the correlators havebeen 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 farlarger complexity of the combinatorial problem posed by the calculation of P a > hints at a strikingdifference 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.19 .2 Height cluster probabilities Cluster height probabilities are a rather obvious generalization of one-site probabilities, by which weask 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 shownbelow.1 2 21 2 31 23 1 2 2 4 3 2 43 2 (5.10)The three clusters on the left belong to the family of weakly allowed subconfigurations , orminimal height clusters, first introduced in [MD91], and which contains the cluster made of asingle height equal to 1. They are minimal subconfigurations in the sense that if one decreasesany of its heights by 1, the clusters become (or contain) forbidden subconfigurations. As was donein 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 finitedeterminants P S ( i ) = det[ I + ∆ − B S ( i )] where the defect matrix B S ( i ) depends on the cluster S considered [MD91, MR01].The three clusters on the right of (5.10) are not minimal and generalize the simple cluster madeof a single height larger or equal to 2. Their level of complexity is comparable to the latter andare best computed using the methods of [KW15]. Explicit calculations become fairly tedious asthe size of the cluster increases. In terms of the subtracted height variables , h a ( i ) ≡ δ h ( i ) ,a − P a , (5.11)the n -point correlation functions are given by σ a ,a ,...,a n ( i , i , . . . , i n ) = E h h a ( i ) h a ( i ) . . . h a n ( i n ) i . (5.12)These are the functions we are primarily interested in for a future comparison with a conformalfield theory. To make the comparison sensible, we have to take the infinite volume limit and thelimit of large separations | i k − i ℓ | → + ∞ . In addition, to avoid the boundary effects –they willbe studied later on–, all insertion points i k are to stay (infinitely) far from the boundaries. Inpractice, one first replaces ∆ Γ by the Laplacian ∆ on Z , 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 newconfigurations ˜ C on a locally modified lattice ˜Γ is simply repeated for the neighbourhood of each In order to make contact with fields, we slightly change the notation h i → h ( i ) for the height at site i . h ( i ) = 1 at site i , a height h ( i ) = 1 at site i , andso on, is equal to P Γ (cid:16) h ( i ) = 1 , h ( i ) = 1 , . . . (cid:17) = det[∆ Γ + B ( i ) + B ( i ) + . . . ]det ∆ Γ = det (cid:20) I +∆ − n B ( i )+ B ( i )+ . . . o(cid:21) . (5.13)The correlators σ , ,..., are obtained by taking appropriate subtractions and the limits discussedabove.The first few n -point correlators can be easily computed for arbitrary configurations of insertionpoints [MR01, PR17]. By construction, the 1-point function vanishes, σ ( i ) = 0 (the relation(4.4) is indeed the main motivation for the subtraction). The 2-point function is found to be( i − i = ~r = r e iϕ ) σ , ( i , i ) = − P r − π − π −
2) cos 4 ϕ ] π r + . . . (5.14)where the dots stand for lower order terms.This first result is instructive for several reasons. First, for large separation distances, thedominant term indicates that the correlation decay is algebraic, which shows that the model iscritical and makes room for a conformal field theoretic description. Second, choosing the scaledimension ∆ = 2, the scaling limit (4.7) indeed retains the first term only, the form of which hasthe expected form (note that the second term, like all other subdominant ones, has only the latticerotation invariance, and is therefore not expected to survive the scaling limit). And third, thedominant term is negative, indicating an anticorrelation between the heights 1. This is consistentwith the fact that the presence of many heights 1 in a configuration makes it more likely to benon-recurrent. Interestingly, the calculation can be carried out in Z d , with the result that thecorrelation decays like r − d , 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 obtainedin [PGPR08, PGPR10] by using classical graph-theoretic techniques, and then reconsidered andextended in [PR17] using the results of [KW15]. Whatever the method used, one has to evaluatethe fractions ˜ X k ( i ) of spanning trees, as defined in Section 5.2, but on a lattice modified aroundthe i where the height 1 is located. This modification affects the toppling (Laplacian) matrix andits inverse, and consequently the whole computation, heavily based on these two matrices. Theresult for a height 1 and a height 2 reads, at dominant order, σ , ( i , i ) = − P r ( log r + (cid:18) γ + 32 log 2 + 16 − π π − (cid:19)) + . . . (5.15)where γ = 0 .
577 216 ... is the Euler constant. The first subdominant correction is of order r − and contains a non-trivial angular dependence, like in (5.14), but also a log r term [PR17]. Theexpressions of σ , and σ , are similar, with different coefficients.The expressions σ a, for a > { h ( z ) , h ( z ) } , one would think that σ , and σ , ought to go over to the 2-point21unctions h h ( z ) h ( z ) i and h h ( z ) h ( z ) i respectively. However, conformal invariance impliesthat the former vanishes identically, whereas the latter is not logarithmic. We could think of com-puting σ , to see what comes out, but large distance correlators with several heights strictly largerthan 1 are far beyond our present computational capabilities. Let us add that the calculationof σ a,b ( i , i ) for a, b > X a − ,b − ( i , i )which would generalize the numbers X a − ( i ) defined earlier and enumerate the spanning trees withfixed numbers of predecessors among the nearest neighbours of i resp. i . Indeed the possibilitythat neighbours of i are predecessors of i , or vice-versa, substantially complicates the matter.Details on how to perform the correct counting have been given in [PR17].To reconcile the previous lattice results and the LCFT predictions, we pause for a while toexamine the effects of a seemingly unrelated observable. In the previous section, the calculation of height probabilities started on a finite grid Γ, where theonly dissipative sites are boundary sites. We did not pay too much attention to exactly whichboundary sites are dissipative; in fact, since the infinite volume limit sends the boundaries off toinfinity, there is no need to know precisely which boundary conditions are used (this is what wemeant when we said that ∆ Γ becomes the Laplacian on Z ). We show now that the situationchanges if we make some of the bulk sites dissipative [PR04]. It is not difficult to understandwhy this is so in terms of spanning trees. We remember that dissipative sites are sites that areconnected 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 fromthe middle of the grid, thereby affecting the macroscopic structure of the spanning trees.To make a bulk site i dissipative, one simply has to connect it to the sink. In the notations ofSection 2, this amounts to increase the value z ⋆i , for instance from z i = 4 (on Z ) to 5 (a highervalue would not make much difference). In turn this changes by 1 the diagonal entry (∆ Γ ) i ,i ofthe toppling matrix, that is, ∆ Γ → ∆ Γ + D i , with ( D i ) i,j = δ i,i δ j,i . More generally, the newtoppling matrix ˜∆ n ≡ ∆ Γ + D i + D i + . . . + D i n defines a new model in which several bulksites i k are dissipative. As a consequence, the height variables at these sites take values in the set { , , , , } .A simple and natural way to evaluate the effect of inserting isolated dissipation is to considerthe change in the number of recurrent configurations, by computing the ratio det ˜∆ n / det ∆ Γ , firstat 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 tocompute since the defect matrix D i has rank 1,det ˜∆ det ∆ Γ = det( I + ∆ − D i ) = 1 + (∆ − ) i,i . (5.16)It is a finite number at finite volume, but diverges in the infinite volume limit, no matter wherethe site i is located. The divergence reflects the fact that the extra value h i = 5 allows enormouslymore recurrent configurations in the modified model . We note that the inverse ratio det ∆ Γ / det ˜∆ is equal to Prob[all h j
4] = Prob[ h i
4] = 1 − Prob[ h i = 5] same divergence is present in the ratios det ˜∆ n / det ∆ Γ , which suggests to change thenormalization and compare the effect of inserting n dissipative sites with respect to the situationwhere there is only one dissipative site, that is, to consider instead the ratios det ˜∆ n / det ˜∆ ,perfectly well-defined. Let us also remark that in the infinite volume limit, the denominator doesnot depend on the location of the (only) dissipative site, so that the ratios are fully symmetric inthe insertion points i k and translation invariant.The first two ratios read, with r = | i − i | for n = 2,det ˜∆ det ˜∆ = 1 , det ˜∆ det ˜∆ = 1 π log r + 2 γ + O ( r − ) , (5.17)where γ = π ( γ + 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 h ω ( z, z ) i = 1 , h ω ( z , z ) ω ( z , z ) i = 1 π log | z − z | + 2 γ . (5.18)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 logarithmicpair ( I , ω ) equal to λ = − π so that L ω = L ω = − π I . The relation h I i = 0, as noted for thefree symplectic fermion theory, is here understood as being given by the inverse of the divergentquantity in (5.16).The lattice calculation of det ˜∆ / det ˜∆ , corresponding to the insertions of three dissipativesites, is not difficult and yields the following 3-point correlation, with z ij ≡ z i − z j , h ω (1) ω (2) ω (3) i = 3 γ + γ π log | z z z | + 116 π " log | z | log (cid:12)(cid:12)(cid:12)(cid:12) z z z (cid:12)(cid:12)(cid:12)(cid:12) + cyclic . (5.19)It is fully consistent with the general 3-point correlators of fields in a logarithmic pair [Fl03]. Manyadditional checks have been carried out [PR04] which all confirm the consistency of the above fieldassignment. It has been shown [JPR06] that the bulk dissipation field can be realized in terms ofsymplectic free fermions as ω ( z, z ) = 12 π θ ˜ θ + γ , (5.20)in the sense that the correlators of this composite field, computed in the symplectic fermion theory,reproduce the above expressions. The multisite height probabilities computed in Section 5.3 were obtained by taking the limitover a sequence of grids of increasing size. Because of the dissipation along the boundaries, theprobabilities 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 impliesthat h i = 5 with probability 1. ω ( ∞ ) in the correlators.Thus the proposal, first made in [JPR06], is that a lattice n -point height correlator is described inthe scaling limit by an ( n + 1)-point field correlator, σ a ,a ,...,a n ( i , i , . . . , i n ) scalim −→ h h a ( z ) h a ( z ) . . . h a n ( z n ) ω ( ∞ ) i . (5.21)It turns out that the proposed field correlations exactly reproduce the form of the lattice resultsobtained in Section 5.3. If { Φ , Ψ } are fields of weights h = h forming a logarithmic pair such that( L − h )Ψ = ( L − h )Ψ = λ Φ, one finds with ∆ = 2 h [JPR06] h Φ( z , z ) Φ( z , z ) ω ( ∞ ) i = A | z | , h Φ( z , z ) Ψ( z , z ) ω ( ∞ ) i = B − λA log | z | | z | , (5.22a) h Ψ( z , z ) Ψ( z , z ) ω ( ∞ ) i = C − λB log | z | + λ A log | z | | z | . (5.22b)Comparing with (4.13), we see that the insertion of dissipation at infinity through ω ( ∞ ) allows anon-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 σ , ( i , i ) and σ , ( i , i ), we infer that the (subtracted) bulk lattice height 1 and height 2 variables convergeto fields, h ( z ) and h ( z ), that form a logarithmic pair of weight ∆ = 2. Moreover if we assignthem the same normalization as their lattice companions ( A = − P ), we find the parameter ofthe logarithmic pair ( h , h ) equal to λ = − . The explicit results for σ , ( i , i ) and σ , ( i , i )[PGPR10] show that the fields h ( z ) and h ( z ) are also logarithmic partners of h ( z ) albeit ofdifferent normalizations and for different values of λ . As noted in Section 4, it means that theycan be written as linear combinations of h ( z ) and h ( z ) with known coefficients, h a ( z ) = α a h ( z ) + β a h ( z ) , α = 0 , α = 1 , α = 8 − π π − , α = − π + 42( π − , (5.23)and other values for the β a [JPR06]. These field assignments predict that the lattice correlation ofheights larger or equal to 2 behave asymptotically as σ a,b ( i , i ) ≃ − α a α b P rr + O (cid:18) log rr (cid:19) , a, b > . (5.24)Because α 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. Numericalsimulations have successfully confirmed the behaviour (5.24) [JPR06]. A lattice proof howeverremains 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 The four height fields h a ( z ) also satisfy the trivial identity h ( z ) + h ( z ) + h ( z ) + h ( z ) = 0. c = − − − σ , , ( i , i , i ) = 0 + . . . (5.25)implying that its scaling limit, corresponding to the CFT 4-point function h h ( z ) h ( z ) h ( z ) ω ( ∞ ) i ,vanishes identically.The lattice 4-point correlator has the expected dominant degree − σ , , , ( i , i , i , i ) = P ( | z z | + 1 | z z | + 1 | z z | − z z z z ) − z z z z ) − z z z z ) + c . c . ) + . . . (5.26)and is much more instructive: it is precisely the expression we obtain for the 5-point CFT cor-relation function h h ( z ) h ( z ) h ( z ) h ( z ) ω ( ∞ ) i if we assume that the height 1 field h ( z ) isa primary field of dimensions ( h, h ) = (1 ,
1) in a CFT with central charge c = −
2, and thatit satisfies a certain degeneracy condition at level 2. This last condition furnishes a differen-tial equation [DFMS97], from which the correlator can be fully determined, the result being ex-actly the function (5.26) ! From this result, one can actually infer that the previous correlator h h ( z ) h ( z ) h ( z ) ω ( ∞ ) i must vanish if it is to be symmetrical in the three insertion points[PR17].The lattice 3-point correlators σ a, , ( i , i , i ), a >
2, have been computed more recently in[PR17]. For simplicity, the three points were assumed to be aligned horizontally in the plane, withreal separations x ij . The following result was obtained to dominant order, σ a, , ( i , i , i ) = α a P x x + . . . , ( a > , (5.27)where the coefficients α a are those given in (5.23). In addition to being very simple, this expressionis surprisingly non-logarithmic. Because its scaling limit should be given by h h a ( z ) h ( z ) h ( z ) ω ( ∞ ) i , it is particularly important to understand it from the CFT point of view. Indeed thecomputation of such a correlator requires additional information on the height fields h a > as log-arithmic partners of h . Since h and h can be regarded as linear combinations of h and h , itis sufficient to consider h .Inspired by the conformal representations appearing in the bosonic sector of the symplectictheory [GK99], the following proposal has been made in [JPR06] regarding the conformal natureof h ( z, z ). A more complete account will be presented in Section 8.25he field h ( z, z ) is not primary since it transforms into h ( z, z ) under dilations. It is alsonot quasi-primary because its L and L transforms generate two new fields, ρ ( z, z ) and ρ ( z, z )respectively, with weights (0 ,
1) and (1 , ρ ( z, z ) is left primary and its L transform is equal to κ I ; likewise ρ ( z, z ) is right primary and its L transform is also equal to κ I . All this results in the following transformation law of h ( z, z ) under a general conformaltransformation z → w ( z ) and z → w ( z ), h ( z, z ) = (cid:12)(cid:12)(cid:12)(cid:12) d w d z (cid:12)(cid:12)(cid:12)(cid:12) (cid:20) h ( w, w ) + log (cid:12)(cid:12)(cid:12)(cid:12) d w d z (cid:12)(cid:12)(cid:12)(cid:12) h ( w, w ) (cid:21) + 12 (cid:18) d w d z (cid:30) d w d z (cid:19) d w d z ρ ( w, w )+ 12 d w d z (cid:18) d w d z (cid:30) d w d z (cid:19) ρ ( w, w ) + κ (cid:18) d w d z (cid:30) d w d z (cid:19)(cid:18) d w d z (cid:30) d w d z (cid:19) , κ = − P . (5.28)This conformal transformation law of h is sufficient to compute correlators involving h , butsubstantially complicates the calculations. Using this transformation and the left and right level2 degeneracy of h ( z, z ), the required correlator can be nonetheless determined. The result reads[PR17] h h ( z ) h ( z ) h ( z ) ω ( ∞ ) i = P
16 1 | z z | (cid:20) z z + 1 z z (cid:21) . (5.29)When the tree points z i are aligned horizontally, it exactly reproduces the lattice result (5.27) for a = 2. The cases a = 3 , α a since h h ( z ) h ( z ) h ( z ) ω ( ∞ ) i = 0. The calculation of occurrence probabilities of minimal subconfigurations have been briefly discussedin Section 5.2. Their correlations can be computed very much like those of heights 1 by using adefect matrix. The calculation of mixed 2-point correlators for about a dozen different minimalsubconfigurations has been reported in [MR01]. It turns out that each such cluster S can bespecified by a triplet ( a, b , b ) of real numbers.We define as before subtracted variables h S ( i ) = δ S ( i ) − P S , (5.30)where δ S ( i ) denotes the event “a minimal subconfiguration S is found around site i ”, and P S ( i ) = E h δ S ( i ) i is the probability of such an event. The mixed correlator of two such variables takes theform σ S,S ′ ( i , i ) = E h h S ( i ) h S ′ ( i ) i = − r (cid:26) aa ′ + ( b b ′ − b b ′ ) cos 4 ϕ (cid:27) + . . . (5.31)We see that the dominant contribution retains an angular dependence, which in this case is notsurprising since the minimal clusters are generally not rotationally invariant. As a matter ofillustration, the cluster reduced to a single height 1 is characterized by the triplet ( a, b , b ) =( P , ,
0) while for the second one in (5.10), one has S = 2 21 : a = (4 − π )(72 − π )(3 π − π , b = 0 , b = − (4 − π )(32 − π )(3 π − π . (5.32)26n the scaling limit, the subtracted cluster variables give rise to the fields h S ( z ), whose mixedcorrelators 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 ofSection 4. Indeed one may check that the explicit fields given by h S ( z, z ) = − (cid:26) a (cid:16) ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ (cid:17) + ( b + i b ) ∂θ ∂ ˜ θ + ( b − i b ) ∂θ ∂ ˜ θ (cid:27) , (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 explainedearlier. In particular, the height 1 field h ( z, z ) is recovered upon setting a = P and b = b = 0.Let us note a generic field h S ( z, z ) is a linear combination of three fields with different conformalweights, namely (1 , ,
0) and (0 , 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 probabilitiesfor minimal clusters can be computed by using defect matrices which implement the local bondmodifications. On the other hand, defect matrices always yield contributions that are given by finitedeterminants of discrete Green matrix entries. In the limit of large separations, the determinantsconverge to polynomial expressions in the Green function and its derivatives. It is therefore nota complete surprise that the associated fields can be constructed out from the symplectic freefermions θ, ˜ θ . Indeed, because the 2-point correlators of θ, ˜ θ are given by the Green function, thecorrelators of any fields that are local in θ, ˜ θ and their derivatives are necessarily polynomials in theGreen function and its derivatives (Wick’s theorem). We expect this observation to extend to allthe observables that correspond to local perturbation of the toppling matrix. Isolated dissipationand minimal cluster variables are among them; the arrow variables discussed in the next sectionare 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 freesymplectic fermion theory ? As explained earlier, these three variables cannot be handled withfinite rank perturbations of the toppling matrix, because they involve non-local constraints on thenearest neighbours (some of them should not be predecessors). Using the technique developed in[KW15], the 1-site probabilities P a > can be efficiently computed. Surprisingly, the details showthat the explicit values are given in terms of a few entries of the lattice Green matrix (at shortdistances), which explains why the values of P a > quoted in (5.9) are not much more complicatedthan for P . This is no longer the case for large distance correlations σ a > , ( i , i ). The analysis[PR17] shows that those correlators are expressed in terms of sums of product of Green matrixentries over a path connecting i to i , and thus in terms of quantities that are not local in theGreen matrix. It supports the view that the height fields h a > do not belong to the free symplecticfermion theory. A detailed analysis of this question has been carried out in [JPR06], and hasreached the same conclusion. Section 8 below summarizes this somewhat strange situation. The qualifier ‘conservative’ means that the defect matrix that implements the bond modifications has zero rowand column sums. The defect matrix used in Section 5.1 to compute the height 1 probability does not have thisproperty. It can however be replaced by another one that does have it [MR01]. An example of a non-conservativebond modification is given in Section 5.7. .7 Spanning tree related variables A recurrent configuration of the sandpile model can be specified as a set of height values oras a spanning tree; the former has the local heights as natural variables, the latter has localconnectivities as natural variables, namely the existence or absence of specific bonds in the spanningtree.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 uniquepath connecting any vertex to the root, and therefore any vertex to any other vertex. A rootedspanning tree can then naturally be oriented by deciding that the edges of the tree all point towardsthe root. As a consequence, in any rooted spanning tree, there is exactly one outgoing edge at eachvertex but the root; there may however be more (or less) ingoing edges (a vertex with no ingoingedge is a leaf). A site j is a then a predecessor of i if the unique path from j to i is consistentlyoriented (equivalently if the unique path form j to i does not pass through the root). As we haveseen, the question of being predecessor is a non-local problem, even if i and j are close to eachother, even nearest neighbours [PP11].Connectivities between neighbouring sites can be handled in much the same way as heights 1 orminimal height cluster variables. To see this, we must first understand why the determinant of thetoppling matrix on a graph counts the number of spanning trees on that graph. In the perspectiveof this section, we can generalize the matrix by assigning arbitrary weights to the oriented edgesof 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 ,∆ i,j = y i = x i⋆ + P j = i x ij for i = j, − x ij for i = j, (5.34)In the context of the sandpile model, the difference x i⋆ = y i − P j = i x ij can be viewed as the weightof the oriented edge from i to the sink ⋆ , so that the conservative vertices have this difference equalto 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 sumover 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,det ∆ = X σ ∈ S N ε σ ∆ ,σ (1) ∆ ,σ (2) . . . ∆ N,σ ( N ) = [ N/ X k =0 ( − k X σ has k proper cycles ∆ +1 ,σ (1) . . . ∆ + N,σ ( N ) (5.35)where the matrix ∆ + is ∆ without the minus signs in the non-diagonal part. The second equalityfollows by combining the signs in the non-diagonal entries of ∆ with the parity of σ : every cycle oflength ℓ > σ brings a sign ( − ℓ − coming from the parity ε σ and another sign( − ℓ from the product of non-diagonal entries of ∆, resulting in an overall sign − ℓ = 1, corresponding to a vertex left invariant by σ , brings no sign.The term k = 0 is simply equal to Q i y i , as the only permutation with no proper cycle is theidentity. In combinatorial terms, the product Q 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 theroot, each configuration being weighted by the product of the weights carried by the arrows. The28eneric term k = 0, apart for the sign ( − k , is a weighted sum of arrow configurations whichcontain 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 formmore loops.By using the inclusion-exclusion principle, one can see that the above alternating sum has theeffect to subtract from the term k = 0 the weights of all the arrow configurations which containat 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 atleast must have its arrow oriented to the sink (a configuration of N arrows on N vertices necessarilycontains a loop). When all weights x ij are equal to 0 or 1, det ∆ is simply the number of spanningtrees.Let us come back to the question of local connectivities on a rectangular grid in Z and computethe probability that the outgoing arrow from the site i is oriented to its right neighbour. Thisamounts to compute the fraction of spanning trees with such an arrow, namely P → ( i ) = P (right arrow at i ) = { i } det ∆ , (5.36)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 neighbourE, one could simply set to 0 the weights between i and its three neighbours S, W and N. It ishowever computationally more efficient to set the weight from i to its right neighbour to x , andtake the limit x → + ∞ so as to give the edges to the other three neighbours a relative weight equalto 0. This implies that we define a new matrix ˜∆ which coincides with ∆ except on two entries,namely ˜∆ i,i +ˆ e = − x and ˜∆ i,i = x + 3. As before we write ˜∆ = ∆ + B → ( i ) for the defect matrix B → ( i ) which is zero everywhere except on the two sites i, i + ˆ e , where it reduces to (cid:0) x −
10 1 − x (cid:1) .Since x is in any case large, we can simply set that part of B → ( i ) to (cid:0) x − x (cid:1) . From Kirchhoff’stheorem, we obtain P → ( i ) = lim x → + ∞ x det[ I + ∆ − B → ( i )] , (5.37)which reduces to a 2-by-2 determinant. In the infinite volume limit, one finds P → ( i ) = , asexpected. If we want to have the arrow at i oriented to its neighbour j other then E, we usesimilar defect matrices B ↑ , B ↓ , B ← , which look the same as B → but with the non-zero 2-by-2block (cid:0) x − x (cid:1) placed on the sites i, j . The four orientations of the arrow at i yield the same result . Multipoint arrow probabilities can be computed in the now usual way, placing appropriatedefect matrices at the different sites. For instance the probability to find a right arrow at two sites i and i is σ → , → ( i , i ) = lim x → + ∞ x det h I + ∆ − { B → ( i ) + B → ( i ) } i . (5.38)In the infinite volume limit and for a large distance, the following two-point probabilities are found29t dominant order, σ → , → ( i , i ) −
116 = 116 π ( z + z ) | z | + . . . (5.39a) σ ↑ , ↑ ( i , i ) −
116 = − π ( z − z ) | z | + . . . (5.39b) σ → , ↑ ( i , i ) −
116 = − i16 π z − z | z | + . . . (5.39c)Looking for symplectic fermion realizations of fields ρ → and ρ ↑ which reproduce these two-pointcorrelators, one quickly sees that these have to include two parts with respective weights (1 ,
0) and(0 , ρ → ( z, z ) = 12 π (cid:16) θ ∂ ˜ θ + θ ∂ ˜ θ (cid:17) , ρ ↑ ( z, z ) = i2 π (cid:16) θ ∂ ˜ θ − θ ∂ ˜ θ (cid:17) , (5.40)in agreement with the fact that ρ ↑ ( z, z ) is formally the rotated form of ρ → ( z, z ) (under z → − i z ).The first two correlators are indeed related by rotation. This observation also suggests that theother 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 oftheir orientation, can be computed. The probability that a single, fixed edge belongs to a tree isthe sum of the probabilities to find it in either of the two possible orientations, is thus equal to .Likewise the probability to find m edges in a tree is the sum of the probabilities to find themin all possible orientations, and so is the sum of 2 m probabilities of m oriented edges. That sumcan however be obtained in one go by replacing the block (cid:0) x − x (cid:1) used for the oriented edgesby (cid:0) x − x − xx (cid:1) and dividing as above the determinant by x m . Because two arrows with oppositeorientations can not occupy the same edge (as they would form a loop), the summation over the2 m terms is correctly realized.Correlations of unoriented edges should decay faster than that of oriented edges, because thesum of the two orientations is zero, in view of the relations ρ ← = − ρ → and ρ ↓ = − ρ ↑ , at least atthe order that was dominant for the oriented edges ( r − ). Indeed explicit calculations yield a r − decay, σ ↔ , ↔ ( i , i ) −
14 = σ l , l ( i , i ) −
14 = − π ( z + z ) | z | + . . . (5.41a) σ ↔ , l ( i , i ) −
14 = 116 π ( z − z ) | z | + . . . (5.41b)From these correlators, the associated fields φ ↔ and φ l must have components with conformalweights (2 , ,
1) and (0 , (cid:0) x − x − xx (cid:1) is conservative.More precisely, the following fermionic expressions reproduce the above correlations, φ ↔ ( z, z ) = 12 π (cid:16) ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ (cid:17) = ( ∂ + ∂ ) ρ → − π (cid:16) θ∂ ˜ θ + θ ∂ ˜ θ (cid:17) , (5.42a) φ l ( z, z ) = 12 π (cid:16) ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ − ∂θ ∂ ˜ θ − ∂θ ∂ ˜ θ (cid:17) = i( ∂ − ∂ ) ρ ↑ − π (cid:16) θ∂ ˜ θ + θ ∂ ˜ θ (cid:17) . (5.42b)30n 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, onewould expect that the fields describing a horizontal resp. vertical unoriented edge are proportionalto the horizontal resp. vertical derivative of the fields describing the oriented edges, namely φ ↔ ∼ ( ∂ + ∂ ) ρ → and φ l ∼ i( ∂ − ∂ ) ρ ↑ . It turns out not to be quite the case. Formulating the sandpile model on a surface with boundaries is important to see how they affectthe statistics of the model. The multisite correlations discussed in the previous section are likely tobe modified by the presence of a boundary, and by the associated boundary conditions. Moreoverthe microscopic variables on a boundary or very close to it will surely have a different behaviourfrom their bulk versions. In the field theoretic description, the boundary fields have to be properlyidentifed, and the way changes of boundary conditions are implemented must be clarified. All thisadds to the known set of bulk fields a number of boundary related fields, and offers the opportunityto further test the consistency of their identification by computing mixed correlations combiningboth types of variables.Surfaces with boundaries arise in the thermodynamic limit when some of the boundaries of thefinite 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 : y > } . There is onlyone 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 twopossible boundary conditions: open and closed. Let us recall that the boundary condition is openresp. closed if the boundary sites are dissipative resp. conservative. As before, a boundary opensite has z ⋆i = 4 while a boundary closed site has z ⋆i = z i = 3. In the scaling limit, it endows H with ahomogeneous boundary condition, open or closed. We also can (and will) consider inhomogeneousboundary conditions, by alternating open and closed stretches on a single boundary.In terms of height variables, the open boundary condition is equivalent to fix all the boundaryheights to 4, whereas the closed condition amounts to constrain the boundary heights not to takethe value 4. The fixed boundary condition with boundary heights equal to 1 is not possible (twoneighbouring 1’s form a forbidden subconfiguration); the fixed boundary conditions with boundaryheights 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 calledwindy boundary conditions will be discussed at the end of this section (as we will see, they are notso far from the possibility just mentioned, namely that of having height variables being equal to 2or 3). No other boundary condition has been considered so far though it would be very surprisingthat no other exist . Indeed the burning algorithm implies that the boundary sites, all with a height equal to 4, will burn at the firststep of the burning process. The sites on the next layer all have z ⋆i = 4, corresponding to the open condition. Of course we talk here of no other universality class of boundary conditions. Many boundary conditions maydiffer in the way they are microscopically defined on the lattice and nevertheless renormalize to the same continuum = 1 a = 2 a = 3 a = 4 c a P = π − π π − π e γ + − π π − π π e γ + − π +2 π π − π +44 π e γ + π − π π d a P = π − π − π π − π +44 π Table 1: Numerical coefficients for one-site height probabilities on the UHP, with e γ = γ + log 2.They satisfy the relations P a c a = P a d a = 0. In this section, we would like to reconsider the multisite height probabilities, but on a domain witha boundary, the upper-half plane (UHP) being the simplest case. The principle underlying thecalculations on the UHP stays the same as on the full plane. The most essential difference is thatthe toppling matrix becomes in the thermodynamic limit the Laplacian matrix on the discrete UHPwith the appropriate boundary condition, open or closed. In this section, we consider homogeneousboundary 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 | y > } . For either boundary condition, theLaplacian matrix ∆ op or ∆ cl is minus the adjacency matrix of the discrete UHP plus a diagonalmatrix, everywhere equal to 4 for the open condition, equal to 4 and 3 respectively for the bulkand boundary sites for the closed condition. By the method of images, the Green matrices G op / cl =(∆ op / cl ) − can be easily computed in terms of that on the full plane, G op( x ,y ) , ( x ,y ) = G ( x ,y ) , ( x ,y ) − G ( x ,y ) , ( x , − y ) , (6.1a) G cl( x ,y ) , ( x ,y ) = G ( x ,y ) , ( x ,y ) + G ( x ,y ) , ( x , − y ) , (6.1b)for y , y >
0. As anticipated, we verify that G op satisfies the Dirichlet condition, namely it isodd under the reflection through the line y = 0 and therefore vanishes on it, and that G cl satisfiesthe Neumann condition, namely it is even under the reflection through the line y = , inducinga vanishing normal derivative in the scaling limit. The calculations of the previous section canthen be generalized to the UHP geometry by using these Green matrices. For the height 1 andfor 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 combinatoricsinvolved is heavier.The simplest case is the 1-site height probability P op / cl1 ( y ) to find a height equal to 1 at adistance y from the boundary. It can be computed by using the formula (5.6) where ∆ − isreplaced, in the infinite volume limit, by one of the two Green matrices given above. This washistorically the first calculations of boundary effects in sandpile models [BIP93], with the result σ op1 ( y ) = P op1 ( y ) − P = P y + . . . , σ cl1 ( y ) = P cl1 ( y ) − P = − P y + . . . (6.2)The analogous results for higher heights were obtained somewhat later [PR05a, JPR06] andwere the first to firmly establish their logarithmic nature. They take the following form, valid for boundary condition in the scaling limit. > σ op a ( y ) = 1 y (cid:18) c a + d a d a log y (cid:19) + . . . , σ cl a ( y ) = − y (cid:18) c a + d a log y (cid:19) + . . . (6.3)up to terms of order O ( y − log y ), which have since been explicitly computed [PR17], as they enterthe calculations of σ op / cl a given below. The coefficients are explicitly known and are collected inTable 1. One may check that the relations (5.23) expressing h and h linearly in terms of h , h are confirmed. The distinctive change of sign between the two boundary conditions, the fact thatfor fixed a , both are controlled by the same constants c a and d a and the equality d = c arestriking. As will be explained below and in one of the next sections, all three features will followfrom the CFT picture.Let us mention that these lattice calculations have been carried out on another lattice realizationof the UHP, namely on the diagonal upper-half plane { ( x, y ) ∈ Z | y > x } , for which the methodof images allows to explicitly compute the Green matrices for the two boundary conditions. Asexpected, the dominant terms are exactly the same as above in terms of the Euclidean distancebetween 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 = ( x , y ) and i = ( x , y ),and which involve the same subtractions as before, σ op / cl a, ( x ; y , y ) = P op / cl a, ( i , i ) − P op / cl a ( i ) P − P a P op / cl1 ( i ) + P a P , (6.4)have been computed in [PR17] when the two sites are far from the boundary and far from eachother, again using the technique developped in [KW15]. They depend on three real variables,the horizontal distance x = x − x between the two sites and their vertical positions y , y . Forsimplicity however, the lattice calculations have been carried out for two vertically aligned sites,that is, for x = 0.Defining the two bivariate functions, P ( u, v ) = 18 u v − u − v ) − u + v ) , Q ( u, v ) = 1( u − v ) − u + v ) , (6.5)the results for a = 1 , σ op1 , (0; y , y ) = σ cl1 , (0; y , y ) = P P ( y , y ) + . . . , (6.6a) σ op / cl2 , (0; y , y ) = P (cid:20) P ( y , y ) (cid:18) log y + γ + 52 log 2 (cid:19) + Q ( y , y ) log (cid:12)(cid:12)(cid:12)(cid:12) y + y y − y (cid:12)(cid:12)(cid:12)(cid:12)(cid:21) + H op / cl ( y , y ) y y ( y − y ) + . . . , (6.6b)where H op ( u, v ) and H cl ( u, v ) are homogeneous polynomials of degree 4 in u, v , with explicitlyknown coefficients. The results for a = 3 , h a ( z, z ). For a homogeneous boundary condition like here, boundary CFT prescribes tocompute bulk correlators on the UHP by viewing a field φ ( z, z ) with conformal weights ( h, h ) as33he 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 2 n chiral fields on the full plane; thecorrelation appropriate for the boundary condition under consideration is accordingly selected inthe solution space of these 2 n -correlators.The above prescription must however be adapted in the case of logarithmic fields because thechiral factorization is not consistent with the non-diagonal action of L . Indeed let us consider alogarithmic 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 , L Ψ = ( L ψ h ) ψ h = ( hψ h + λ φ h ) ψ h = h Ψ + λ φ h ψ h , (6.7)that the chiral factorization of the primary partner is Φ = φ h ψ h . The same argument with L 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 h a ( z, z ) i op . Note that, unlike the correlations on the plane discussedin Section 5, we do not insert dissipation at infinity since the whole boundary is open and thereforedissipative. From the prescription recalled above, these 1-point functions should have the generalform of chiral 2-point functions. If ψ and φ denote chiral versions of the height 2 and height 1fields respectively, with h = h = 1, the chiral factorizations of the height fields h and h read h ( z, z ) = ψ ( z ) φ ( z ) and h ( z, z ) = ψ ( z ) ψ ( z ). The CFT formalism gives the general forms (4.14), h h ( z, z ) i op = h φ ( z ) ψ ( z ) i = B ( z − z ) , h h ( z, z ) i op = h ψ ( z ) ψ ( z ) i = C − λB log ( z − z )( z − z ) . (6.8)With the value λ = − noted in Section 5, these forms exactly reproduce the lattice results (6.3),including the relation d = c = B .For σ cl a ( y ) and since the closed boundary is not dissipative, we insert by hand dissipation atinfinity, so that σ op a ( y ) should be given by h h a ( z, z ) ω ( ∞ ) i cl . Using the same chiral factorization asabove leads to a 3-point function. The selection of a physically sensible solution leads to the samegeneral form as for the open boundary condition [JPR06].The conformal calculations required to account for σ op a, are only technically more involved.The needed chiral correlators are h φ ( z ) ψ ( z ) φ ( z ) ψ ( z ) i for a = 1 and h ψ ( z ) ψ ( z ) φ ( z ) ψ ( z ) i for a = 2. Both can be computed from the primary nature of the chiral field φ , as established inSection 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 ofthe two heights, z = ( x , y ) and z = ( x , y ), h h ( z , z ) h ( z , z ) i op = P (cid:26) z − z ) ( z − z ) − | z − z | − | z − z | (cid:27) , (6.9)and h h ( z , z ) h ( z , z ) i op = P y y t − t + 4 t − − t ) π − π − P (cid:18) log y + γ + 52 log 2 (cid:19) + P y y t ( t − − t ) (cid:18) log (1 − t ) + y y (cid:19) − t − t , t = − y y | z − z | . (6.10)34ne can check that setting x = x exactly reproduces the lattice results σ op1 , (0; y , y ) and σ op2 , (0; y , y ) reported above.The analogous calculation for the closed boundary has been carried in the case a = 1, yieldingthe 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 ω mustbe added).Similar calculations with isolated bulk dissipation instead of height variables have been con-sidered in [PR04]; it was found in all cases that the CFT predictions compare successfully withthe lattice results. 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 boundarycondition, both on the lattice and in the emerging field theory. Like in the previous section, weconsider 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 anopen boundary site has ∆ i,i = z ⋆i = 4 and a closed one has ∆ i,i = z ⋆i = 3, we simply lower by 1the diagonal entry ∆ i,i to close an open site, and we increase it by 1 to open a closed site (as wedid in Section 5.4 to introduce dissipation at bulk sites). We do it either way for n consecutiveboundary sites to change the boundary condition on a interval I of length n , that is, we do thefollowing change on the toppling matrix ∆ → ∆ ± D I , where D I implements the diagonal shiftsdescribed above.Let us examine the effect of closing n consecutive sites in an otherwise open boundary. Wedecide to measure this effect as in Section 5.4, namely by comparing the number of recurrentconfigurations 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 correspondingtoppling 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 theinfinite volume limit of the ratio as Z op ( n ) Z op = det ∆ op ( n )det ∆ op = det[∆ op − D I n ]det ∆ op = det[ I − G op D I n ] = det[ I − G op ] i,j ∈ I n , (6.11)where ( D I n ) i,j = δ i,j for i, j ∈ I n , is zero elsewhere, and I n = { ( ℓ,
1) : 1 ℓ n } is the set ofsites being closed. Using the relation (6.1a) expressing G op in terms of the Green matrix on thefull plane Z , the matrix in the determinant reads (cid:18) I − G op (cid:19) i,j ∈ I = (cid:18) δ ℓ,ℓ ′ − G ( ℓ, , ( ℓ ′ , + G ( ℓ, , ( ℓ ′ , − (cid:19) ℓ,ℓ ′ n . (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 theFourier coefficients of the following symbol, which has a so-called Fisher-Hartwig singularity, σ op ( k ) = √ − cos k · n √ − cos k − √ − cos k o . (6.13)35or large n , the asymptotics of such Toeplitz determinants is well-known (see f.i. [DIK13]), andleads to the following result [Ru02], Z op ( n ) Z op ≃ A n / e − π n , n ≫ , (6.14)with G= 0 . ... the Calatan constant. The proportionality constant A is explicitly known butis unimportant here.What if we consider the opposite situation, in which we open n consecutive sites of a closedboundary ? Reasoning as above, we quickly get the corresponding ratio, Z cl ( n ) Z cl = det ∆ cl ( n )det ∆ cl = det[∆ cl + D I n ]det ∆ cl = det[ I + G cl ] i,j ∈ I n , (6.15)which is also a Toeplitz determinant. However this one is infinite –each entry is infinite– for thesame reason we have pointed out in Section 5.4. Adopting the same point of view, we similarlyevaluate the effect of opening n sites with respect to the situation where only one site is open. Onetherefore considers instead the ratio Z cl ( n ) Z cl (1) , which one can write as Z cl ( n ) Z cl (1) = 1 b det (cid:16) b ℓ − ℓ ′ (cid:17) ℓ,ℓ ′ n , (6.16)where the entries b m are the Fourier coefficients of a symbol σ cl ( k ) given by σ cl ( k ) = 12 (1 − cos k ) α · n √ − cos k + √ − cos k o , α = − . (6.17)Its Fourier coefficients are well-defined for α > − , diverge in the limit α → − but nonethelesskeep the ratio (6.16) finite. Remarkably, for large n , it takes the form [Ru02] Z cl ( n ) Z cl (1) ≃ A n / e π ( n − , n ≫ , (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 π N 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 Nf b in case the boundary condition b is the same at all boundarysites. For the partition functions considered above, the boundary conditions only differ on thelower edge of the grid, so that for large N ≫ n ≫ Z op ( n ) Z op ≃ e ( N − n ) f op + nf cl e Nf op = e − n ( f op − f cl ) , Z cl ( n ) Z cl (1) ≃ e ( N − n ) f cl + nf op e ( N − f cl + f op = e ( n − f op − f cl ) . (6.19)The free energies f op and f cl represent (the logarithm of) the effective number of values taken bythe boundary heights in the set of recurrent configurations. The number of possible values takenby the height at an open boundary site is 4, and is 3 at a closed site. If these numbers values geteffectively reduced in the set of recurrent configurations, one should expect that the number of36alues at an open boundary site remains larger than that at a closed site, implying f op − f cl > f op − f cl = π , in agreement with the aboveresults. To fix the ideas, the effective number of values taken by a boundary height is e f op = 3 . 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, isimplemented by the insertion in the correlators of a specific field φ a , b ( x ). Such boundary conditionchanging fields are usually expected to be chiral primary fields, and satisfy φ a , b ( x ) = φ b , a ( x )when the boundary conditions a and b do not carry an intrinsic orientation (see Section 6.4 forcounterexamples). The insertion of the product φ a , b ( x ) φ b , a ( x ) accounts for the change at x fromcondition a to condition b, and then back from b to a at x , but does not account for the exponentialterms related to the difference of boundary free energies of condition a versus condition b, namelythe 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 underlyingCFT, which itself applies to all the models in the universality class to which the sandpile modelbelongs.It follows that the effect of changing the boundary condition given above, in which we omit theexponential terms, should correspond to the 2-point function h φ op , cl (0) φ cl , op ( n ) i = h φ cl , op (0) φ op , cl ( n ) i .The two are indeed equal on the lattice and asymptotic to n , and from this, we infer that theboundary condition changing field φ op , cl ( x ) = φ cl , op ( x ) is a chiral conformal field of weight h = − ,with a correlator given by h φ op , cl ( x ) φ cl , op ( x ) i op = A | x − x | / . (6.20)For physical reasons, we might worry about having a correlator that actually increases with thedistance, suggesting somehow the existence of a strange interaction that would get stronger atlarger distances. There is nothing strange however, as it does not really correspond to the physicalcorrelation 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 inthe lattice properties of the model.One of the strongest consequence of the primary nature of φ op , cl and the assumed structureof the conformal module that contains it, is that any correlator where this field appears mustsatisfy a second-order partial differential equation , the precise form of which depends on theother fields involved. A first and simple test is to look at a 4-point function , for instance h φ op , cl ( x ) φ cl , op ( x ) φ op , cl ( x ) φ cl , op ( x ) i , which should describe the effect of closing the sites ontwo disjoint interval [ x , x ] and [ x , x ] in the otherwise open boundary of the UHP. Using theglobal conformal invariance, one can reduce the partial differential equation to a second-orderordinary differential equation. In the two-dimensional solution space, we select the only solution In the correspondence between statistical system and field theory, the boundary condition changing fields aresomehow special. They describe the effects of a change of boundary condition but are not associated to a latticeobservable, unlike the height fields h a for instance. Again the technical assumption is that the field φ op , cl is degenerate at level 2, similarly to the height 1 field h ,see Section 5.5. The CFT is really defined on the UHP plus the point at infinity. The boundary must therefore be thoughtof as the real line plus the two points ±∞ identified, and forming a loop closing at infinity. Any change ofboundary condition thus involves an even number of insertions of φ op , cl . For instance h φ op , cl (0) φ cl , op ( ∞ ) i changesthe boundary condition from open to closed on the positive real axis. h φ op , cl ( x ) φ cl , op ( x ) ih φ op , cl ( x ) φ cl , op ( x ) i when the two intervals areinfinitely distant. This unique solution solution reads (with x ij = x i − x j ) h φ op , cl ( x ) φ cl , op ( x ) φ op , cl ( x ) φ cl , op ( x ) i op = 2 A π ( x x ) / (1 − t ) / K ( t ) , t ≡ x x x x , (6.21)where K ( t ) = ´ π θ √ − t sin θ is the complete elliptic integral.To compare with a lattice calculation, we take x i integers, with x , x and x all large, andtry to compute the determinant in (6.11) with I = I ∪ I the union of the two intervals [ x , x ] and[ x , x ]. This determinant is no longer Toeplitz, which makes it difficult to compute its asymptoticsanalytically. Dividing it by the prefactor ( x x ) / , it can however be evaluated numerically as afunction of t by varying the lengths of the intervals and their separation distance. The agreementwith (6.21) is more than satisfactory [Ru02].The opposite situation –two open intervals in a closed boundary– has also been considered. Theappropriate 4-point function can be obtained from (6.21) by making a simple cyclic permutation( x , x , x , x ) → ( x , x , x , x ), with the result that K ( x ) gets replaced by K (1 − x ). An equallysuccessful agreement was observed [PR04]. Many other crosschecks have been done, confirmingthat the open/closed boundary condition changing field is indeed a primary field with conformalweight h = − . One of them, particularly convincing, is presented in the next section. In Section 6.1, we have computed the lattice 1- and 2-site height probabilities on the UHP, witheither the open or the closed boundary condition. Here we would like to revisit these results inthe light of what we have learned of the boundary condition changing field, in terms of which oneshould be able to relate the probabilities for the two boundary conditions. In particular, we wouldlike to understand the 1-site probabilities σ op a ( y ) and σ cl a ( y ), σ op a ( y ) = 1 y (cid:18) c a + d a d a log y (cid:19) + . . . , σ cl a ( y ) = − y (cid:18) c a + d a log y (cid:19) + . . . (6.22)One can do this by computing, on the CFT side, a more general probability. Namely we look forthe probability to find a height equal to a at a distance y from the boundary, when the boundarycondition is mixed, namely open everywhere except on the interval [ x , x ] where the condition isclosed. The two homogeneous open and closed conditions can be recovered in the limits x → x and x → −∞ , x → ∞ . Let us denote by h h a ( z, z ) i mix the corresponding quantity in the CFT,given by h h a ( z, z ) i mix = h φ op , cl ( x ) φ cl , op ( x ) h a ( z, z ) i op h φ op , cl ( x ) φ cl , op ( x ) i op , (6.23)where the division by h φ op , cl ( x ) φ cl , op ( x ) i op comes from the fact that we want to evaluate theprobability to have a height 1 in front of a mixed boundary condition, and not the combinedeffects of having a height 1 and the closing the boundary between x and x . The denominator isknown from (6.20).To compute the numerator, we represent the height fields in terms of the chiral fields as h ( z, z ) = φ ( z ) ψ ( z ) and h ( z, z ) = ψ ( z ) ψ ( z ) (as usually, considering the heights a = 1 , φ op , cl . Because ψ is the chiral logarithmic part-ner of φ , the general solution for h φ op , cl ( x ) φ cl , op ( x ) ψ ( z ) ψ ( z ) i op in fact depends on that of h φ op , cl ( x ) φ cl , op ( x ) φ ( z ) ψ ( z ) i op .All calculations done, one finds that they depend on two integration constants c and d insuch a way that the ratios (6.23) take the following forms where y = Re z [JPR06], h h ( z, z ) i mix = d y t √ t , t = ( x − z )( x − z )( x − z )( x − z ) , (6.24a) h h ( z, z ) i mix = 12 y t √ t ( c + d √ t ) √ t + d log y − d i y (1 − t )( x − x ) (cid:20)
11 + t − √ t (cid:21)) . (6.24b)Although t is complex, both expressions are real on account of t ∗ = 1 /t .Let us now discuss the above two limits x → x and x → −∞ , x → ∞ . For convenience,we set x = − x and examine the limits x → + and x → + ∞ . To compute the two limits, theimportant thing to notice is that the complex variable t , now equal to t = ( x + z )( x − z )( x + z )( x − z ) , (6.25)has complex norm equal to 1 and loops anticlockwise around the origin as x varies from 0 + to+ ∞ , starting from 1 + 0i to 1 − t itself goes to 1 in both limits but √ t goes to+1 when x → + and goes to − x → + ∞ . The actual limits yield h h ( z, z ) i op = lim x → + h h ( z, z ) i mix = d y , h h ( z, z ) i cl = lim x → + ∞ h h ( z, z ) i mix = − d y , (6.26a) h h ( z, z ) i op = 1 y (cid:18) c + d d log y (cid:19) , h h ( z, z ) i cl = − y (cid:18) c + d log y (cid:19) , (6.26b)in complete agreement with the lattice results: the change of the overall sign between the openand closed boundary conditions, the specific dependence on the two coefficients c and d , andthe equality c = d are all accounted for ! The conformal approach however cannot fix the twocoefficients c and d ; these must be determined by lattice calculations.The expressions (6.24) can also be tested in situations where the boundary condition along thereal axis is no longer homogeneous. A particularly instructive case is when the boundary conditionis closed on the negative part of the real axis, and open on the positive part, corresponding to thelimits x → −∞ and x →
0. The conformal transformation w = Lπ log z can be used to map theUHP onto an infinite strip of width L , with open boundary condition on the left side, closed on theright 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] thatthe conformal predictions and the numerical plots match remarkably well, thereby confirming oncemore all the field identifications made so far. 39 .4 Wind on the boundary The open and closed boundary conditions are very natural as the very definition of the sandpilemodel uses dissipative and conservative sites. One may wonder what other type of boundarycondition could be thought of. Perhaps we could think of alternating open and closed boundarysites; we expect however that such a boundary condition would flow to the open condition in thescaling limit, as numerical experiments confirm. We have already commented on the possibilityto uniformly fix the boundary heights. Fixing the boundary heights to 2 or to 3, or even to 2 or3, seems difficult. The two boundary conditions, different from open and closed, which have beenconsidered 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 boundarysites to be uniformly left or uniformly right . In terms of height values, either condition meansthat 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, towardsthe root). The converse is however not true: recurrent configurations with height values equal to2 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. Thisallows to evaluate the effects of inserting a stretch of left or right arrows into an open or a closedboundary, similarly to what we did in Section 6.2. We refer the reader to [Ru07] for the details ofthe analysis, and restrict here to a summary of the results.An obvious but unusual feature of the boundary conditions ← and → is that they are intrins-ically oriented. It implies that the boundary condition changing field φ a , → turning the boundarycondition 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 arejust the identity). We already know φ op , cl = φ cl , op , and likewise, if a ∈ { op , cl } is unoriented andb ∈ {← , →} is oriented, we expect the identifications φ a , → = φ ← , a and φ a , ← = φ → , a on the basis ofa left-right reflection symmetry. These identifications have been confirmed and reduce the numberof distinct fields to seven.There is an additional subtlety for the field that changes the orientation from → to ← . Indeedthe right and left arrows ◦→ • ←◦ point to the same boundary site (in black), and whether thatsite is open or closed may be relevant. Indeed if it is open, the flow of arrows, which eventuallyterminates at the root, can go directly to the root; if it is closed, it must necessarily go upwardsinto 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 confirmthat they are distinct fields as their conformal weights are different.We therefore have eight distinct boundary condition changing fields. A mix of analyticalcalculations and numerical simulations has been used to determine the conformal weights of theseeight fields. The results are given in Table 2.The more delicate question of the exact nature of all these fields has been addressed by con-sidering the fusion of the representations to which they belong. Loosely speaking, the fusion rules In contrast, the outgoing arrow of a closed boundary site of the UHP can point left, up or right, while that ofan open boundary site can point in any of the four directions, a down arrow pointing to the root. [ φ a , b ] open closed → ← open − − −
18 38 → ← − φ a , b which implement a change of boundary conditionfrom a (row label) to b (column).implement the composition law φ a , b ⋆ φ b , c ≃ φ a , c of boundary condition changing fields in the limitwhere the insertion points coincide. The ensuing consistency conditions suggest that all of themare primary fields, except two, which could belong to logarithmic representations (i.e. reducibleindecomposable with Jordan cells). Also the fields of weight 0 are non-trivial, that is, not equalto the identity (they are found to be degenerate at level 3). Relying on these proposals, various4-point correlators have been computed and successfully compared with numerical simulations.We refer to [Ru07] for more details on these specific points. The boundary condition changing fields are not the only ones to live on a boundary. The latticemodel includes observables in the bulk as well as on the boundaries. Those in the bulk have beendiscussed at length and give rise in the scaling limit to non-chiral fields Φ( z, z ), characterized bya 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 boundaryfield associated wih a boundary observable and its conformal weight depend on the boundarycondition.In the Abelian sandpile model, only the boundary fields arising from the height variables andfrom the insertion of isolated dissipation have been studied on the UHP. In both case, only openand closed boundaries have been considered.The case of isolated dissipation is simpler and has been examined in details in [PR04], whereisolated dissipation has been considered on a closed boundary only. The calculation proceeds muchlike those for the bulk, reviewed in Section 5.4, for which the same regularization is used. Theresults 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 variouscombinations of logarithms like for its bulk version. On an open boundary, already dissipative, thedissipation field ω op ( x ) is expected to be a descendant of the identity. Isolated dissipation is thesimplest 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 can41e 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 proportionalto the boundary stress-energy tensor T ( x ) of the symplectic theory, a non-primary chiral field ofweight h op = 2, and a descendant of the identity since T ( x ) ∼ ( L − I )( x ).Boundary height variables are more complicated than dissipation but simpler than the bulkheight variables. The first results have been derived by Ivashkevich in [Iv94], where the one-and two-site height probabilities on open and closed boundaries were obtained. The probabilitiesinvolving heights 1 only are no more complicated than in the bulk and can be easily obtained byusing 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 termsof the number of predecessors among their nearest neighbours. So it leads essentially to the sameproblems of computing non-local contributions. Both in the bulk and on a boundary, one canwrite 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 tocalculate probabilities. In the bulk, the linear system is underdetermined and cannot be invertedto provide the required non-local contributions, and then the probabilities themselves. The mainobservation made in [Iv94] was that on a boundary, the linear system can be inverted, and thereforeallows 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 werecomputed exactly. For comparison purposes, we reproduce here their numerical values (the exactvalues can be found in [Iv94]) and recall those in the bulk, as given in Section 5.1, P = 0 .
073 63 , P = 0 .
173 90 , P = 0 .
306 29 , P = 0 .
446 17 , (6.27) P op1 = 0 .
103 82 , P op2 = 0 .
216 57 , P op3 = 0 .
316 23 , P op4 = 0 .
363 38 , (6.28) P cl1 = 0 .
113 38 , P cl2 = 0 .
216 571 , P cl3 = 0 .
316 225 . (6.29)On the open boundary, for which the comparison makes more sense, lower heights are thus morelikely.Mixed 2-site correlators τ op a,b ( x , x ) and τ cl a,b ( x , x ) on an open or closed boundary were alsocomputed in [Iv94]; all of them were found to decay like | x − x | − . Although logarithmic conformalfield 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 abovethat 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, whateverthe boundary condition, have a conformal dimension equal to h op = h cl = 2. But like for the otherobservables 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 contri-butions using defect matrices, it suggests again to look for field candidates constructed out fromthe symplectic free fermions θ, ˜ θ . This was done independently in [Je05b] and in [PR05b], followinghowever two different approaches: the former computed various 3-point correlators whereas thelatter considered 2-point correlators only but in the massive extension of the sandpile model (seeSection 7.1). The massive extension indeed allows to distinguish more efficiently different fieldswhich would otherwise have the same 2-point correlators in the non-massive (critical) limit.42he results are as follows. The four height fields on an open boundary are all proportional toa single field, h op a ( x ) = O a ∂θ ∂ ˜ θ, a , (6.30)with explicit normalization constants O a and where the θ, ˜ θ fields satisfy the Dirichlet boundarycondition. Thus on an open boundary, the four height fields and the dissipation field turn out tobe all proportional to each other. On a closed boundary, the three height fields are distinct andgiven by h cl1 ( x ) = C ∂θ ∂ ˜ θ, h cl2 ( x ) = C ∂θ ∂ ˜ θ + 12 π θ ∂∂ ˜ θ, h cl3 ( x ) = C ∂θ ∂ ˜ θ − π θ ∂∂ ˜ θ, (6.31)where the θ, ˜ θ fields now satisfy the Neumann boundary condition. In both cases, the boundarycondition means that the correlators are computed using the Wick theorem with the Wick con-tractions 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 theirbulk version.The question of the nature of the height fields on the windy boundary conditions discussed inthe previous section is definitely interesting, but has not been considered so far. This long section on boundaries has been largely devoted to a discussion of the open and closedboundary conditions, the best known and most studied ones. To finish, it is worth pointing outthat a duality exists between these two boundary conditions, which has not been fully investigatednor exploited. This duality follows from a duality relation for planar graphs, well-known in graphtheory, and acquires in the framework of the sandpile model an interesting flavour. It has beenconsidered and discussed in [IPRH05, IPR07] in the dimer model, intrinsically related, like thesandpile model, to spanning trees.Let us consider a rectangular portion of Z , that is, the graph Γ made of a rectangular arrayof vertices, in which two adjacent vertices are linked by a single edge. The boundary conditionschosen for the boundary vertices determine the extended graph Γ ⋆ , obtained from Γ by addingthe sink vertex and the edges connecting the open boundary sites to the sink. The graph Γ ⋆ corresponding to a 3 × ⋆ is embedded in the plane , the faces of Γ ⋆ are the connected componentsof its complementary in the plane (for a finite graph, there is thus a large outer face, encircling thegraph). The definition of the dual graph (Γ ⋆ ) ∗ is standard: the vertices of (Γ ⋆ ) ∗ are associated tothe faces of Γ ⋆ , and two such vertices are connected if their corresponding two faces are separatedfrom 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 ex-changed: if a boundary is homogeneously open resp. closed in Γ ⋆ , it becomes homogeneouslyclosed resp. open in (Γ ⋆ ) ∗ . In addition, the dual graph (Γ ⋆ ) ∗ is the extension (Γ ∗ ) ⋆ by a sink of adual rectangular grid Γ ∗ , of size slightly different from the original grid Γ. This requires Γ ⋆ to be planar, and therefore excludes that some of the bulk vertices and some of the boundaryvertices be open (dissipative) at the same time, except in a few very special cases. × ⋆ is equal to the number ofspanning trees on its dual (Γ ⋆ ) ∗ . In fact, for every spanning tree T on Γ ⋆ , there is a unique dualspanning tree T ∗ on (Γ ⋆ ) ∗ such that the two are perfectly interdigitating: the edges of T ∗ areexactly those of (Γ ⋆ ) ∗ which cross the edges of Γ ⋆ not used in T , and vice-versa. An example ofthis 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 equalityof their partition functions can be explicitly written for rectangular grids. If Γ is an L × L rectangular grid with k of its four boundaries being open, the other 4 − k being closed, the dualΓ ∗ is an L ′ × L ′ rectangular grid with swapped boundary conditions, and the following identityholds, Z [ k op | (4 − k ) cl] ( L , L ) = Z [(4 − k ) op | k cl] ( L ′ , L ′ ) . (6.32)The dimensions are related as follows: L ′ i = L i + 1 resp L i − L i of Γ are both open resp. closed, and L ′ i = L i otherwise. If k = 4, the dual rectangle has all itsboundaries closed with a single boundary site open.The isomorphism of the two descriptions may be hard to formulate in concrete terms for theheight variables as it is defined for the associated trees. Its practical utility remains to be seen.44 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.
In the standard sandpile we have studied so far, the sites in the bulk of grid, that is, the vastmajority of sites, are conservative. This meant that when such a site topples, it loses a certainnumber of sand grains which are all redistributed to its nearest neighbours. Sand moves in thegrid but remains conserved. Dissipative sites must be present for the dynamics of the model to bewell-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 macroscopicparts of it, inducing strong correlations between distant heights. In the long run, the system entersa critical state described by the invariant measure P , characterized by infinite correlation lengthsin the infinite volume limit, and algebraic decays of the correlation functions. The field theoryemerging 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 tointroduce a fair amount of dissipation so as to make the range of the avalanches shorter. It is notcompletely clear what a fair amount means, as there are several ways to introduce dissipation. Inthe most common version, every site is made dissipative, with a dissipation rate controlled by anexternal parameter. In this case, it has been argued that indeed criticality is broken, resulting inan exponential decay of the correlations [GLJ97, TK00, MR01]. A mathematically rigorous proofthat all correlations decay exponentially has been provided in [MRS04]. Presumably a non-zerodensity of dissipative sites could be a sufficient to break criticality, but to our knowledge, thispossibility has not been investigated. In any case, the field theory emerging from the dissipativesandpile model must be massive, with mass(es) inversely proportional to the lattice correlationlength(s).To make all sites dissipative, one can simply add to the toppling matrix of the standard modelan integer multiple of the identity matrix, ∆ → ∆( t ) = ∆ + t I with t an integer, while leavingall non-diagonal entries unchanged. According to the update of the heights after the toppling ofsite j , namely h i → h i − ∆ j,i ( t ), a toppled site loses t sand grains more than what it used tolose (whether or not the toppled site is on a boundary). That this change makes the correlationfunctions 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 inverseLaplacian ∆ − ( t ) has a kernel given at large distances by G i ,i ( t ) ≃ π K (cid:16) | i − i |√ t (cid:17) + . . . , anddecays exponentially like e − r √ t at large distances ( K is the modified Bessel function). Thus allmultisite probabilities examined in the earlier sections, for observables like minimal cluster vari-ables, 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 thereason explained above: there is a loss of sand each time a site topples, which makes the typicalavalanches 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 dissipative45odel. The technique explained in Section 5.3 in terms of defect matrices goes through. Atdominant order, the result, which is the off-critical extension of (5.14), reads [MR01] σ , ( i , i ; t ) = − t P ( K ′′ − K K ′′ + 1 π K ′ + 1 + π π K ) + . . . (7.1)where the argument of the Bessel functions is r √ t and P on the r.h.s. is the critical probability;the dots stand for higher orders in t . We see that the correlation decays exponentially, with acorrelation length proportional to ξ ∼ t − / .How do we compute the scaling limit in the massive model ? The general discussion in Section 3suggested that setting i = ~xε in the lattice corrrelator and taking the limit over ε (after multiplyingthe correlator by a suitable power of ε ) yields the field theoretic correlator. This cannot be theright way to proceed in the dissipative model. Because the correlators decay exponentially, thelimit for ε going to zero of exp ( −| ~x − ~x |√ t/ε ) vanishes whathever the power of ε it is multipliedby. The only way to get a non-trivial limit is to take a double limit: as we take the large distancelimit by setting i = ~xε , we simultaneously take the large correlation length limit by accordinglyadjusting the dissipation rate. In the present case, we should take the latter proportional to ε :we therefore set t = M ε , with M playing the role of a mass (inverse correlation length in thecontinuum field theory).Looking at the lattice correlator (7.1), we see that the factor t carries the overall dimension ofthe fields involved: t is proportional to M , and thus inversely proportional to a distance to thefourth power. It replaces the explicit dependence in r in the non-dissipative model. Eventuallywe find that the scaling limit of the correlator (7.1) islim ε → ε − σ , (cid:16) zε , wε ; ε M (cid:17) = − M P ( K ′′ − K K ′′ + 1 π K ′ + 1 + π π K ) , (7.2)where the argument of the Bessel function is now M | z − w | . It is straightforward to check thatthe M → − P / | z − w | , obtained in Section 5.The last question is: the expression above is the correlator of what field and in what fieldtheory ? The most obvious guess turns out to be correct: let us look in the massive extension ofthe free symplectic fermion theory. It contains the same two fields as before, which simply acquirea mass through a mass term in the action, S = 1 π ˆ d z d z (cid:16) ∂θ∂ ˜ θ + M θ ˜ θ (cid:17) . (7.3)The 2-point correlators of the two fundamental fields are now given by h θ ( z, z ) θ ( w, w ) i = h ˜ θ ( z, z )˜ θ ( w, w ) i = 0 , h θ ( z, z )˜ θ ( w, w ) i = K ( M | z − w | ) . (7.4)Using Wick’s theorem, it is a simple matter to check that the following local field, h ( z, z ; M ) = − P (cid:20) ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ + M π θ ˜ θ (cid:21) , (7.5)46as a 2-point correlator in the massive fermionic theory that is precisely given by (7.2). The 3-and 4-point correlators of the same field have been checked to reproduce the corresponding latticeresults. The field h ( z, z ; M ) is therefore what the height 1 variable in the dissipative sandpilemodel converges to in the scaling limit.Similar correlators have been computed for many minimal clusters in [MR01], with an unex-pectedly simple result. The field describing the minimal cluster variable S in the dissipative modelappears to be simply given by h S ( z, z ; M ) = h S ( z, z ) − P S N S M π θ ˜ θ, (7.6)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 boundarieshave been explicitly evaluated in the dissipative model [PR05b]. For them too, it is found that theboundary fields given in Section 6.5 get additional terms proportional to M θ ˜ θ .The nature of the higher height fields remains elusive but is definitely worth investigating as itwould add a most valuable and crucial element of understanding of the sandpile model. Universality is the statement that the large distance properties of statistical models should onlydepend on some gross features of the way they are defined; microscopic details which becomeinvisible 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. Insandpile models, these would include the precise way sand is deterministically redistributed amongneighbours (provided some form of isotropy is preserved), or, to a certain extent, the specific graphor lattice on which the model is defined. Features that do matter are a substantial introductionof dissipation, as we have seen in the previous section, a directed redistribution of sand aftertoppling [DR89], a dynamics with stochastic toppling rules [Ma91], the formulation of the modelon a hierarchical geometric structure like the Bethe lattice [DM90], and of course a change ofdimensionality of the underlying lattice.Very early on, universality with respect to the planar lattice on which the sandpile is beingformulated has been tested via a renormalization group approach [PP97, LH02] and numericalsimulations [HL03]. More recently, exact calculations of height correlations have been carried outon the honeycomb and triangular lattices.In [ADMR10], all calculations of height 1 correlations presented in the previous sections havebeen 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 closedboundaries, 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 squarelattice, the dominant terms are exactly identical, up to normalizations. The same distinctivefeatures are found, like the fact that the 3-site bulk correlation vanishes in the scaling limit (the The insertion by hand of the dissipation field ω ( ∞ ) at infinity in the field theoretic correlator is not requiredin the dissipative model, as dissipation is present everywhere in the bulk. − − P = 112 , P = 724 , P = 58 , (7.7)while those on the infinite triangular lattice are somewhat more complicated, like P = 1175864 − √ π − π + 30 √ π + 45 π − √ π ≃ . , (7.8)and very similar expressions for P a .Concerning the nature of the height variables in the scaling limit, the results confirm whatthe reader has probably already suspected: far from boundaries, the height 1 variable becomes aprimary field with conformal weights ( h, h ) = (1 , exhibit the same bulk andboundary behaviours as on the square lattice. Thus for what concerns the type of the underlyinglattice, universality has been explicitly and successfully verified. 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 someof the most peculiar and not so well understood issues. We will almost exclusively discuss thenon-chiral bulk fields, but before coming to those, we briefly comment on the chiral boundaryfields encountered so far.The boundary fields have been somewhat less investigated than the bulk fields. We have en-countered two types of boundary fields, those arising from boundary observables and the boundarycondition changing fields. In the first class, we have considered the height fields on open and closedboundaries and the dissipation field. Except for the dissipation on a closed boundary, none of themis logarithmic and no evidence of a logarithmic partner has been found. All can be expressed aslocal fields in the symplectic fermions. Some of the calculations done on the square lattice could not be worked out. For instance, we could not finda proper method of images to compute the Green matrix on the triangular half-plane with the closed boundarycondition, and therefore could not investigate that boundary condition.
48n the second class, we found primary fields of weight − and , which are both standard fieldsin a c = − in the symplectic fermion theory [GK99]. The status ofthe 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 fieldswhich fit well within the symplectic theory. From this point of view the bulk fields are somehowmore intriguing.Most of the bulk fields we have encountered seem to have a realization in terms of symplecticfermions, by which we mean that the fermionic expressions reproduce the known correlators. Afew have not been realized in this way so far, namely the height variables h a > not equal to 1,logarithmic partners of the height 1 field h , as well as the two fields ρ and ρ , to which theytransform under L and L respectively.Although we have not given any physical interpretation of ρ and ρ , they appear to be relatedto the derivatives of the dissipation field ω [PR17], ρ = δ L − ω, ρ = δ L − ω, (8.1)where δ is a constant which may depend on the lattice considered, and equal to δ = π P on thesquare lattice. In addition, the primary field h may be consistently identified as being proportionalto the derivatives of ρ and ρ , L − ρ = L − ρ = βλh , β = 12 , (8.2)where λ is defined from L h = L h = h + λh and depends on the normalizations of h , h .Combining these relations with the previous ones yields the somewhat surprising result that theheight 1 field is proportional to the Laplacian of the dissipation field, h ∼ ∂∂ω . The correl-ator (5.19) confirms this: applying ∂ ∂ ∂ ∂ on it indeed yields a multiple of 1 / | z − z | , itselfproportional to h h ( z , z ) h ( z , z ) ω ( ∞ ) i .From these observations, it follows that all bulk fields encountered so far, namely h a> , h , ρ, ρ, ρ → , ρ ↑ , φ S , φ ↔ , φ l , ω, I , (8.3)belong to the same conformal representation, as they are all related to each other by the action ofVirasoro modes L n or L n . Indeed ρ → and ρ ↑ are not quasi-primary and transform to a multiple of I under L or L , while φ S , φ ↔ and φ l are a linear combinations of h and the chiral and antichiralstress-energy tensors T and T . In fact in terms of fermions, all these fields, except h a> , areproportional to or are linear combinations of I , θ ˜ θ , θ∂ ˜ θ , θ∂ ˜ θ , ∂θ∂ ˜ θ , ∂θ∂ ˜ θ , ∂θ∂ ˜ θ and ∂θ∂ ˜ θ . Clearlythe main question is: do the fields h a> also have a realization in terms of symplectic fermions ?In the symplectic fermion theory, the conformal representation which contains the fields quotedabove is larger because it contains many other fields, like θ∂θ or ∂θ∂θ , which have not yet beenfound in the sandpile model. Among other peculiarities, the fermionic theory also contains fourlogarithmic pairs ( φ αβ , ψ αβ ) of weight (1 , φ αβ = ∂θ α ∂θ β for the primary fields, and Very much like the spin field of the Ising model belongs naturally to the free Majorana fermionc theory with c = , 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 firststudied in [Ro96] in their chiral version. As far as we know, it has been first noticed in [GK96] forthe case of modules containing rank 2 Jordan blocks, that these representations are characterizedby an intrinsic complex parameter β , known as a logarithmic coupling, an indecomposabilityparameter or a beta-invariant. The parameter β is crucial because it specifies the equivalenceclass of such representations, whose general structure was further studied in [KR09] in the rank 2case. The non-chiral staggered modules are far less understood and documented, and reflect thedifficulty to formulate a consistent and local logarithmic CFT; see however [DF08] and [Ri12]. Itis nonetheless believed that the parameter β present in the chiral representations plays the samerole of equivalence class label in the non-chiral ones, even if the latter may have more than onesuch label.Concentrating on the action of the chiral Virasoro modes, the parameter β arises when weconsider the triangular relations satisfied by a generic logarithmic pair ( φ, ψ ) of weights (1 ,
1) andthe associated ρ , φ ψρ The arrows coming out of ψ indicate the actions ( L − ψ = λφ and L ψ = ρ . It is importantto note that if the normalization of ψ is fixed, those of λφ and ρ are fixed as well (the value of λ depends on the way φ is normalized). The vertical arrow indicates that L − ρ is proportional to λφ , L − ρ = β ( λφ ) , (8.4)the proportionality factor β being intrinsic to the representation as all normalizations have alreadybeen fixed. In addition, these relations are invariant under the change ψ → ψ + αφ because thefield φ is primary ( L φ = 0), and so they do not depend on which logarithmic partner is considered.To answer the above question thus amounts to check whether the sandpile representation andthe symplectic representation have the same value of β . The value of β in the sandpile model hasbeen given above: if the pair ( φ, ψ ) is chosen to be ( h , h ), with the same normalization as theheight variables on the square lattice, in which case λ = − , then one finds β = [JPR06].The field h has been already identified in terms of fermions, and yields a natural choice forthe primary field φ on the symplectic side, φ θ = − P (cid:16) ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ (cid:17) . (8.5)As mentioned earlier, the lattice results in the scaling limit are consistent with h being degenerateat level 2, namely ( L − − L − ) h = 0, see Section 5.5. The same equation is satisfied by φ θ .The only candidate for the logarithmic partner of φ θ is proportional to θ ˜ θ ( ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ ) upto an irrelevant multiple of φ θ . By computing its conformal transformations via its OPE with thechiral stress-energy tensor, one finds that the following normalization, ψ θ = − P θ ˜ θ ( ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ ) , (8.6)50atisfies L ψ θ = ψ θ − φ θ , for the same value λ = − . The same OPE reveal in addition that ρ θ = L ψ θ = − P (cid:16) θ ∂ ˜ θ + ∂θ ˜ θ (cid:17) , (8.7)from which, upon using ∂∂θ = ∂∂ ˜ θ = 0, one obtains L − ρ θ = ∂ρ θ = − P (cid:16) ∂θ ∂ ˜ θ + ∂θ ∂ ˜ θ (cid:17) = 12 φ θ . (8.8)Comparing with (8.4), the value of the logarithmic coupling is found to be β θ = − h a> , and therefore does not appear to be the correct CFT to describe the scaling limit of thesandpile model .As one might suspect, the value of β has strong consequences on correlation functions involving ψ . A detailed comparison between β = versus the fermionic realization β θ = − h corresponding to avalue β = − . On general grounds, this can also be understoodfrom 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 presentcase, the singular logarithmic field is a combination of a descendant of ψ at level 5 and a descendantof ρ at level 6, with the following explicit dependence on β [KR09], ξ = (cid:16) L − − L − L − + 12 L − (cid:17)(cid:16) L − − L − (cid:17) ψ − β (cid:20) − ( β + 1) L − L − + (14 β + 5) L − L − L − − βL − − β − L − L − + 8 βL − L − − (5 β + 2) L − L − + 4 βL − (cid:21) ρ. (8.9)Using the relations L ψ = ρ, ( L − ψ = λφ , as well as the degeneracy condition ( L − − L − ) φ = 0(and the value c = − ξ satisfies L ξ = L ξ = 0 provided the identity L − ρ = βλφ holds. A rather convincing confirmation for the value of β = in the sandpile modelis therefore to check that the various correlators involving h indeed satisfy the condition ξ = 0for β = . It has been done for the correlator (6.24b).The situation seems therefore to be the following. The sandpile model contains a conformallogarithmic representation whose structure is very similar to the one appearing in the symplecticfermion theory, but which is nevertheless inequivalent to it. As far as the logarithmic partner ψ isnot brought in, the two representations look the same; this explains why some of the fields can berealized in terms of symplectic fermions. However the fermionic theory does not contain the β = representation found in the sandpile model, from which one concludes that it does not describe itsscaling limit.To characterize the CFT that does describe the sandpile model, even if a Lagrangian realizationof it cannot be found, remains an enormous challenge. At the moment, this looks to be an extremelyambitious question in view of the (very) small number of fields which have been successfullyidentified. As an example, the correlator h h ( z, z ) i mix displayed in (6.24b), and corresponding to the four-point function h φ op , cl ( x ) φ cl , op ( x ) h ( z, z ) i , can be computed upon assuming that h is a logarithmic partner of h carrying ageneric value of β = 0. Its general form is given in [Ru13]. eferences [AD95] A.A. Ali and D. Dhar, Breakdown of simple scaling in Abelian sandpile models in onedimension , Phys. Rev. E51 (1995) R2705.[ADMR10] N. Azimi-Tafreshi, H. Dashti-Naserabadi, S. Moghimi-Araghi and P. Ruelle,
TheAbelian sandpile model on the honeycomb lattice , J. Stat. Mech. (2010) P02004.[Ba96] P. Bak,
How Nature Works: The Science of Self-Organised Criticality , Springer-Verlag NewYork 1996.[BIP93] J. Brankov, E. Ivashkevich and V.B. Priezzhev,
Boundary effects in a two-dimensionalAbelian sandpile , J. Phys. I 3 (1993) 1729.[BTK87] P. Bak, C. Tang and K. Wiesenfeld,
Self-organized criticality: an explanation of the 1/fnoise , Phys. Rev. Lett. 59 (1987) 381.[Ca84] J. Cardy,
Conformal invariance and surface critical behavious , Nucl. Phys. B240 (1984)514.[Ca96] J. Cardy,
Scaling and Renormalization in Statistical Physics , Cambridge Lecture Notes inPhysics, Cambridge: Cambridge University Press (1996).[CPS08] S. Caracciolo, G. Paoletti and A. Sportiello,
Explicit characterization of the identityconfiguration in an Abelian sandpile Model , J. Phys. A41 (2008) 495003.[Cr91] M. Creutz,
Abelian sandpiles , Comput. Phys. 5 (1991) 198.[CS12] S. Caracciolo and A. Sportiello,
Exact integration of height probabilities in the Abeliansandpile model , J. Stat. Mech. (2012) P09013.[DF08] A.-L. Do and M. Flohr,
Towards the construction of local logarithmic conformal field the-ories , Nucl. Phys. B802 (2008) 475.[DFMS97] P. Di Francesco, P. Mathieu and D. Sénéchal,
Conformal field theory , Springer (1997).[Dh90] D. Dhar,
Self-organized critical state of sandpile automaton models , Phys. Rev. Lett. 64(1990) 1613.[Dh06] D. Dhar,
Theoretical studies of self-organized criticality , Physica A369 (2006) 29.[DIK13] P. Deift, A. Its and I. Krasovsky,
Toeplitz matrices and Toeplitz determinants under theimpetus of the Ising model. Some history and some recent results , Commun. Pure Appl. Math.66 (2013) 1360.[DM90] D. Dhar and S.N. Majumdar,
Abelian sandpile model on the Bethe lattice , J. Phys. A 23(1990) 4333.[DR89] D. Dhar and R. Ramaswamy,
Exactly solved model of self-organized critical phenomena ,Phys. Rev. Lett. 63 (1989) 1659. 52DRSV95] D. Dhar, P. Ruelle, S. Sen and D.-N. Verma,
Algebraic aspects of Abelian sandpiles , J.Phys. A: Math. Gen. 28 (1995) 805.[Fl03] M. Flohr,
Bits and Pieces in Logarithmic Conformal Field Theory , Int. J. Mod. Phys. A18(2003) 4497.[GK96] M.R. Gaberdiel and H.G. Kausch,
Indecomposable fusion products , Nucl. Phys. B477(1996) 293.[GK99] M.R. Gaberdiel and H.G. Kausch,
A local logarithmic conformal field theory , Nucl. Phys.B538 (1999) 631.[GLJ97] P. Ghaffari, S. Lise and H.J. Jensen,
Nonconservative sandpile models , Phys. Rev. E56(1997) 6702.[GR06] M.R.Gaberdiel and I. Runkel,
The logarithmic triplet theory with boundary , J. Phys. A:Math. Gen. 39 (2006) 14745.[GRR13] A. Gainutdinov, D. Ridout and I. Runkel eds,
Logarithmic conformal field theory , specialissue of J. Phys. A: Math. Theor. 46 (2013).[Gu93] V. Gurarie,
Logarithmic operators in conformal field theory , Nucl. Phys. B410 (1993) 535.[He99] M. Henkel,
Conformal Invariance and Critical Phenomena , Springer (1999).[HL03] C.-K. Hu and C.Yu Lin,
Universality in critical exponents for toppling waves of the BTWsandpile model on two-dimensional lattices
Physica A318 (2003) 92.[IPR07] N.Sh. Izmailian, V.B. Priezzehv and P. Ruelle,
Non-Local Finite-Size Effects in the DimerModel , SIGMA 3 (2007) 001.[IPRH05] N.Sh. Izmailian, V.B. Priezzehv, P. Ruelle and C.-K. Hu,
Logarithmic Conformal FieldTheory and Boundary Effects in the Dimer Model , Phys. Rev. Lett. 95 (2005) 260602.[Iv94] E.V. Ivashkevich,
Boundary height correlations in a two-dimensional Abelian sandpile , J.Phys. A: Math. Gen. 27 (1994) 3643.[Je05a] M. Jeng,
Conformal field theory correlations in the Abelian sandpile model , Phys. Rev. E71(2005) 016140.[Je05b] M. Jeng,
The four height variables, boundary correlations, and dissipative defects in theAbelian sandpile model , Phys. Rev. E71 (2005) 036153.[JPR06] M. Jeng, G. Piroux and P. Ruelle,
Height variables in the Abelian sandpile model: scalingfields and correlations , J. Stat. Mech. (2006) P10015.[Je98] H.J. Jensen,
Self-Organized Criticality , Cambridge University Press 1998.[KR09] K. Kytölä and D. Ridout,
On staggered indecomposable Virasoro modules , J. Math. Phys.50 (2009) 123503. 53KW15] R.W. Kenyon and D.B. Wilson,
Spanning trees of graphs on surfaces and the intensity ofloop-erased random walk on planar graphs , J. Amer. Math. Soc. 28 (2015), 985.[KW16] A. Kassel and D.B. Wilson,
The looping rate and sandpile density of planar graphs , Amer.Math. Monthly 123 (2016) 19.[LBR02] Y. Le Borgne and D. Rossin,
On the identity of the sandpile group , Discrete Math. 256(2002) 775.[LH02] C.-Yu Lin and C.-K. Hu,
Renormalization-group approach to an Abelian sandpile model onplanar lattices , Phys. Rev. E66 (2002) 021307.[Ma91] S.S. Manna,
Two-state model of self-organized criticality , J. Phys A 24 (1991) L363.[MD91] S.N. Majumdar and D. Dhar,
Height correlations in the Abelian sandpile model , J. Phys.A: Math. Gen. 24 (1991) L357.[MD92] S.N. Majumdar and D. Dhar,
Equivalence between the Abelian sandpile model and the q → limit of the Potts model , Physica A185 (1992) 129.[MR01] S. Mahieu and P. Ruelle, Scaling fields in the two-dimensional Abelian sandpile model ,Phys. Rev. E64 (2001) 066130.[MRS04] C. Maes, F. Redig and E. Saada,
The infinite volume limit of dissipative Abelian sand-piles , Commun. Math. Phys. 244 (2004) 395.[PGPR08] V.S. Poghosyan, S.Y. Grigorev, V.B. Priezzhev and P. Ruelle,
Pair correlations insandpile model: a check of logarithmic conforml field theory , Phys. Lett. B659 (2008) 768.[PGPR10] V.S. Poghosyan, S.Y. Grigorev, V.B. Priezzhev and P. Ruelle,
Logarithmic two-pointcorrelators in the Abelian sandpile model , J. Stat. Mech. (2010) P07025.[PP97] Vl.V. Papoyan and A.M. Povolotsky,
Renormalization group study of sandpile on the tri-angular lattice , Physica A246 (1997) 241.[PP11] V.S. Poghosyan and V.B. Priezzhev,
The problem of predecessors on spanning trees , Act.Polytech. 51 (2011) 59.[PPR11] V.S. Poghosyan, V.B. Priezzhev and P. Ruelle,
Return probability for the loop-erasedrandom walk and mean height in the Abelian sandpile model: a proof , J. Stat. Mech. (2011)P10004.[Pr94] V.B. Priezzhev,
Structure of Two-Dimensional sandpile I. Height Probabilities , J. Stat.Phys. 74 (1994) 955.[Pr12] G. Pruessner,
Self-Organised Criticality: Theory, Models and Characterisation , CambridgeUniversity Press 2012.[PR04] G. Piroux and P. Ruelle,
Pre-logarithmic and logarithmic fields in a sandpile model , J.Stat. Mech. (2004) P10005. 54PR05a] G. Piroux and P. Ruelle,
Logarithmic scaling for height variables in the Abelian sandpilemodel , Phys. Lett. B607 (2005) 188.[PR05b] G. Piroux and P. Ruelle,
Boundary height fields in the Abelian sandpile model , J. Phys.A: Math. Gen. 38 (2005) 1451.[PR17] A. Poncelet and P. Ruelle,
Multipoint correlators in the Abelian sandpile model , J. Stat.Mech. (2017) 123102.[PR18] A. Poncelet and P. Ruelle,
Sandpile probabilities on triangular and hexagonal lattices , J.Phys. A: Math. Theor. 51 (2018) 015002.[PRZ06] P.A. Pearce, J. Rasmussen and J.-B. Zuber,
Logarithmic minimal models , J. Stat. Mech.(2006) P11017.[Ri12] D. Ridout,
Non-chiral logarithmic couplings for the Virasoro algebra , J. Phys. A: Math.Theor. 45 (2012) 255203.[Ro96] F. Rohsiepe,
On reducible but indecomposable representations of the Virasoro algebra ,arXiv:hep-th/9611160.[RS92] P. Ruelle and S. Sen,
Toppling distributions in one-dimensional Abelian sandpiles , J. Phys.A: Math. Gen. 25 (1992) L1257.[Ru02] P. Ruelle, A c = − boundary changing operator for the Abelian sandpile model , Phys.Lett. B539 (2002) 172.[Ru07] P. Ruelle, Wind on the boundary for the Abelian sandpile model , J. Stat. Mech. (2007)P09013.[Ru13] P. Ruelle,
Logarithmic conformal invariance in the Abelian sandpile model , J. Phys. A:Math. Theor. 46 (2013) 494014.[TK00] T. Tsuchiya and M. Katori,