STOCHASTIC HOMOGENIZATION OF GAUSSIAN FIELDS ON RANDOM MEDIA

In this article, we study stochastic homogenization of non-homogeneous Gaussian free fields Ξ D and bi-Laplacian fields Ξ b,a D . They can be characterized as follows: for f = δ the solution u of ∇ · a∇u = f , a is a uniformly elliptic random environment, is the covariance of Ξ D . When f is the white noise, the field Ξ b,a D can be viewed as the distributional solution of the same elliptic equation. Furthermore, we consider fields on domains D ⊂ R such that the boundary is smooth enough and on the discrete torus. Based on stochastic homogenization techniques applied to the eigenfunction basis of the Laplace operator ∆ in a smooth domain D, we will show that such families of fields converge to an appropriate multiple of the GFF resp. bi-Laplacian. The limiting fields are determined by their respective homogenized operator ā∆, where ā is a deterministic constant depending on the law of a. The proofs are based on the techniques introduced by [AKM19] and [GNO14].


Introduction
Heat flow or conductance of charge through materials with randomly placed impurities (porous media) are typically modelled using parabolic or elliptic equations in divergence form and diffusion matrices a modelling the environment, [Vaz07,CW84,PSQ16]. The coefficients modelling the porous media vary on a microscopic scale whereas on a macroscale the media shows effective behaviourā. Stochastic homogenization aims at identifyinḡ a and quantifying the convergence of the solution of the heterogeneous equation to the solution of an appropriate deterministic constant-coefficient equation. There is a rich mathematical theory developed in the last years, [AKM19, Tar09, TKO12, BD98, Jur80, Koz79, PV79, GNO14, BFFO17] just to mention a few.
In this article, we want to study two problems of stochastic homogenization of nonhomogeneous Gaussian fields, namely the non-homogeneous Gaussian free field (GFF) and bi-Laplacian field, in a domain D ⊂ R d and on the discrete torus. The (homogeneous) GFF Ξ g D and bi-Laplacian field Ξ b D are the two most prominent examples of the family of fractional Gaussian fields, [LSSW16]. The GFF has a wide range of connections including scaling limits of observables of interacting particle systems [Ken01], quantum field theory [Ber16] and invariant measures of Allan-Cahn type stochastic partial differential equations [Hai15]. The bi-Laplacian field or membrane model is related to scaling limits of uniform spanning trees [SW13], odometer functions of sandpile models [LMPU16,CHR18] and general interface models [Fun05].
For the first problem, consider the steady state solution u of the following parabolic equation ∂ t u = ∇ · a∇u + f in some domain D ⊂ R d , d ≥ 2, with regular boundary ∂D, f ∈ L 2 (R d ), and the randomly uniformly elliptic matrix-valued function a modelling the environment. The steady state solution u satisfies the elliptic equation (1.1) on D and u = 0 on ∂D. If f = δ is the Dirac delta function, the Gaussian field whose covariance function u satisfying (1.1) will be identified as the non-homogeneous GFF, denoted by Ξ g,a D . For f = ξ white noise the distributional solution will be referred to as the non-homogeneous bi-Laplacian field Ξ b,a D . In order to deal with the non-homogeneity of the problem (1.1), the idea of stochastic homogenization is to introduce a parameter 0 < ε 1 in order to represent the ratio of the microscopic and macroscopic scales. The equation (1.1) is then rescaled as − ∇ · a ε ∇u ε = f, where a ε (·) = a(·/ε). Furthermore, there will be an effective stationary fieldā, which depends on the law P of a, quantifying the convergence of u ε towards u as ε → 0, solution of − ∇ ·ā ∇u = f (1.3) in D ⊂ R d . The second model concerns a non-homogeneous GFF Ξ g,a D,N and bi-Laplacian field Ξ b,a D,N defined on the discrete torus T d N = [− 1 2 , 1 2 ) d ∩ 1 N Z d with random conductances a. The gradient and divergence operators in (1.2) and (1.3) are replaced respectively by their discrete versions and some averaging condition needs to be imposed for the model to be well-defined on the torus which has no boundary.
In particular, we will prove in Theorem 3.1 that Ξ g,aε as ε → 0, respectively in some appropriate Sobolev space. In the second result, Theorem 3.2, we will prove (1.4) with Ξ g,aε D replaced by Ξ g,a D,N and statement (1.5) with Ξ b,aε D replaced by Ξ b,a D,N as N → ∞. To the knowledge of the authors, these results are new and do not exist in the literature yet. By analyzing the stochastic homogenization of the bi-Laplacian, we are extending the techniques given in [AKM19,GNO14] not only to incorporate the randomness coming from f = ξ, but also the fact that the white-noise is not defined pointwise. Indeed, instead of ξ, we could take f to be a fixed distribution with a given negative Sobolev regularity.
Due to the influence of the random environment a, the convergence of u ε → u as ε → 0 will be proved on a Sobolev space of regularity 1/2 less, H −β−1/2 0 (D), than the optimal space for which u is well-defined, namely H −β 0 (D) for some appropriate β > 0. The main novelty of our proof is to write G ε , the Green's function, the solution of (1.2) for f = δ, in terms of {φ k } k≥0 , the orthonormal eigenfunctions basis of ∆. Using the fact that the operator ∇ · a ε ∇ in (1.2) is self-adjoint, we have that where ∇ · a ε ∇ ϕ ε k = −λ kā φ k , λ k is the eigenvalue of ∆ associated with φ k , andā is a constant that only depends on the law of a. With this simple expression, we can reduce most of our difficulties to applications of known results in stochastic homogenization of PDEs.
The structure of the manuscript is as follows. In Section 2 we will introduce the notation and define non-homogeneous Gaussian free fields and bi-Laplacian fields in regular domains and on the torus. Our main results are presented in Section 3, and their proofs are postponed to Section 4. Finally, in Section 5 we discuss some possible extensions.
2.1.1. Continuous spaces. Let D ⊂ R d , d ≥ 2, be either a convex bounded domain or a bounded domain such that ∂D ∈ C 1,1 (R d ). We will call such domains regular domains. Consider {φ k } k∈N be an orthonormal basis of L 2 (D) composed by eigenfunctions which are enumerated so that λ k is non-decreasing in k . Notice that by classical regularity theory, we have that φ k are all C ∞ c (D). We will denote by f, g the L 2 (D) inner product if f, g ∈ L 2 (D) and (by abuse of notation) the pairing between a distribution f and a smooth function g.
Notice that because we are restricting ourselves to mean-zero functions, the equation (2.3) indeed defines a norm.
We can also easily characterise the dual space of H β 0 (D). For β > 0, the space H −β 0 (D) is given by the continuous linear functionals on H β 0 (D). The weak- * topology in this space induces the norm given by (2.3), but with β substituted by −β. To highlight whether we are talking about a space of functions or a space of distributions, unless stated otherwise, we take β to be positive.
In the case D = T d , it will be more convenient to index the eigenfunctions/eigenvalues by k ∈ Z d . This is because the Fourier basis given by {φ k } k∈Z d , where φ k := exp(2πιk · x) is also an orthonormal basis of eigenvectors of the Laplacian in T d with periodic boundary conditions.
For strictly technical reasons, we will need to introduce another type of fractional Sobolev spaces denoted by W α,p (D). We will use these spaces in the proof of Lemma 4.2 where we rely on the bounds provided by Theorem 4.1 given in [AKM19]. We refer to Remark 1 for a comment on the relation between the two types of Sobolev spaces.
Definition 2.2. The fractional Sobolev space W β,p (D) for any β > 0, p ≥ 1 and d ≥ 2 is defined by · stands for the floor function and is called the Gagliardo semi-norm for δ ∈ (0, 1).
Remark 1. Notice that it only makes sense to compare the two types of spaces in the case p = 2, as this is the only case for which W β,p (D) is a Hilbert space. Observe that in the proof of Lemma 4.2 we only consider p > 2. The W β,2 (D)-norm of a mean zero function f can be written in terms of the L 2 (D) norm of the singular integral definition of the fractional Laplacian operator restricted to D. On the other hand, the H β 0 (D)-norm of f is the L 2 (D) norm of the eigenvalue/spectral definition of fractional Laplacian. These spaces do not coincide in general, unless β ∈ N or if we take the domain D to be either R d or T d . In the former case, the fractional part of the W β,2 (D)-norm is omitted, whereas in the latter we can express singular integrals in terms of Fourier transforms.
Finally, for two non-negative functions f, g we write if there exists a positive constant C = C(γ 1 , . . . , γ k ) > 0 such that for all x to be understood as a graph with a periodic boundary. In this context, we use E d N to denote the edges of the graph Z d N . For a function f : T d N −→ C, we will denote its restriction by T d N as is an orthonormal basis of the space 2 (T d N ) := C T d N with respect to the inner product When discussing scaling limits of discrete models, we will always use x, y to denote points of Z d or Z d N , and usex,ŷ to denote points in T d N . 2.2. Gaussian fields. In this section we will introduce the Gaussian free field and bi-Laplacian field in continuous domains, D ⊂ R d . First we will consider the homogeneous and then the random media setting. Throughout this article, we distinguish between two scenarios for D, in which we either take D to be a regular domain (as defined at the beginning of the previous Section 2.1), or D = T d .
2.2.1. Gaussian fields in the continuum.

