Smoothness of correlation functions in Liouville Conformal Field Theory

In this article we prove smoothness of the correlation functions in probabilistic Liouville Conformal Field Theory. Our result is a step towards proving that the correlation functions satisfy the higher Ward identities and the higher BPZ equations, predicted by the Conformal Bootstrap approach to Conformal Field Theory.

Quantization of the Liouville theory then amounts to defining the measure e −S L (X,g) DX on a space of generalized functions from C to R (a negative order Sobolev space) so that the observables F of the random field X are given by the path integral where Z is a normalization constant and DX denotes a formal infinite dimensional Lebesgue measure on the chosen space of generalized functions. The resulting theory exhibits conformal symmetry and is called the Liouville Conformal Field Theory (LCFT). The motivation for studying this theory comes from the hope to understand conformal metrics of the form e γX(z) g(z)|dz| 2 on C, where X is the random field with law (1.2).
The rigorous definition of the path integral (1.2) was given by David-Kupianen-Rhodes-Vargas in [5] by using Gaussian Multiplicative Chaos (GMC) methods. For us the relevant observables of the field will be the correlations of V α (z) := e αφ(z) , where φ(z) = X(z) + Q 2 ln g(z). Thus we consider where α i are real numbers (the Liouville momenta) with certain restrictions and z i ∈ C (the insertions). These are called vertex operators in the physics literature, and they are relevant for understanding the conformal metrics e γX(z) g(z)|dz| 2 . The correlations (1.3) with real weights α i are relevant for many conjectures related to scaling limits of random planar maps coupled to certain statistical physics models [5,14,20]. From the CFT point of view one would also want to define these for complex α i , but at this point it remains a challenge to give a rigorous probabilistic definition in the complex case. In the CFT language the fields e αφ(z) are supposed to be the primary fields of the LCFT when α belongs to the spectrum, which is supposed to be Q + iR. Many other quantities of the theory are then supposed to be expressible in terms of the correlation functions of the primary fields.
In the physics literature the quantum Liouville theory was first considered in the context of String Theory by Polyakov [16] as a building block for Liouville Quantum Gravity. Physics reviews of Liouville theory can be found in [15,22,21].
In addition to the path integral formulation, Liouville theory has also been studied (in the physics literature) by using the Conformal Bootstrap method, developed by Belavin-Polyakov-Zamolodchikov in [3]. One goal of the recent mathematical study of LCFT is to unify the path integral and the Conformal Bootstrap approaches, since the equivalence of these two formulations has been controversial even for physicists. Rigorous results in this direction can be found in [11,12,13]. For references on the Bootstrap see e.g. [21,22]. One mathematical consistency check of the equivalence of the Path Integral and the Conformal Bootstrap was done in [1], where it was shown that the one-point function of LCFT on the torus agrees with the predictions of the Conformal Bootstrap in the large moduli limit. In [2] the authors derived fusion estimates for the four point function of LCFT on the Riemann sphere and showed that they agree with the predictions of the Conformal Bootstrap.
In this article we will establish the smoothness of the correlation functions (1.3) with respect to the insertions (z i ) N i=1 , which is required for rigorously proving the CFT structure of Liouville theory, predicted by the Conformal Bootstrap approach (see Section 1.2). We now state the exact form of our theorem.
and that g is any diagonal Riemannian metric g(z)|dz| 2 . Then the functions The correlation functions were shown to be C 2 on this domain in [11] and our smoothness proof will be partially based on an iteration of their C 1 -argument.
Remark 1.2. The generalization of Theorem 1.1 for arbitrary Riemannian metrics is simple since an arbitrary Riemannian metric g ′ can be written as g ′ = ψ * g where ψ * is a pullback of a diffeomorphism ψ : C → C and g is a diagonal metric. The correlation functions are supposed to satisfy diffeomorphism covariance from which the generalization follows. Of course this requires defining the Liouville theory for arbitrary metrics. On surfaces with genus 2 or higher this has already been done in [8].
1.2. Perspectives. The smoothness of the correlation functions is needed for the program of deriving the Conformal Bootstrap postulates from the path integral. The Conformal Bootstrap approach predicts that the correlation functions (1.3) appear in certain partial differential equations of arbitrarily high order.
The first set of equations are the Conformal Ward identities. These are supposed to emerge from a variation of the background metric g. More precisely, let g = 2 i,j=1 g ij dx i ⊗ dx j be a Riemannian metric and fix some nice functions (f ij ) 2 i,j=1 . We define smooth variations of this metric by g ij ε = g ij + εf ij where g ij are the components of the inverse matrix of (g ij ) 2 i,j=1 . Then we expect where vol g (d 2 z) is the volume form of g and T ij is called the stress-energy tensor. In a CFT two of the components of T are nontrivial, see [7]. In the (z,z) coordinates they are T (z) := T zz (z) andT (z) := Tzz(z). Then according to  the Conformal Ward identities for any non-negative integers M and N are Here ∆ α = α 2 (Q − α 2 ) and c L = 1 + 6Q 2 is the central charge of LCFT. The other nontrivial componentT is supposed to satisfy similar identities. In [11] these were proven for M ∈ {0, 1}. In the proof the authors defined via a regularization procedure and then computed (1.5) by Gaussian integration by parts. The computations get quite lengthy and thus to prove the identities for all M one should take the variational relation (1.4) as the definition of T .
In the variational computation one should use the fact that two smooth metrics g ′ and g on C are related by where ψ : C → C is a diffeomorphism, ψ * is the associated pullback and ϕ : C → R is a smooth function. This means that on the Riemann sphere two metrics are equivalent modulo a diffeomorphism and a conformal factor e ϕ . On higher genus surfaces one has to also take into account the moduli space.
The dependency of N i=1 V α i (z i ) g on ϕ is explicitly given by the Weyl anomaly [5] and thus the differentiation with respect to this factor is easy. The only thing left to do is to investigate the ψ dependency and for this the smoothness of the correlation functions is needed. In [7] this computation is done in a general axiomatic CFT setting where the author assumes the Weyl anomaly, diffeomorphism covariance and some regularity for the correlation functions.
The Ward identities are needed for the construction of representations of the Virasoro algebra. For this the canonical construction of the Hilbert space associated to the LCFT should be carried out and then the generators of the Virasoro algebra should act on some dense subspace of this space, see [14]. This will be carried out in a future work.
The other set of partial differential equations that the correlation functions are supposed to satisfy are the Belavin-Polyakov-Zamolodchikov equations (BPZ equations). More precisely the (r, 1)-order degenerate field where the differential operator D r is given by D r = r k=1 (n 1 ,...,n k )∈N k n 1 +...+n k =r The degenerate field of order (1, r), given by V − 2(r−1) , is supposed to satisfy a similar equation where γ 2 gets replaced by 2 γ . In [11] the BPZ equations were proven for the (2, 1) and (1, 2) degenerate fields by using Gaussian integration by parts. The BPZ equations are essential for proving integrability of LCFT. They were used in the proof of the DOZZ-formula [12,13] for the 3-point function of LCFT on the sphere, and after this similar methods were used for obtaining integrability results for one dimensional GMC measures on the unit circle [17] and on the unit interval [18]. The unit circle computation was based on a boundary LCFT, which is defined in [9]. The connection between the unit interval computation and LCFT is not clear, although the methods used are very similar to the methods used in [12].
Acknowledgements. I would like to thank Antti Kupiainen and Yichao Huang for many fruitful discussions and for giving comments on the manuscript. This project has received funding from the ERC Advanced Grant 741487 (QFPROBA).

