characterisation the continuum Gaussian in

. — We prove that under certain mild moment and continuity assumptions, the d -dimensional continuum Gaussian free ﬁeld is the only stochastic process satisfying the usual domain Markov property and a scaling assumption. Our proof is based on a decomposition of the underlying functional space in terms of radial processes and spherical harmonics. Résumé (Une caractérisation du champ libre gaussien dans le continu en toute dimension) Nous montrons que, sous de faibles hypothèses de moment et de continuité, le champ libre gaussien dans le continu à d dimensions est le seul processus stochastique satisfaisant à la propriété habituelle de Markov sur le domaine et une propriété d’échelle. Notre preuve est basée sur une décomposition de l’espace fonctionnel sous-jacent en termes de processus radiaux et d’harmoniques sphériques.


Introduction
The continuum Gaussian free field (GFF) is a generalisation of Brownian motion, taking the index set to higher dimensions, which is defined on the unit ball as follows.
where G B denotes the zero boundary Green's function for the Laplacian in B. Note in particular that G B is identically zero in R d B, so that (h B , f ) is almost surely zero as soon as f has support outside of B.
In dimension d = 1, this corresponds to a slightly awkward definition of the Brownian bridge on the interval (−1, 1): in this concrete one-dimensional case the GFF can actually be defined as an a.s.continuous function, i.e., indexed by points instead of test functions.In d 2 the GFF does not make sense as a pointwise defined function, thus the more general definition.However, one can still heuristically see the continuum GFF as a Gaussian height function parameterised by a domain D ⊆ R d [She07].The Gaussian free field has played many roles in probability and mathematical physics since the 1970s: it is the stationary solution of the stochastic heat equation; is used to describe both free and interactive Euclidean Quantum field theories; it appears as a corrector in stochastic homogenisation; to name just a few.Recently, the 2d GFF has played a crucial role in the probabilistic study of statistical physics models, Schramm-Loewner evolution, Liouville quantum gravity and Liouville field theory [Dub09,MS16,DS11,DKRV16].It is both the proven and the conjectured scaling limit of several natural discrete height functions [NS97, Ken01, RV07,BLR20] There are many characterisations of Brownian motion that single it out among 1d stochastic processes and it is similarly natural to also ask how the GFF in higher dimensions can be characterised.Recently, the 2d continuum GFF was characterised as the only conformally invariant field, satisfying a domain Markov property and certain minimal moment and continuity assumptions [BPR20,BPR21].
This short article provides a first characterisation of the d-dimensional continuum Gaussian free field for d 3, answering Question 6.1 in [BPR20].It also provides a new characterisation of the GFF in d = 2 that relaxes the assumption of conformal invariance and thereby generalises the main result of [BPR20] via a new more elementary proof.Our characterisation is based only on scaling and the domain Markov property and in particular, we do not even need to assume rotational invariance.
More precisely, but still somewhat informally, we prove that, under certain mild moment and continuity assumptions, the d-dimensional Gaussian free field is the only random Schwartz distribution in d 2 that satisfies the following domain Markov property -for each ball, the field inside can be written as a sum of an independent scaled and translated copy of the original field plus the "harmonic extension" of its behaviour on the boundary.We prove the theorem by identifying the covariance structure as the Green's kernel (lighter step), and then proving Gaussianity (more involved).The second step in particular uses a completely different argument to that in [BPR20] -instead of working with circle averages, we use a decomposition of L 2 (B) based on spherical harmonics.The two steps in our proof are quite independent of each other (although we use the knowledge of the covariance kernel to simplify some arguments later on), but both make strong use of the domain Markov property.Our proof, which does not even assume rotational invariance, stresses the nice interplay between the GFF and harmonic functions.Recently, there has been a lot of interest in proving scaling limit results for random surface models, e.g., [DCHL + 19, BLR19,GM21].The characterisation given in this article and its mild generalisations below could potentially be helpful for obtaining such results, by giving a new way to identify the GFF as a continuum scaling limit.Importantly, one that does not rely on conformal nor rotational invariance -conditions that are both very difficult to prove from lattice considerations (though there is important recent progress [DCKK + 20]).
As mentioned, essentially the only property that is needed to characterise the continuum Gaussian free field is its domain Markov property.It should be emphasised that, as stated, this domain Markov property also encapsulates a scaling property for the field (that could however be disposed of to some extent, see comments below).
) and zero as soon as f has support outside of B. We define the random Schwartz distribution h a+rB on domains a+rB as the rescaled image of h B under translation by a and scaling (1) by r: We then say that h satisfies the domain Markov property if the following holds: -Suppose that a + rB ⊆ B. Then we can write where the two summands are independent, ϕ a+rB B is a stochastic process that a.s.corresponds to integrating against a harmonic function when restricted to a + rB and h a+rB B is equal in law to a translated and scaled copy h a+rB .
We can now state the characterisation: ) and zero as soon as f has support outside of B.
If h satisfies the following conditions, then for some constant c the stochastic process We have not striven for most general technical assumptions, but rather have tried to keep the proofs light and self-contained.However, we believe that several assumptions can most likely be relaxed.Let us remark on those and the case d = 1.
-In fact, a similar characterisation holds also in d = 1 with basically the same proof.In this case one would just work directly in the space of continuous functions, and state the zero boundary condition pointwise.
-We believe that the moment assumption can be probably relaxed using methods of [BPR21].
-Although using spherical harmonics is key in our proofs, one can also adapt the method to work for other reasonable domains, e.g. when one considers boxes instead of spheres both as the original domain and for the DMP.Indeed, in this case one can still make sense of spherical harmonics defined on the boundary of the box, defined with respect to harmonic measure on the boundary as seen from the centre of the box.A characterisation of harmonic functions using such "box averages" will then give the required harmonicity for the covariance in the first part of the proof and everything in the second part will work the same way.More generally, a similar strategy ought to work for smooth convex domains, where there is a point of homothety inside the domain.The proof would, however, become more technical, and hence is not pursued here.
-We use the same form of domain Markov property (DMP) as in [BPR20], but only for balls.The DMP plays a key role in our argument -it is used both in showing that the covariance kernel is harmonic off the diagonal, and in proving Gaussianity.With the current approach it seems difficult to relax the harmonicity condition for the extension of the boundary data, but it would be very interesting to determine if this is possible.
-On the other hand, one can easily envisage replacing the condition of having an exact (scaled) copy of the field in the DMP with a more relaxed, martingale type of condition.Namely, we could ask that h B inside any sub-ball a + rB ⊂ B, can be decomposed as the independent sum of an a.s.harmonic function, and a random Schwartz distribution that has mean zero when tested against any test function, with some conditions on the second moment.Indeed, this is known to suffice in the case of Brownian motion.We believe that our proof would adapt, as long as we imposed certain growth conditions for the variance of averages around a point, and certain homogeneity conditions.For the sake of keeping things simple and short, we will not pursue this direction here.
-It might seem a bit surprising that we are not using rotational invariance.However, the reason is that the DMP already implies rather easily that the covariance is a harmonic function, and this is a very strong property.
-Rather than assuming that h is almost surely a Schwartz distribution, we could instead assume that: (1) . Indeed, (2) with respect to the usual topology on C ∞ c (R d ) in which fn → f iff all derivatives of fn converge to the corresponding derivative of f , uniformly on the closure of supp(f ).
J.É.P. -M., 2022, tome 9 we will use linearity of the process (h B , f ) repeatedly, but will only use the fact that f → (h B , f ) is a Schwarz distribution for justifying (2.10) below, which would actually follow more directly had we assumed (2).
-Finally, let us comment on the exact form of the zero boundary condition.In [BPR20] this convergence was only asked for rotationally symmetric functions.We need the slightly generalised form (that is implied for example by conformal invariance and the rotationally invariant form) in only two places: to determine that the covariance is the Green's kernel and to obtain the uniqueness of the domain Markov property.Notice that for pointwise defined functions our condition is weaker than the usual pointwise zero boundary condition.
Remark 4. -When d = 2 this is indeed a generalisation of the result in [BPR20].Namely, in [BPR20] the authors assume that the law of a field h D is given for every simply connected domain D, and the family of laws satisfy conformal invariance (if ϕ : D → D is conformal then the law of the pushforward of h D by ϕ, as a generalised function, is equal to that of h D ) and a more general domain Markov property (whenever D ⊂ D we can write h D = h D D + ϕ D D as an independent sum, with ϕ D D harmonic in D and h D D equal in law to h D ).In particular, under these assumptions, h B does satisfy (A) of Theorem 3.
Therefore, if we take the assumptions of [BPR20] as an input, (3) our result gives that h B is equal in law to some multiple of the zero boundary Gaussian free field in B. Then the assumption of conformal invariance identifies the law of h D for arbitrary simply connected D, and we reach the same conclusion as [BPR20].
We will now present the argument in three sections.First, we discuss some immediate consequences of the assumptions (along similar lines to [BPR20], so we will keep this brief).In Section 3, using the DMP, scaling and translation invariance, we show how to deduce that the covariance kernel is the Green's kernel (Proposition 6).Finally, the main part of the paper is Section 4, where we prove Gaussianity (Proposition 14) -we do this using solely the DMP, and a decomposition of the underlying functional space using spherical harmonics.Theorem 3 is an immediate consequence of Proposition 6 and Proposition 14.Sections 3 and 4 can be read quite independently of each other, although we use the identification of the covariance kernel to simplify some arguments in the latter.

