Properties of the Gradient Squared of the Discrete Gaussian Free Field

In this paper we study the properties of the centered (norm of the) gradient squared of the discrete Gaussian free field in Uε=U/ε∩Zd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U_{\varepsilon }=U/\varepsilon \cap \mathbb {Z}^d$$\end{document}, U⊂Rd\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U\subset \mathbb {R}^d$$\end{document} and d≥2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d\ge 2$$\end{document}. The covariance structure of the field is a function of the transfer current matrix and this relates the model to a class of systems (e.g. height-one field of the Abelian sandpile model or pattern fields in dimer models) that have a Gaussian limit due to the rapid decay of the transfer current. Indeed, we prove that the properly rescaled field converges to white noise in an appropriate local Besov-Hölder space. Moreover, under a different rescaling, we determine the k-point correlation function and joint cumulants on Uε\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$U_{\varepsilon }$$\end{document} and in the continuum limit as ε→0\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\varepsilon \rightarrow 0$$\end{document}. This result is related to the analogue limit for the height-one field of the Abelian sandpile (Dürre in Stoch Process Appl 119(9):2725–2743, 2009), with the same conformally covariant property in d=2\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d=2$$\end{document}.


Introduction
The Gaussian free field (GFF) is one of the most prominent models for random surfaces.It appears as scaling limit of observables in many interacting particle systems, see for example Jerison et al. (2014), Kenyon (2001), Sheffield (2007), Wilson (2011).It serves as a building block for defining the Liouville measure in Liouville quantum gravity (see Ding et al. (2021) and references therein for a list of works on the topic).
Its discrete counterpart, the discrete Gaussian free field (DGFF), is also very well-known among random interface models on graphs.Given a graph Λ, a (random) interface model is defined as a collection of (random) real heights Γ = (Γ (x)) x∈Λ , measuring the vertical distance between the interface and the set of points of Λ (Funaki (2005), Velenik (2006)).The discrete Gaussian free field has attracted a lot of attention due to its links to random walks, cover times of graphs, and conformally invariant processes (see Barlow and Slade (2019), Ding et al. (2012), Glimm and Jaffe (1987), Schramm and Sheffield (2009), Sheffield (2007), among others).In the present paper, we will consider the DGFF on the square lattice, that is, we will focus on Λ ⊆ Z d , in which case the probability measure of the DGFF is a Gibbs measure with formal Hamiltonian given by H(Γ ) = 1 2d x, y: x−y =1 V Γ (x) − Γ (y) , (1.1) where V(ϕ) = ϕ 2 /2.We will always work with 0-boundary conditions, meaning that we will set Γ (x) to be zero almost surely outside Λ.For general potentials V(•) the Hamiltonian (1.1) defines a broad class of gradient interfaces which have been widely studied in terms of decay of correlations and scaling limits (Biskup and Spohn (2011), Cotar et al. (2009), Nadaf and Spencer (1997)), among others.
The gradient Gaussian free field ∇Γ is defined as the gradient of the DGFF Γ along edges of the square lattice.This field is a centered Gaussian process whose correlation structure can be written in terms of T (•, •), the transfer current (or transfer impedance) matrix (Kassel and Wu  (2015)).Namely, if we consider the gradient ∇ i Γ (•) := Γ (• + e i ) − Γ (•) in the i-th coordinate direction of R d , we have, for x, y ∈ Z d , 1 ≤ i, j ≤ d, that E ∇ i Γ (x)∇ j Γ (y) = G Λ (x, y) − G Λ (x + e i , y) − G Λ (x, y + e j ) + G Λ (x + e i , y + e j ) = T (e, f) where e = (x, x + e i ) and f = (y, y + e j ) are directed edges of the grid and G Λ (•, •) is the discrete harmonic Green's function on Λ with 0-boundary conditions outside Λ.Here T (e, f) describes a current flow between e and f.
The main object we will study in our article is the following.Take U to be a smooth, connected, bounded subset of R d .Consider the recentered squared norm of the gradient DGFF, formally denoted by on the discretized domain U ε = U/ε ∩ Z d , ε > 0, d ≥ 2, with Γ a 0-boundary DGFF on U ε .The colon : (•) : denotes the Wick centering of the random variables.In the rest of the paper we will simply call Φ ε the gradient squared of the DGFF.Let us remark that we do not consider d = 1 here since in one dimension the gradient of the DGFF is a collection of i.i.d.Gaussian variables.

