Fractional Gaussian fields: a survey

We discuss a family of random fields indexed by a parameter $s\in \mathbb{R}$ which we call the fractional Gaussian fields, given by \[ \mathrm{FGF}_s(\mathbb{R}^d)=(-\Delta)^{-s/2} W, \] where $W$ is a white noise on $\mathbb{R}^d$ and $(-\Delta)^{-s/2}$ is the fractional Laplacian. These fields can also be parameterized by their Hurst parameter $H = s-d/2$. In one dimension, examples of $\mathrm{FGF}_s$ processes include Brownian motion ($s = 1$) and fractional Brownian motion ($1/2<s<3/2$). Examples in arbitrary dimension include white noise ($s = 0$), the Gaussian free field ($s = 1$), the bi-Laplacian Gaussian field ($s = 2$), the log-correlated Gaussian field ($s = d/2$), L\'evy's Brownian motion ($s = d/2 + 1/2$), and multidimensional fractional Brownian motion ($d/2<s<d/2 + 1$). These fields have applications to statistical physics, early-universe cosmology, finance, quantum field theory, image processing, and other disciplines. We present an overview of fractional Gaussian fields including covariance formulas, Gibbs properties, spherical coordinate decompositions, restrictions to linear subspaces, local set theorems, and other basic results. We also define a discrete fractional Gaussian field and explain how the $\mathrm{FGF}_s$ with $s \in (0,1)$ can be understood as a long range Gaussian free field in which the potential theory of Brownian motion is replaced by that of an isotropic $2s$-stable L\'evy process.


Introduction Introduction
The d-dimensional fractional Gaussian field h on R d with index s ∈ R (abbreviated as FGF s (R d )) is given by where W is a real white noise on R d and (−∆) −s/2 is the fractional Laplacian on R d . In Sections 2 and 3, we will review classical and recent literature on the fractional Laplacian (see, e.g., [LD72,Sil07,CSS08,CG11]) and show how to assign rigorous meaning to (1.1).
Our goal is to provide a mathematically rigorous, unified, and accessible account of the FGF s (R d ) processes, treating the full range of values s ∈ R and d ∈ N. This paper is fundamentally a survey, but we also present several basic facts that we have not found articulated elsewhere in the literature. Many of these are generalizations of classical results that had previously only been formulated for specific d and s values.
We hope that this survey will increase the circulation of basic information about fractional Gaussian fields in the mathematical community. For example, the vocabulary and content of the following statements should arguably be well known to probabilists, but the authors were unaware of much of it until recently: • In dimension 3, the Gaussian field with logarithmic correlations has been used as an approximate model for the gravitational potential of the early universe; its Laplacian is a FGF −1/2 (R 3 ) and has been used to model the perturbation from uniformity of the mass/energy density of the early universe 1 . • In dimension 4, the so-called bi-Laplacian field has logarithmic correlations, and its Laplacian is white noise. • In any dimension, Lévy Brownian motion can be defined as a random continuous function whose restriction to any line has the law of a Brownian motion (modulo additive constant). In dimension 5, the Laplacian of Lévy Brownian motion is the Gaussian free field. When H < 0 (grey shaded region) the FGF is defined as a random tempered distribution, not a function. When H ∈ (0, 1), the FGF is defined as a random continuous function modulo a global additive constant. Generally, for integers k > 0 and H ∈ (k, k + 1), the FGF is a translation invariant random k-times-differentiable function defined modulo polynomials of degree k. For integer H = k ≥ 0, the FGF is a random (k − 1)-times-differentiable function (or distribution if k = 0) defined modulo polynomials of degree k.
We also hope that this text will be a useful reference for experts in the study of Gaussian fields; to this end, we provide a robust account of the regularity of FGF fields, the long and short range correlation formulae, conditional expectations given field values outside of fixed domains, the Fourier transforms and spherical coordinate decompositions of the FGF, and various bounded-domain definitions of the FGF.
The family of fractional Gaussian fields includes several well-known Gaussian fields such as Brownian motion (d = 1 and s = 1), white noise (s = 0), the Gaussian free field (s = 1), and the log-correlated Gaussian field (s = d 2 ). Given s ∈ R and d ≥ 1, the Hurst parameter H is defined by If h is an FGF s (R d ), the Hurst parameter describes a scaling relation satis-fied by h. For a > 0, the field x → h(ax) has the same law as x → a H h(x). 2 When d = 1 and H ∈ (0, 1), the FGF s (R d ) process is commonly known as fractional Brownian motion with Hurst parameter H, and is the subject of an extensive literature (see the survey [CI13]). Brownian motion itself corresponds to H = 1/2 and s = 1.
The law of a Brownian motion or fractional Brownian motion B t , indexed by t ∈ R and defined so that B 0 = 0, is not translation invariant. However, the law of Brownian motion is translation invariant if we consider Brownian motion as a random process defined only modulo a global additive constant. In other words, Brownian motion has stationary increments. Similarly, the indefinite integral of a Brownian motion can be interpreted, in a translation invariant way, as a random function defined modulo the space of linear functions. We generally interpret all of the FGF s processes as translation invariant random distributions, but in some cases they are defined modulo a space of polynomials. More precisely, when H < 0, FGF s (R d ) is a translation invariant random tempered distribution (that is, a generalized function) on R d . When H > 0, FGF s (R d ) is a translation invariant random element of the space C H −1 (R d ) modulo the space of polynomials on R d of degree no greater than H . This means that h is defined as a linear functional on the subspace of test functions φ satisfying R d φ(x)L(x)dx = 0 for all polynomials L of degree H . Alternatively, at the cost of breaking translation invariance, we may define FGF s (R d ) as a random element of C H −1 (R d ) by fixing the derivatives of h at 0 up to order H − 1. The FGF covariance structure is described by the Hurst parameter H. When H is a positive non-integer, we have Cov[(h, φ 1 ), (h, φ 2 )] = C(s, d) Note that H is an affine function of s and can be used instead of s to parameterize the family of FGFs. We use the parameter s in part to highlight the connection to the fractional Laplacian and white noise. With our convention, white noise is FGF 0 (R d ) and the Gaussian free field is FGF 1 (R d ).
However, in many of our formulas and theorems H will be the more natural parameter to use; thus, we fix the relationship (1.2) and reference both H and s throughout the paper.