Mathematical background
In this section we quickly review the rigorous definitions behind the probabilistic approach to quantum Liouville theory. Similar discussions can be found in [14,23,5,11,12].
2.1. Gaussian Free Field and Gaussian Multiplicative Chaos. In this paper we will work with the round metric given by For this choice the scalar curvature is constant R g (z) = 2 for all z ∈ C.
The usual starting point for defining the measure e −S L (X) DX is to separate the free field part |∂ z X| 2 dz and to think of the measure as Now the e − 1 π |∂zX(z)| 2 d 2 z DX factor can be naturally thought of as the (non-normalized) distribution of the Gaussian Free Field (GFF), which formally is a Gaussian process (X(z)) z∈C with covariance and with vanishing mean over the Riemann Sphere In other words C g is the zero mean Green function of the Laplace-Beltrami operator ∆ g . From (2.3) we see that the variance EX(z 1 ) 2 is infinite. This means that in reality the GFF is a random generalized function rather than a function. More precisely, the probability law of X lives on the negative order Sobolev space H −1 ( C, g) which is the continuous dual of H 1 ( C, g) (see Appendix A for definitions). Thus X is a Gaussian process (X(f )) f ∈H 1 ( C,g) with covariance given by The fact that X is a random generalized function poses a problem in defining the measure (2.2) since we would like to think of it as the GFF measure multiplied by some Radon-Nikodym derivative, but now the term e γX(z) becomes ill-defined since the exponential of a generalized function is not defined. This is where the theory of Gaussian Multiplicative Chaos steps in, since it provides a framework for defining exponentials of logarithmically correlated Gaussian fields. This work goes back to Kahane [10]. For a more recent review see [19].
We define e γX(z) g(z)d 2 z to be the ε → 0 limit of the measures where d 2 z is the Lebesgue measure on C and X ε (z) denotes a regularization which we choose to be a smooth mollification (another common regularization is the circle average). More precisely let ρ be a non-negative C ∞ (R) function with compact support and define ρ ε (z) = ε −2 ρ(|z| 2 /ε 2 ). Then the regularization of X is defined by X ε = ρ ε * X. We also adapt the notation In [4] it was shown that for γ ∈ (0, 2) the measures (2.5) converge weakly in probability as ε → 0 and we denote the limit by M γ (d 2 z). This measure is called the GMC associated to X with respect to the measure g(z)d 2 z. The following result goes back to Kahane [10].
Proposition 2.1. (Kahane Convexity Inequality) Let X and Y be two continuous Gaussian fields on C such that for all x, y ∈ C Then for all convex F : R + → R with at most polynomial growth at infinity and f : When applying Kahane convexity to the GFF one has to use the regularized field X ε because of the continuity assumption, but usually this is not a problem.