Definition 5 (Scaling function).
-In what follows we write s(r) for the function on (0, ∞) defined by − log r when d = 2 and r 2−d when d 3.
Acknowledgements.-We would like to thank here H. Duminil-Copin and V. Tassion for asking us whether only rotational invariance should suffice for the characterisation of the 2d continuum GFF.This discussion happened in a nice and friendly workshop in Fribourg, organised by I. Manolescu, who we would hereby like to thank too.We would (3) also using the assumptions of [BPR20] on K 2 , which would equivalently work for our proof (as explained in the bullet point above).
finally like to thank N. Berestycki, J.C. Mourrat, G. Ray and W. Werner for many interesting and fruitful discussions on this topic, and two anonymous referees who pointed out several places where more clarity was required.

Immediate consequences
Here we discuss some immediate properties of our assumptions.The section is selfcontained, but we remain brief, as similar properties have been shown in detail in the 2d case in [BPR20].
The domain Markov decomposition is unique.-Indeed, suppose that for some a, r we had two decompositions as in (A) of Theorem 3: Then for any z ∈ a + rB, by harmonicity, there exists a sequence (f n ) n 0 of functions satisfying the zero boundary condition assumption (C) from Theorem 3, and such that (ϕ, f n (a + r•)) = ϕ(z) for all n and any harmonic ϕ in a + rB.The zero boundary condition then implies that (h a+rB and h i is independent of {ϕ, (h j ) j =i } for each i.Indeed, for any fixed i we can write ϕ = ϕ i − j =i h j and since h j is zero in B i for all j = i, ϕ = ϕ i is harmonic in B i .But also for every j = i, h i is zero in B j and therefore h j is measurable with respect to ϕ i .So the collection {ϕ, (h j ) according to the domain Markov decomposition of h B in z 1 + εB, z 2 + εB respectively.Then by (2.1) and (2.2), the quantity is well-defined, i.e., it does not depend on ε > 0 satisfying the above conditions.Note that, by harmonicity, we can alternatively write J.É.P. -M., 2022, tome 9 for any smooth functions η z1 = η ε z1 , η z2 = η ε z2 of mass one, that are supported in z 1 + εB, z 2 + εB and rotationally symmetric about z 1 , z 2 respectively.
For every δ > 0, the following crude upper bound for To justify this, observe that by Cauchy-Schwarz and the domain Markov property, it suffices to show that for any z ∈ (1 − δ)B and ε < δ Indeed, as long as 1) , then the domain Markov property implies that E(ϕ εB Further, by the domain Markov property again and also scaling, the right hand side can be written as the sum ) is well-defined, i.e., it does not depend on the choice of a i satisfying the above conditions.Again we can write this as E 4 i=1 (h B , η zi ) for any smooth functions (η zi ) 1 i 4 , with η zi having mass one, being rotationally symmetric about z i and supported in z i + a i B for i = 1, 2, 3, 4.
Using the same argument as in the 2-point case, one can now bound |k 4 (z 1 , . . ., z 4 )|.Indeed, by Hölder's inequality it suffices this time to obtain bounds on the fourth moments of ϕ z+εB B (z) for ε > 0. This can be done very similarly to the 2-point bound.One first treats z = 0 using the i.i.d nature of increments -the only change here is that when expanding E ( ) -and then transports this bound to general z using translation invariance and scaling.The result is that for any given δ > 0, there exists for any z ∈ (1−δ)B and ε < δ (we omit the details).This leads to the final conclusion: there are constants C(δ) > 0 such that for any distinct The diagonal contribution of the 2-point function is negligible.-We can now show that k 2 is the covariance function of the field, in the sense that for any where the right-hand side is well defined as the limit (see just below for a proof).In other words, we show that K 2 is an integral kernel with density k 2 , and that there is no "diagonal contribution" of K 2 , so it does not matter than k 2 is only defined off the diagonal.This already rules out, for example, the possibility that the field is a white noise or one of its derivatives. (4) Let us now justify (2).For each w ∈ B and 0 < ε < d(w, ∂B), let η ε w be a unit-mass radially symmetric mollifier in the ε ball around w, as in (2.3).Then if f 1 , f 2 ∈ C ∞ c (B) (we fix these from now on), the functions The assumption that h B is almost surely a Schwartz distribution therefore implies that (2.9) where the limit inside the expectation is almost sure.Moreover, (h where by (2.6), the domain Markov property and scaling, we can bound We would like to thank A. Sepúlveda for raising the need for extra clarity on this point.
J.É.P. -M., 2022, tome 9 for C uniform over ε (small enough) and z in supp(f 1 ) ∪ supp(f 2 ).Thus the family is uniformly integrable and we can take the limit outside of the expectation in (2.9), i.e., we have (2.10) We are now in a good position to prove (2), but still need to deal with the contribution to K 2 from points "near the diagonal".For this, we break up the right-hand side of (2.10) as where by (2.3), we have for each ε by Cauchy-Schwarz, where by (2.5), the domain Markov property and scaling, we can bound E((h for some constant C. Note that this constant C does not depend on ε, and the bound holds for all z in the compact supports of f 1 and f 2 simultaneously.This implies that lim sup and hence that (z) to be the ε-spherical average of h B around z.The spherical average admits the following natural approximations.Write ρ ε n for a sequence of smooth test functions with total mass one that are rotationally symmetric about z and supported in the annular region for each n by the Markov property, where the second term on the right-hand side goes to 0 as n → ∞ by the zero boundary condition assumption.It therefore follows that and moreover this convergence is uniform for, say, ε ∈ (δ, 1).
Using the same construction as above with z = 0 and ε = r n converging to 1 as n → ∞, we also see that lim n→∞ E(h rn (0) 2 ) = lim n→∞ E(h B , ρ rn n ) 2 ), and this final limit is 0, again by the zero boundary condition.Since E(h r (0) 2 ) is decreasing as r ↑ 1, this implies that (2.11) E(h r (0) 2 ) 0 as r ↑ 1.

Covariance is the Green function
Let us start by showing that scaling and translation invariance together with the domain Markov property already imply that the covariance kernel is the Green's function: Proposition 6. -The function k 2 (x, y) (defined for x = y) is a positive multiple of the zero boundary Green's function.In particular, by (2), this implies that the bilinear form K 2 is a multiple of the map for each w ∈ ∂(x + ηB).The conclusion is that for any given η > 0, and any x ∈ B such that |x − y| ∧ d(x, ∂B) > η: Note that the support of ρ η x does not include y by assumption, so the integral against k 2 (•, y) is perfectly well defined.(The proof above is verbatim that given in [BPR20, Lem.2.9], which applies directly to all d 2.) (3.1) implies that k y is indeed harmonic in B {y}.Note that in particular we have continuity in this region, which we did not assume a priori.
(5) Note that for any fixed y, the support of fn will not intersect y and so (ky, fn) makes perfect sense as the integral of ky against fn in B. The proof in [WP22] works exactly the same if we use this "zero boundary condition" for ky rather than a pointwise zero boundary condition.
We now check the boundary condition for k y .If f n are a sequence of functions as in (C) of Theorem 3, we have (k y , f n ) = k 2 (y, x)f n (x) dx by definition of k y .Then, by the harmonicity shown above (and because f n has support outside of some fixed neighbourhood B δ (y) of y for n large enough), this integral is equal to k 2 (w, x)f n (x)F (w) dxdw where we can take F to be a fixed smooth function supported in B δ (y) that is radially symmetric about y.Due to (2), we see that this is equal to , which is in turn bounded (using Cauchy-Schwarz) by the square root of E((h B , F ) 2 )E((h B , f n ) 2 ).Applying the boundary condition (C), we see that this indeed converges to 0 as n → ∞.Now, notice that by (2.4), for any δ > 0 there are some constants Thus by Bôcher's theorem [ABR01, Chap.III], we conclude that there is some harmonic function ν(w) In particular, there is some b > 0 such that k 2 (0, w) − bs(|w|) is harmonic and bounded in (1 − 2δ)B {0} and thus can be extended to a harmonic function on (1 − 2δ)B.Note that b must be positive, since k 2 (0, z)ρ |w| 0 (dz) = E(h |w| (0) 2 ) is positive and increasing to ∞ as w ↓ 0 by the domain Markov property.
From here, we return to a fixed y ∈ B. We choose ε > 0 such that d(y, ∂B) > 2ε and write where η x and η y are smooth functions with mass one, radially symmetric around x, y and supported in small non-intersecting balls around x and y respectively.Using the decomposition and harmonicity of ϕ y+2εB B we have where the second term on the right-hand side is bounded by Cauchy-Schwarz and (2.4).The first term is equal to bs(|x − y|) + O(1) in some neighbourhood of y by the definition of the distribution of h y+2εB

B
, and the previous paragraph.We have therefore shown that x → k 2 (x, y) satisfies the condition ( ) and is therefore equal to a multiple (which must actually be b) of G B (x, y).Since y ∈ B was arbitrary we are done.
Then we may define (h B , f ) to be the L 2 limit of (h B , f n ) as n → ∞.The limit does not depend on the approximation.
In what follows, we will therefore use the notation (h B , f ) for f ∈ H −1 (B) without further justification.

Gaussianity
It now remains to argue that the field is Gaussian.We do this via a decomposition of the field into a sequence of radial processes.We show using the domain Markov property that each of these processes is Gaussian, and that moreover, the whole sequence is jointly Gaussian.
First, we will see a bound for the 4-point function that is reminiscent of Wick's theorem: this will help us deduce continuity of our processes.In the second subsection, we do the basic case -the case of spherical averages.Finally, we extend this to a wider range of processes, obtained from so called spherical harmonics.The Gaussianity of the field is proved in this final subsection.4.1.Weak Gaussianity in terms of the 4-point function.-To prove continuity of our radial processes, we need the following bound that is implied by our assumption on existence of fourth moments.Notice that this can be seen as establishing a very weak form of Wick's theorem, i.e., getting us closer to Gaussianity.Recall that for r > 0, h r (0) = ϕ rB B (0) is the spherical average at radius r around the origin, and converges to 0 in probability as r ↑ 1.
Before proving the lemma, let us show that for any 0 < r < 1 − η, and for some constant C(η) depending only on η, where A a = {z 1 , z 2 , z 3 , z 4 ∈ B 4 : d(z i , z j ) a for some i = j}.Indeed, by symmetry it suffices to bound the integral over the region where min( However, on this region, by (2.7), we can bound |k 4 | above by an η-dependent constant times Then expanding the product of the two sums, we can show the desired bound for the integral of each term separately (uniformly in a), using Cauchy-Schwarz in the integral over z 3 and z 4 for all terms, except s(|z 3 − z 4 |)s(|z 1 − z 4 |) which can be integrated directly.Here we are using the fact that for any fixed z 1 ∈ ∂(rB), we have that ∂(rB) s(|z 1 −z|)ρ r 0 (dz) = O(s(r)).In particular, applying dominated convergence, this shows that the integral in (4.1) is well-defined and finite.
J.É.P. -M., 2022, tome 9 The expression for the 4-point function now also follows from dominated convergence.
Proof of (4.1).-Fix r ∈ (0, 1) and for each n 1 partition the sphere ∂((r − 2 −n )B) into regions each having diameter no larger than 2 −n .For each of these regions we can then choose a ball of maximal radius centred at a point in the region, so that the ball intersected with the sphere lies inside the region, but within distance 2 −2n from the boundary of the region.This produces a sequence {z i,n } 1 i mn and radii {r i,n } 1 i mn (all less than 2 −(n+1) ) such that the balls z i,n + r i,n B do not intersect.
We can now set ν i,n for each i, n to be a smooth mollifier supported on z i,n + r i,n B, that is radially symmetric about z i,n and has total mass one.Using that for any distinct 1 i, j, k, m n and every n.Dominated convergence using (4.3) then allows us to conclude.
The Wick-type of bound is slightly trickier: Proof of (4.2). -Fix z 1 , z 2 , z 3 , z 4 distinct and for each j write a j = min i =j d(z i , z j )/2.Denote B j := (z j + a j B) ∩ B for each j so that the B j do not intersect.Write ρ j for uniform measure on ∂(z j + a j B) and ρ j for harmonic measure on ∂B j ∂B, seen from z j .
Claim 9. -We can write This claim could be checked straightforwardly if we could apply the domain Markov property inside the regions B 1 to B 4 , so it should not be too surprising.However, since the Markov property has only been assumed for balls, we will have to do a little bit of work to show this.This will be carried out shortly, but let us first see how (4.2) follows.
First, notice that using the domain Markov property we can write ) for each j .By the same dominated convergence argument as for (4.2), and writing k 2B 4 for the four-point function of h 2B , we see that the latter is equal to ∂(zj +aj B) |k 2B 4 (x, y, z, w)| ρ 1 (dx) . . .ρ 4 (dw).But then by (4.3) and using that each z j + a j B is far from the boundary of 2B we see that this is less than some constant, not depending on δ, times (4.5) In the case that a j > δ, sup ∂Bj ∂B d ρ j /dρ j can be bounded above by the probability that d-dimensional Brownian motion on the half space {(x 1 , . . ., x d ) : x 1 > 0}, started from (δ/a j ) reaches the boundary of the unit sphere before hitting the hyperplane {x 1 = 0}.This is bounded above by a constant times (δ/a j ) [Bur86,Th. 4.4].When δ > a j we have sup ∂Bj ∂B d ρ j /dρ j = 1.So overall, we obtain the bound (where from now on a b means a Cb for some constant C depending on the dimension): Thus applying Cauchy-Schwarz to (4.4) we see that On the other hand, we can lower bound g(z 1 , z 2 , z 3 , z 4 ) using the explicit expression In particular, when |x| ).This implies that, on the region |x − y| > δ, we have Without loss of generality we can now assume that a 1 a 2 a 3 a 4 .Combining the above lower bounds for G B and the definition of g with the upper bounds for k 4 , we obtain (4.2) under the condition that a 2 > a 1 /10.Note that the δ −η correction is only actually needed when the dimension d = 2.
When a 2 a 1 /10, i.e., one point is considerably further than the rest, our bounds do not suffice -this is because we haven't taken properly care of cancellations occurring in the third moment, when three points are together.Let us do that now.Notice that when a 2 a 1 /10 we have B 2 , B 3 , B 4 ⊂ z 2 + a 1 B; we write ρ i for harmonic measure seen from z 2 , z 3 , z 4 on ∂(z 2 + a 1 B) ∂B.We need an extension of Claim 9, that first separates the three points, and looks at the occurring cancellations: Claim 10. -In the case that a 2 < a 1 /10 we can further write k 4 (z 1 , z 2 , z 3 , z 4 ) as Again, we postpone the proof of the claim and first see how it implies (4.2).From the same Cauchy-Schwarz argument and bounds as in a 2 > a 1 /10 case (noting also that d(z j , ∂(z 2 + a 1 B)) > a 1 /2 for j = 3, 4), we can deduce that the first term in (4.6) is δ −η g(z 1 , z 2 , z 3 , z 4 ) for some η ∈ [0, 1).Again the δ −η correction is only needed when d = 2.
To deal with the latter terms in (4.6), we use the fact that the covariance of h B is equal to bG B .Indeed, using the explicit expression for G B we see that the latter term is given by 4b 2 σ∈S(2,3,4) and since G D G D for D ⊂ D, (4.2) then follows in the case a 2 a 1 /10 too.
J.É.P. -M., 2022, tome 9 It now remains to prove the claims.
Proof of Claims.-We start with (4.4).As mentioned above, we cannot apply the Markov property directly inside the regions B j .Instead we will work in 2B and use the equivalent of (4.4) for the field h 2B , together with the Markov decomposition In this decomposition, the two summands are independent, h B 2B has the law of h B and ϕ B 2B is harmonic inside B. For each i = 1 . . .4, let η i be a smooth function supported in B i , rotationally symmetric about z i , such that k 4 (z 1 , z 2 , z 3 , z 4 ) = E( i (h B , η i )).Let us also denote by ρ i the harmonic measure on ∂B i (now including the part on ∂B) seen from z i .
By the Markov decomposition for the field h 2B inside {z 1 + a 1 B, . . ., z 4 + a 4 B}, as in (2.2) with n = 4, we have that Opening the brackets and using the independence of h B 2B and ϕ B 2B we can write this further as As h B 2B has the law of h B , to prove that i=1 (h B , ρ i ) it suffices to show that all terms with |S| = 4 cancel out.Since (h B , ρ i ) = (h B , ρ i ) a.s.for each i, this proves the claim.
To show the cancellation, first observe that as ϕ B 2B is harmonic inside each B j , we have that (ϕ B 2B , η j ) = (ϕ B 2B , ρ j ) for every j = 1, . . ., 4 and thus the terms with |S| = 4 cancel out.Also, both h B 2B and ϕ B 2B are mean zero, so all terms with |S| = 1 or |S| = 1 also cancel out.We are left to consider the cases when |S| = 2 and |S | = 2.

But we already know that
To prove this, we will use the fact that any continuous stochastic process (indexed by positive time) with independent increments, is Gaussian.The fact that h r (0) has a continuous modification in r comes from Lemma 8 of the previous subsection.
Proof of Lemma 11. -Since the variance of this process increases as r ↓ 0, it is natural to parameterise the time so that the process starts at r = 1.Thus we set X t = h 1−r (0).Then X t → 0 in probability as t ↓ 0 and X has independent increments by the domain Markov property (more precisely (2.1)).Thus, to show that it is a Gaussian process, by [Kal97,Th. 11.4], it suffices to prove that it admits a continuous modification in t.
We apply Kolmogorov's continuity criterion to show this.By the domain Markov property and scaling assumption it suffices to control the behaviour at time 0, i.e., it is enough to show that for some C ∈ (0, ∞), some η < 1 and all δ ∈ (0, 1): We will generalise the case of spherical averages to a wider class of processes, stemming from so called spherical harmonics.The interest comes from the following classical theorem (see e.g.[SW71, Chap.IV]).
J.É.P. -M., 2022, tome 9 Remark 13. -In fact one can write out a specific collection of such functions using Legendre polynomials and Bessel functions and choose e n,j,i to be eigenfunctions of ∆.This is not necessary here.
For example, in d = 2, one has M 0 = 1 and M n = 2 for n 1, with the ψ-s given by the usual Fourier series on the circle.That is, form an orthonormal basis of L 2 (B), where (J n ) n 0 are the Bessel functions and α n,k are the zeroes of J n for each n.
Using these notations, the main result of this section can be stated as follows: Proposition 14 (Gaussianity).-The random variables (h B , e n,j,i ) n∈N,j∈{1,...,Mn},i∈N are jointly Gaussian.In particular To prove Proposition 14 we will choose appropriate radial functions g n,j for which we can verify that h B tested against g n,j (r)ψ n,j (θ) on each sphere at radius r is a Gaussian process in r ∈ (0, 1).The key observation is the following.For r ∈ (0, 1) and a smooth function ψ : ∂B → R let ν ψ r be the signed measure defined by the condition that for all functions φ such that φ(r, θ) ∈ L 1 (∂B, ρ 1 0 (dθ)), where as before ρ 1 0 is uniform measure on ∂B.
Lemma 15. -Suppose that ϕ is a harmonic function in B. Suppose also that ψ(θ) is a smooth function defined on ∂B such that (r, θ) → ψ(θ)r n is harmonic in B. Then r −n ν ψ r (ϕ) is constant as a function of r on (0, 1).
Finally, if the r i 's are not distinct, the same argument still works: the "pairs" just mentioned will simply be larger tuples, concluding the proof of the lemma.
We are now ready to prove Proposition 14.
The fact that e n,j,i form a basis of L 2 (B), together with linearity of (h B , f ) and that h B is zero outside of B, now implies that (h B , f ) f ∈C ∞ c (R d ) is a Gaussian process.
J.É.P. -M., 2022, tome 9 has the law of c times a d-dimensional zero boundary Gaussian free field in B. (A) Domain Markov property with scaling for balls, as in Definition 2. (B) Moments.We have E

B
− h a+rB B , f ) → 0 a.s.along a subsequence as n → ∞.Thus, it must be the case that ϕ a+rB B (z) = ϕ a+rB B (z) a.s.Applying this for a dense collection of z and using harmonicity of ϕ a+rB B , ϕ a+rB B proves the uniqueness.We will use the following consequences of this uniqueness repeatedly: -if a + rB ⊂ a + r B ⊂ B then (2.1) ϕ a+rB B − ϕ a +r B B is independent of ϕ a +r B B and equal in law to (r/r ) (2−d)/2 ϕ (a−a )+(r/r )B B ; -we can apply the domain Markov property in several balls at once.More precisely, if B 1 , . . ., B n are n balls, and h with the X m i.i.d. each having the law of ϕ 0.5B B (0). Adding up the variances gives the desired bound.When z = 0 with |z| = r < 1, let r = 1/(1 + r) ∈ (1/2, 1), ε = ε/(1 + r) and z = −r z.By applying the Markov property for h z +r B B in ε B we can write ϕ ε B B (0) = ϕ z +r B B (0) + ϕ(0) where the summands are independent and by translation invariance and scaling, ϕ(0) has the law of (r ) 2−d times ϕ z+εB B (z).But since the variance of the summands add up, the variance of ϕ(0) is no greater than the variance of ϕ ε B B (0). Hence the claim follows from the case z = 0.The 4-point function.-Similarly, due to (2.2), for z 1 , . . ., z 4 ∈ B with min j =i as required.A similar argument appears in [BPR20, Proof of Lem.2.18].Spherical averages.-For z ∈ B and ε < d(z, ∂B) we define h ε (z) = ϕ z+εB B for the zero boundary Green's function in B. We are going to use the following characterisation of G B (see for example [WP22, Lem.3.7]).( ) Suppose that for y ∈ B, k y (x) is a harmonic function defined in in B {y}, such that k y (x) − bs(|x − y|) is bounded in a neighbourhood of y for some b > 0 and such that (k y , f n ) → 0 as n → ∞ for any sequence of functions f n as in our zero boundary condition (C) of Theorem 3. (5)(6) Then k y (x) = bG B (x, y) for all x = y; x, y ∈ B. Let us first see that for fixed y ∈ B the function k y (•) := k 2 (y, •) is harmonic away from {y}.For this, suppose that η > 0 and x ∈ B is such that |x − y| ∧ d(x, ∂B) > η.Then if we choose a such that |x − y| > η + 2a, we have k 2 (x, y) = E(ϕ x+(η+a)B B (x)ϕ y+aB B (y)) by definition, where ϕ x+(η+a)B B is almost surely harmonic in x + ηB.This means that ϕ x+(η+a)B B satisfies ϕ x+(η+a)B B (x) = ϕ x+(η+a)B B (w)ρ η x (dw) (the mean value property) with probability one, where ρ η x is uniform measure on ∂(x + ηB).Moreover, the domain Markov decomposition (and its uniqueness) mean that ϕ w+(a/2)B B (w) = ϕ x+(η+a)B B (w) + Z with Z centred and independent of ϕ y+aB B (y), so again using the definition of k 2 , we have that k 2 (w, y) = E(ϕ x+(η+a)B B (w)ϕ y+aB B (y))