Examples
The simplest example of a fractional Gaussian field is FGF 0 (R d ), which is white noise. We denote by S(R d ) the space of Schwartz functions on R d , and we let S (R d ) be its dual, the space of tempered distributions (see Section 2 for details). If h ∈ S (R d ) and φ ∈ S(R d ), we use the notation (h, φ) for h evaluated at φ. White noise (surveyed in [Kuo96]) is a random element of S (R d ) with the property that for φ 1 , φ 2 ∈ S(R d ), the random variables (h, φ 1 ) and (h, φ 2 ) are centered Gaussians with covariance Taking d = 1 and s = 1, we see that (−∆) −s/2 is the antiderivative operator. It follows that FGF 1 (R) is the antiderivative of one dimensional white noise, which is a Brownian motion interpreted as a real-valued function modulo constant. If we fix the constant by setting the value at 0, we get ordinary Brownian motion.
If s = 1 and d ∈ N, then FGF 1 (R d ) is a d-dimensional generalization of Brownian motion called the Gaussian free field (GFF). As surveyed in [She07], the GFF is a random tempered distribution on R d (defined modulo additive constant if d = 2) with covariance given by where Φ is the fundamental solution of the Laplace equation in R d . The two-dimensional GFF (which is the same as FGF 1 (R 2 )) has been studied in a wide range of contexts in recent years. It can be obtained as a scaling limit of random discrete models, such as domino tilings [Ken01], as well as continuum models, such as those arising in random matrix theory [RV06]. It is central to conformal field theory and Liouville quantum gravity [She10,DS11] and has many connections to the Schramm Loewner evolution [Dub09,SS10,MS12a,MS12b,MS,MS13]. For all d ∈ N, the d-dimensional GFF exhibits a certain Markov property: For each fixed domain D ⊂ R d , if we are given the restriction a GFF h to R d \ D, then the conditional law of h restricted to D is given by a conditionally deterministic function (the harmonic extension of the field from ∂D to D) 3 plus an independent zero-boundary GFF defined on D. In Section 5 we will establish an analogous property that applies when h is an FGF s (R d ) with s ≥ 0. Namely, if we are given the restriction of h to R d \ D, then the conditional law of h restricted to D is given by a conditionally deterministic function (the so-called s-harmonic extension of the field from R d \ D to D) plus a random function (the so-called zero-boundarycondition FGF s on D). If s ∈ N, then the conditionally deterministic function depends on the restriction to ∂D of h and its derivatives up to a certain order. This follows from the fact that (−∆) s is a local operator when s ∈ N.
As previously mentioned, another generalization of Brownian motion is the fractional Brownian motion (FBM). Fractional Brownian motion appears to have been first introduced by Kolmogorov in 1940 [Kol40], and the term "fractional Brownian motion" was introduced by Mandelbrot and Van Ness in 1968 [MVN68]. As motivation, Mandelbrot and Van Ness discuss various empirical studies of real world processes (the price of wheat, water flowing through the Nile, etc.) that had been made by Hurst, who found different scaling exponents in different settings.
The definition of fractional Brownian motion can be extended to describe a random function modulo additive constant on R d when d > 1. Given H ∈ (0, 1) we define the FBM (also called the fractional Brownian field) on R d as a mean-zero Gaussian process (B H t ) t∈R with covariance where H is the Hurst parameter of the field. We will prove in Section 6 that the multidimensional fractional Brownian motion defined this way is equivalent to FGF s (R d ), where H = s − d 2 ∈ (0, 1). In the case H = 1/2, this multidimensional process was introduced by Lévy in 1940 and is known as Lévy Brownian motion [Lév40]. General processes including multidimensional fractional Brownian motion are discussed in Yaglom in 1957 and by Gangolli in 1967 [Yag57, Gan67]. (Gangolli gives general analytic arguments for positive definiteness of covariance kernels that apply in this case.) Fractional Brownian motion is studied in more detail in works of Mandelbrot, as referenced in [Man75]. More detailed and modern discussions of fractional Brownian motion (including topics such as excursion set theory, Hausdorff dimension, Hölder regularity, etc.) can be found in [AT07,Adl10].
The log-correlated Gaussian field (LGF) is a random element h of the space of tempered distributions modulo constants and has covariance given by In two dimensions, the LGF coincides with the GFF (up to a constant factor). We will see in Section 3 that the d-dimensional LGF is a multiple of FGF d/2 (R d ). In recent years the log-correlated Gaussian field has enjoyed renewed interest because of its relationship to Gaussian multiplicative chaos. For a survey article of Gaussian multiplicative chaos see [RV13]. Furthermore, the LGF in R 3 plays an important role in early universe cosmology, where it approximately describes the gravitational potential function of the universe at a fixed time shortly after the big bang; see [DRSV] for more discussion and references.
Another noteworthy subclass of the fractional Gaussian fields is FGF 2 (R d ), which is known as the bi-Laplacian Gaussian field. The discrete counterpart of the bi-Laplacian Gaussian field is called the membrane model in physics literature; for a mathematical point of view see [Sak03], [Kur07], [Kur09], and [Sak12]. In dimension at least five, there is a natural discrete field associated with the uniform spanning forest on Z d whose scaling limit is FGF 2 (R d ) [SW13].

Fractional Gaussian fields in one dimension
The FGF s (R d ) processes are easiest to classify and explain when d = 1. We first consider H = s − d 2 ∈ (0, 1) (so that s ∈ (1/2, 3/2)), in which case the FGF s (R) is a Gaussian random function h : R → R which we interpret as being defined modulo an additive constant. This means that while the quantity h(t) is not a well-defined random variable for t ∈ R, the quantity h(t 1 ) − h(t 2 ) is a well-defined random variable for t 1 , t 2 ∈ R. When H ∈ (0, 1), the FGF s (R) is the stationary-increment form of the fractional Brownian motion with Hurst parameter H. The law of the fractional Brownian motion is determined by the variance formula When H = 0, so that s = 1/2, the FGF s (R) is the log-correlated Gaussian field (LGF), which is defined as a random tempered distribution modulo additive constant.
When d = 1 the weak derivative of an FGF s (R) is an FGF s−1 (R). Thus all FGF s (R) processes may be obtained by either integrating or differentiating fractional Brownian motion (with s ∈ (1/2, 3/2)) or the LGF (s = 1/2) an integral number of times. From this, it is clear that if an FGF s (R), for s ∈ (1/2, 3/2], is defined modulo additive constant in a translation invariant way, then the distributional derivatives FGF s−1 (R), FGF s−2 (R), etc. are defined without an additive constant. Thus the FGF s (R) is defined as a random tempered distribution without an additive constant when s ≤ 1/2. Similarly, if the FGF s (R), for s ∈ (1/2, 3/2] is defined modulo additive constant (in a translation invariant way), then the indefinite integrals FGF s+1 (R), FGF s+2 (R), etc. are respectively defined modulo linear polynomials, quadratic polynomials, etc.
The following proposition, rephrased and proved as Theorem 7.1 in Section 7, is one reason that the one-dimensional case is significant.

Interpretation as a long range GFF
The Gaussian free field FGF 1 (R d ) can be approximated by the discrete Gaussian free field, which only has nearest neighbor interactions. This discrete Markov property gives rise to the domain Markov property of the Gaussian free field in the limit [She07]. In Section 12, we construct a discrete version of FGF s for s ∈ (0, 1) by introducing a discrete fractional gradient to play the role of the discrete gradient in the definition of the discrete GFF. The fractional gradient involves long range interactions, which may be viewed as the reason that the Markov property fails for FGF s when s is not an integer.
The comparison between the short range FGF s (R d ) (when s ∈ Z) and the long range FGF s (R d ) (when s / ∈ Z) may also be seen from the point of view of the corresponding potential theories. As an illustration, consider GFF and FGF s for 0 < s < 1. The covariance kernel for the Gaussian free field is given by the solution of the ordinary Laplace equation −∆ f = φ. As we will see, the counterpart for FGF s with 0 < s < 1 is the fractional Laplacian equation (−∆) s f = φ. The Laplacian is a local differential operator, while (−∆) s for s ∈ (0, 1) is a non-local pseudo-differential operator and (−∆) s f (x) depends on the values of f (x) for all x ∈ R d . Another way to see the distinction between the s = 1 and s ∈ (0, 1) cases is to recall that the Green's function for the Dirichlet Laplacian is given by the density of the occupation measure of a Brownian motion (see [MP10], for example), which is continuous. The corresponding process when s ∈ (0, 1) is an isotropic 2s-stable Lévy motion, which is a jump process.

Preliminaries
In this section we remind the reader of some definitions and facts regarding tempered distributions and homogeneous Sobolev spaces. Some of the following notation and ideas are from [Tri83], to which we refer the reader for more discussion on homogeneous spaces. We will introduce and construct several linear spaces. To aid the reader in keeping track of the various definitions, we include a glossary of these definitions in the appendix on page 67.