Liouville correlation functions.
We start by defining the path integral (1.2) by setting Here E is the expectation with respect to the GFF X and the integral over c corresponds to a zero mode 1 . The observable F : H −1 ( C, g) → R is arbitrary as long as the integral converges. This definition is the same as the one given in [12], and differs slightly from the original definition in [5]. 1 Now c + X is a field where X is the GFF with zero mean and c is distributed according to the Lebesgue measure. Thus the law of c + X is the pushforward of dc ⊗ dµX under the map (c, X) → c + X where dµX is the law of X. Note that this is not a finite measure.
Another way to view this is to recall that the GFF has the series representation The formula comes essentially from plugging in the field φ = c + X + Q 2 ln g into a path integral with the regularized action ( 1 π |∂ z φ| 2 + µε γ 2 2 e γφε ) d 2 z with the Euclidean metric g ≡ 1 and adding an integration over the zero mode c and taking the ε → 0 limit. Absorbing the g dependency to the field via adding the Q 2 ln g term is common, see for example the discussion in Section 3.4 of [20].
Recall that V α (z) = e αφ(z) , so the vertex operators correspond to the choice F (X) = e αX(z) in (2.6). To define F rigorously for this choice of F one has to regularize the exponential in a similar manner as in the GMC definition above. Thus we define .
From now on we denote By G(z) and G(x; z) we denote the corresponding ε → 0 limits. In [5] it was shown that G(z) exists and is non-zero if and only if the momenta α i satisfy the Seiberg bounds In particular this implies that we need N 3. By performing the c integral and using the Cameron-Martin theorem we arrive at Thus the correlation functions can be expressed as integrals of explicit functions against the GMC measure. This was initially shown in [5] and by using this formula it is possible to derive fusion estimates that tell us the singular behaviour of G when two of the points z i merge [11]. In the smoothness proof we will use a slight generalization of the fusion estimate from [11] (see Section 3). The correlation functions satisfy the following useful integral identity.
where Xn are i.i.d. standard Gaussians, and en and λn are the eigenfunctions and eigenvalues of −∆g, respectively. Then the zero mode corresponds to adding the constant eigenfunction e0 into this series. This has eigenvalue 0 so the term in the series would be X 0 √ 0 e0. If we interpret X 0 √ 0 as a Gaussian with infinite variance, that is, the Lebesgue measure, then we end up with the random field c + X.
The same identity holds when G ε is replaced by G.

