Scaling limits in divisible sandpiles: a Fourier multiplier approach

In this paper we complete the investigation of scaling limits of the odometer in divisible sandpiles on $d$-dimensional tori generalising the works Chiarini et al. (2018), Cipriani et al. (2017, 2018). Relaxing the assumption of independence of the weights of the divisible sandpile, we generate generalised Gaussian fields in the limit by specifying the Fourier multiplier of their covariance kernel. In particular, using a Fourier multiplier approach, we can recover fractional Gaussian fields of the form $(-\Delta)^{-(1+s)} W$ for $s>0$ and $W$ a spatial white noise on the $d$-dimensional unit torus.


Introduction and main results
Gaussian random fields arise naturally in the study of many statistical physical models. In particular fractional Gaussian fields FGF s (D) := (−∆) −s/2 W, where W denotes a spatial white noise, s ∈ R and D ⊂ R d , typically arise in the context of random phenomena with long-range dependence and are closely related to renormalization. Examples of fractional Gaussian fields include the Gaussian free field and the continuum bi-Laplacian model. We refer the reader to Lodhia et al. (2014) and references therein for a complete survey on fractional Gaussian fields. In this paper we study a class of divisible sandpile models and show that the scaling limit of its odometer functions converges to a Gaussian limiting field indexed by a Fourier multiplier.
The divisible sandpile was introduced by Peres (2009, 2010) and it is defined as follows. A divisible sandpile configuration on the discrete torus Z d n of side-length n is a function s : Z d n → R, where s(x) indicates a mass of particles or a hole at site x. Note that here, unlike the classical Abelian sandpile model (Bak et al., 1987, Járai, 2018, s(x) is a real-valued number. Given (σ(x)) x∈Z d n a sequence of centered (possibly correlated) multivariate Gaussian random variables, we choose s to be equal to s(x) = 1 + σ(x) − 1 n d z∈Z d n σ(z). (1.1) If a vertex x ∈ Z d n is unstable, i.e. s(x) > 1, it topples by keeping mass 1 for itself and distributing the excess s(x) − 1 uniformly among its neighbours. At each discrete time step, all unstable vertices topple simultaneously. The configuration s defined as (1.1) will stabilize to the all 1 configuration. The odometer u : Z d n → R ≥0 collects the information about all mass which was emitted from each vertex in Z d n during stabilization. Our main theorem states that u, properly rescaled, converges to a Gaussian random field in some appropriate Sobolev space.
In Cipriani et al. (2017) the authors consider divisible sandpiles with nearest-neighbor mass distribution, and show that for any configuration s given by (1.1) where the σ's are i.i.d. with finite variance the limiting odometer is a bi-Laplacian Gaussian field (−∆) −1 W on the unit torus T d (or FGF 2 (T d ) in the notation of Lodhia et al. (2014)). Relaxing the second moment assumption on σ leads to limiting fields which are no longer Gaussian, but alpha-stable random fields, see Cipriani et al. (2018c). On the other side if one keeps the second moment assumption and instead redistributes the mass upon toppling to all neighbours following the jump distribution of a long-range random walk one can construct fractional Gaussian fields with 0 < s ≤ 2 (Chiarini et al. (2018)). To summarize, Gaussian fields appear under the assumption of finite second moments in the initial configuration, while tuning the redistribution of the mass leads to limiting interfaces with smoothness which is at most the one of a bi-Laplacian field.
One natural question arises: what kind of sandpile models give rise to odometer interfaces which are smoother that the bi-Laplacian? It turns out that to obtain limiting fields of the form (−∆) −s/2 W such that s > 2 the long-range dependence must show up in the initial Gaussian multivariate variables σ rather than in the redistribution result. The novelty of the present article is that it complements Chiarini et al. (2018), Cipriani et al. (2017Cipriani et al. ( , 2018c by removing the assumption of independence of the weights (in the Gaussian case), and in addition provides an example of a model defined on a discrete space which scales to a limiting field (−∆) −s/2 W such that s > 2 . In the proof, we start by defining a sequence of covariance matrices (K n ) n∈N for the weights of the sandpile. Their Fourier transform K n is assumed to have a pointwise limit K as n goes to infinity. Under suitable regularity assumptions, K defines the Fourier multiplier of the covariance kernel of the limiting field. A key idea in the proof is to absorb the multiplier K into the definition of the abstract Wiener space. This defines a new Hilbert space where we will construct the limiting field. Note that this approach is different from the one in Chiarini et al. (2018), Cipriani et al. (2017), where instead the covariance structure of the odometer was given. Furthermore we would like to stress that the scaling factor a n ∼ n −2 used for convergence (see Theorem 1) is dimension-independent in contrast to the above mentioned works. This follows from the fact that, being the σ's correlated, the dimensional scaling is absorbed in the covariance structure of the odometer (see Lemma 7).
Let us finally remark that depending on the parameters s, d the limiting field will be either a random distribution or a random continuous function. More precisely, if the Hurst parameter H of the FGF s field is strictly negative, then the limit is a random distribution while for H ∈ (k, k + 1), k ∈ N ∪ {0}, the field is a (k − 1)-differentiable function (Lodhia et al. (2014), with the caveat that the results presented therein are worked out for R d or domains with zero boundary conditions). In the case of H ≥ 0, a stronger result could be pursued, namely an invariance principle à-la Donsker (as for example in Cipriani et al. (2018a, Theorem 2.1), Cipriani et al. (2018b, Theorem 3)). To keep the same outline for all proofs we will treat the limiting field a priori as a random distribution and thus prove finite dimensional distribution convergence by testing the rescaled odometer against suitable test functions.

Main result.
Notation. In all that follows, we will consider d ≥ 1. We are going to work with the spaces Z d n := [−n/2, n/2] d ∩ Z d , the discrete torus of side-length n, and T d := [−1/2, 1/2] d , the ddimensional torus. Moreover let B(z, ρ) a ball centered at z of radius ρ > 0 in the ∞ -metric. We will use throughout the notation z·w for the Euclidean scalar product between z, w ∈ R d . We will let C, C , c . . . be positive constants which may change from line to line within the same equation. We define the Fourier transform of a function f ∈ L 1 (T d ) as f (y) := T d f (z) exp −2πıy · z dz for y ∈ Z d . We will use the symbol · to denote also Fourier transforms on Z d n and R d (cf. Subsection 2.1 for the precise definitions).
We can now state our main result. We consider the piecewise interpolation of the odometer on small boxes of radius 1/2n and show convergence to the limiting Gaussian field Ξ K depending on the covariance K n of the initial sandpile configuration. This field can be represented in several ways: a convenient one is to let it be the Gaussian field with characteristic functional where f belongs to the Sobolev space We will give the analytical background to this definition in Subsection 2.2.
Theorem 1. For n ∈ N, let K n : Z d n → R >0 be such that Furthermore assume it is even and that exists. Consider an initial sandpile configuration defined by (1.1) where the collection of centered with K n the inverse Fourier transform of K n on Z d n . Let (u n (x)) x∈Z d n be the odometer associated with the collection (σ(x)) x∈Z d n via (1.1). Let furthermore a n := 4π 2 n −2 . We define the formal field on T d by Then a n Ξ K n (x) converges in law as n → ∞ to Ξ K in the topology of the Sobolev space H (for the analytical specification see Subsection 2.2). The field Ξ K and the space H −ε K (T d ) depend on K, the inverse Fourier transform of K as in (1.3).
Structure of the paper. We will give an overview of the needed results on divisible sandpiles and Fourier analysis on the torus in Section 2. The proof of the main result will be shown in Section 3. In Section 4 we discuss two classes of examples. In the first class, we consider weights with summable covariances, leading to a bi-Laplacian scaling limit. In the second class the limiting odometer is a fractional field of the form (−∆) −(1+s) W, s > 0.

Preliminaries
2.1. Fourier analysis on the torus. We will use the following inner product for 2 (Z d n ): Let ∆ g denote the graph Laplacian defined by Consider the Fourier basis of the same space given by the eigenfunctions of the Laplacian (2.1) The corresponding eigenvalues {λ w } w∈Z d n are given by Given f ∈ 2 (Z d n ), we define its discrete Fourier transform by

Consider the Fourier basis {φ
and denote Finally, we write C ∞ (T d )/ ∼ for the space of smooth functions with zero mean, that is, the space of smooth functions modulo the equivalence relation of differing by a constant.

Abstract Wiener Spaces and continuum fractional Laplacians.
In this Subsection our aim is to define the appropriate Sobolev space in which the convergence of Theorem 1 occurs. To do so, we repeat the classical construction of abstract Wiener spaces as done in Cipriani et al. (2017), Silvestri (2015).
be the limiting Fourier multiplier obtained by means of (1.3). For a ∈ R the following Proof. The linearity and conjugate symmetry are immediate. Furthermore, since (1.3) holds it follows that where the sum converges because a < 0 and f ∈ ( Define H a K (T d ) to be the Hilbert space completion of C ∞ (T d )/ ∼ with respect to the norm · K,a . Our goal is to define a Gaussian random variable . We do this by constructing an appropriate abstract Wiener space for Ξ K . We first of all recall the definition of such a space (see Stroock (2008, §8.2)).
In order to construct a measurable norm · B as above it is sufficient to construct a Hilbert-Schmidt operator on H and set · B := T · H . For a ∈ R define the continuum fractional Laplace Let us remark that we need mean zero test functions in order to cancel the atom at ν = 0 which arises from taking an inverse (a < 0) of the Laplacian in the definition above.
In the following we construct an orthonormal basis on H a K (T d ) . We set Note that f ξ is indeed an orthonormal basis of L 2 (T d )/ ∼ as we show now.
is an orthonormal basis of H a K (T d ) under the norm · K,a . Proof. First we observe that the f ξ 's are orthogonal: Next we show that all g ∈ H a K (T d ) have a Fourier expansion in the f ξ 's. Indeed, choose any g ∈ H a K (T d ); then by definition there exists a Cauchy sequence So in fact for ξ fixed, { g n (ξ)} n≥1 is a Cauchy sequence and thus has a limit g(ξ). We define h := ξ∈Z d \{0} g(ξ) f ξ . Note we have for all ∈ N: . This can be seen by applying Fatou's lemma: In this way we see that we must have g = h = ξ∈Z d \{0} g(ξ) f ξ , so g has a f ξ -Fourier expansion.
Next we will define a Hilbert-Schmidt operator on H a K (T d ). Lemma 4. Let ε > d/4 − a. The norm defined by We have the following for all ν ∈ Z d \ {0}: In this way, we see that Definition 2.2 (Definition of the limit field). Our AWS is the triple where ε is as in Lemma 4. We will choose from now on a := −1 and denote · K, −1 simply as · K . The measure µ − is the unique Gaussian law on H −ε K whose characteristic functional is given in (1.2). The field associated to Φ will be called Ξ K .
2.3. Covariance kernels. We are going to show that positive, real Fourier coefficients on Z d n correspond to a positive definite function (x, y) → K n (x − y) on Z d n × Z d n by proving the analog of Bochner's theorem on the discrete torus.
Lemma 5. The function (x, y) → K n (x, y) = K n (x − y) on (Z d n ) 2 is positive definite, and thus a well-defined covariance function, if and only if the Fourier coefficients K n are real-valued, even and positive.
Proof. Assume first that K n is positive definite. Then we have for any function c : Z d n → C that is not the zero function that 0 < x,y∈Z d n K n (x − y)c(x)c(y).
We then find As this needs to hold for all functions c : Z d n → C, we necessarily have K n (ξ) ∈ R >0 . Since K is even, we also have thus K n is even on Z d n . The other direction of the proof can be obtained in a similar way hence we will omit the proof here.

Divisible sandpile model and odometer function.
A divisible sandpile configuration s = (s(x)) x∈Z d n is a map s : Z d n → R where s(x) can be interpreted as the mass or hole at vertex x ∈ Z d n . It is known (Levine et al., 2015, Lemma 7.1) that for all initial configurations s such that x∈Z d n s(x) = n d the model will stabilize to the all 1 configuration. We consider in this paper initial divisible sandpile configurations of the form (1.1) where (σ(x)) x∈Z d n are multivariate Gaussians with mean 0 and stationary covariance matrix K n given by E(σ(x)σ(y)) = K n (x, y) (2.5) where K n is as in Theorem 1. Let us remark that putting K(z) := 1l z=0 retrieves the case when the σ's are i.i.d. We will study the following quantity. Let u n = (u n (x)) x∈Z d n denote the odometer (Levine et al., 2015, Section 1) corresponding to the divisible sandpile model specified by the initial configuration (1.1) on Z d n . In words, u n (x) denotes the amount of mass exiting from x during stabilization. Let with g z (x, y) the expected amount of visits of a simple random walk on Z d n starting from x and visiting y before being killed at z. The following characterization of u n is similar to Levine et al. (2015, Proposition 1.3) and follows a close proof strategy. Lemma 6. Let (σ(x)) x∈Z d n be a collection of centered Gaussian random variables with covariance given in (2.5) and consider the divisible sandpile s on Z d n given by (1.1). Then the sandpile stabilizes to the all 1 configuration and the distribution of the odometer u n (x) is given by Here (η(x)) x∈Z d n is a collection of centered Gaussian random variables with covariance Proof. By Lemma 7.1 in Levine et al. (2015) the sandpile stabilizes and the odometer u n satisfies = v + c for some constant c ∈ R, but since min u = 0, it must hold that c = 0 a.s. Now, since each v(x) is a linear combination of Gaussian random variables, v is again Gaussian with covariance ( 2.6) Observe that The expectation in the summation (2.6) can be calculated in the following way: If we now plug this into (2.6), we obtain where (η(x)) x∈Z d n is a collection of centered Gaussians with Now since u n − v is constant and min u n = 0, we conclude from (2.7) the desired statement: The next Lemma is concerned with yet another decomposition of the odometer function, namely it allows us to express its covariance in terms of Fourier coordinates. It is the analog of Cipriani et al. (2017, Proposition 4).
Lemma 7. Let u n : Z d n → R ≥0 be the odometer function as in Lemma 6. Then Here (χ z ) z∈Z d n is a collection of centered Gaussians with covariance Proof. We denote g x (·) := g(·, x). Subsequently we utilize Plancherel's theorem on the covariance matrix of η (see Lemma 6) to find We find g x (0) = z∈Z d n g(z, x), which does not depend on x. In this way, we see that the first term is constant, and thus vanishes in a similar way as in the proof of Proposition 1.3 of Levine et al. (2015) and the proof of Proposition 4 in Cipriani et al. (2017). Considering now the second summand above, we recall Equation (20) in Levine et al. (2015), which states that for ξ 0, We obtain To show the positive-definiteness of H n (x, y) we will show that for any function c : Z d n → C, such that c is not the zero function, x,y∈Z d n H n (x, y)c(x)c(y) > 0. First of all, since K n (z − z ) is positive definite, we have (2.8) Using Plancherel, we obtain 0 < (2.8) = n 2d Since this needs to hold for all such functions c(x), we conclude that K n (ξ) > 0 for all ξ ∈ Z d n . Next we find which concludes the proof.

Proof of Theorem 1
In this section we will prove Theorem 1 using the fact that convergence in distribution for the fields Ξ K n is equivalent to showing that for all mean-zero f ∈ C ∞ (T d )/ ∼ we have, as n goes to infinity, Observe first that is a linear combination of Gaussians, and thus Gaussian itself for each n. In order to now prove the convergence a n Ξ K n , f d → Ξ K , f it is enough to show convergence of the first and second moment for any mean-zero f ∈ C ∞ (T d )/ ∼. To this end, note that by Lemma 7 so in fact a n Ξ K n , f = 4π 2 n −2 z∈T d n u n (nz) Because f is a mean-zero function, we can neglect the random constant C and with a slight abuse of notation we write a n Ξ K n , f = 4π 2 n −2 As we have E[χ nz ] = 0 for all z ∈ T d n , it follows that E[ a n Ξ K n , f ] = 0 for all n. For the second moment note first that a n Ξ K n , f 2 = 16π 4 n −4 Here we have defined and essentially just approximated the integral with the value at the centre of the box B(z, 1 2n ). The proof will now proceed in two steps: we will first show that the first term of (3.2) goes to the desired limit variance (Proposition 8), and after that we will show that the other terms in (3.2) converge to zero (Proposition 9).

Proposition 8. We have
where · K is defined in (2.3).
Proposition 9. We have that First we will give the proof to Proposition 8. As a preliminary tool, we need to recall the following bound on the eigenvalues λ ξ .
Lemma 10 (Cipriani et al. (2017, Lemma 7)). There exists c > 0 such that for all n ∈ N and w ∈ Z d n \ {0} we have Proof of Proposition 8. First of all, If we now use the notation f n : Z d n → R for f n (·) = f (·/n) then we find that Furthermore, the left-hand side term above can be seen as an approximation of a Riemann integral, so we have the convergence f n (ξ) → f (ξ) for all ξ ∈ Z d as n → ∞, where f is the standard Fourier transform on T d . We find that (4.2) then simplifies to For the next step in the proof we use Lemma 10. Since K n is non-negative, we have that on the one hand and on the other that (3.6) We will show in the following that A, which is the right-hand side of (3.4), exhibits the following limit: and B, C vanish as n → ∞. As in Cipriani et al. (2017) we split the proof into two cases. First we give a more direct proof for d ≤ 3 and then consider d ≥ 4. We have that | f n (ξ)| 2 is uniformly bounded in ξ, so We can now use the dominated convergence theorem to obtain that In d ≥ 4 we use a mollifying procedure. Take any φ ∈ S(R d ), the space of Schwartz functions, such that φ is compactly supported on the unit cube [−1/2, 1/2) d with integral 1. We write φ κ (·) := κ −d φ (·/κ) for κ > 0. In order to show the convergence of (3.7), we split the terms: Next, we take from Cipriani et al. (2017) the following bound, where C > 0 is some constant: Since the right-hand side above converges to uniformly in n. Now taking the limit κ → 0 in (3.8) gives that lim κ→0 lim sup For the other term, we observe that since φ κ is smooth, φ κ (ξ) decays rapidly in ξ, and since the f n (ξ) are uniformly bounded, we can use the dominated convergence theorem to obtain We have that | φ κ (ξ)| ≤ 1 and φ κ (ξ) → 1 as κ → 0. Applying the dominated convergence theorem once more, we see that To conclude it remains to show that B, C as defined in (3.6) vanish. This can be achieved in a similar way to Cipriani et al. (2017, proof of Proposition 5). The key point is to observe that ξ∈Z d n \{0} K n (ξ)| f n (ξ)| 2 ≤ C is uniformly bounded in n, thanks to the uniform upper bound on K n (·) and (3.9).
Remark 11. Let us remark that the sum (3.3) may be absolutely converging for dimensions higher than 4 (for example with the choice K(ξ) = ξ −s for some s > 0 in d ≤ 3 + s ). This essentially indicates in which dimensions the limiting field is a proper Gaussian field. In the cases when the sum is not finite per se we need to use the mollifying procedure.
We continue and prove Proposition 9.
Proof of Proposition 9. Let us calculate E[R 2 n (u)] in the following way: because the K n (ξ) are uniformly bounded and ξ ≥ 1. Now write H n (x) := H n (x/n). Then the above becomes The last inequality relies on the bound given in Cipriani et al. (2017, Lemma 8): This then concludes the proof to Proposition 9.
3.1. Tightness in H −ε K . To complete the proof, we show that the convergence in law a n Ξ K n d → Ξ K holds in the Sobolev space H −ε K (T d ) for any ε > max{1 + d/4, d/2}. We state the following Theorem: Theorem 12. The sequence (a n Ξ K n ) n∈N is tight in H −ε K (T d ), in other words, for all δ > 0 there exists R δ > 0 such that Proof. The proof of this Theorem is analogous to the proof of tightness in Cipriani et al. (2017, Section 4.2). We first apply Markov's inequality and see the assertion follows as we can choose R δ such that P a n Ξ K n H We calculate the expectation and obtain E a n Ξ K We have that both 1l B(x, 1 2n ) and φ ξ ∈ L 2 (T d ) so by Cauchy-Schwarz inequality F n,ξ ∈ L 1 (T d ). Next we claim that, for some C > 0, sup ξ∈Z d sup n∈N x,y∈T d n n −4 H(nx, ny)F n,ξ (x)F n,ξ (y) ≤ C . (3.10) Remark that similarly to Cipriani et al. (2017, Equation (4.5)), we have for some C > 0. We write G n,ξ : Z d n → R for G n,ξ (·) := F n,ξ (·/n). Using this, we find x,y∈T d n n −4 H(nx, ny)F n,ξ (x)F n,ξ (y) n −4 λ −2 ξ K n (ξ)n 2d | G n,ξ (ξ)| 2 (3.11) ≤ Cn 2d Here we have exploited the fact that sup ξ∈Z d n K n (ξ) < ∞ by (1.3). Now by the triangle inequality, We then use this bound to obtain This is a constant that does not depend on n or ν, so the claim (3.10) is proven. Using it, we have by the Euler-Maclaurin formula and the boundedness of K(·) E a n Ξ K The last estimate is due to the fact that −ε < −d/2.

Some examples
In this Section we want to give some concrete examples of initial distributions and Gaussian fields that can be generated via scaling limits of the odometer. 4.1. Initial Gaussian distributions with power-law covariance. We would like to look at the case in which the initial distribution of the σ's is (σ(x)) x∈Z d n ∼ N(0, K n ) when the covariance matrix K n is polynomially decaying. As an example consider The corresponding realizations of the odometer function are indicated in Figures 1-2 and superposed in Figure 3.

4.2.
A bi-Laplacian field in the limit. The next Proposition shows that, even if (1.3) does not hold, but instead lim n→∞ n d K n (ξ) exists and is finite for all ξ, one can rescale the weights σ to go back to the setting of Theorem 1.
Proposition 13. Consider the divisible sandpile configuration defined in (1.1) with (σ(x)) x∈Z d n a sequence of centered multivariate Gaussians with covariance K n satisfying • lim n→∞ K n (w) = K(w) exists for all w ∈ Z d , • K ∈ 1 (Z d ), Figure 3. Superposition of the odometers associated to K + n and K − n .
Let (u n (x)) x∈Z d n be the associated odometer and furthermore b n := 4π 2 n (d−4)/2 C −1/2 K . We define the formal field on T d by Then b n Ξ K n converges in law as n → ∞ to Ξ K = FGF 2 (T d ) in the topology of the Sobolev space H −ε K (T d ), where ε > max {d/2, d/4 + 1}.
Proof. The basic idea is to use, rather than the weights σ as in the assumptions, the rescaled weights σ (x) := n d/2 σ(x), x ∈ Z d n . The Fourier transform of the associated covariance kernel is now K n = n d K n . Now observe that K n (ξ) = w∈Z d K(w)1l w∈Z d n e −2πıw·ξ/n .
By dominated convergence and the fact that K ∈ 1 (Z d ) we deduce lim n→∞ K n (ξ) = C K , ξ ∈ Z d . (4.1) We repeat the computation of (3.2) for b n Ξ K n , f 2 and we obtain as leading term To compute the variance E[ b n Ξ K n , f 2 ], we have From this point onwards the proof of Proposition 8 applies verbatim. In view of (4.1), the rescaling C −1 K in the variance is done to obtain a limiting field with characteristic functional which identifies the FGF 2 (T d ).
4.3. Fractional limiting fields. We can now use Theorem 1 to construct (1+s)-Laplacian limiting fields for arbitrary s ∈ (0, ∞). Define an initial sandpile configuration such that lim n→∞ K n (ξ) = K(ξ) := ξ −4s , ξ ∈ Z d \ {0}. (4.3) We do need K(0) > 0 to ensure positive definiteness of the kernel, so one can choose any arbitrary constant to satisfy this constraint. Using the above Fourier multiplier (observe that it is indeed uniformly bounded for s ∈ (0, ∞) and hence satisfies (1.3)) and applying Theorem 1, we find the following limiting distribution: for all f ∈ C ∞ (T d )/ ∼ we have that Ξ K , f is a centered Gaussian with variance

Initial Gaussian distributions with fractional Laplacian covariance.
We have seen in the previous example that if we define an initial configuration on Z d n with covariance given by K n (·) = · −4s ∨ , where the inverse Fourier transform is on Z d n , the limiting field is Gaussian with variance (for every f ∈ C ∞ (T d ) with mean zero): However, the covariance K n does not necessarily agree with that of the discrete s-Laplacian field on Z d n . Namely, we recall the definition of (minus) the discrete fractional Laplacian (−∆ g ) s on Z where the weights p