Tempered Distributions and Sobolev spaces
Fix a positive integer d, and denote by S(R d ) the real Schwartz space, defined to be the set of real-valued functions on R d whose derivatives of all orders exist and decay faster than any polynomial at infinity. A multiindex β = (β 1 , . . . , β d ) is an ordered d-tuple of nonnegative integers, and the order of β is defined to be |β| := ∑ d j=1 β j . We equip S(R d ) with the topology generated by the family of seminorms f n,β := sup The space S (R d ) of tempered distributions is defined to be the space of continuous linear functionals on S(R d ) .
We take the convention that the Fourier transform F acting on a Schwartz function φ on R d is the function which we will often abbreviate as φ(ξ). The complex Schwartz space (the space of functions whose real and imaginary parts are in S(R d )) is closed under the operation of taking the Fourier transform [Tao10, Section 1.13], so the inverse Fourier transform F −1 is well-defined on the complex Schwartz space and satisfies the formula We define the Fourier transform f of a tempered distribution f by setting ( f , φ) := ( f , φ), so that F and F −1 may be interpreted as operators For the fundamentals of the theory of distributions, we refer the reader to [Lax02,Appendix B] or [Tao10]. For a more detailed introduction to distribution theory we refer to [FJ98] and [Hör03].
which is equal to the image of S r (R d ) under the inverse Fourier transform operator. We define the Fourier transform of an element of S r (R d ) as an Define the space and equipH s (R d ) with the inner product .
We define the Sobolev spaceḢ s (R d ) to be the Hilbert space completion of H s (R d ), which we continuously embed in S H (R d ) as follows. If { f n } n≥1 is a Cauchy sequence inH s (R d ) and φ ∈ S H (R d ), then by Plancherel and The first factor on the right-hand side tends to 0 as min(m, n) → ∞ and the second factor is finite Cauchy in R, which implies that we can define a linear map f : we conclude that f is a continuous functional on S H (R d ). Therefore, we may realizeḢ s (R d ) as a subset of S H (R d ) by identifying each Cauchy se- We can characterizeḢ s (R d ) in another way which will be useful for the following section. Note that if (φ n ) n∈N is anḢ s (R d )-Cauchy sequence of Schwartz functions converging to f in S H (R d ), then ( φ n ) n∈N is Cauchy in L 2 (R d , |ξ| 2s dξ), where |ξ| 2s dξ denotes the measure whose density with respect to Lebesgue measure is ξ → |ξ| 2s . Therefore, there exists g ∈ L 2 (R d , |ξ| 2s dξ) to which φ n converges with respect to the L 2 (R d , |ξ| 2s dξ) norm. Furthermore, it is straightforward to verify using the Cauchy-Schwarz inequality that g = f ∈ S H (R d ). Therefore, where f ∈ L 2 (|ξ| 2s dξ) means that there exists g ∈ L 2 (|ξ| 2s dξ) such that

The Fractional Laplacian
The fractional Laplacian generalizes the notion of a power (−∆) s of the Laplacian from nonnegative integer values of s, for which it is defined as a local operator by iterating the Laplacian, to all real values of s. A standard reference for the fractional Laplacian is [LD72]. Here we use ideas from Section 2 of [Sil07]. Let k ∈ {−1, 0, 1, 2, . . .}, and let φ ∈ S k (R d ).
If s > − 1 2 (d + k + 1), then we set agrees with the local definition of −∆ when s = 1. Because of the singularity at the origin in its Fourier transform, (−∆) s φ is not necessarily Schwartz. However, it is real-valued, smooth, and has polynomial decay at infinity: (2.3) Furthermore, (−∆) s φ is real-valued and smooth.
Proof. The proof of smoothness is routine: we write the inverse Fourier transform using its definition and differentiate under the integral sign. To show that (−∆) s φ is real-valued, we note that a function ψ ∈ L 1 (R d ) is the Fourier transform of a real-valued function in L 1 (R d ) if and only if ψ(−ξ) and ψ(ξ) are complex conjugates for all ξ ∈ R d . The function ξ → |ξ| −2s φ(ξ) has this property whenever φ does, so we may conclude that (−∆) s φ is real-valued.
where C is some constant. To obtain the desired bound for the first integral, we write the integral in spherical form and apply integration-byparts with respect to the radial coordinate. For the second integral, we bound φ 2 (x) by a constant times sup |β|=k+1 |∂ β φ(0)||ξ| −k−1 near the origin, using Taylor's theorem.
is bounded for all multi-indices α. These spaces interpolate between C ∞ (R d ) and S(R d ), as the derivatives of their elements decay polynomially at a rate indexed by s. In particular, implies that φ vanishes except possibly at the origin. This implies that φ is a polynomial, which in turn implies that It is straightforward to verify that this definition agrees with (2.2) when f ∈ S(R d ). Observe that (−∆) s 1 (−∆) s 2 = (−∆) s 1 +s 2 for all s 1 , s 2 ∈ R. We will consider two important examples of elements of ((−∆) s S k (R d )) : (i) Elements of homogeneous Sobolev spaces. Let s ∈ R and H = s − d/2. It is straightforward to verify that f ∈Ḣ s (R d ) determines an el- Interpreting f as a linear functional on (−∆) s S k (R d ) by integration against a test function, the continuity of f with respect the U s+(k+1)/2 (R d ) topology follows from Proposition 2.1.
The following proposition gives an alternative representation of the fractional Laplacian in the case 0 < s < 1.
When s ∈ {0, 1, 2, . . .}, the fractional Laplacian (−∆) s coincides with the poly-Laplacian, a fundamental example of a higher order elliptic operator, obtained by iterating the Laplacian operator. For s ∈ (0, 1), (−∆) s is a classical example of a non-local pseudo-differential operator. These two classes generate all the operators of the form (−∆) s for s ≥ 0 in the sense that (−∆) s can be written as a composition of (−∆) s− s and (−∆) s .
For the properties of the poly-Laplacian, we refer the reader to [GGS10] and references therein. For more properties of (−∆) s where s ∈ (0, 1), see [Sil07] and reference therein.

White Noise
On a finite dimensional Hilbert space H with inner product (·, ·) H , one White noise on R d can be regarded as a standard Gaussian on L 2 (R d ). We will define W to be a random generalized function such that However, it is not obvious that there exists a measure on S (R d ) satisfying these conditions. Since we will rigorously construct the FGF in Section 3.1 in the same manner, we will review a construction of white noise following [Sim79].
Proof. We briefly sketch the proof given in [Sim79, Theorem 2.3] for the case d = 1; the case d > 1 may be proved similarly. We introduce coordinates to the space S(R) by writing each function φ ∈ S(R) and the topology of S(R) is equivalent to the one induced by the family of seminorms · m . Furthermore, S (R d ) is isomorphic to x n y n . Bochner's theorem states that characteristic functions of R n -valued random variables are in one-to-one correspondence with normalized, continuous, positive definite functions on R n . Using Bochner's theorem, we conclude for all n ∈ N 0 , there is a measure µ n on span(φ 1 , . . . , φ n ) such that By the uniqueness part of Bochner's theorem, these measure are consistent. By the Kolmogorov extension theorem, there exists a measure µ on R N 0 such that (2.5) holds for all φ in the linear span of {φ n } ∞ n=0 (that is, when at most finitely many of φ's coordinates are nonzero). It may be shown using the continuity of Φ that µ(s ) = 1 (see [Sim79] for details), which allows us to restrict µ to obtain a probability measure on s and conclude that (2.5) holds for all φ ∈ S(R).
We will use Theorem 2.3 in conjunction with the following proposition, which gives sufficient conditions for a functional to be positive definite. Proposition 2.4. Let (S, (·, ·)) be an inner product space. Then the functional Proof. Let v 1 , . . . v n be elements of S, and choose an orthonormal basis e 1 , . . . , e m of the span of {v 1 , . . . v n }. Let Z = (Z 1 , . . . , Z m ) be a vector of independent standard normal real random variables, and note that for all We will apply Theorem 2.3 to construct a measure µ on S (R d ) which we will refer to as white noise W. Recall that S(R d ) is a nuclear space and let us define the functional By Proposition 2.4, this functional is positive definite. Since it is also continuous and satisfies Φ 0 (0) = 1, Theorem 2.3 implies that there is a unique probability measure on S (R) having Φ 0 as its characteristic function, which we define as white noise W. In particular we have the relation which implies for every f ∈ S(R d ) the random variable (W, f ) is a centered Gaussian with variance f 2 L 2 (R d ) . An alternative to the preceding view of white noise as a random tempered distribution is to regard white noise as a collection of random variables {(W, f ) : f ∈ S(R d )}. The advantage of this perspective is that we may extend this collection so that (W, f ) is a well-defined random variable for all f ∈ L 2 (R d ). However, in this construction f → (W, f ) is no longer almost surely continuous. Recall the following definition from [Jan97] or [She07].
Definition 2.5. A Gaussian Hilbert space is a collection of Gaussian random variables on a common probability space (Ω, F , µ) which is equipped with the L 2 (Ω, F , µ) inner product and is closed with respect to the norm of the L 2 (Ω, F , µ) inner product.
To define a Gaussian Hilbert space , this map is an isometry. Since L 2 (Ω) is complete, we may extend this isometry to an operator from so if f and g are orthogonal with respect to the L 2 (R d ) inner product, then (W, f ) and (W, g) are independent. We may rewrite the above expression as Cov and say that W has covariance kernel δ(x − y) (here δ(x) dx is notation for the Dirac measure which assigns unit mass to the origin). In Section 3.2, we will compute the covariance kernel of the FGF s (R d ) for general s and d.