Homogeneous fields.
Definition 2.3. The Gaussian free field Ξ g D on D is a multivariate Gaussian field defined by where ∆ is the Laplacian operator and δ is the Dirac delta function. For D = T d we substitute the assumption on the boundary condition in (2.5) by T d G T d (x, y)dy = 0. A common representation of G D in spectral terms is the following: (2.6) We will use a similar representation for the non-homogeneous case, as it will allow to quantify the distance between the fields we are interested in. Note that for any given test function f we can write with δ k,k representing the Kronecker delta function, and k, k ∈ N. Let ξ denote the (spatial) white noise in L 2 (R d ): In the following definition we follow the approach in [LSSW16].
where ξ denotes the white noise defined in (2.8).
Let us remark that, for ξ(k) := ξ, φ k , we have that satisfies the distributional equation (2.9). Moreover, if D ⊂ R d is a domain due to Weyl's law which states that λ k D,d k 2/d (see [Lax02, Appendix B] for a proof), we have that Ξ b D belongs to H β 0 (D) for all β < − d 4 + 1. Similarly, we can prove that Ξ g D belongs to H β 0 (D) for all β < − d 4 + 1 2 . Notice that for sufficiently small dimensions, the fields can be interpreted as functions in positive Sobolev spaces.
Gaussian fields in random media. Let us introduce the random environment a and present some assumptions. We will follow the approach in [AKM19]. The probability space of the environment will be denoted by (Ω(Λ),F, P), which we will define in the sequel.
Call R d×d sym the set of real valued symmetric matrices. We will be interested in maps We then define the space Ω(Λ) := {a ∈ R d×d sym : a Lebesgue measurable and satisfying (2.11)}. (2.12) We will consider F U the σ-algebra generated by the maps . Let P be a probability distribution on (Ω(Λ),F) with the following properties: is the Hausdorff distance. (A3) Isotropic in law for each I : R d −→ R d isometry such that I maps each of the coordinate axes to another, we have that P • I −1 = P.
Let us remark that Assumptions (A1) and (A2) are standard in [AKM19], but (A2) can be relaxed to sufficiently fast decaying covariances. Assumption (A3) is a matter of convenience, in this case, the homogenized operator is given by a multiple of the Laplacian (instead of ∇ ·ā ∇ whereā is a matrix). In this case, the asymptotic behaviour of the L p (R d ) norms of the φ k (the eigenvalues of the limiting operator) are well understood.
Definition 2.5. Let (Ω(Λ) ⊗ Ω,F ⊗ F, P ⊗ P), D ⊂ R d a regular domain, and a ∈ Ω(Λ) fixed. Furthermore, let G a D (·, ·) be the solution of We define Ξ g,a D a non-homogeneous Gaussian free field as the Gaussian random field with mean 0 and covariance given by Definition 2.6. Let (Ω(Λ) ⊗ Ω,F ⊗ F, P ⊗ P), D ⊂ R d a regular domain, and a ∈ Ω(Λ) fixed. We define Ξ b,a D a non-homogeneous bi-Laplacian as the distributional solution to (2.14) For a given random environment a ∈ Ω(Λ), stochastic homogenization techniques will allow determining limiting fields of Ξ g,ε D := Ξ g,aε The main idea of this article is to obtain a good representation for the non-homogeneous fields not in terms of the eigenfunctions of ∇ · a ε ∇ but instead, in terms of the eigenfunctions of the liming operator. For this, we will need to introduce a sequence of functions that we call the pseudo-eigenfunctions The deterministic constantā is positive and has a special meaning in stochastic homogenization where it is referred to as effective/homogenized coefficient, see Remark 2 for more details.
In particular, we will show that the non-homogeneous GFF Ξ g,ε D can be written as a mean zero infinite Gaussian vector Similarly, for the non-homogeneous bi-Laplacian field whereξ ε (k) = ξ, ϕ ε k and ϕ ε k is defined in (2.15). Remark 2 (The effective coeffientā ). The effective coefficientā only depends on the law of the environment P. In particular in [AKM19,Equation 3.92], the authors provide an explicit formula forā in terms of the expected energy and flux of the first-order correctors Φ. That is, for any unit vector e ∈ R d : and the first-order corrector Φ e is the difference between the solution of the Dirichlet problem and an affine function defined in [AKM19, Section 3.4]. A similar characterisation holds in the discrete setting (which we introduce in the next section). For more information see for instance [GO11].
2.2.2. Gaussian fields in discrete spaces.
wherex ∼ŷ denotes thatx andŷ are nearest neighbours in T d N . Note that we will use the same definition of the normalized discrete Laplacian as [GNO14]. In this case, we take the Green's function as the solution of It will be particularly useful for us that Moreover, we have that the Green's function of the discrete torus T d N can be written as is the eigenvalue of ∆ N associated with φ N k . Definition 2.7. The discrete Gaussian free field on the torus is a multivariate Gaussian vector N (0, 1).
is the solution of the finite difference equation Gaussian fields with random conductances. In this context, we will use the stochastic homogenization bounds from [GNO14] instead of [AKM19]. We will keep some of the same notation used in the continuous case in order to make a clear analogy between the two. The underlying probability space will be denoted by (Ω dis (Λ),F, P). Similarly to the infinite case, we define (Ω dis N (Λ),F N ) by substituting Z d for Z d N and E d by E N . Let P be any product probability measure on Ω dis (Λ).
We define a projection operator Π N : where (e i ) d i=1 is the canonical basis of R d , a ∈ Ω dis (Λ) and e ∈ E d N . On the other hand, we define an extension operator from Ω dis where e * is the only edge in E d N such that e * = {x * , y * } with x * ≡ x mod Z d N and y * ≡ y mod Z d N . For a ∈ Ω dis (Λ), we write a N := Ext N •Π N (a). Notice that, Π N • Ext N is the identity in Ω dis N (Λ) and Ext N •Π N (a) e = a e for all a ∈ Ω dis (Λ) whenever e = {x, y} with x, y ∈ Z d N . Let P N be short for P(Π −1 N ·), which is a product measure in Z d N . Finally, for any f : T d N −→ R and any a ∈ Ω(Λ), we can write Definition 2.9. Let (Ω dis (Λ) ⊗ Ω,F ⊗ F, P ⊗ P). We define the discrete non-homogeneous Gaussian free field to be the Gaussian vector (Ξ g,a N (x))x ∈T d N with mean 0 and covariance and δ N x,ŷ = z∈Z d δ Nx,Nŷ+z , with δx ,ŷ being the standard Kronecker delta function. Call (ξ N (x))x ∈T d N the i.i.d collection of standard normal random variables.
is called discrete non-homogeneous bi-Laplacian field. Again, we are using the notation (ξ) N := 1 The formal field on D = T d as for i ∈ {b, g} is defined by where Ξ g,a N is defined in Definition 2.9 and Ξ b,a N in Definition 2.10 and c g = (2d) −1/2 resp. c b = (2d) −1 .
Remark 3. In (2.26), there is a clear abuse of notation as we are using Ξ i,a N (·) to denote both a random vector (which has well-defined values for each choice ofẑ) and for a distribution given by the sum of delta functions (and therefore has no well-defined notion of value at a given point). We chose to do so to simplify the notation.