Now recalling that g(z 1
, z 2 , z 3 , z 4 ) is the four-point function for the zero boundary GFF in B, we see that the integral on the right hand side is the fourth moment of the spherical average of the GFF at radius 1 − δ.But this spherical average is a centred Gaussian with variance − log(1 − δ) when d = 2 and 1 − (1 − δ) 2−d when d > 2 -see [WP22, Eq. (13)].Since these are both of order δ as δ → 0 the fourth moment is of order δ 2 by Gaussianity of the GFF itself.4.3.Gaussianity in the general case.-In what follows, we will often use the coordinates z = (r, θ) = (|z|, z/|z|) for a point in R d , d 2.
r n ψ] − r n ψ∆ϕ),where we are integrating against the standard volume measure on r 0 B on the righthand side, and the induced measure on its boundary on the left, which is a multiple of uniform measure.Now the right-hand side is zero as both r n ψ and ϕ are harmonic by assumption.Thus we deduce that n ∂B ϕ(r 0 , θ)ψ(θ) ρ 1 0 (dθ) = r 0 ∂B ψ(θ) d dr ϕ(r, θ) r=r0 ρ 1 0 (dθ).
Definition 2 (Domain Markov Property (DMP) with scaling for balls) Let d 2 and B ⊆ R d denote the open unit ball.Suppose that h is a random Schwartz distribution with support on B, i.e., a random continuous functional