Proof. After the change of variables
On the other hand, by differentiating (2.7) with respect to µ we get The identity for the ε → 0 limits follows from Dominated Convergence, since in [11] it was shown that G ε has an integrable dominant uniformly in ε.
2.3. The first derivative of the correlation functions. Throughout this section we assume that the insertion points are distinct, that is, z ∈ U N = {(z 1 , . . . , z N ) ∈ C N : z i = z j , ∀ i = j}. In [11] it was shown by using Gaussian integration by parts that A priori we do not know if the ε → 0 limit of the integral exists. In [11] the convergence was shown by studying the integral transform

Indeed by integration by parts we have
If we contour integrate this relation along the contour ∂B(z i , r), we get where we used the Residue Theorem. On the other hand, the left-hand side can be written as We end up with the fundamental identity The last term we can write as ε,ε is antisymmetric in x and y. The crucial observation is that the ε → 0 limit of the integrand will turn out to be absolutely integrable. Now (2.11) leads to where the last term comes from integrating by parts the left-hand side of (2.11). From this we can take the ε → 0 limit and we see that z → G(z) is C 1 (checking the ∂z i -derivative works the same way). In the end we obtain is convergent by the two-point fusion estimate from [11] and the following lemma.
Lemma 2.3. Let K be a bounded set containing the origin and B = B(0, r) for some r > 0. Then Proof. Denote x = (x 1 , x 2 ) and y = (y 1 , y 2 ). The integral is convergent if and only if the integral is convergent (or we could use squares instead of the balls in our smoothness proof). We can also use the norm |x−y| = |x 1 −y 1 |+|x 2 −y 2 | since all norms are equivalent. After a change of variables The singular part of this integral is around the origin, and computing the part over a small square of the first term yields a < 3.
The fusion estimate (Proposition 5.1 from [11]) tells us that when x and y fuse and the z i 's are away from x and y, we have where ζ > 0 depends on γ. For higher derivatives we need a generalization of this estimate for a fusion of n separated pairs of points.
Finally, the boundary integral term in (2.12) is convergent since the map x → G(x; z) is continuous when x stays away from the z i 's. The conclusion is that the limit of conditionally convergent integrals on the left-hand side of (2.12) is expressible as a sum of absolutely convergent integrals. Iteration of a process like this is our strategy for the smoothness proof: we will differentiate the resulting absolutely integrable terms and the derivatives are given by conditionally convergent integrals. Then we simplify these conditionally convergent integrals into sums of absolutely convergent integrals and repeat.