The FGF on R d
We provide the construction of FGF s (R d ) following the same procedure as used in Section 2.3 for white noise. We also compute the covariance kernel for the FGF s (R d ).

Definition of FGF s (R d )
We begin with some heuristic motivation for the rigorous construction that follows. We want to define h to be a standard Gaussian onḢ s (R d ). As a first guess, we might try to define a random element h ofḢ However, sinceḢ s (R d ) is infinite dimensional, no such random element exists [Jan97,She07]. However, we note that when h, Therefore, substituting (3.2) into (3.1) and defining φ := (−∆) s f , we find that it is reasonable to change the desired relation from (3.1) to The advantage of this formulation is that we may reinterpret it by replac- The norm on the right-hand side can be rewritten as then we say that h is a fractional Gaussian field with parameter s on R d and write h ∼ FGF s (R d ); note that by abuse of notation we refer to either h or its law as FGF s (R d ). We note that when h ∼ FGF s (R d ) and a > 0, the scaling relation We now provide a construction establishing the existence of fractional Gaussian fields. We would like to apply the Bochner-Minlos theorem with , not for all φ ∈ S(R d ). Therefore, we define a functional (3.5) which is finite for all Schwartz functions and which reduces to Let {φ α : α is a multi-index} be a collection Schwartz functions such that (3.5) By Proposition 2.4, C s is positive definite. Since C s is also continuous and satisfies C s (0) = 1, we may apply the Bochner-Minlos theorem to conclude that there is a random we obtain a random element of S H (R d ) which satisfies (3.4) (note that this restriction is necessary so that the definition does not depend on the arbitrary choice of functions φ α ).
As we did for white noise (see page 20), we may define a Gaussian Hilbert to the random variable (h, φ); we extend this isometry to an operator from T s (R d ) to L 2 (Ω). Writing φ ∈ T s (R d ) as a limit of functions in S H (R d ) and considering the limit of the corresponding characteristic functions, we conclude that We now make sense of the expression In this way, we have constructed a coupling between In this sense we can say that h = (−∆) −s/2 W. For a coupling in which this equation holds almost surely, see Proposition 6.3. Remark 3.1. Computing ||φ|| 2Ḣ −s (R d ) amounts to computing the covariance kernel of the FGF s (R d ), which will be done in Section 3.2.

The FGF covariance kernel
(3.6) We call G s (x, y) a covariance kernel of the FGF s (R d ). We point out that there can be more than one function G s satisfying (3.6). For example, if H ≥ 0 and G s (x, y) satisfies (3.6), then so does G s (x, y) + g(x, y) for any polynomial g in x or in y of degree no greater than H .
In this section we compute covariance kernels for the fractional Gaussian fields on R d . For most positive values of s, we find that G s (x, y) = C(s, d)|x − y| 2H for some constant C(s, d). When s < 0 the formula is similar but involves some derivatives of the delta function, and when H is a nonnegative integer there is a logarithmic correction. The constant C(s, d), and therefore also the correlation of FGF s (R d ), is positive when s ∈ (0, d/2), is (−1) s when s is a negative non-integer, and is (−1) 1+ H when H is a positive non-integer. The statement and proof of the following theorem are adapted from [LD72, Chapter 1, §1]. Theorem 3.3. Each of the following holds.
and H is not a nonnegative integer, then and .
Remark 3.4. In case (ii) above, we can also write Similarly, in case (iii), Then we may compute the covariance: where in the third line we used the Plancherel theorem, and in the last line we used the following Fourier transform formula given in [LD72, Chapter 1, §1]: It is important to note that (3.7) is only valid for 0 < s < d/2 when the class of test functions is taken to be S(R d ). Indeed, |ξ| −2s is not a tempered distribution when s ≥ d/2 (due to the singularity at the origin), and C(s, d)|x| 2H is not a tempered distribution when s ≤ 0. Therefore, we extend the Fourier transform formula (3.7) outside of the region where ψ(x, s) is an analytic continuation from 0 < s < d/2 to s ∈ (−k − 1, −k]. The result for (iii) follows from the equality Finally, to obtain (iv) we will take a limit as t → s of both sides of (3.8) see [LD72,p. 50] for more details. Since φ 1 and φ 2 are in S k (R d ), we have We use Taylor's theorem to write and substitute into (3.8). Taking t → s and using lim t→s (t − s)C(t, d) = Res t=s C(t, d), we obtain (iv). The formula for c We may therefore define the setḢ s 0 (D) to be the closure of C ∞ c (D) iṅ H s (R d ) and equip it with theḢ s (R d ) inner product. 27 We will construct a fractional Gaussian field FGF s (D) for all allowable domains D ⊂ R d (see Remark 4.3). The following lemma gives sufficient conditions for a domain to be allowable. Proof. Let D ⊂ R d be a domain, and let φ ∈ S(R d ) and g ∈ C ∞ c (D). We have by the Plancherel formula and Cauchy-Schwarz. If 0 ≤ s < d/2, we conclude that φ Ḣ−s (R d ) is finite and therefore that D is allowable.
x α dx for every multi-index α satisfying |α| ≤ H (such a function may be constructed via a Gram-Schmidt procedure). Since ηg = 0, we have By our choice of η, the Fourier transform of φ − η vanishes to order H at the origin, so φ − η Ḣ−s (R d ) is finite. Therefore D is allowable.
Suppose that s − d/2 > 0 is not an integer and that D R d . Without loss of generality, we may assume D does not contain the origin. Let P φ be the unique polynomial of degree H such that all the derivatives up to order H of F (φ − P φ (D)δ) are zero, where P(D) denotes the differential operator corresponding to a polynomial P and δ denotes a unit Dirac mass at 0. Then by a constant times |ξ| H +1 near the origin and by a constant times |ξ| H as ξ → ∞.
Let φ ∈ S(R d ), and let D be an allowable domain. By the definition of allowability, (φ, ·) L 2 (R d ) is a continuous linear functional onḢ s 0 (D). Therefore, by the Riesz representation theorem for Hilbert spaces, there exists a unique f ∈Ḣ s . Writing out the definition of ( f , g)Ḣ s (R d ) and using the Plancherel formula, we see that this implies that f is the unique solution of the distributional equation

The zero-boundary FGF in a domain
Let D R d be an allowable domain, let s ≥ 0, and define the functional for φ ∈ S(R d ). Since C D,s is continuous by the definition of allowability, we may use Proposition 2.4 and the Bochner-Minlos Theorem to C D, s to conclude that there is a unique random element the support of h D is almost surely contained in D. We call 5 h D the zero-boundary FGF on D, abbreviated as FGF s (D). Remark 4.3. We construct h D ∼ FGF s (D) only when D is allowable because we want to ensure that h D is a tempered distribution (rather than a tempered distribution modulo a space of polynomials).
We can also define a Gaussian Hilbert space version of FGF s (D), following the corresponding discussion FGF s (R d ) in Section 3.1. In this way we obtain a collection of random variables If s is an even positive integer, then f Ḣ s . If s is an odd positive integer, then for all f ∈ C ∞ 0 (D). Therefore, if s = 0 then h D is white noise on D, and if s = 1 then h D is the GFF on D. Thus FGF s (D) generalizes the domain versions of white noise and the Gaussian free field.