k-point correlation functions
Our first main result determines the k-point correlation functions for the field Φ ε on the discretized domain U ε and in the scaling limit as ε → 0. We defer the precise statement to Theorem 1 in Section 3, which we will now expose in a more informal way.Let ε > 0 and k ∈ N and let the points x (1) , . . ., x (k) in U ⊂ R d , d ≥ 2, be given.Define x ε to be a discrete approximation of x (j) in U ε , for j = 1, . . ., k.Let Π([k]) be the set of partitions of k objects and S 0 cycl (B) be the set of cyclic permutations of a set B without fixed points.Finally let E be the set of coordinate vectors of R d .Then the k-point correlation function at fixed "level" ε is equal to Moreover if x (i) = x (j) for all i = j, the scaling limit of the above expression is η(σ(j)) G U x (j) , x (σ(j)) (1.3) where G U (•, •) is the continuum Dirichlet harmonic Green's function on U.As a corollary (Corollary 1) we also determine the corresponding cumulants on U ε and in the scaling limit.
Let us discuss some interesting observations in the sequel.The k-point correlation function of (1.3) has similarities to the k-point correlation that arises in permanental processes, see Eisenbaum and Kaspi (2009), Hough et al. (2009), Last and Penrose (2017) for relevant literature.In fact, in d = 1 one can show that the gradient squared is exactly a permanental process with kernel given by the diagonal matrix whose non-zero entries are the double derivatives of G U (McCullagh and Møller, 2006, Theorem 1).In higher dimensions however we cannot identify a permanental process arising from the scaling limit, since the directions of derivations of the DGFF at each point are not independent.Nevertheless the 2-point correlation functions of Φ ε are positive (see Equation (4.20) in Section 4), which is consistent with attractiveness of permanental processes (Last and Penrose, 2017, Remark on p. 139), and the overall structure resembles closely that of permanental processes marginals.
In d = 2, the limiting cumulants κ of our field are interestingly connected to the cumulants of the height-one field h ε (x (Dürre, 2009a,  Theorem 2).Theorem 1 will imply that for every set of ℓ ≥ 2 pairwise distinct points in d = 2 one has see Dürre (2009b, Theorem 6).
We would also like to point out that the apparently intricate structure of Equations (1.2)-(1.3)and of Dürre's Theorem 2 can be unfolded as soon as one recognizes therein the structure of a Fock space.We will discuss this point in more detail in Subsection 3.1, where in particular in Corollary 2 we will derive a Fock space representation of the k-point function for the height-one field.We will pose further questions on this matter in the Discussion Section 5.
Due to the similar nature of the cumulants in the height-one field of the sandpile and our field, we show in Proposition 1 that in d = 2 the k-point correlation functions are conformally covariant (compare Dürre (2009a, Theorem 1), Kassel and Wu (2015, Theorem 2)).This hints at Theorems 2 and 3 of Kassel and Wu (2015), in which the authors prove that for finite weighted graphs the rescaled correlations of the spanning tree model and minimal subconfigurations of the Abelian sandpile have a universal and conformally covariant limit.

Scaling limit
The second main result of our paper is the scaling limit of the field towards white noise in some appropriate local Besov-Hölder space.As we will show in Theorem 2, Section 3, as ε → 0 the gradient squared of the discrete Gaussian free field Φ ε converges as a random distribution to spatial white noise W: for some explicit constant 0 < χ < ∞.The result is sharp in the sense that we obtain convergence in the smallest Hölder space where white noise lives.The constant χ, defined explicitly in (3.6), is the analogue of the susceptibility for the Ising model, in that it is a sum of all the covariances between the origin and any other lattice point.We will prove that this constant is finite and the field Φ ε has a Gaussian limit.Note that Newman (1980) proves the same result for translation-invariant fields with finite susceptibility satisfying the FKG inequality.In our case we do not have translation invariance since we work on a domain, so we are not able to apply directly this criterion.From a broader perspective there are several other results in the literature that obtain white noise in the limit due to an algebraic decay of the correlations, see for example Bauerschmidt et al. (2014).Note that our field can be understood in a wider class of models having correlations which depend on the transfer current matrix T (•, •).An interesting point mentioned in Kassel and Wu (2015) is that pattern fields of determinantal processes closely connected to the spanning tree measure and T (•, •) (for example the spanning unicycle, the Abelian sandpile model (Dürre (2009a)) and the dimer model (Boutillier (2007))) have a universal Gaussian limit when viewed as random distributions.Correlations of those pattern fields can be expressed in terms of transfer current matrices which decay sufficiently fast and assure the central limit-type behaviour which we also obtain.
Let us comment finally on the differences between expressions (1.3) and (1.6).The scaling factors are different, and this reflects two viewpoints one can have on Φ ε : the one of (1.3) is that of correlation functionals in a Fock space, while in (1.6) we are looking at it as a Gaussian distributional field (compare also Theorems 2 and 3 in Dürre (2009a)).This is compatible, as there are examples of trivial correlation functionals which are non-zero as random distributions (Kang and Makarov (2013)).
The novelty of the paper lies in the fact that we construct the gradient squared of the Gaussian free field on a grid, determine its k-point correlation function and scaling limits.We determine tightness in optimal Besov-Hölder spaces (optimal in the sense that we cannot achieve a better regularity for the scaling limit to hold).Furthermore we show the "dual" behavior in the scaling limit of the gradient squared of the DGFF as a Fock space field and as a random distribution.As mentioned before we recognize a similarity to permanental processes, and it is worthwhile noticing that for general point processes there is a Fock space structure, see e.g.Last and Penrose (2017, Section 18).Since there is a close connection to the height-one field via correlation structures, we also unveil a Fock space structure in the Abelian sandpile model.

Proof ideas
The main idea for the proof of results (1.2)-(1.3) is to decompose the k-point correlation function in terms of complete Feynman diagrams (Janson (1997)).Then we can use the Gaussian nature of the field to expand the products of covariances as transfer currents.To determine the scaling limit we will use developments from Funaki (2005) and Kassel and Wu (2015).Let us stress that the proof of the scaling limit of cumulants differs from the one of Dürre (2009a, Theorem 2) who instead uses the correspondence between the height-one field and the spanning tree explicitly to determine the limiting observables.The proof of the scaling limit (1.6) is divided into two parts.In a first step (Proposition 2) we prove that the family of fields under consideration is tight in an appropriate local Besov-Hölder space by using a tightness criterion of Furlan and Mourrat (2017).The proof requires a precise control of the summability of k-point functions, which is provided by Theorem 1 and explicit estimates for double derivatives of the Green's function in a domain.Observe that, even if the proof relies on the knowledge of the joint moments of the family of fields, we only use asymptotic bounds derived from them.More specifically, we need to control the rate of growth of sums of moments at different points.The second step (Proposition 3) consists in determining the finite-dimensional distributions and identifying the limiting field.We will first show that the limiting distribution, when tested against test functions, has vanishing cumulants of order higher or equal to three, and secondly that the limiting covariance structure is the L 2 (U) inner product of the test functions.This will imply that the finite-dimensional distributions of converge to those corresponding to d-dimensional white noise.For this we rely on generalized bounds on double gradients of the Green's function from Lawler and Limic (2010) and Dürre (2009a).

Structure of the paper
The structure of the paper is as follows.In Section 2 we fix notation, introduce the fields that we study and provide the definition of the local Besov-Hölder spaces where convergence takes place.Section 3 is devoted to stating the main results in a more precise manner.The subsequent Section 4 contains all proofs and finally in Section 5 we discuss possible generalizations and pose open questions.
We will write |A| for the cardinality of a set A. For any finite set A we define Π(A) as the set of all partitions of A. Let Perm(A) denote the set of all possible permutations of the set A (that is, bijections of A onto itself).When A = [k] for some k ∈ N, we might also refer to its set of permutations as S k .If we restrict S k to those permutations without fixed points, we denote them as S 0 k .Call S cycl (A) the set of the full cyclic permutations of A, possibly with fixed points.More explicitly, any σ : When this condition is relaxed to all A ′ with A ′ > 0 we obtain the set of all cyclic permutations without fixed points which is called S 0 cycl (A).
Definition 1 (Cumulants of a random variable).Let X be a random variable in R d with all moments finite.The cumulants (κ n (X)) n∈N are obtained from the power series expansion of the cumulant generating function as Let V ⊂ R d be a finite set of points, and (X v ) v∈V be a family of random variables in R d with all finite moments.Definition 2 (Cumulants of a random vector).The joint cumulants κ(X v : v ∈ V) of the random vector (X v ) v∈V are defined as

Functions of the Gaussian free field and white noise
Let U ⊂ R d , d ≥ 2, be a non-empty bounded connected open set with smooth boundary.Denote by (U ε , E ε ) the graph with vertex set U ε := U/ε ∩ Z d and edge set E ε .For an (oriented) edge e ∈ E ε of the graph, we denote by e + its tip and e − its tail, and write the edge as e = (e − , e + ).Consider E := {e i } 1≤i≤d , the canonical basis of R d .Since we will use approximations via grid points, we need to introduce, for any t ∈ R d , its floor function as For any finite set V ⊂ Z d define the discrete (normalized) Laplacian by where δ is the Dirac delta function.
Denote by G 0 (•, •) the Green's function for the whole grid Z d when d ≥ 3, or with a slight abuse of notation the potential kernel for d = 2.This abuse of notation is motivated by the fact that we will only be interested in the discrete differences of G 0 , which exist for the infinite-volume grid in any dimension.Notice that for y ∈ U, where ∆ denotes the continuum Laplacian.
Definition 6 (Gradient squared of the DGFF).The discrete stochastic field Φ ε given by : is called the gradient squared of the DGFF, where : : denotes the Wick product; that is, : The family of random fields (Φ ε ) ε>0 is a family of distributions, which acts on a given where •, • denotes the scalar product in L 2 (U).
When no ambiguities appear, we will write For any given function f : Z d × Z d → R, we define the discrete gradient in the first argument and direction e i ∈ E as

and analogously
for the second argument.Once again, when no ambiguities arise, we will write ∇ (1) i for ∇ (1) e i , and analogously for the second argument.
For a continuum function g : e i g(x, y) denotes the partial derivative of g with respect to the first argument in the direction of e i , while ∂ (2) e i g(x, y) corresponds to the second argument, also in the direction e i .The same abuse of notation on the subindex e i applies here.

Besov-Hölder spaces
In this Subsection we will define the functional space on which convergence will take place.We will use Furlan and Mourrat (2017) as a main reference.Local Hölder and Besov spaces of negative regularity on general domains are natural functional spaces when considering scaling limits of certain random distributions or in the context of non-linear stochastic PDE's, see e.g.Furlan and Mourrat (2017), Hairer (2014) especially when those objects are well-defined on a domain U ⊂ R d but not necessarily on the full space R d .They are particularly suited for fields which show bad behaviour near the boundary ∂U.
for some positive integer r ∈ N, that is, they belong to the set of r times continuously differentiable functions on R d with compact support.
Define moreover which makes (φ n,x ) x∈Λn an orthonormal basis of V n resp.(ψ The local Besov space B α,loc p,q (U) is the completion of C ∞ (U) with respect to the family of semi-norms We will use the following embedding property of Besov spaces in the tightness argument.
Lemma 1 (Furlan and Mourrat (2017, Remark 2.12)).For any Finally let us define the functional space where convergence will take place, the space of distributions with locally α-Hölder regularity.
Definition 9 (Hölder spaces).Let α < 0, r 0 = −⌊α⌋.The space C α loc (U) is called the locally Hölder space with regularity α ∈ R on the domain U.It is the completion of C ∞ c (U) with respect to the family of semi-norms Note that by Furlan and Mourrat (2017, Remark 2.18) one has

Main results
The first result we would like to present is an explicit computation of the k-point correlation function of the gradient squared of the DGFF field Φ ε defined in Definition 6.
Remark 1.It will sometimes be useful to write (3.1) as the equivalent expression 3) where the condition of σ belonging to full cycles of B without fixed points is inserted in the no-singleton condition of the permutations π.
Remark 2. From the above expression it is immediate to see that the 2-point function is given by which will be useful later on.
The following Corollary is a direct consequence of Theorem 1. j) for all i = j, then As already mentioned in the introduction, comparing our result with Dürre (2009a, Theorem 2) we obtain (1.4).
The following proposition states that in d = 2 the limit of the field Φ ε is conformally covariant with scale dimension 2. This result can also be deduced for the height-one field for the sandpile model, see Dürre (2009a, Theorem 1).
, and x (j) ε j∈[k] be as in Theorem 1. Furthermore let h : U → U ′ be a conformal mapping and call h ε x (j) := ⌊h x (j) /ε⌋.Then Finally we will show that the rescaled gradient squared of the discrete Gaussian free field will converge to white noise in some appropriate locally Hölder space with negative regularity α in d ≥ 2 dimensions.This space is denoted as C α loc (U) (see Definition 9).
Theorem 2. Let U ⊂ R d with d ≥ 2. The gradient squared of the discrete Gaussian free field Φ ε converges in the following sense as ε → 0: where the white noise W is defined in Definition 7.This convergence takes place in C α loc (U) for any α < −d/2, and the constant χ defined as Remark 3. Let us remind the reader that C α loc (U) with α < −d/2 are the optimal spaces in which the white noise lives.See for example Armstrong et al. (2017, Proposition 5.9).