n-pair fusion estimate
In this section we prove a result concerning the singular behaviour of G(z) when multiple pairs of the points z i merge.
Proof. Our proof will be based on the n = 1 proof which was done in [11]. From the proof of Proposition 5.1 in [11] we get (see Appendix B for details) and D j is an annulus around x j with radii R j and |x j − y j | (we can assume these radii to be so small that the distance between D j and ∂A j is positive). The moment q is given by where α i are assumed to satisfy the Seiberg bounds (2.8). We define the notation (1) First assume that the (W j ) j are completely independent. Then we estimate (using first the estimate ( n j=1 W j ) q n j=1 W q/n j and then independence and finally the proof of Proposition 5.1 in [11], summarized in Appendix B) The series k is convergent since The inequality follows from the Seiberg bound (2.8). The claimed estimate follows by simplifying the exponent of |x j − y j |. (2) In reality W j are not independent. However, we can reduce everything to the independent case by using Proposition 2.1 and the spatial Markov property of the GFF from Appendix A. Decompose X = n j=1 X j + ϕ where X j is a Dirichlet GFF in a simply connected neighbourhood U j of D j , independent of the other X k 's, and ϕ is smooth in ∪ j U j (the harmonic extension of X from the boundaries ∂U j into the domains U j ). We can choose the U j 's to be such that they are disjoint and ∂D j ⊂ U j . Inside each U j the component X j is independent of ϕ. Thus we are using Proposition A.1 multiple times, each time using Harm(U j , g) ⊕ H 1 0 (U j , g) as the Hilbert space decomposition. For x, y ∈ D 1 ∪ D 2 we have EX(x)X(y) = EX 1 (x)X 1 (y) + EX 2 (x)X 2 (y) + EX 1 (x)ϕ(y) + Eϕ(x)X 2 (y) + Eϕ(x)ϕ(y) c δ + E(X 1 (x) + X 2 (x))(X 1 (y) + X 2 (y)) , for some constant c δ , since EX 1 (x)ϕ(y) = 0 for x, y ∈ U j and bounded for x ∈ D 1 , y ∈ D 2 . Also sup x,y∈D j |Eϕ(x)ϕ(y)| < ∞ since ϕ is bounded inside D j . Let N be an independent centered Gaussian with variance c δ . Then and M γ is the chaos measure of X = n j=1 X j +ϕ and M γ is the chaos measure of n j=1 X j . Since N is Gaussian, the exponential moment exists and thus the factor E[e −qγN ] is just a finite constant. Now we can estimate as in the independent case (3.2).

Sketch of the proof
In this section we sketch the smoothness proof. The actual detailed proof is given in the next section. Fix some insertion z i ∈ {z 1 , . . . , z N } and define B j = B(z i , r/j). By A j we denote a closed annulus containing ∂B j such that all the A j 's are disjoint. We define the functions Recall the derivative formula 4.1. Second derivative. Next we compute the second derivative of the second term, which is the most problematic term. We have We have to simplify the terms with 1 z i −x k , k = 1, 2, 3, since they are not absolutely convergent. In all the cases the simplification will work in essentially the same way so we focus on the case k = 3. The simplification follows from an analogue of the identity (2.11): We want to solve for the j = i term in the first sum. The other terms in the first sum are automatically convergent. In the last term we insert 1 = 1 B 2 (x 4 ) + 1 (B 2 ) c (x 4 ) and it becomes since the integral over 1 B 2 (x 4 ) vanishes by antisymmetry. Next we deal with the 1 x 3 −x j , j = 1, 2 terms. We proceed similarly as above, that is, we insert 1 = 1 B 2 (x 1 ) + 1 (B 2 ) c (x 1 ). We symmetrize the integral over the first factor and get Next we want to argue that when x 1 , x 3 ∈ B 2 , the factor F 1 (x 1 , x 2 ) − F 1 (x 3 , x 2 ) behaves like (x 3 − x 1 ) multiplied by something that is integrable against G. Indeed, we have and both of these factors are smooth (and bounded) when x 1 , x 3 ∈ B 2 . This fact is quite intuitive since the singular behaviour of the function F 1 (x, y) happens on the circle ∂B 1 and thus the singularity is gone when x is restricted to the smaller ball B 2 . From this we infer that where H is bounded. When we insert this to (4.3) we get a nice integrable term. The 1 (B 2 ) c (x 1 )term is Terms like this are shown to be integrable by using the estimate because after inserting this, we get terms where one variable has singularity only on one of the circles ∂B j (by disjointness of the annuli A j ) and since the radius j is different for each F -factor, none of these singularities "stack". This is proven in Proposition 5.6.
Next we integrate by parts the left-hand side of (4.2). We get To show that this is convergent we again split the integral into 1 A 1 ∩B 1 (x 1 )1 A 1 ∩(B 1 ) c (x 2 ) and the complement. In the complement F is bounded so the integral is clearly absolutely integrable. For the other part our fusion estimate 3.1 implies where C depends on δ = min N i,j=1;i =j |z i − z j | and ζ > 0. This is integrable.