Covariance kernel for the FGF on the unit ball
Let s ≥ 0, and let D be an allowable domain. As usual, we say that a function G s D : for h D ∼ FGF s (D) and for all φ 1 , φ 2 ∈ C ∞ c (D). We treat each of the cases (i) s is an integer, (ii) s ∈ (0, 1), and (iii) s is a non-integer greater than 1.
Suppose that s is a positive integer, and let φ ∈ S(B). By (2.65) in Chapter 2 of [GGS10], the unique solution of (4.1) is and dy, (4.4) which shows that G s B is the FGF s (D) covariance kernel. Suppose that 0 < s < 1. Let X t denote a 2s-stable symmetric Lévy process, and let τ B be the first time X exits B. Recall the definition of the constant C d,s in Theorem 3.3, and define u(x, y) = (2/π) 2s C d,s |x − y| 2H . By the potential theory of 2s-symmetric stable processes, (see, for example, [CS98]), the function

It follows that for all
] is the FGF s (D) covariance kernel. The following explicit formula for G B s is given as Corollary 4 in [BGR61]: wherek s,d = Γ(d/2) 4 s π d/2 Γ(s) 2 .
Suppose that s > 1 is not an integer and φ ∈ S(R d ). We claim that G s B (x, y) = B G s (x, u)G s− s (u, y)du is the covariance kernel for FGF s (D). Indeed, we may write (−∆) s = (−∆) s (−∆) s− s and calculate which implies that G s B is the FGF s (D) covariance kernel by (4.4). Remark 4.4. Similar results may be obtained for a more general class of domains D. The ingredients are the corresponding potential theory of the poly-Laplacian and fractional Laplacian for s ∈ (0, 1). Proof. If f ∈ Har s (D) and g ∈Ḣ s 0 (D), then ( f , g)Ḣ s (R d ) = ((−∆) s f , g) = 0. Therefore, Har s (D) andḢ s 0 (D) are orthogonal subspaces ofḢ s (R d ). Let f ∈Ḣ s (R d ). Since D is allowable, ((−∆) s f , ·) L 2 (R d ) is a continuous functional onḢ s 0 (D). Therefore, there exists f D ∈Ḣ s 0 (D) such that for all g ∈Ḣ s 0 (D), we have ( f , g)Ḣ s (R d ) = ( f D , g)Ḣ s (R d ) . In particular, this implies that ((−∆) s ( f − f D ), g) = 0 for all g ∈ C ∞ c (D), which means that

Projections of the FGF
Thus we can write f as a sum of elements of Har(D) andḢ s Observe that Proposition 5.1 implies that Har s (D) is a closed subspace oḟ H s (R d ). We define the projection operators P D f = f D and P Har D f = f − f D . We will make sense of P D h and P Har D h almost surely, although these are defined a priori only for h ∈Ḣ s (R d ) and not for arbitrary elements of S H (R d ).
We begin by observing that the solution f of (4.1) is given by . Therefore, we may apply the Bochner-Minlos theorem to the functional for P = P D and for P = P Har D to obtain random tempered distributions h D and h Har D , respectively. We call h Har D the s-harmonic extension of h restricted to R d \ D. In Section 8, we will show that h Har D is smooth in D almost surely.

By the definition of h Har
When s is a positive integer, the operator (−∆) s is local, in that case we have a stronger result: h Har D D is measurable with respect to the σ-algebra generated by the intersection of the value of h on every neighborhood of the boundary (that is, the action of h on test functions supported on a neighborhood of the boundary). This is a generalization of the corresponding Markov property for the Gaussian free field [She07].