Fock space structure
Let us discuss in the following the connection to Fock spaces.We start by reminding the reader of the definition of the continuum Gaussian free field (GFF).
Definition 10 (Continuum Gaussian free field, Berestycki (2015, Section 1.5)).The continuum Gaussian free field Γ with 0-boundary (or Dirichlet) conditions outside U is the unique centered Gaussian process indexed by where G U (•, •) was defined in Definition 4.
We can think of it as an isometry Γ : H → L 2 (Ω, P), for some Hilbert space H and some probability space (Ω, F , P).To fix ideas, throughout this Section let us fix H := H 1 0 (U), the order one Sobolev space with Dirichlet inner product (see Berestycki (2015, Section 1.6)).Note that, even if the GFF is not a proper random variable, we can define its derivative as a Gaussian distributional field.
Definition 11 (Derivatives of the GFF, Kang and Makarov (2013, pg. 4)).The derivative of Γ is defined as the Gaussian distributional field ∂ i Γ , 1 ≤ i ≤ d, in the following sense: There is however another viewpoint that one can take on the GFF and its derivatives, and is that of viewing them as Fock space fields.This approach will be used to reinterpret the meaning of Theorem 1.For the reader's convenience we now recall here some basic facts about Fock spaces and their fields.Our presentation is drawn from Janson (1997, Section 3.1) and Kang and Makarov (2013, Sec. 1.2-1.4).
For n ≥ 0, we denote H ⊙n as the n-th symmetric tensor power of H; in other words, H ⊙n is the completion of linear combinations of elements f 1 ⊙ • • • ⊙ f n with respect to the inner product The symmetric Fock space over H is Fock(H) := n≥0 H ⊙n .
We now introduce elements in Fock(H) called Fock space fields.We call basic correlation functionals the formal expressions of the form for p ∈ N, x 1 , . . ., x p ∈ U, and X 1 , . . ., X p derivatives of Γ .The set S(X p ) := {x 1 , . . ., x p } is called the set of nodes of X p .Basic Fock space fields are formal expressions written as products of derivatives of the Gaussian free field Γ , for example 1 ⊙ Γ , ∂Γ ⊙ Γ ⊙ Γ etc.A general Fock space field X is a linear combination of basic fields.We think of any such X as a map u → X(u), u ∈ U, where the values X = X(u) are correlation functionals with S(X ) = {u}.Thus Fock space fields are functional-valued functions.Observe that Fock space fields may or may not be distributional random fields, but in any case we can think of them as functions in U whose values are correlation functionals.
Our goal is to define now tensor products.We will restrict our attention to tensor products over an even number of correlation functionals, even if the definition can be given for an arbitrary number of them.The reason behind this presentation is due to the set-up we will be working with.
Definition 12 (Tensor products in Fock spaces).Let m ∈ 2N.Given a collection of correlation functionals with pairwise disjoint S(X j )'s, the tensor product of the elements X 1 , . . ., X m is defined as where the sum is taken over Feynman diagrams γ with vertices u labeled by functionals X pq in such a way that there are no contractions of vertices in the same S(X p ). E γ denotes the set of edges of γ.One extends the definition of tensor product to general correlation functionals by linearity.
The reader may have noticed that (3.7) is simply one version of Wick's theorem.It is indeed this formula that will allow us in Subsection 4.2 to prove Theorem 1, and that enables one to bridge Fock spaces and our cumulants in the following way.For any j ∈ [k], k ∈ N, i j ∈ [d] and x (j) ∈ U one can define the basic Fock space field X i j x (j) := ∂ i j Γ x (j) .Introduce the correlation functional j) . (3.8) We obtain now the statement of the next Lemma.
Lemma 2 (k-point correlation functions as Fock space fields).Under the assumptions of Theorem 1, where Here |π| stands for the number of blocks of the partition π.
The Fock space structure is more evident from the Gaussian perspective of the DGFF, but (1.4) together with Dürre's theorem entail a corollary which we would like to highlight.We remind the reader of the definition of the constant C in (1.5).
Corollary 2 (Height-one field k-point functions, d = 2).With the same notation of Theorem 1 one has in d = 2 that where As before, |π| stands for the number of blocks of the partition π.