Main results
Theorem 3.1. Let D ⊂ R d be a regular domain, (Ω(Λ) ⊗ Ω,F ⊗ F, P ⊗ P) the underlying probability space, where P satisfies Assumptions (A1)-(A3). Then there exists a positive deterministic constantā such that the following statements hold.
(1) The non-homogeneous GFF defined in (2.16) satisfies (2) The non-homogeneous bi-Laplacian field defined in (2.19) satisfies where P is a product measure. Then there exists a positive deterministic constantā such that (1) The discrete non-homogeneous GFF defined in (2.26) satisfies (2) There exists a coupling between ξ and the sequence ξ N , such that the discrete nonhomogeneous bi-Laplacian defined in (2.26) satisfies as N → ∞. The convergence holds in H −β 0 (D), β > d 4 − 1 2 . For more details onā in both contexts, see Remark 2. Note that in Theorem 3.2 we only take the limit in N → ∞ to obtain the homogenized limiting field. This is due to the translation invariant nature of the environment a.
Finally, let us mention that one could easily prove that non-homogeneous bi-Laplacian fields converge almost surely (in more irregular Sobolev spaces than the ones stated above) to their homogeneous counterpart by a simple Borel-Cantelli argument.

Proofs
4.1. Proof of Theorem 3.1. We start this section by stating some necessary results regarding stochastic homogenization. Then we proceed to prove some useful lemmata necessary to recover our representation of the Green's functions before finally proving the main theorems.
Under the Assumptions (A1)-(A3), there are good bounds for stochastic homogenization of the solutions of elliptic partial differential equations. Indeed, one can prove that there is a positive constantā (see Remark 2), such that for u, any sufficiently regular function, we can estimate its L 2 (D) distance to the function u ε given by − ∇ · a ε ∇ u ε =ā ∆u in D, u ε = 0 in ∂D. (4.1) We will state a simplified version of a result from [AKM19] which will be one of the main ingredients for our proofs. For this, we use the stochastic integrability notation O s (·) according to the law P. Given a random variable X and parameters s, θ ∈ (0, ∞), we say that where a + := max{0, a}.
Proof. Fix α ∈ (0, 1], κ > 0 and let p = 2+δ with δ > 0 to be defined later to be sufficiently small. By applying Theorem 4.1 to u = φ k , where φ k is the eigenfunction associated to λ k the k-th eigenvalue of (−∆) we have that for p ∈ (2, ∞), (4.6) By [BL12, Theorem 6.4.5], for any given a function f ∈ C ∞ c (R d ), s 1 , s 2 ∈ R, p 1 , p 2 > 1 and θ ∈ (0, 1), we have as long as This can also be extended to any such function f ∈ W s θ ,p θ (D). Now, if we choose s θ = 1 + α, p θ = 2 + δ, s 1 = p 1 = 2 and s 2 = 0, we have that θ = (1 − α)/2 and Notice that as δ −→ 0 + , we have that p 2 −→ 2 + for any α ∈ (0, 1]. Furthermore, by classical results any f ∈ C ∞ c (D), such that f | ∂D ≡ 0 can be naturally extended to so that Ext f L p 2 (R d ) = f L p 2 (D) and Ext f W 2,2 (R d ) = f W 2,2 (D) . As each of the φ k are in C ∞ c (D)∩H 2,2 0 (D), we can use (4.7), and the isometry properties above to get where in the last inequality we estimated the W 2,2 (D)-norm of φ k by the L 2 (D)-norm of ∆φ k = −λ k φ k and that ||φ k || L 2 (D) = 1. Now, we use [Gri02, Theorem 1], which states that Using the interpolation version of Hölder inequality, the bound (4.9), and the fact that φ k L 2 (D) = 1, we have that (4.10) The proof is completed by choosing δ small enough so that Consider G ε D , the fundamental solution of ∇ · a ε ∇, that is Classical results of partial differential equations provide some useful properties for G ε , for instance the symmetry of Green's function, that is G D (x, y) = G D (y, x), see [GW82, Theorem 1.3] for the case d ≥ 3 and [DM95, Section 6] for the case d = 2.
The following lemma provides a representation for G ε D . It introduces a natural way of comparing it to the Green's function ofā ∆, and gives the main intuition behind our results. Remark that this lemma justifies the representation (2.19) of the non-homogeneous bi-Laplacian fields.
Proof. Remember that for any a ∈ Ω(Λ), the Green's function is pointwise well-defined in D := D × D \ {(x, x) : x ∈ D}. Call for simplicity G ε x := G ε D (x, ·) for some fixed x ∈ D. We have that for d ≥ 3, G ε x ∈ W 1,1 0 (D) (see [GW82]) and for d = 2 we have that G ε x ∈ L 1 (D) (see [TKB13]). In either case, we can see G ε x as a tempered distribution and since φ k ∈ C ∞ c (D), we get that G ε x , φ k is well-defined.
In order to calculateĜ ε x (k) := G ε x , φ k , note that as G ε is the Green's function of ∇ · a ε ∇ and the operator ∇ · a ε ∇ is self-adjoint. Therefore we have that Hence, Finally, by integrating in x, we have that as long as β > d 4 − 1 2 , using that α and κ can be made arbitrarily small and Weyl's law. This also implies that P-a.s and almost all x according to the Lebesgue, G ε Stochastic homogenization of the GFF. For this proof we will use a truncation argument and apply the following theorem. for all τ > 0, and X n,K We will apply the theorem by choosing S = H −β 0 (D), X = Ξ g D , X n = Ξ g,ε D and µ = P⊗P. Define the following truncations X n,K := Ξ g,ε K and Z K := Ξ g K , centered Gaussian fields such that By truncating the covariance of Z K and X n,K we essentially reduce the problem to convergence of finite dimensional Gaussian vectors. Since Ξ g,ε K and Ξ g K are determined by their covariance structure, the convergence in distribution result will follow from proving that their covariance matrices converge. This will follow from Lemma 4.2 since we truncated for the k's to be in the set {1, . . . , K}.
In the remainder we demonstrate that (4.13) holds. Define the error field by (4.14) Using the definition of the norm H −β 0 (D) and the monotone convergence theorem we have Using the triangular inequality first and then Lemma 4.2, we get that the sum above is bounded above (up to a multiplicative constant) by which is finite as long as β > d 4 , as we can choose κ, α to be arbitrarily small together with Weyl's law. We can then use Markov's inequality with the product measure P ⊗ P, to get that Applying Theorem 4.4, we conclude the proof.
4.2. Proof of Theorem 3.2. The strategy to prove the convergence result on the discrete torus will follow similar ideas as in the continuum. Indeed, we can expand the Green's function for ∇ N · a N ∇ N in terms of the eigenfunctions of the discrete Laplacian. The advantage in this case is that the eigenfunctions are precisely the Fourier basis of 2 (T d N ). Discrete Fourier analysis was also the main tool of previous works on scaling limits of odometer fields, see [CJR21,CdGR19]. In the case of random environments note that the Fourier basis ceases to be a basis of eigenfunctions for the operator ∇ N · a N ∇ N . Again, the key idea is to show that, for large N , the Fourier basis will be close to the basis of eigenvectors. Let The next theorem states good bounds for the difference between f N andf N , the function in 2 (T d N ) solving the finite difference equation whereā is a positive constant that only depends on the law of a. The next theorem is going to serve as the equivalent of Theorem 4.1 for the discrete context and is a simple adaptation of [GNO14, Corollary 1.2].
Theorem 4.5 (Corollary 1.2, [GNO14]). There exists a deterministic positive constant a only depending on the law of a and d with the following property. Given N ≥ 1 with f ∈ C ∞ (T d ), let f N be the solution of (4.18) andf N be the solution of (4.19).
We will be interested in the case g N = −ā λ are the eigenvalues of the normalized discrete Laplacian operator. Very conveniently, in this casef N = φ N k . In the following lemma, we will denote by ϕ N k the solution to (4.18). Now, let us state the representation of the Green's functions G N,a defined in (2.24).
Lemma 4.6. For all a ∈ Ω(Λ) and any N ≥ 1, we have (4.21) for allx,ŷ ∈ T d N . This proof is similar to the proof of Lemma 4.3 and will be therefore omitted. It is much simpler as the sum only has a finite number of terms. The next lemma is just an application of Theorem 4.5.
Lemma 4.7. Stochastic homogenization of the discrete non-homogeneous GFF. The proof is very similar to the proof of Theorem 3.1 (1) hence we will point out the main differences. We use Theorem 4.4 with S = H −β 0 (T d ), with β > d/4 and X := Ξ g D , X N,K := Ξ g,a N,K and Z K :=ā −1/2 Ξ g K , where Ξ g,a N,K := and a ∧ b := min{a, b}, Ξ g,a D,N was defined in (2.26) and (4.23) The proof of the argument follows similarly but using Lemma 4.7 instead of Lemma 4.2 to prove the convergence of Ξ g,a N,K toā −1/2 Ξ g D,K .
Stochastic Homogenisation of the discrete non-homogeneous bi-Laplacian. Let τ > 0, we want to prove that lim N →∞ . For the discrete case, we compare the non-homogeneous field to a discrete homogeneous one. That is, we will first show that Ξ b D,N (the homogeneous formal field) converges to Ξ b D in probability according to an appropriate Sobolev norm, as long as we choose a suitable coupling between ξ and ξ N in (2.23). For this, we will take ξ N ( where ξ is the same sample of the noise used in the definition of (2.9), 1l A denotes the indicator function of the set A and B N ( in this same Sobolev space using estimates given in [GNO14].
We will prove the first point described above in the next proposition. (4.24) Proof. Again, it will be convenient to study the action of the field on the basis of eigenfunctions of the Laplacian. Remember that, as we are in the torus (both in the discrete and the continuous cases), such basis of eigenfunctions is given by the Fourier basis φ k := exp(2πιk · x). Let us start by calculating explicitly Ξ b N (φ k ) for k ∈ Z d N , where in the first identity we used that ẑ∈T d N G N (ẑ,ŷ) = 0 for anyŷ. A simple computation shows that (see [CHR18,Lemma 7 On the other hand, using the expansion of the Green's function G T d on the torus, we we get that for any β, we have proving (4.24).
We now proceed to show that Ξ b,N Err converges in probability to 0 in the space and remark that here we are committing the same abuse of notation as described in Remark 3. Let us start by computing Ξ b,N Err (φ k ), for fixed k ∈ Z d .
Expanding the previous expression we get By using Lemma 4.7, we get Finally, let β > d 4 − 1 2 and choose κ > 0 sufficiently small so that β > d 4 − 1 2 + κ 2 , then we have which vanishes as N goes to infinity.

Discussion
In this section we provide a few remarks regarding possible generalization.
Convergence in probability of the non-homogeneous GFF. A natural question is whether the convergence of the (non-homogeneous) GFF can be improved to a convergence in probability. For this, we need a coupling in terms of a elliptic PDE which allows us to employ stochastic homogenization techniques. Indeed, there is a natural coupling by writing the desired fields as the solution to PDEs similar to the one found in [GM17], in which the authors define the notion of generalized Gaussian free field.
Unfortunately, to be able to extract results from such coupling, we would need to obtain a bound similar to the one found in Lemma 4.2 but for the quantity b ε ∇ϕ ε k −b∇φ k L 2 (D) where b ε := √ a ε andb := √ā . This is not expected to converge as ε vanishes. Indeed, one can see in [AKM19,GNO14] that ϕ ε k needs the so called first-order correctors in order to converge to its homogenized counter part in H 1 0 (D). Discretized domains. In this article, we derive scaling limits in the continuous setting in a domain and for the discrete setting on a torus. However, we do believe that such results could be extended to discretized domains. The proof would require both adapting the techniques of quantitative discrete stochastic homogenization near the boundary. One would also need to account for the fact that the discretized eigenfunctions of the Laplacian cease to be the same as the eigenfunctions of the discretized Laplacian. One should be able to bound such deterministic errors in such context, similarly to [BIK15] where they bound such error for manifolds without boundary.
General fractional fields. We focused in the case of the GFF and bi-Laplacian fields, one could wonder what happens for general fractional Gaussian fields, see [LSSW16]. This seem much harder, as one would need to deal with stochastic homogenization for non-local operators, which seems far from the scope of the current available theory of quantitative stochastic homogenization.