4.2.
Higher derivatives. Denote x = (x 1 , . . . , x 2n ). When we start to compute higher order derivatives, integrals of the form will start to appear. For each new derivative ∂ z i we add new γ-insertions (insertions with Liouville momentum α = γ) to G and a factor F j (x a , x b ), with some indices a = b, to F , where the index j tells the radius of the ball appearing in the definition of F j (we increment j after each differentiation so all the balls have different radii). The most singular term appearing in the formula for ∂ n z i G(z) will then be where x = (x 1 , . . . , x n ) and y = (y 1 , . . . , y n ). The proof of convergence of this integral is essentially the same as in the C 2 -case we did above: just use the estimate (4.4) and the n-pair fusion estimate 3.1. Then to show differentiability of this integral in the z i 's, we have to take the derivative of G in the above integral and this leads us to investigate the integrals where k ∈ {1, . . . , 2n + 1}. This we simplify by using the same integration by parts argument as in the C 2 -case. Thus we write the integral in two different ways (by using the derivative formula ∂ x k G and by integration by parts). Then we can solve for (4.5). In the resulting expression we have to simplify the integrals . The latter part produces the integral which converges by the fusion estimate. The remaining part we symmetrize (as before). Note that the factors F j that depend on x k are smooth and bounded in B n+1 and the same holds for x l .
Thus after symmetrizing, the parts which depend on these variables produce a O(x k − x l )-term and the parts that do not depend on these variables are invariant. Thus we get where ϕ(x) is bounded. Integrability of this follows from the fusion estimate.

Proof of smoothness
In this section we give the detailed proof of Theorem 1.1. We want to iterate the computation we did in the C 1 proof. What we will observe later is that the derivatives of G(z) can be expressed as integrals of singular functions against G with additional γ-insertions.
Let N 3 denote the amount of insertion points in the correlation function and n ∈ Z + the amount of partial derivatives we are taking. Fix a sequence of numbers i(j) ∈ {1, . . . N } where j ∈ {1, . . . , n}. We study the nth order partial derivative where z = (z 1 , . . . , z N ) ∈ U N . Let B j = B(z i(j) , r/j), A j be the corresponding closed annulus containing ∂B j so small that (A j ) n j=1 are disjoint (see Figure 1, which describes the notation for the partial derivative ∂ z 1 ∂ z 2 ∂ z 1 G(z 1 , z 2 , z 3 , z 4 )), and To set up a suitable induction, we define the following function classes.
We will need the following two properties of this function class.
Proof. Follows from the definition of F j and property 4 in Definition 5.1.

Simplification.
Remark 5.3. Strictly speaking we should do the following computations with the regularized functions G ε and 1 (xa−x b )ε,ε , but this does not affect any of the algebraic manipulations we do, and interchanging limits and integrals works in the end easily by our integrability result (Proposition 5.6) and the estimates for sup ε>0 G ε (x, y; z) in [11]. Thus we choose to not write all the epsilons in our computations.
We symmetrize the first term to get By definition of F n+1 the function F is smooth in x a and x b in the integration domain and thus the difference is not a function of x a and x b , and ϕ is bounded everywhere and smooth outside of the circles ∂B j . Indeed, the function H will consists of the factors F j in F that are not functions of x a and and clearly the function F n+1 (x a , x b )F (x) belongs to F n+1 .
Next we show that the conditionally convergent integrals that pop up when we compute derivatives of G can be expressed in terms of absolutely convergent integrals.
We want to solve for the j = i term in the first sum. The j = i terms need no simplification. The rest of the terms simplify correctly by Lemma 5.4. When we integrate by parts the left-hand side we get where H ∈ F n+1 does not depend on x k and ϕ is bounded everywhere and smooth outside the circles ∂B j . Thus The boundary integral is convergent because by Lemma 5.7.

5.2.
Integrability. Next we show that the integrals appearing in Proposition 5.5 are absolutely convergent.
Proposition 5.6. Let F ∈ F n . Then Proof. Let a(j), b(j) be as in Definition 5.1. We may assume that where J ⊂ {1, . . . , n}. Since ϕ is globally bounded, we may also assume that ϕ ≡ 1. Next, we can apply KPZ-formula 2.2 to obtain i.e. we integrate away the γ-insertions which do not appear in the function F . We split the integrals and estimate as follows Each of the x-variables appear in only one of the F j -factors since otherwise the integrand vanishes (because A j are disjoint). Thus we can use the n-pair fusion estimate (Lemma 3.1) to get This converges by Lemma 2.3.