Previous results from literature
Let us now expose some important results that we will refer to throughout the proofs.They refer to partially known results and partially consist of straightforward generalizations of previous results.
Our computations will rely on the fact that the distribution of the gradient field ∇ i Γ , i ∈ [d], is well-known.The following result is quoted from Funaki (2005, Lemma 3.6).
Lemma 3. Let Λ ⊂ Z d be finite, and let (Γ x ) x∈Λ be a 0-boundary conditions DGFF on Λ (see Definition 5).Then Consequently, we can directly link the gradient DGFF to so-called transfer current matrix where e, f are oriented edges of Λ (see Kassel and Wu (2015, Section 2)).Equivalently we can write From Lemma 3, it is clear that we need to control the behaviour of double derivatives of discrete Green's function in the limit ε → 0. In order to find the limiting joint moments of the point-wise field Φ ε (x) we will need the following result about the convergence of the discrete difference of the Green's function on U ε (see Definition 3) to the double derivative of the continuum Green's function G U (•, •) on a set U (see Equation (2.1)).This result follows from Theorem 1 of Kassel and Wu (2015).
Lemma 4 (Convergence of the Green's function differences).Let v, w be points in the set U, with v = w.Then for all a, b ∈ E, The next lemma is a generalization of Dürre (2009b, Lemma 31) for general dimensions d ≥ 2. The proof is straightforward and will be omitted.It provides an error estimate when replacing the double difference of G Uε (•, •) on the finite set by that of G 0 (•, •) defined on the whole lattice.
Lemma 5. Let D ⊂ U be such that the distance between D and U is non-vanishing, that is, dist (D, ∂U) := inf (x,y)∈D×∂U |x − y| > 0. There exist c D > 0 and ε D > 0 such that, for all ε ∈ (0, ε D ], for all v, w ∈ D and also An immediate consequence of (4.4) and the expression (3.1) in Theorem 1 for two points gives us the following bound on the covariance of the field: Corollary 3. Let D, v and w be as in Lemma 5. Then On the other hand, we will also make use of a straightforward extension of Lawler and Limic (2010, Corollary 4.4.5)for d = 2 and Lawler and Limic (2010, Corollary 4.3.3)for d ≥ 3, yielding the following Lemma.
Lemma 6 (Asymptotic expansion of the Green's function differences).As |v| → +∞, for all i, j ∈ [d] ∇ (1) i ∇ (2) The following technical combinatorial estimate, which is an immediate extension of a corollary of Dürre (2009b, Lemma 37), will be important when proving tightness of the family (Φ ε ) ε , in order to bound the rate of growth of the moments of Φ ε , f for some test function f: ,

Proof of Theorem 1
The strategy to prove the first theorem is based on decomposing the k-point functions into combinatorial expressions that involve basically covariances of Gaussian random variables.This is made possible by our explicit knowledge of the Gaussian field which underlies Φ ε .These covariances can be estimated using the transfer matrix (Equation (4.2)), whose scaling limit is well-known: it is the differential of the Laplacian Green's function (cf.Kassel and Wu (2015, Theorem 1)).
The symbol ν(γ) is the so-called value of the Feynman diagram γ.When γ is complete, for a collection of vertices {ξ α } 1≤α≤m , m ∈ 2N 0 , it is defined as and α s = β s for all s.We refer to Janson (1997, Chapter 1.5) for an exposition of Feynman diagrams in this context.

Proof of Theorem 1. Let us compute the function
From Definition 6 of Φ ε x (j) ε we know that with E the canonical basis of R d .In our case we have k products of the Wick product 2 : (indexed by j, not i j ).So we can identify ξ ij in Theorem 3 with ∇ i j Γ x (j) ε . Since these terms are squared, we need to consider each of them in the product twice.
Let us denote x (j) (we drop the dependence on ε to ease notation).Also to make notation lighter we fix the labels i j for the moment and keep them implicit.We then define U := x (1) , x (1) , . . ., x (k) , x (k) , where each copy is considered distinguishable.We also define FD 0 as the set of complete Feynman diagrams on U such that no edge joins x (i) with (the other copy of) x (i) .That is, a typical edge b in a Feynman diagram γ in FD 0 is of the form x (j) , x (m) , with j = m and j, m ∈ [k].Thus by (4.6) we have , where E γ are the edges of γ (note that the edges of γ connect edges of Now we would like to express Feynman diagrams in terms of permutations.We first note that any given γ ∈ FD 0 cannot join x (i) with itself (neither the same nor the other copy of itself).So instead of considering permutations σ ∈ Perm(U) we consider permutations σ ′ ∈ S k , being S k the group of permutations of the set [k].Any γ ∈ FD 0 is a permutation σ ∈ Perm(U), but given the constraints just mentioned, we can think of them as permutations σ ′ ∈ S k without fixed points; that is, with c(σ ′ ) a constant that takes into account the multiplicity of different permutations σ that give rise to the same σ ′ depending on its number of subcycles, which we will calculate later on.
Let us disassemble this expression even more.In general σ ′ can be decomposed in q cycles.Since σ ′ ∈ S 0 k (in particular, it has no fixed points), there are at most ⌊k/2⌋ cycles in a given , where the notation j ∈ σ ′ h means that j belongs to the domain where σ ′ h acts (non trivially).Now we note that a cyclic decomposition of a permutation of the set [k] determines a partition π ∈ Π [k] (although not injectively).This way, a sum over the number of partitions q and σ ′ ∈ S 0 k with q cycles can be written as a sum over partitions π with no singletons, and a sum over full cycles in each block B (that is, those permutations consisting of only one cycle).Hence where we also made the switch between B∈π and σ∈Scycl(B) by grouping by factors.Alternatively, we can express this average in terms of S 0 cycl (B), the set of full cycles without fixed points, as As for the multiplicity factor c(σ), a straightforward counting calculation shows that it only depends on the size of the partitions B, and is equal to Finally, we need to put back the subscript i j in the elements x (j) and sum over i 1 , . . ., i k ∈ E. Note that for any function f : , and grouping the η(j)'s according to each block B ∈ π we get .
Regarding the transfer matrix T , using Equation (4.2) we can write the above expression as obtaining the first result of the theorem.Finally, using Lemma 4 we obtain the second statement.

Proof of Corollary 1 and Proposition 1
Proof of Corollary 1. Recall that Definition 2 yields The result then follows from expressions (3.1) and (3.2) in Theorem 1.
The equality in absolute value between our cumulants and those of Dürre (2009a, Theorem 1) allow us to adapt his proof and conclude that, in the case of d = 2, our field is conformally covariant with scale dimension 2.
Proof of Proposition 1.It is known (Berestycki, 2015, Proposition 1.9) that the continuum Green's function G U (•, •), defined in Equation (2.1), is conformally invariant against a confor- mal mapping h : Recalling expression (3.5) for the limiting cumulants we see that, for any integer ℓ ≥ 2, From the cumulants expression we deduce that, for a given permutation σ and assignment η, each point x (j) will appear exactly twice in the arguments of the product of differences of G U ′ .Thus, using the chain rule and the Cauchy-Riemann equations, for a fixed σ we obtain an overall factor ℓ j=1 h ′ x (j) 2 after summing over all η.We then obtain The result follows plugging this expression into the moments.

Proof of Theorem 2
The proof of this Theorem will be split into two parts.First we will show that the family (Φ ε ) ε>0 is tight in some appropriate Besov space and then we will show convergence of finite-dimensional distributions ( Φ ε , f i ) i∈[m] and identify the limit.

Tightness
Recall that the local Besov space B α,loc p,q (U) was defined in Definition 8 and the local Hölder space C α loc (U) in Definition 9.

Finite-dimensional distributions
Proposition 3. Let U ⊂ R d and d ≥ 2. There exists a normalization constant χ > 0 such that, for any set of test functions , m ∈ N , the random elements Φ ε , f i converge in the following sense: as ε → 0.

Proof of Proposition 2
We will use the tightness criterion given in Theorem 2.30 in Furlan and Mourrat (2017).First we need to introduce some notation.Let f and (g Let K ⊂ U be compact and k ∈ N. We say that the pair We say that the set K is a spanning sequence if it can be written as where (K n ) is an increasing sequence of compact subsets of U such that n K n = U, and for every n the pair (K n , k n ) is adapted.
Theorem 4 (Tightness criterion, Furlan and Mourrat (2017, Theorem 2.30)).Let f, (g ) with the support properties mentioned above, and fix p ∈ [1, ∞), q ∈ [1, ∞] and α, β ∈ R satisfying |α| , |β| < r, α < β.Let (Φ m ) m∈N be a family of random linear forms on C r c (U), and let K be a spanning sequence.Assume that for every (K, k) ∈ K , there exists a constant c = c(K, k) < ∞ such that for every m ∈ N, and Then the family (Φ m ) m is tight in B α,loc p,q (U).If moreover α < β − d/p, then the family is also tight in C α loc (U).Proof of Proposition 2. We will consider an arbitrary scaling ε γ , γ ∈ R, and then choose an optimal one to make the fields tight.We define Φ ε as the scaled version of Φ ε , that is, The family of random linear forms (Φ m ) m∈N in Theorem 4 is to be identified with the fields ( Φ ε ) ε>0 taking for example ε decreasing to zero along a dyadic sequence.Now let us expand the expressions (4.9) and (4.10) in Theorem 4. To simplify notation, let us define , and analogously for g (i) .
In the proof we will set p ∈ 2N.This will not affect the generality of our results because of the embedding of local Besov spaces described in Lemma 1.This means that we can read (4.9) forgetting the absolute value in the left-hand side.Let us rewrite the p-th moment of (4.11)We will seek for a more convenient expression to work with.If we allow ourselves to slightly abuse the notation for Φ ε , then we can express it in a piece-wise continuous fashion as where S a (y) is the square of side-length a centered at y.Under a change of variables, if we define U ε := U ∩ εZ 2 (mind the superscript and the definition which is different from that of This way, expression (4.11) now reads Therefore the left-hand side of expression (4.9) from Theorem 4 is upper-bounded by Analogously, expression (4.10) from Theorem 4 reads for some δ > 0 and R such that (4.8) holds.Let us first consider (4.13).Given that supp g (i) 2 n (• − x) ⊂ B x (R2 −n ) we can restrict the sum over y j to the set We now bound (4.13) separately for the cases 2 n ≥ Rε −1 and 2 n < Rε −1 .
If 2 n ≥ Rε −1 , we have The sum over Ω n,x can be bounded by a sum over a finite amount of points independent of n, since under the condition 2 n ≥ Rε −1 the set Ω x,n has at most 3 d points for any x and n.
Using the fact that Sε(y j ) which gives the bound Then, for any γ ≤ 0 we can bound the above expression by a constant multiple of 2 −γn .On the other hand, if 2 n < Rε −1 , we have Sε(y j ) We also note that Using this and calling N := ⌊2R2 −n ε −1 ⌋, we obtain y 1 ,...,yp∈Ωn,x Let us first study the behaviour of this expression for p = 2.By Corollary 3 we get Observe that Theorem 4 allows the constant c to depend on k and K so the symbol K is not an issue.Let us now analyze E Φ ε (y 1 ) . . .Φ ε (y p ) for an arbitrary p.Looking at expression (3.3) we observe the following: any given partition π ∈ Π([p]) with no singletons can be expressed as π = {B 1 , . . ., B k }, with 1 ≤ k ≤ p such that 1≤i≤k n i = p, with n i := |B i |.Then the cumulant corresponding to any given B i (see Corollary 1) is proportional to a sum (over σ ∈ S cycl (B i ) and η : B i → E) of terms of the form η(σ(j)) G Uε y j , y σ(j) .
Using (4.4) we can bound this expression (up to a constant) by where the minimum takes care of the case in which B i has repeated values.Lemma 7 implies that E Φ ε (y 1 ) . . .Φ ε (y p ) can be bounded (up to a constant depending on K) by a sum of terms of magnitude Since the sum takes place over partitions of the set [p] with no singletons, putting everything back into (3.3)we see that the term with the largest value of |π| will dominate for large N.For p even this happens when π is composed of cycles of two elements, in which case |π| = p/2.Hence, for p even.Finally, If γ ≥ −d/2 then we can bound the above expression by a constant multiple of 2 dn 2 2 −(γ+ d 2 )n = 2 −γn .Otherwise, we cannot bound it uniformly in ε, as the bound depends increasingly on ε as it approaches 0. Now we need to obtain similar bounds for (4.9), which applied to our case takes the expression given in (4.12).For the case 2 k ≥ Rε −1 we have As before, only if γ ≥ −d/2 we have the required bound.
Theorem 4 now implies that under scaling ε −d/2 the family (Φ ε ) ε>0 is tight in B α,loc p,q (U) for any α < −d/2, any q ∈ [1, ∞] and any p ≥ 2 and even.Using Lemma 1 this holds for any p ∈ [1, ∞].This way, the family is also tight in C α loc (U) for every α < −d/2.Remark 4. Observe that the scaling ε −d (the one used for the joint moments in Theorem 1) is outside the range of γ required for the tightness bounds, and therefore it will give a trivial scaling.

Proof of Proposition 3
The proof of this proposition will be divided into three parts.Firstly, we will determine the normalizing constant χ and show that it is well-defined, in the sense that it is a strictly positive finite constant.Secondly, recalling Definition 1, we will demonstrate that the n-th cumulant κ n ( Φ ε , f ) of each random variable Φ ε , f , f ∈ L 2 (U), vanishes for n ≥ 3 and finally that the second cumulant κ 2 ( Φ ε , f ), which is equal to the covariance, converges to the appropriate one corresponding to that of white noise.The ideas are partially inspired from Dürre (2009b, Section 3.6).
For the rest of this Subsection we will work with test functions f ∈ C ∞ c (U).The lifting of the results to every f ∈ L 2 (U) follows by a standard density argument (Janson, 1997, Chapter 1, Section 3).Let us first derive a convenient representation of the action Φ ε , f defined in Equation (2.3).More precisely, defining Φ ε , f S as where R ε (f) denotes the reminder term that goes to 0 in L 2 , as we show in the next lemma.
where A x := a ∈ U : ⌊a/ε⌋ = x .It is easy to see that |A x | ≤ ε d , and given that the support of f is compact and strictly contained in U, for ε sufficiently small (depending on f), the distance between this support and the boundary ∂U will be larger than √ dε.So there is no loss of generality if we assume that |A x | = ε d .Now, we can rewrite (4.14) as Let us call I(x) the term The set A x is not a Euclidean ball, but it has bounded eccentricity (see Stein and Shakarchi (2009, Corollary 1.7)).Therefore we can apply the Lebesgue differentiation theorem to claim that I(x) will be of order o(1), where the rate of convergence possibly depends on x and f.To see statement (4.14), we square the expression in (4.15) and take its expectation, obtaining where we used the Cauchy-Schwarz inequality.By Corollary 3 the expectation on the righthand side can be bounded as while the second term in (4.16) is of order o(ε −d ).With the outer factor ε 2d (4.16) goes to 0, as we wanted to show.
Let us remark that, by the previous lemma, proving finite-dimensional convergence of will be equivalent to proving finite-dimensional convergence of
Proof.Recall that by Corollary 1, for n ≥ 2, the n-th cumulant satisfies with D := supp f, which is compact inside U.The goal now is to show that First, we note from the cumulants expression (3.4) and bound (4.4) in Lemma 5 that, for any set V of (possibly repeated) points of D ε , with |V| = n, we have Using the above expression and Lemma 7, it is immediate to see .
We observe in particular that for d ≥ 2 this expression goes to 0 for any n ≥ 3. Furthermore, going back to (4.18), since f is uniformly bounded this shows that for n ≥ 3 the cumulants κ n go to 0 as ε → 0, as we wanted to show.

Discussion and open questions
In this paper we studied properties of the gradient squared of the discrete Gaussian free field on Z d such as k-point correlation functions, cumulants, conformal covariance in d = 2 and the scaling limit on a domain U ⊂ R d .One of the most striking result we have obtained is the "almost" permanental structure of our field contrasting the block determinantal structure of the height-one field of the Abelian sandpile studied in Dürre (2009a,b), Kassel and Wu (2015).We plan to investigate implications of these structures further in the future.
The other direction of research we would like to put forward is the existence of an anomalous scaling for "squared gradient fields".It is apparent from the combination of Theorem 1 and Theorem 2 that in order to obtain a limit for Φ ε , f one has to ignore the way correlations decay in the gradient squared of the DGFF (namely, like ε 2d /|v − w| 2d for v, w in the interior of U) and rather scale by the sum of all covariances in U ε , which is a factor of order ε d .The tightness result implies that this is the only possible non-trivial renormalization of this object.However, we conjecture that by inserting new parameters that make covariances decay more mildly we would obtain a non-trivial limiting object.
In fact, the idea of the proof for tightness in Proposition 2 is based on the application of a criterion by Furlan and Mourrat (2017) for local Hölder and Besov spaces.The proof requires a precise control of the summability of k-point functions, which is provided by Theorem 1 and explicit estimates for double derivatives of the Green's function in a domain.Observe that the proof is based only on the growth of sums of moments at different points.Thus this technique can be generalized to prove tightness of other fields just by having information on these bounds, which is usually easier to obtain than the whole expression on the joint moments.
Regarding the convergence of finite-dimensional distributions, Proposition 3, note that this strategy can be generalized to prove convergence to white noise of other families of fields, given the relatively mild conditions that we used from the field in question.Among them, one only requires knowledge on bounds of sums of joint cumulants, the existence of an infinite volume measure, and the finiteness of the susceptibility constant.Note that similar scaling results were given for random fields on the lattice satisfying the FKG inequality in Newman (1980).