Fractional Brownian motion and the FGF
The d-dimensional fractional Brownian motion B with Hurst parameter H > 0 is defined to be the centered Gaussian process on R d with (6.1) The existence of such a process is guaranteed by the general theory of Gaussian processes (for example, see Theorem 12.1.3 in [Dud02]). The special case H = 1 2 is called Lévy Brownian motion [Lév40], [Lév45].
Proof. Let x ∈ R d . Since the Fourier transform of δ x is ξ → e 2πix·ξ , one may verify from the definition of theḢ −s (R d ) norm that δ x − δ 0 is an element ofḢ −s (R d ) and therefore an element of T s where G s (x, y) = C(s, d)|x − y| 2H . Combining (6.1) and (6.2), we see that C(s, d)B and h have the same covariance structure. Since both are centered Gaussian processes, this implies that they have the same law.
Since (6.1) and (6.2) show that h(0) = B(0) = 0 almost surely, Proposition 6.1 establishes that the FGF s (R d ) can be identified as (a constant multiple of) the fractional Brownian motion by fixing its value to be zero at the origin.
Denote by C k,α (R d ) the space of functions on R d all of whose derivatives of order up to k exist and are α-Hölder continuous. Note that the differentiability and Hölder continuity of a function-modulo-polynomials is well-defined, because adding a polynomial to a function does not affect its regularity properties. (ii) Suppose that 1 < H < 2, and let s = d/2 + H. As in the case H ∈ (0, 1), it is straightforward to verify that ∂ α δ x − ∂ α δ 0 ∈ T s (R d ) when |α| ≤ 1 and x ∈ R d . Therefore, if h ∼ FGF s (R d ), we may fix all derivatives of h of order up to 1 to vanish at the origin. In this way we obtain a scale-invariant function h 0 whose restriction to S 1 (R d ) coincides with h. Since |h 0 (x)| has the same law as |x| H h 0 (1) by scale invariance, we have E|h 0 ( which implies that h 0 satisfies condition (2.4) with s = 1/2 almost surely (see Section 2.2). Therefore, h := (−∆) 1/2 h 0 is well-defined as a random element of S 0 (R d ). Furthermore, since Thus h is α-Hölder continuous for all α < s − 1 by the preceding case. By the proof 6 of [Sil07, Proposition 2.8], h is almost surely in C 1,α (R d ).
(iv) If H = 2, then we may apply the same reasoning we applied in case (ii), leveraging the H = 1 case.
As an application of the ideas presented in this section, we construct a coupling of all the fractional Gaussian fields on R d . Proof. We will start with an FGF with Hurst parameter 2 and apply the fractional Laplacian to obtain FGFs with Hurst parameters in (0, 2). The remaining FGFs are then obtained by applying integer powers of the Laplacian to FGFs with Hurst parameter in (0, 2]. Let h 2+d/2 ∼ FGF 2+d/2 (R d ). As discussed in the proof of Proposition 6.2 case (ii), we can fix the values and first-order derivatives of h to vanish at the origin to obtain a scale-invariant random function h 0 whose restriction whenever s ∈ (0, 1/2] and k = 1 or when s ∈ (1/2, 1) and k = 0. Therefore, we may define h s +d/2 = (−∆) 1−s /2 h 2+d/2 for all s ∈ (0, 2). If

Restricting FGFs
In this section we study how fractional Gaussian fields behave when restricted to a lower dimensional subspace.
We regard R d−1 as a subspace of R d by associating (x 1 , . . . , , which means that (h, φ ↑ ) is a well-defined random variable almost surely (see Section 3.1). Moreover, h d almost surely determines a ran- holds almost surely, where C is a constant depending only on d and s.
We refer to h d−1 as the restriction of h d to R d−1 .
Proof. Let {η k } k∈N be an approximation to the identity, which means that , because applying the definition of a convolution and making a substitution w = x − y yields satisfies the formula in Theorem 3.3 where we replace R d by R d−1 and set φ 1 = φ 2 = φ. This is the covariance structure of FGF on R d−1 with the same Hurst parameter as h d , up multiplicative constant. In other words, there is a constant C so that if we define , then h d−1 has the law of an FGF s−1/2 (R d−1 ) restricted to Φ. Therefore, h d−1 extends uniquely to a tempered distribution on R d−1 , and it satisfies (7.1) Since h d−1 is a function of h d almost surely, we can define a measurable function R on S such that h d−1 = Rh d . We call R the restriction operator.
We can see that R maps an FGF to a lower dimensional FGF with the same Hurst parameter. By applying R repeatedly, we can restrict an FGF(R d ) to an FGF(R d ) with the same Hurst parameter, as long as d > −2H.
When the Hurst parameter is positive, FGF s (R d ) is a pointwise-defined random function, so R agrees with the usual restriction of functions. In particular, we note that the restriction of a multidimensional fractional Brownian motion with Hurst parameter H to a line through the origin is a linear fractional Brownian motion with Hurst parameter H.

Regularity of FGF(D)
Let s ≥ 0, and let h ∼ FGF s (R d ) be coupled with h D ∼ FGF s (D) and h Har D as in Section 5, so that h = h D + h Har D . In this section, we will show that h Har D is smooth in D almost surely. First, we record some results on the fractional Laplacian following Section 2 of [Sil07].
Lemma 8.1. If s is a positive integer and g is a distribution on D such that (−∆) s g is smooth on D, then g is smooth on D.
Proof. Let s = 1; the case s > 1 follows by induction. If ∆g = 0, the desired result is Weyl's lemma (see Appendix B of [Lax02]). If ∆g is not zero, suppose that U ⊂ D is an arbitrary ball, and let g 1 be a function which is smooth on U such that (−∆)g 1 | U = (−∆)g| U [Fol99, Corollary 2.20]. Applying the result for the case ∆g = 0, we see that g − g 1 is smooth on U, and hence so is g. Since U was arbitrary, g is smooth on D.
Since P is smooth, we see that g is smooth in B whenever (−∆) s g = 0 in B. We can now prove the main result in this section. (ii) Suppose that d is even, 1 < H < 2 and D is a ball whose closure does not contain the origin. By the scale invariance of h, there exists c > 0 so that have E|∇h(x)| = c|x| H−1 for all x ∈ R d . Thus which implies that |∇h| satisfies condition (2.4) with s = H − 1 almost surely.
Therefore, |∇h Har D | satisfies (2.4) almost surely. By the same argument as in Case (i) above, (iii) If d is even, H ∈ {0, 1} and D is a ball, then s is an integer. So h Har D is smooth by Lemma 8.1.
(iv) If D R d is allowable, suppose that U ⊂ D is an arbitrary ball. Since (E|h(x)|) 2 , the arguments for the preceding cases imply that h Har U,D is smooth on U. By the formula we see that h Har D is also smooth on U. Since U is arbitrary ball contained in D, this implies that h Har D is smooth in D. If U is an arbitrary allowable domain in D, again by h Har U = h Har U,D + h Har D , we see that h Har U,D is smooth on U.
(v) Suppose that d is even and s > 0. In the preceding cases we have (vi) Suppose that d is odd, H > 0, and H is not an even integer. Suppose D is an allowable domain in R d and regard R d as a subspace of R d+1 by mapping x ∈ R d−1 to (x 1 , x 2 , . . . , x d−1 , 0). Since h = h D + h Har D where h D and h Har D are independent, h Har D is the conditional expectation of h given h on R d \ D. So if we regard R d \ D as a closed set in R d+1 , the restriction of h Har R d+1 \(R d \D) has the same law as h Har D on R d . Since the restriction of a smooth function is smooth, we conclude that h Har D is a smooth function in D almost surely.
(vii) If d is odd and s > 0, we apply the argument in Case (v) to the result from Case (vi).
Since h = h D + h Har D and h Har D is smooth in D, the regularity of FGF s (D) is the same as the regularity of FGF s (R d ). In other words, h D is has α-Hölder derivatives of order up to k, where k = H − 1 and α < H − H (Proposition 6.2).

The eigenfunction FGF
Let D ⊂ R d be a bounded domain, and let s ∈ (0, 1). In this section we discuss an different notion of a fractional Gaussian field on D, which we call the eigenfunction FGF and denote EFGF s (D).
The eigenfunction FGF is based on the following definition of a fractional Laplacian operator on D. Following Section 2.3 in [She07], we let { f n } n∈N be an orthonormal basis of eigenfunctions of the Dirichlet Laplacian on D, arranged in increasing order of their corresponding eigenvalues λ n > 0. We define for [She07]. We call (−∆) s D the eigenfunction fractional Laplacian operator on D. This fractional Laplacian operator determines a Hilbert space, analogous toḢ s 0 (D), with inner product given by Note that {λ −s/2 n f n } n∈N defines an orthonormal basis with respect to this inner product. We define EFGF s (D) to be a standard Gaussian on this space; more precisely, let {Z n } n∈N be an i.i.d. sequence of standard normal random variables and set for all φ ∈ C ∞ c (D), By Weyl's law, λ n = Θ(n 2/d ) as n → ∞, so the sum on the right-hand side converges almost surely for each φ. Furthermore, the functional h defined this way is a continuous functional by the same argument given for the GFF case in [She07]. We define EFGF s (D) to be the law of h.
Both the fractional Laplacian and the eigenfunction fractional Laplacian can be understood in terms of a local operator in d + 1 dimensions. In [CS07], the fractional Laplacian is realized as a boundary derivative for an extension problem in R d × [0, ∞). A corresponding analysis for the eigenfunction Laplacian is developed in [CT10] and [CDDS11] by considering a similar extension problem in D × [0, ∞). We carry out an analogous comparison between FGF s (R d ), FGF s (D), and EFGF s (D) by realizing each as a restriction of a higher-dimensional random field that can be understood as a Gaussian free field with spatially varying resistance (see Propositions 9.1, 9.2 and 9.3 below).
Let s ∈ (0, 1), and define α = 1−2s 1−s ∈ (−∞, 1). For simplicity, we will assume d ≥ 2. We introduce the coordinates (x 1 , . . . , x d , z) for R d+1 , and we define the following variant of the gradient operator. For φ ∈ S(R d+1 ), we set We will use as a space of test functions. Integrating by parts (see [Bas98, Chapter 7] for more details), we find that for all φ ∈ S sym (R d ), we have where the operator L α is defined by By the Bochner-Minlos theorem, we can define a random tempered distribution h α for which and (−L α ) −1 φ satisfies −L α (−L α ) −1 φ = φ and vanishes at infinity-see the proof of Proposition 9.1 for the existence of such a function. We then restrict the domain of h α to S sym (R d+1 ), so that (9.1) holds with φ in place of φ.
Since the right-hand side of (9.1) reduces when α = 0 to the GFF characteristic function evaluated at φ, we may think of h as a symmetrized and re-weighted 7 version of the Gaussian free field on R d+1 . We define restric- where δ 0 denotes the unit Dirac mass at z = 0. See the proof of Theorem 7.1 for an explanation of why the random variable on the right-hand side is well-defined. More precisely, we will show that the covariance kernel of h α | R d ×{0} is that of an FGF s (R d ). It follows from continuity of FGF s (R d ) (as a functional on S(R d )) that this restriction can be defined on a countable dense subset of S(R d ) and continuously extended to obtain a random tempered distribution.
Proposition 9.1. The restriction of h α to R d × {0} is an FGF s (R d ), up to multiplicative constant.
Proof. Let (X t ) t≥0 be a diffusion in R d+1 with a standard Brownian motion B in the first d coordinates and the process Z t := √ 2 δ Y t δ in the last coordinate, where δ = 2(1 − s) and Y is δ-dimensional Bessel process reflected symmetrically at 0. An application of Itō's formula reveals that L α is the generator of X. We define the Green's function where Q(x, ) := {y ∈ R d+1 : |x k − y k | < for all 1 ≤ k ≤ d + 1}. Since d + 1 ≥ 3 implies that X is transient, the limit on the right-hand side is well-defined-the proof is similar to the proof for the case α = 0 [MP10, Section 3.3]. Since L α is the generator of X, we have for all φ ∈ S sym (R d+1 ) (see Chapter II in [Bas98]). Therefore, We denote by ( s ) s≥0 the local time of Z at 0 and by (τ t ) t≥0 the inverse function of [RY99]. Then (X τ t ) t≥0 is a 2s-stable Lévy process in R d , and the integral s 0 1 2 1 {Z u ∈(− , )} du converges almost surely to s as → 0 [MO69]. Therefore, the Green's function of the a 2s-stable Lévy process in In other words, the restriction to {z = 0} of the Green's function of X is equal to the Green's function of a 2s-stable Lévy process in R d . The latter is proportional to |x 1 − x 2 | 2s−d [CS98, (1.1)], and the covariance kernel of FGF s (R d ) is also proportional to |x 1 − x 2 | 2s−d by Theorem 3.3. Since the law of a centered Gaussian process is determined by its covariance kernel, this concludes the proof.
In Propositions 9.2 and 9.3 below, we discuss projections of h α onto certain subdomains of R d+1 . These projections are analogous to those discussed for the FGF in Section 5. We state these propositions using terminology described in Remark 5.3, to which we refer the reader for a rigorous interpretation.
(We will use C to denote a generic constant whose value may change throughout the proof.) It is straightforward to verify that h α minus the conditional expectation of h α given its values on ∂D × R is equal in law to the field h cyl α whose covariance kernel is given by the Green's function of the diffusion X defined in the proof of Proposition 9.1 stopped upon hitting the cylinder ∂D × R. Equivalently, the covariances of h cyl α are given in terms of the inverse L −1 α of the operator L α with zero boundary conditions on and w λ (z) = w λ (−z) for all z ∈ R. A symbolic ODE solver may be used to express w λ in terms of the modified Bessel function of the second kind K s as w λ (z) = Cλ s/2 z s/(2−2s) K s 2(1 − s)z We define the operator L α,λ = −λ + ∂ ∂z z α ∂ ∂z . Integration by parts reveals that for for some constant C, where in the last step we have used the expansion as t → 0 + . We may restate (9.2) by writing δ 0 = Cλ −s (−L α,λ )w λ , where δ 0 denotes the unit Dirac mass at the origin. Therefore, using the relation where {η k } n∈N is an approximation to the identity, we have

Local Sets of the FGF on a Bounded Domain
The concept of a local set of the Gaussian free field is developed in [SS10]. It turns out to be an important concept and tool in the study of couplings between the GFF and random closed sets such as SLE ( [SS10], [MS12a], [MS12b], [MS], [MS13]). The theory of local sets of the Gaussian free field carries over to the FGF setting with minimal modification.
Let Γ D be the space of all closed non-empty subsets of D. We endow Γ with the Hausdorff metric induced by Euclidean distance: the distance between sets S 1 , S 2 ∈ Γ is d Haus (S 1 , S 2 ) := max sup where dist(x, S) := inf y∈S |x − y|. Note that Γ is naturally equipped with the Borel σ-algebra induced by this metric. Furthermore, Γ D is a compact metric space [Mun99,. Note that the elements of Γ are themselves compact.
Given A ⊂ Γ, let A δ denote the closed set containing all points in Γ whose distance from A is at most δ. Let A δ be the smallest σ-algebra in which A and the restriction of h (as a distribution) to the interior of A δ are measurable. Let A = δ∈Q,δ>0 A δ . Intuitively, this is the smallest σ-algebra Lemma 10.2 may be proved by making minor modifications to the proof of Lemma 3.9 in [SS10] to generalize from the setting s = 1, d = 2 to arbitrary s ∈ R and d ≥ 1.
We say a random closed set A coupled with an instance h of the FGF, is local if one of the equivalent conditions in Lemma 10.2 holds. For any coupling of A and h, we use the notation C A to describe the conditional expectation of the distribution h given A. When A is local, C A is the distribution h 1 described in (iii) above.
Given two distinct random sets A 1 and A 2 (each coupled with a FGF h), we can construct a coupling (h, A 1 , A 2 ) such that the marginal law of (h, A i ) (for i ∈ {1, 2}) is the given one, and conditioned on h, the sets A 1 and A 2 are independent of one another. This can be done by first sampling h and then sampling A 1 and A 2 independently from the regular conditional probabilities. The union of A 1 and A 2 is then a new random set coupled with h. We denote this new random set by A 1∪ A 2 and refer to it as the conditionally independent union of A 1 and A 2 . The following lemma is analogous to [SS10, Lemma 3.6], Lemma 10.3. If A 1 and A 2 are local sets coupled with the GFF h on D, then their conditionally independent union A = A 1∪ A 2 is also local. Moreover, given A and the pair (A 1 , A 2 ), the conditional law of h is given by C A plus an instance of the FGF on D \ A.

An example of a local set
Certain level lines of the Gaussian free field are studied in [SS10] and shown to be local sets. We will show that certain level sets of fractional Gaussian fields with positive Hurst parameter are also local sets.
Let c 1 , c 2 > 0, let s > d/2 and let h be the FGF s on the unit ball B in R d with boundary values c 1 on the upper hemisphere, −c 2 on the lower hemisphere, and zero outside a compact set. Then there is a unique surface whose boundary equals between the boundary of the upper hemisphere and on which h = 0. This surface separates a region where h is positive and a region where h is negative. We call this interface the level set of h and denote it by L. To see that L is a local set, fix δ > 0 and let L δ be the intersection of D with the union of all closed boxes of the grid δZ d that intersect L. For each fixed closed set C, the event {L δ = C} is determined by h| C . Given a deterministic open set U ∩ C = ∅, the projection of h tȯ H s 0 (U) is independent of h| C . Thus L δ is local. Letting δ → 0, we see that L is local.
When r 1 and r 2 are far apart, G s (r 1 , r 2 ) is approximately a constant times (r 1 + r 2 ) 2H since 2 F 1 (a, b, c; z) approaches a constant as z → 0. So we see that long-range covariances of h are determined by H, while local covariances are dictated by the parameter s − 1/2. In the following proposition, we show that in fact s − 1/2 also governs the almost-sure regularity of sample paths of h. We prove such a statement only for s − 1 2 ∈ (0, 1), but we remark that in general the spherical average process is differentiable s − 1 2 − 1 times, and those derivatives are α-Hölder continuous for all α less than the fractional part of s − 1 2 .
where C 1 is some constant. Since h(r 1 ) − h(r 2 ) is Gaussian, (11.2) holds for all m ∈ N, for some constants C m . Applying the Kolmogorov-Chentsov continuity theorem with suitably large m, we conclude that h, and therefore also h, has a version which is almost surely α-Hölder continuous for all α < s − 1/2.

Background on spherical harmonic functions
We write the Laplacian in spherical coordinates as x d ] is said to be harmonic if ∆φ = 0. Suppose that φ is harmonic and homogeneous of degree k. Let f = φ| S d−1 , and note that we have φ(ru) = f (u)r k for all u ∈ S d−1 and r ≥ 0. Writing ∆φ = 0, using (11.3), and setting r = 1 yields In other words, f is an eigenfunction of We mention a few basic results about spherical harmonics that appear, for example, in [SW71, Chapter IV, §2]. Assume d ≥ 2, let A k be the set of homogeneous degree k harmonic polynomials on R d and let H k be the space of functions on S d−1 obtained by restricting functions in A k . An important property is that the spaces H k are pairwise orthogonal (for the L 2 (S d−1 ) inner product) and their union is dense in L 2 (S d−1 ). This means that we can define, for each fixed k, an orthonormal basis {φ k,j : 1 ≤ j ≤ dim(H k )} of H k which is the restriction of the harmonic polynomials {P k,j : 1 ≤ j ≤ dim(H k )} ⊂ A k , so that the collection of all φ k,j is an orthonormal basis of L 2 (S d−1 ) .
We will need the following important theorem concerning the behaviour of harmonic polynomials under the Fourier transform [Ste70, pg. 72]. We say that a function f : R d → C is radial if f (x) = f (y) whenever |x| = |y|. We occasionally abuse notation and write f (r) where f is radial and r ≥ 0, with the understanding that we mean f ((r, 0, . . . , 0)).
Theorem 11.2. Let P k (x) be a homogeneous harmonic polynomial of degree k in R d . Suppose that f is radial and that P k f ∈ L 2 (R d ). Then the Fourier transform of P k f is of the form P k g, where g is a radial function. Moreover, the induced transform F d,k ( f ) := g depends only on d + 2k. More precisely, we have F d,k = i k F d+2k,0 .
Applying the Fourier transform on both sides (which is the inverse Fourier transform evaluated at −x) and using the theorem again, we obtain where (−∆) s/2 R d+2k f is the fractional Laplacian on R d+2k acting on f interpreted as a function on R d+2k (that is, we define f (x) for x ∈ R d+2k to be Remark 11.4. Let P k,j f 1 and P k ,j f 2 ∈Ḣ s (R d ). Then We see that the right hand side of (11.5) (for (k, j) = (k , j )) can be rewritten as f 1 , where the radial functions f i are treated as functions defined on R d+2k (as described in the remark above). We thus have a unitary correspondence between elements For k ∈ N and 1 ≤ j ≤ dim(H k ), we define the Hilbert spaceḢ s k,j (R d ) to be the space of all functions of the form P k,j f where x → f (|x| R d+2k ) ∈ H s (R d+2k ) is radial and P k,j f ∈Ḣ s (R d ) . By 11.5, we see thatḢ s k,j (R d ) are orthogonal. In fact, they also spanḢ s (R d ): Lemma 11.5.Ḣ s k,j (R d ) are orthogonal subspaces spanningḢ s (R d ).
Proof. We only need to check the spanning condition. Since S(R d ) is dense inḢ s (R d ), it suffices to show that all g ∈ S(R d ) can be written as a linear combination of terms inḢ s k,j (R d ). To do this, we use the stated fact that {ω → φ k,j (ω) : k ∈ N, 1 ≤ j ≤ dim H k } a basis for L 2 (S d−1 ). We compute for every sphere of radius |x|: and see that g(x) = ∑ k,j ρ k,j (|x|)φ k,j (x/|x|) = ∑ k,j |x| −k ρ k,j (|x|)P k,j (x). Define g k,j (x) = |x| −k ρ k,j (|x|)P k,j (x) and let χ R (x) be the characteristic function of an annulus of radii 1/R and R, where R > 1. It is clear (by orthogonality of φ k,j (x/|x|)), thus by Fatou's Lemma g k,j ∈ L 2 (R d ). Hence, it follows that the Fourier transform g k,j exists and is in L 2 (R d ). Following the same reasoning as above with ξ → |ξ| 2s g k,j (ξ), we have that x → ρ k,j (|x|)P k,j (x) ∈Ḣ s k,j (R d ) as required.
Proof. Since each side of (12.3) is a bilinear form in f and g, it suffices to show that the formula holds with f = g. We simplify the left-hand side of (12.3) to obtain Changing variables for x + y in the second square-bracketed expression shows that the left-hand side of (12.3) is equal to which equals the right-hand side of (12.3) by Proposition 2.2.

The discrete fractional Gaussian field
In this section we define a sequence of discrete random distributions converging in law to the fractional Gaussian field FGF s (D), where s ∈ (0, 1) and D ⊂ R d is a sufficiently regular bounded domain. We follow the strategy of [Cap00] and prove convergence using a random walk representation of the field covariances. This method was introduced by Dynkin [Dyn80].
Suppose that D ⊂ R d is a bounded domain and s ∈ (0, 1). For δ > 0, define V δ := δZ d ∩ D. Recall that the zero-boundary discrete Gaussian free field (DGFF) is defined to be the mean-zero Gaussian field with density at f ∈ R V δ proportional to where C d is a constant and where we interpret the expression in parentheses as a quadratic form in the variables { f (x) : x ∈ δZ d ∩ D} by substituting zero for each instance of the variable f (x) for all x / ∈ D. Observing that the sum in (12.4) is a rescaled discretized version of the L 2 norm of the gradient of f , we define the zero-boundary discrete fractional Gaussian field DFGF s (D) by replacing this expression with a rescaled discretized L 2 norm of the fractional gradient of f . More precisely, we let , where x = (x 1 , . . . , x d ), and define h δ ∼ DFGF δ s (D) to be a Gaussian function h δ with density at where we interpret the expression in parentheses as a quadratic form in the variables { f (x) : x ∈ δZ d ∩ D} (as we did for the DGFF). Observe that this quadratic form includes long-range interactions, unlike the quadratic form for the GFF which includes only nearest-neighbor interactions. The constant C d,s is chosen so that the discrete FGF converges to the FGF with no further normalization-see (12.9) below to understand the role that this constant plays in the calculation.
We interpret h δ as a linear functional on C ∞ c (D) by setting (h δ , φ) := ∑ x∈V δ h δ (x)φ(x)δ d , for all φ ∈ C ∞ c (D). (12.5) To motivate (12.5), we note that the right-hand side is approximately the same as the integral of an interpolation of h δ against φ. The following theorem is a rigorous formulation of the idea that the DFGF converges to the FGF as δ → 0 when D is sufficiently regular. The idea of its proof is to compare a random walk describing the covariance structure of the DFGF to the 2s-stable Lévy process describing FGF covariances. Recall that D is said to be C 1,1 if for every z ∈ ∂D, there exists r > 0 such that B(z, r) ∩ ∂D is the graph of a function whose first derivatives are Lipschitz [CS98]. T 0 1 X t ∈A dt] for all Borel sets A ⊂ D. We have T < ∞ almost surely, and µ X,x is a finite measure-see the proof of Lemma 12.3 where a stronger statement is proved. By Lemma 12.3, µ X n ,x → µ X,x weakly. Since weak convergence implies convergence of integrals against bounded continuous functions, we have where the quantity denoted o(1) is uniformly bounded as x varies over the support of φ and tends to 0 as δ → 0 for each fixed x. Substituting (12.10) into (12.7) and using the convergence of the Riemann integral (as well as dominated convergence to handle the o(1) term), we obtain G(x, y) dx dy (12.11) as δ → 0, where G is the density of the occupation measure (that is, the Green's function) of Y. This Green's function is in turn equal to G s D (x, y) (see (4.2)), the Green's function of the fractional Laplacian [CS98]. Therefore, the right-hand side of (12.11) is equal to E[(h, f ) 2 ], as desired.
Lemma 12.3. Let (X n ) n≥1 be a sequence of processes in R d converging in law with respect to the J 1 metric to a symmetric α-stable process X. Let D ⊂ R d be a C 1,1 domain. If T is the hitting time of R d \ D, then the occupation measure of X T n converges weakly to the occupation measure of X T .
Proof. For n ≥ 1, denote by µ n the occupation measure of X T n : Recall that for probability measures, convergence with respect to π is equivalent to weak convergence [Bil99]. Since weak convergence of a sequence of finite measures (µ n ) n≥1 to a nonzero measure µ is equivalent to weak convergence of the normalized measures µ n /µ n (R d ) → µ/µ(R d ) along with convergence of the total mass (that is, µ n (R d ) → µ(R d )), we see that convergence with respect to π is equivalent to weak convergence for finite measures too. Therefore, it suffices to show that for all > 0 and A ⊂ R d , we have µ n (A) ≤ µ(A ) + and µ(A) ≤ µ n (A ) + . Since µ n (R d \ D) = µ(R d \ D) = 0, it suffices to consider A ⊂ D. For η > 0, define D η = {x ∈ D : dist(x, ∂D) > η}. For η > 0, define B η to be the event that X stopped upon exiting D η is contained in D η . By integrating the upper bound in Theorem 1.5 in [CS98], we conclude that B η has probability tending to 0 as η → 0. Furthermore, for each positive integer n, the event E n that |X n+1 − X n | is larger than the diameter of D has probability

Notation
We fix the relation H = s − d/2 for the definitions of the following spaces. We refer the reader to the referenced page numbers for the spaces' topologies.
The Schwartz space of real-valued functions on R d whose derivatives of all orders exist and decay faster than any polynomial at infinity.

13
S r (R d ) For r ∈ R, denotes S max(−1, r ) (R d ) 13 S k (R d ) For k ∈ {−1, 0, 1, . . .}, denotes the space of continuous linear functionals on S k (R d ). Equivalently, S k (R d ) may defined to be the space S (R d ) of tempered distributions modulo polynomials of degree less than or equal to k.

13
H s (R d ) The subspace of S H (R d ) consisting of functions whose Fourier transform ξ → f (ξ) is in L 2 (|ξ| 2s dξ)

T s (R d )
The closure of S H (R d ) inḢ −s (R d ). This space serves as a test function space for FGF s (R d ).
The closure of the space of restrictions to D of Schwartz functions under the metric d(φ, ψ) = φ − ψ Ḣ−s (D) . This space serves as a test function space for FGF s (R d ).