Boundary terms.
We still have to show that the boundary term appearing in Proposition 5.5 is integrable and that the boundary terms are differentiable and satisfy the analogues of Propositions 5.5 and 5.6.
Proof. This follows from the fact that in the proof of Lemma 3.1 the δ-dependent constant C δ satisfies sup δ∈K |C δ | < ∞ whenever K is a compact set disjoint from the origin. This is easy to see since by taking a bit more care of C δ one sees that it can be chosen to be Cδ a for some constant C and a < 0.
The above lemma says that the boundary integrals that we see are integrals of bounded functions over compact sets so they converge. In addition to this we need that the boundary integral terms appearing in our iteration are differentiable with respect to the insertions z i . This follows from the following lemma. Then for k ∈ L ∪ {2n + 1} we have For k ∈ L c we have the same formula without the boundary integral term since Proof. The case k ∈ L ∪ {2n + 1} is exactly the same as for Proposition 5.5 and the case k ∈ L c follows from the definition of F n .

5.4.
Proof of Theorem 1.1. By combining the derivative formulas (2.10) and (2.12) with Proposition 5.5 and the corresponding results for the boundary terms in Section 5.3, we see that where C J,k are some constants and F J,k ∈ F k is a linear combination of functions of the form where ϕ is bounded and a(j) = b(j) are some arbitrary choice of indices. Combining this with Proposition 5.6 together with Section 5.3 we see that all these integrals are absolutely convergent.
Taking ∂z i -derivatives works the same way since in the derivation of the derivative formula (2.10) one uses Gaussian integration by parts, which leads to terms containing derivatives of the form ∂ x C g (x, y) of the correlation of the GFF. When computing ∂z i instead of ∂ z i the derivative ∂ x C g (x, y) gets replaced by ∂xC g (x, y) and the essential term in C g (x, y) is ln 1 |x−y| which is symmetric in x − y andx −ȳ. So in the end everything works the same way in the ∂z i case.
We have established the smoothness in the case of the round metric (2.1). The generalization for any diagonal metric follows from the Weyl anomaly (Theorem 3.11 in [5] Now clearly if we have smoothness in the metric g, then we get smoothness for any metric e ϕ g in the same conformal class, that is, for any diagonal metric.

Appendix A. Decompositions of the Gaussian Free Field
Let g be a Riemannian metric on the Riemann Sphere C. The associated Laplace-Beltrami operator is given by where (g ij ) 2 i,j=1 are the components of the inverse of g. The operator −∆ g has eigenfunctions (ϕ i ) ∞ i=0 with non-negative eigenvalues (λ i ) ∞ i=0 . The eigenvalues satisfy 0 = λ 0 < λ 1 < λ 2 < . . . and the eigenfunctions (ϕ i ) ∞ i=0 form an orthonormal basis of L 2 (g). We define the Sobolev spaces H s ( C, g) on the Riemann sphere by We denote the subspace of zero mean elements of H 1 ( C, g) by H 1 0 ( C, g) where vol g is the volume form of g. For s < 0 this means the elements satisfying f, 1 = 0 where ·, · is the dual bracket. From (A.1) we see that this subspace corresponds to the elements satisfying f 0 = 0. From the sequence representation it is easy to see that the continuous dual of H s ( C, g) is H −s ( C, g) and the continuous dual of H s 0 ( C, g) is H −s 0 ( C, g). The zero mean spaces become Hilbert spaces when endowed with the inner product Note that in the case s = 1 the Sobolev norm agrees with the Dirichlet energy , where ∇ g is the g-gradient. If g ′ = e ϕ g is another metric in the conformal class of g, then H 1 0 ( C, g ′ ) = H 1 0 ( C, g) by conformal invariance of the Dirichlet energy. Furthermore, the covariance of the zero mean GFF satisfies for any f, h ∈ H 1 0 ( C, g). The duality between H 1 0 ( C, g) and H −1 0 ( C, g) and the above equality then lead to the following decomposition result.
Proposition A.1. Let p : H 1 0 ( C, g) → H 1 0 ( C, g) be an orthogonal projection. Then the zero mean GFF X can be decomposed as where f ∈ H 1 0 ( C, g) and X 1 and X 2 are independent Gaussian processes satisfying The covariances of X 1 and X 2 are given by Proof. Define X 1 (f ) := X(pf ) and X 2 (f ) := X(f − pf ). Then it is clear that X = X 1 + X 2 . Independency follows from the following computation The projection p : H 1 0 ( C, g) → H 1 0 ( C, g) induces an orthogonal projection p * on H −1 0 ( C, g) (the adjoint) by g)) .
The formulas for the covariances of X 1 and X 2 follow directly from the definitions.
If the projection p has an integral kernel K, that is, then the covariance kernels of X 1 and X 2 are explicitly given by . The corresponding decomposition of the GFF is X = ϕ + X D where ϕ is harmonic inside D (the harmonic extension of X from ∂D into D) and X D is the GFF in D with the boundary condition X| ∂D = 0 (the Dirichlet GFF). In Section 3 we will use this decomposition. It is often called the spatial Markov property of the GFF. The decomposition is useful when estimating the correlation E[X(x)X(y)] when x and y are taken from two sets separated by a positive distance, see the end of proof of Proposition 3.1 for the concrete computation. The corresponding decomposition of the GFF is X(z) = X r (|z|) + Y (z) where t → X r (e −t ) − X r (1) is the Brownian motion and Y is a Gaussian process called the lateral noise, see [6].
Appendix B. Lemma for the fusion estimate In this section we work with the GFF with zero mean over the unit circle, which we denote by X 0 . It has the covariance C 0 (x, y) = ln 1 |x − y| + ln |x|1 {|x| 1} + ln |y|1 {|y| 1} . (B.1) Changing the zero mean GFF to the zero circle average GFF corresponds to shifting the constant c in c + X. Indeed if X is the GFF with zero mean over the whole complex plane, then The term X − 2π 0 X(e iθ ) dθ 2π can be identified as the zero circle average GFF X 0 . We use the radial decomposition of the GFF X 0 (x) = X r (|x|) + Y (x) from Example 2. Plugging this into the chaos measure M 0 γ of X 0 we get where M γ (d 2 x, Y ) is the GMC measure of the Gaussian field Y . Inside the unit disk (the purely log-correlated region of X 0 ) integrals of the GMC measure can now be written as where B s = X r (e −s ) is a Brownian motion and µ Y is independent of (B s ) s 0 .
Recall the formula , B(α) is a constant depending on γ and the α i 's and We want to derive a fusion estimate for G(x, y; z) in the case when x j merges with y j while the pair (x j , y j ) stays away from all the other insertions. Using (B.3) we get G ε (x, y; z) n j=1 1 A j ∩B j (x j )1 A j ∩(B j ) c (y j ) C(δ, α, g) where A j and B j are as in the beginning of Section 5 and and D j is an annulus with center at x j , inner radius |x j − y j | and outer radius R j . We want to choose such R j that for x ∈ D j we have We can assume that |x j − y j | is really small for every j ∈ {1, . . . , n}. Then we choose any R j such that D j ⊂ A j and the above estimate holds. We can also estimate By combining these two estimates we get (B.4).
Next we derive another estimate that we will need. Inside D j we have |x − y j | ε C|x − x j | ε . Also, the correlation of the regularized field X ε satisfies where C is uniform in ε, and thus we can use the Kahane Convexity Inequality 2.1 to pass to the non-regularized measure M 0 γ (d 2 x). Without loss of generality we can assume that the points x j and y j fuse at the origin. By using the radial decomposition (B.2) of the GFF about the origin we get where P s = B s + (2γ − Q)s. We split the integral by using the following events We estimate the resulting integrals using the following lemma which is a special case of Lemma 6.5 in [11].