Recursive construction of the operator product expansion in curved space

I derive a formula for the coupling-constant derivative of the coefficients of the operator product expansion (Wilson OPE coefficients) in an arbitrary curved space, as the natural extension of the quantum action principle. Expanding the coefficients themselves in powers of the coupling constants, this formula allows to compute them recursively to arbitrary order. As input, only the OPE coefficients in the free theory are needed, which are easily obtained using Wick's theorem. I illustrate the method by computing the OPE of two scalars $\phi$ in hyperbolic space (Euclidean Anti-de Sitter space) up to terms vanishing faster than the square of their separation to first order in the quartic interaction $g \phi^4$, as well as the OPE coefficient $\mathcal{C}^{\mathbb{1}}_{\phi \phi}$ at second order in $g$.


Introduction
The operator product expansion, first conceived by Wilson [1], has found many applications in quantum field theory, especially conformal theories and in the context of the AdS/CFT correspondence [2,3]. The OPE is the statement that [1,4,5] where O A denotes a composite operator, the expectation value · Ω is taken in any suitable (interacting) state Ω which may include other (spectator) fields, where the sum on the righthand side runs over all composite operators in the theory, and where the OPE coefficients C B A 1 ···A k (x 1 , . . . , x k ; y) are independent of the state Ω, such that one really has an "expansion of operators". The concrete meaning of ∼ depends on the theory: while for conformal field theories it has been shown that the OPE is convergent [6][7][8][9], such that one can replace ∼ by =, in general it is only supposed to hold as an asymptotic expansion. That is, if one scales the insertion points x i to the expansion point y and includes in the sum on the right-hand side of (1.1) all operators O B up to a fixed dimension [O B ] ≤ ∆, the difference between left-and right-hand side vanishes like where d is the largest distance between y and any of the x i .
On the other hand, in recent work by Holland, Hollands and Kopper [10,11] it has been proven that the OPE (1.1) actually also converges for (Euclidean) scalar λφ 4 theory, to all orders in perturbation theory. Since the rate of convergence depends on the order of perturbation theory (and becomes worse at higher orders), this result is not as far-reaching as the non-perturbative convergence in conformal field theories, but shows that the OPE is much better behaved as originally expected. It seems thus possible to define a quantum field theory not by its correlation functions (obtained, e.g., from a classical Lagrangian in perturbation theory), but instead by the OPE coefficients which encode the algebraic properties of the theory and the one-point expectation values (form factors) O B (y) Ω , which determine the quantum state Ω. For this, one needs in addition certain factorisation and associativity conditions on the OPE coefficients. The factorisation conditions are obtained by performing the OPE (1.1) for a subset of the operators first, and then a second one for the remaining operators, which leads to (1.2) while the associativity conditions come from demanding that two different ways of performing such a factorisation agree. These conditions hold again as equalities in conformal field theories [7,9] for suitable configurations of the points x i , y and z (where the associativity conditions lead to crossing symmetry constraints), and have also been proven for scalar λφ 4 theory to all orders in perturbation theory [12,13]. It remains to give an algorithm to determine the OPE coefficients themselves. In free theories, they can be obtained directly by defining composite operators as the normalordered ones, using Wick's theorem and expanding the normal-ordered products around the expansion point [12,14]. For example, in flat space one has φ(x)φ(y) = :φ(x)φ(y): G + G(x, y)1 = G(x, y)1 + :φ 2 (z): G + [(x − z) µ + (y − z) µ ]:(φ∂ µ φ)(z): G + · · · , (1.3) where G(x, y) = φ(x)φ(y) 0 is the free two-point function and normal-ordering is performed with respect to G, and can read off the two-point OPE coefficients In general interacting theories, the textbook way to compute the coefficients is to take certain high-momenta limits of correlation functions [14], which then gives the most singular parts of the OPE coefficients, the ones that contribute most to the expansion (1.1) for small separations. Obviously, this contradicts the idea to determine quantum field theories by OPE coefficients and form factors without taking recourse to correlation functions, and other constructions are needed. For conformal field theories, conformal symmetry and associativity place severe restrictions on the OPE coefficients. This is especially fruitful in two dimensions (where the above conditions can be formulated using vertex algebras [15][16][17]), and one can nonperturbatively construct a large class of models (see, e.g., [18][19][20][21][22] for very general results). In higher dimensions, these conditions are less restrictive, but can still be used to constrain and even numerically construct conformal theories, an approch known as conformal bootstrap [6,18,[23][24][25][26][27][28][29]. For non-conformal field theories, no such strong restrictions on the OPE coefficients exist, and other ideas are needed.
One natural condition that can be imposed in a general field theory on the OPE are the (interacting) equations of motion, e.g., . . , x k ; y) (1.5) for a massless scalar λφ 4 interaction in flat space. Using associativity and this condition, an explicit recursive construction of the OPE coefficients as a power series in λ was given by Holland [30], which involves an inversion of ∂ 2 at each step and is technically involved. Much more useful would be a formula for the dependence of OPE coefficients on coupling constants that only depends on the coefficients themselves, involving the composite operator conjugate to the coupling constant. Namely, expanding such a formula in perturbation theory one could compute the OPE coefficients recursively in perturbation theory, using as input only the OPE coefficients of lower orders. Such formulas are closely related to Schwinger's quantum action principle expressing the change in correlation functions due to a change in couplings, which was shown to hold in renormalised perturbation theory (with an appropriate choice of renormalisation conditions) for both massive and massless fields in flat space by Zimmermann, Lowenstein, Schroer, Gomes and Lam [31][32][33][34][35]. The possibility of extending the quantum action principle to OPE coefficients was mentioned by Wilson [1]. That the counterterms in a renormalized action principle involve the OPE coefficients was shown by Sonoda [36], who then derived from this a variational formula for the dependence of the OPE coefficients on coupling constants [37]. A similar formula was also obtained by Guida and Magnoli [38] and Bochicchio and Becchetti [39,40], but all these formulas either still need to be renormalised or depend on correlation functions (in addition to OPE coefficients). The first fully renormalised variational formula that only involves the OPE coefficients themselves was derived (in perturbation theory) for scalar field theory in flat Euclidean space by Holland and Hollands [12,13] and generalised to gauge theories by Holland and the author [41].
In this paper, I derive an analogous formula valid also for the case of more than one coupling, and in curved space. Consider an interaction Lagrangian where dx ≡ √ det g d n x in an n-dimensional Riemannian space with metric g µν , and the couplings g A are in general position-dependent. The main result of this paper is the fol-lowing: there exists a renormalisation scheme such that the OPE coefficients fulfill (3.38) δ δg D (z) C B A 1 ···A k (x 1 , . . . , x k ; y) = −C B DA 1 ···A k (z, x 1 , . . . , x k ; y) . . . , x k ; y) C B DC (z, y; y) C C A 1 ···A k (x 1 , . . . , x k ; y) (1.7) In this formula, the first term involving the OPE coefficient C B DA 1 ···A k with an additional insertion of the interaction operator O D corresponding to the coupling g D is the naively expected one which would also follow from a formal path-integral treatment, with the minus sign arising because I work in a Riemannian setting where the action enters in the form exp(−S) instead of exp(iS). The terms in the second line are UV counterterms: if z is close to one of the x i , the OPE coefficient C B DA 1 ···A k (z, x 1 , . . . , x k ; y) factorises according to the condition (1.2) and one obtains , this is less singular than d(x i , z) −n and thus integrable over z. Similarly, the terms in the third line of (1.7) are IR subtractions, ensuring the decay of the right-hand side if z is far from the expansion point y (and thus all points). Finally, the terms in the last line arise from the specific choice of renormalisation scheme, and in special cases can be absorbed into redefinitions of the couplings (e.g., in flat space and if there is only one coupling as in λφ 4 theory [12,13]). I derive formula (1.7) in the algebraic approch to quantum field theory on curved space(-time)s [42][43][44], which has the advantage of cleanly separating the algebraic issues in quantum field theory (including renormalisation and possible renormalisation amibiguities) from the choice of states. Section 2 gives a short overview of the algebraic approach, including the exact definition of composite operators in curved space that is used and the definition of (physically acceptable) Hadamard quantum states. In section 3, I give the definition of the OPE coefficients, derive formula (1.7), and show that with this definition the OPE exists in the asymptotic sense (1.1) and fulfills the factorisation and associativity conditions (1.2) in any Hadamard state. Finally, in section 4 I use the formula (1.7) to compute the OPE of φ with itself in hyperbolic space (Euclidean anti-de Sitter space) to first order in g for a gφ 4 interaction, and conclude in section 5. As usual, = c = 1.

Quantum field theory in curved space from an algebraic point of view
Consider a complete smooth Riemannian manifold M of dimension dim M = n ≥ 2, with smooth metric g µν and covariant derivative ∇. An important concept in Riemannian geometry is the exponential map exp x around a point x, which determines points y in a neighbourhood of x as endpoints of geodesics passing through x. That is, given the geodesic from x to y and the tangent vector ξ = ξ(x, y) along this geodesic at the point x, we have exp x (ξ) = y. This is well-defined in a sufficiently small neighbourhood of x, called normal geodesic neighbourhood, where a unique geodesic between x and y exists, and the length of the tangent vector |ξ| ≡ g µν (x)ξ µ ξ ν is equal to the geodesic length d(x, y). In the following, we work exclusively in such a normal geodesic neighbourhood, which is not a restriction since the OPE is anyway a local expansion. Nevertheless, if M has non-positive curvature (e.g., flat or hyperbolic space), the exponential map is even globally well-defined and the normal geodesic neighbourhood of any point is equal to the whole M [45]. A useful coordinate system are Riemann normal coordinates (RNC), where geodesics are straight lines and consequently y = exp x (ξ) = x + ξ. Since RNC are uniquely defined, partial derivatives in RNC can be expressed using covariant derivatives, for which we use the short-hand notations x and symmetrised over all indices is equal to ∇ k T , symmetrised over all indices, which is shown as follows: In particular, since ξ is an arbitrary vector, it follows for τ = 0 (thus z = x) that Consider now in RNC the totally symmetrised expression for some tensor T . Replacing a partial by a covariant derivative, this is equal to Distributing the partial derivatives in the second term, one obtains a sum of products of totally symmetrised derivatives acting on the Christoffel symbol, and derivatives acting on T . Since however the totally symmetrised derivatives acting on the Christoffel symbol vanish at x (2.2), all these terms vanish. Inductively one thus obtains In particular, for a scalar φ one has ∂ k ξ φ = ∇ k φ.

The free theory
For simplicity, we restrict to the case of a real scalar field φ. The free theory depends on the equation-of-motion operator P ≡ − + m 2 + ξR, where m is the mass of the scalar field and ξ denotes a non-minimal coupling to the scalar curvature R. We assume that there exists a Green's function G(x, y) y) is the Hadamard parametrix that is uniquely determined by the local geometry and singular as x → y, while W (x, y) is a smooth function depending on the boundary conditions. Analogously to the Lorentzian case [46,47], the parametrix has the expansion where c n is a normalisation constant depending on the dimension n, σ(x, y) = 1 2 d(x, y) 2 is Synge's world function equal to one half of the geodesic distance squared between the points x and y, and U (x, y) and V (x, y) are smooth symmetric functions with a computable expansion in powers of σ, see, e.g. [48][49][50] and references therein. In even dimensions, the expansion of U is finite and V has an infinite expansion (convergent in analytic spaces), while for odd dimensions V = 0 and U has an infinite expansion (again convergent in analytic spaces). The expansion coefficients fulfill recursion relations that are obtained by imposing the equation of motion for G(x, y) and comparing powers of σ (expanding also W (x, y)), but their concrete form is irrelevant for the purposes of this paper, where it is only important that they are completely determined by the local geometry, concretely by the metric and curvature tensors along the geodesic from x to y. On the other hand, W (x, y) determines the state-dependence of the correlation functions and depends on global boundary conditions.
In the Lorentzian case, the parametrix H is constrained by the fact that one needs to obtain the correct commutator [φ(x), φ(y)] = i[G adv (x, y) − G ret (x, y)], where G adv/ret are the advanced and retarded propagators for P . This is ensured by a proper i prescription in the parametrix (2.6), which leads to retarded, advanced, Feynman, etc. parametrices. Furthermore, W (x, y) is constrained by the requirement of positivity, i.e., G(x, y) must be a distribution whose Fourier transform is positive, which is necessary for the probability interpretation. Contrary to the Lorentzian case, in the Euclidean case the fields φ commute everywhere, and the positivity requirement is replaced by reflection positivity, which ensures that correlation functions can be Wick-rotated back to Lorentzian signature and can then be interpreted as the correlation functions of a unitary theory [51,52]. However, for a general Riemannian manifold M there is no corresponding Lorentzian manifold M obtained by any kind of Wick rotation, such that in the following I do not impose any condition on W (x, y) (other than it being smooth).
In the algebraic approach [53][54][55], one first constructs the algebra of non-interacting fields A 0 (M, g), which in the Riemannian case is just the commutative algebra generated by the φ's and the identity 1. Non-interacting quantum states · ω are then linear functionals on A 0 , and if ω is a free state with two-point function G one has the usual Wick theorem As is well known, one has to perform normal ordering to define powers of φ because of the singularities of G(x, y) as y → x. There are essentially two ways: one is to perform normal ordering with respect to the full two-point function, which gives and :φ(x 1 ) · · · φ(x k ): G ω = 0. However, since the state-dependent part W (x, y) of G depends on global boundary conditions, the so-defined normal-ordered products do not transform locally and covariantly under diffeomorphisms. That is, under a change of coordinates that leaves x 1 , . . . , x k (and connecting geodesics) invariant, the normal-ordered product :φ(x 1 ) · · · φ(x k ): G might nevertheless change because the value of W at the points x i also depends on the geometry at other far-away points, and in particular on boundary conditions.
To remedy this situation, one has to perform normal ordering with respect to the Hadamard parametrix only [55,56]: The Wick theorem for the Hadamard-normal-ordered products then reads and we see that the expectation values of normal-ordered products are smooth functions, such that the limit of coinciding points x i → x can be taken. Since by construction the Hadamard parametrix only depends on the local geometry, it follows that :φ(x 1 ) · · · φ(x k ): H does not change under diffeomorphisms that leave x 1 , . . . , x k and the connecting geodesics invariant. Composite operators that transform covariantly are then obtained by taking the limit x i → x, see the next subsection. For the product of two normal-ordered quantities, one obtains the formula where the variational derivatives formally act to the left and right like on classical fields, as indicated by the arrows. This formula follows directly from the definition (2.9), and an analogous formula with all H's replaced by G's holds for products normal ordered with respect to the full two-point function.
Lastly, we note that one can also consider more general (non-free) Hadamard states, which are linear functionals on A 0 whose connected two-point function has the Hadamard form H(x, y) + W (x, y) as above, and whose connected n-point functions for n = 2 are smooth functions W (x 1 , . . . , x n ). 1 This is especially important because unlike in the Lorentzian case where there are plenty of free Hadamard states [59,60], on a large class of Riemannian manifolds there is a unique fundamental solution for the Green's function [61][62][63], and thus a unique free ("ground" or "vacuum") state. However, from a given Hadamard state one can easily construct a new one by adding a spectator field and defining for any smooth function f with φ(f ) = φ(x)f (x) dx and φ(f * )φ(f ) ω = 0, which gives another (non-free) Hadamard state.

Composite operators
A composite operator is given by products of derivatives of φ, and I use an abstract multiindex notation, where a label A = (a 1 , . . . , a n ) with a i ∈ N 0 is associated to the operator such that ∂ A O A = 1. Two operators are identified if their labels only differ by a permutation, which corresponds to the fact that the fields φ commute; the number of possible permutations for a given A is given by where the multi-index C = A · B is given by C = (a 1 , . . . , a |A| , b 1 , . . . , b |B| ), or any permutation thereof. In particular, when in later formulas one sums over multi-indices A and B with the restriction that A · B = C for some given multi-index C, this sum extends over multi-indices A and B corresponding to distinct composite operators O A and O B (i.e., counting A and A that differ by a permutation only once), and such that A · B is some permutation of the given multi-index C. In other words, the sum extends over composite : H are defined in the obvious way by normal-ordering ∇ a 1 φ(x) · · · ∇ an φ(x). In the following, I concentrate on the Hadamard-normal-ordered quantities, but analogous formulas (with H replaced everywhere by G) hold for normal ordering with respect to the full two-point function. Since as explained above in normal-ordered expressions one can take the limit of coinciding points, and since symmetrised covariant derivatives are equal to partial derivatives in RNC, we have and similar expressions for the product of more than one composite operator. Since expectation values of normal-ordered quantities (2.10) are smooth functions of the points, it follows that also the expectation value of composite operators : H ω is a smooth function of the points x i . For the OPE, we need to extract the term with a given composite operator O A from a sum of normal-ordered expressions, which I do as in [64][65][66][67]  holds, which the exact definition of δ B A given later in equation (2.21). D B x is defined to act on a normal-ordered product of m fields φ by where the ξ i are the tangent vectors along the geodesic from x to x i , and where S denotes the idempotent symmetrisation with respect to all possible permutations of the b i . This ensures that the definition does not depend on the choice of an ordering among the b i nor among the x i (since the normal-ordered product is symmetric under an exchange of its arguments). This definition is extended to all normal-ordered products by linearity and by postulating that D commutes with taking derivatives in RNC and taking limits. Taking derivatives one obtains and taking the limit This corresponds to a formal Taylor expansion of the φ(x i ) around x, picking the coefficient which has the correct number of derivatives to correspond to O B . As one checks easily, the inclusion of the combinatorial factor P B in the definition without any combinatorial factors when summing only over the multi-indices B corresponding to distinct operators O B .

Renormalised products
In the Riemannian setting, (free) renormalised products (called E-products in [68]) are the analogue of the (free) time-ordered products in the Lorentzian case. In the algebraic approach, they are constructed recursively in the number of arguments, imposing certain conditions that are detailed in the following [48, 53-56, 68, 69]. The renormalised product of a single argument is just the Hadamard-normal-ordered product which was defined in the previous subsections. The recursive construction then proceeds by imposing the conditions of symmetry, factorisation, scaling, field independence, and local covariance, together with some other technical conditions which are not needed for the present paper. Since all fields commute, the symmetry condition is just the obvious one where the dots stand for an arbitrary number of other operators. The Riemannian analogue of the well-known causal factorisation of time-ordered products in the Lorentzian case (for two arguments, and similar expressions in the general case) is given by if all the points x i are distinct from all the y j . That is, since in the Riemannian case all points are space-like related, we have factorisation in either order if not all points coincide.
Field independence is the condition and is clearly satisfied for the renormalised products of a single argument (2.23). Local covariance requires that the renormalised products transform locally and covariantly under diffeomorphisms, which is the reason that one needs to use the Hadamard-normal-ordered product in (2.23). Field independence and local covariance ensure the existence of a local Wick expansion (shown in appendix A) where the r are numerical (locally covariant) distributions which are smooth functions for non-coinciding points. Using the local Wick expansion for each of the renormalised products in the factorisation condition (2.26), and formula (2.11) for the product of two Hadamardnormal-ordered products, the distributions r are already fully determined except for points which all coincide. The extension of r to those points corresponds to renormalisation, and the renormalisation freedom is given by local terms (δ distributions and their derivatives, which correspond to polynomials in momenta in flat space). Furthermore, the condition of local covariance severely restricts the coefficients of these local terms, which must be polynomials of curvature tensors and their covariant derivatives of the correct dimension, together with powers of mass and other parameters of the theory [48,55,56]. Lastly, the scaling condition is a bound on the degree of divergence as the points are scaled together, and can be given both for the renormalised products R and the numerical distributions r, and the latter will be important for this paper. The scaling of points is defined using RNC: writing y = exp x (ξ) with the tangent vector ξ along the geodesic from x to y, the rescaled point is given by y τ ≡ exp x (τ ξ) for 0 ≤ τ ≤ 1. The scaling condition then reads where I write In the Lorentzian case, one then includes the interaction L by the Bogoliubov formula [70] and defines the interacting time-ordered products T L as Since the time-ordered products factorise only if its arguments are causally separated, this gives a non-trivial result. For example, to first order one obtains for the interacting observables/composite operators the expression 32) and the second term is non-vanishing for (the part of) L in the past light cone of x.
In the Riemannian case, this does not work for the simple reason that the renormalised products factorise everywhere (except for coinciding points). Instead, one has to define expectation values in the interacting theory (the Schwinger functions) by the Gell-Mann-Low formula [71] for some state ω of the non-interacting theory. 2

The operator product expansion
One would like to define the OPE coefficients C B A 1 ···A k (x 1 , . . . , x k ; y) such that the OPE holds for all Schwinger functions (2.33) as an asymptotic equality: In perturbation theory, multiplying both sides with R[exp(−L)] ω this is equivalent to 2) which is easier to work with. Instead of performing the OPE for all operators appearing in the expectation value ] ω , one can also perform an expansion of the first m operators, and then an expansion of the remaining ones. Demanding that this agrees with the expansion for all operators at once, one obtains the factorisation conditioñ By demanding that two different factorisations agree, the associativity condition follows automatically: and other associativity conditions are obtained by permuting the operators O A i before performing the OPE, using the symmetry condition (2.24) for the renormalised products. In subsections 3.1 to 3.3 I will first define the OPE coefficients and derive a recursion formula for them, which then allows a quite short proof of the OPE (3.2), as well as the factorisation (3.3) and associativity (3.4) conditions for the coefficients as asymptotic equalities.

Definition of OPE coefficients
Assuming that an OPE of the form (3.2) exists, one can obtain the OPE coefficients by extracting the coefficient of the corresponding operator. That is, one would like to define where D B L,y is a linear functional fulfilling with δ B A defined by (2.21). This "projection" of the renormalized product onto a given operator was first proposed by Gorishny, Larin and Tkachov [64,65], and independently by Bostelmann [66] (in the algebraic approach to QFT); the generalisation to curved spacetimes is due to Hollands [67]. In the free theory where from which an expression for the difference D B L,x −D B x can be obtained to arbitrary order in perturbation theory. Namely, consider the Neumann series formed by the right-hand side interpreted as multiplication operators from the set of composite operators or multi-indices to itself. The formal sum Since the sum over C, and also the sums implicit in the definition (3.9) ofD B x C are infinite, this is only a formal solution. Nevertheless, I will now proceed and see later on how to turn this into a well-defined solution. Using which is a recursive definition for the OPE coefficients in perturbation theory since D B y E is at least of first order. By restricting the sum over E to operators with [O E ] < ∆ for some ∆ > 0, this becomes well-defined in perturbation theory. That is, the OPE coefficients with cutoff ∆ and expansion point y are defined as a formal power series in perturbation theory by the recursion (3.14) One checks that the same recursion formula can be obtained from the definition of Hollands [67].
To show existence of the OPE in any form it is obviously important that no artificial cutoffs are in place. In the free theory, the definition 3.14 of the OPE coefficients does not depend on ∆, but including interactions one has to perform a redefinition of composite operators to be able to take the limit ∆ → ∞. In general, such a redefinition is of the form where the Z A B (x) are in general functions of x. In order for this redefinition to correspond to an allowed mixing of composite operators, Z must be lower triangular, i.e.
, which ensures that a composite operator only mixes with operators of equal or lower dimension such that , and that the sum over B is finite. For the purposes of removing the cutoff ∆, the redefinition will be at least of first order in L and depend on ∆, such that Under such a redefinition, the OPE coefficients for O B , where the inverse is given as a formal power series in perturbation theory by the geometric series (3.17) and is again lower triangular. The OPE coefficients should also transform covariantly, which means that (in the sense of formal perturbation theory) (3.18) However, since there is a cutoff ∆ on operator dimensions it is only necessary to demand this equality for operators of low dimension; concretely I set The mixing matrix Z ∆ should be determined such that the recursion formula for the redefined coefficients reads since then the limit ∆ → ∞ can be taken without problems. Combining all of the above, Z ∆ needs to be determined such that By "cancelling" the OPE coefficients on both sides, one thus obtains an explicit recursion formula ) and the mixing matrix is lower triangular, the very last term not contribute, and it follows that for each ∆ > 0, there exists an admissible redefinition of composite operatorsÕ A , given by equation (3.23), such that the redefined OPE coefficients C B,∆ A 1 ···A k satisfy the recursion (3.20) without the very last term. Furthermore, from the explicit formula (3.23) for the mixing matrix Z ∆ , it follows that the redefinition is stable in the sense that for all composite operators (3.20), and the limit ∆ → ∞ can be taken. The resulting OPE coefficientsC B A 1 ···A k (x 1 , . . . , x k ; y) are then obtained recursively as (3.24)

Variational formula for the original OPE coefficients
I start with the simpler case of the original definition (3.14) with cutoff ∆. Since the renormalised products are defined in the free theory, the renormalisation procedure is in particular independent of any coupling constant. Since also the definition of composite operators (2.13) and the Hadamard-normal ordering (2.17), (2.9) are independent of the coupling, one simply has This formula is basically the naive one that one expects from a path integral definition of correlation functions, with the small (but important) difference that both sides are fully renormalised. Since also the free functional D B y (2.19) is independent of the coupling, it is straightforward to take a derivative of the recursion (3.14) and use equation (3.25) to obtain where I used equation (3.14) to express D B y acting on a renormalised product again by OPE coefficients. Renaming C ↔ D in the last line, it is seen that the three last lines are of higher order in perturbation theory, and that the unique solution of this equation reads That is, the variational derivative of an OPE coefficient with respect to a coupling constant is given by the (naively expected) coefficient with an additional insertion of the corresponding interaction operator O I , and a correction term that involves a product of two OPE coefficients. Since a variational derivative was taken, to obtain the correction to the OPE coefficient itself one has to integrate this equation over z; the correction term ensures (at least formally) that this integral is well defined in the IR. Namely, if the point z is far away from the x i , the OPE coefficient C B,∆ A 1 ···A k I (x 1 , . . . , x k , z; y) should factorise exactly in the way given by the sum, such that the whole expression vanishes at least in the formal limit ∆ → ∞. On the other hand, since the OPE coefficients are by definition well-defined (renormalised) distributions, there is no UV problem when z is close to any of the x i . In the next subsection, we will see how this works when the limit ∆ → ∞ can actually be taken.

Variational formula for the redefined OPE coefficients
I now assume that the interaction L has the generic form Since the inverse mixing matrix is lower triangular, if there is a finite number of couplings g A before the redefinition, there is also have a finite number of redefined couplingsg A after it. Consider the quantity It is a long but straightforward computation to check that is a solution of this recursive equation (uniquely defined in perturbation theory), using the recursive definition of the OPE coefficients (3.14) in the form (3.34) Since the limit ∆ → ∞ exists for the redefined OPE coefficients, the mixing matrix and Γ ∆ EA B (from the definition (3.30)), and the sum over E is finite if the number of couplings g E is finite, it is possible to take the limit ∆ → ∞ of this last equation.
In this limit, I now compute and thus for the redefined OPE coefficients (3.24) and using the result (3.34) (in the limit ∆ → ∞) it follows that (3.36) The unique solution (in perturbation theory) is obviously given by the first three lines, but one still has to express the variational derivative using the new couplingsg E (3.29). Using the definition of Γ EA C (3.30) and the sum (3.34) (in the limit ∆ → ∞), I compute and thus This equation agrees with the one derived in flat space [12,13,41] except for the addition of the last term. However, this term is of higher order in couplings, and expanding equation (3.38) in powers of the couplings the right-hand side is of lower order, and the OPE coefficients can be constructed recursively. Expressing also the variational derivatives of renormalised products (3.35) in terms of the redefined coupling, using the relation (3.34) (in the limit ∆ → ∞) I also obtain which is now a recursive equation that can again be solved order by order in couplings. In principle, one would like to perform also a redefinition of the couplingsg A →ĝ A to remove the last term in equations (3.38) and (3.39), and simply have instead of equation (3.37). However, this is impossible in general, as can be seen as follows: for such a redefinition to exist, one must satisfy the integrability conditions which using equations (3.30) and (3.34) (in the limit ∆ → ∞) reduces to The result (3.42) then tells us that the commutator of two covariant variational derivatives is given by else (3.45) is the torsion of the covariant variational derivative. A change of couplings is a change of coordinates (a diffeomorphism) in this manifold, but since torsion is a tensor it cannot vanish in any coordinate system if it doesn't vanish in a given one. Note that the problem doesn't arise from the fact that the couplings are position-dependent: even if one takes the adiabatic limit where the couplings become constant, and consequently has to integrate the torsion (3.45) over x and y, it can still be non-vanishing if O A = O B . Accordingly, if there is only one coupling and the adiabatic limit is taken where this coupling becomes constant, an additional redefinition of this coupling exists such that the last line of formula (3.38) is absent.

Existence, factorisation and associativity of the OPE
With the recursive formula (3.38) for the OPE coefficients, it is now easy to show that the OPE is an asymptotic expansion, and that the coefficients factorise and fulfill the associativity condition. Consider thus the remainder of the OPE (3.2) where the index ∆ > 0 determines how many terms of the OPE are subtracted, and the remainder of the factorisation (3.3) where ∆ > 0 has the same meaning as before, while the index m denotes that factorisation occurs between the first m and the last k − m operators. Using the variational formula (3.38), equation (3.39) and that the (non-interacting) expectation value · ω is independent of the couplings, a long but straightforward computation yields . . . , x k ; y, z) .

(3.49)
Together with the variational formula (3.38) for the OPE coefficients, with these equations it is now easy to show inductively that as the insertion points x i scale towards the expansion point y, the OPE coefficients scale in the appropriate way, that the OPE is an asymptotic expansion, and that the OPE coefficients factorise properly. The scaling of the points is defined in RNC as for the renormalised products: given x = exp y (ξ) with ξ the tangent vector at the geodesic from x to y, the rescaled point reads x τ ≡ exp y (τ ξ), such that x 1 = x and x τ → y as τ → 0. Furthermore, instead of the fixed expansion points y and z one also needs rescaled expansion pointsỹ τ = exp y (τξ) andz τ = exp z (τζ) for some tangent vectors ξ andζ (which may be zero, in which case we have again fixed expansion points). At order in perturbation theory, I then want to show inductively that C B,( ) For (super-)renormalisable interactions, both the scaling of the OPE coefficients themselves and the speed at which the OPE remainder vanishes (or the factorisation holds) are thus independent of the order of perturbation theory, while for non-renormalisable interactions the OPE coefficients become more singular at each order, such that one needs to include more terms in the OPE for a given accuracy. 3 To show the scalings (3.50), I start with the free theory, where the OPE can be thought of as a generalised Taylor expansion (with D B as the derivative operator), as will be clear from the following. Using Wick's theorem (2.10), one computes for even n where the sum is over all (unordered) pairings p of {1, . . . , n}, and :φ(x 1 ) · · · φ(x n ): H ω = 0 for odd n. Since W (x, y) is a smooth function of the points x and y, one can expand this in a Taylor series with remainder around a point y, assuming all points x i = exp y (ξ i ) (with the tangent vectors ξ i along the geodesic from x i to y) lie in a common normal geodesic neighbourhood. The Taylor expansion reads where ξ = (ξ 1 , . . . , ξ n ), k = (k 1 , . . . , k n ) is a standard multi-index of length n, and where I used that derivatives of fields in composite operators are exactly defined using derivatives in Riemann normal coordinates (2.17). On the other hand, we have (2.19) and thus Combining this with the fact that since O k does not depend on the order of the k i in the multi-index k, and that the sum over all multi-indices k can then be replaced by a sum over all operators B by including an additional factor P B that counts how many permutations of indices give the same operator, it follows that (3.56) The analoguous expression for composite operators is obtained by acting with derivatives (using that D B commutes with derivatives), and setting some of the x i equal according to equation (2.17). In this way, one obtains (3.57) where the remainder has a lower order than before because derivatives have been taken. Using now the local Wick expansion (2.28) of the renormalised products, it follows that (3.58) Since the sums over the B i are finite and D C is linear and only acts on normal-ordered products, the sums over C and the B i can be interchanged and the local Wick expansion used in reverse, such that the first term reads , and the scaling (3.50b) for the remainder of the OPE holds in the free theory (also for a rescaled expansion pointỹ τ , since the difference between y andỹ τ is of higher order in τ ).
Note that this result is really a statement about the smoothness of the state-dependent part W (x, y) of the two-point function. In fact, the same proof applies to any linear functional that is smooth and commutes with derivatives, and is not specific to the expectation value · ω . That is, given any linear functional F such that F (:O A 1 (x 1 ) · · · O An (x n ): H ) is smooth in the positions x i and commutes with derivatives with respect to the x i , we have by the above argument that which is the scaling (3.50a) of the OPE coefficients. On the other hand, taking (3.62) one computes under the assumption that Assume now that inductively the scaling (3.50a) has been shown for all orders up to order , such that it can be used for the terms on the right-hand side. Since the couplings is of perturbative order − − 1 because of the additional functional derivative with respect to the coupling, but has an additional insertion of [O C ] according to equation (3.38). Finally, in RNC centered at y and with n = dim M we have [45] dz = det g(z) d n exp y (ξ) = d n ξ + O ξ n+1 (3.66) and thus dz τ = O(τ n ), such that with the estimate for the left-hand side follows. In exactly the same way, the other scalings (3.50b) and (3.50c) follow from equations ( in Poincaré coordinates, which (in contrast to Lorentzian AdS) cover the whole space. Hyperbolic space is maximally symmetric, and the geodesic distance d(x, y) between two points x and y can be expressed as in terms of the chordal distance In a maximally symmetric space, all tensors depending on two or more points and transforming covariantly under the isometries of that space can be expressed using the geodesic distance between the given points, the tangent vectors along the geodesic connecting them and the tensor of parallel transport along the geodesic [78,79]. In turn, these can be expressed using the chordal distance and derivatives thereof. In particular, the tangent vector ξ(x, y) µ at x along the geodesic from x to y is given by a multiple of ∇ µ x u(x, y), which is determined by imposing the normalisation ξ(x, y) µ ξ(x, y) ν g µν (x) = d(x, y) 2 . Using (which is easily obtained from the above expressions), the tangent vector thus reads which has the right flat- Similarly, the tensor of parallel transport g µν (x, y) can be defined by g µν (x, y)ξ(y, x) ν = −ξ(x, y) µ and the coincidence limit lim y→x g µν (x, y) = g µν (x) [79]. From the second condition, it follows that g µν (x, y)g µν (x, y) = d + 1, and one computes using also that and that covariant derivatives with respect to different points commute. In the flat-space limit we have g µν (x, y) → − 1 2 ∂ x µ ∂ y ν d(x, y) 2 = η µν for d(x, y) → (x − y) 2 as required. The free propagator of a scalar of mass m satisfies the equation of motion Conformal coupling is obtained for ∆ = (d + 1)/2, and since the maximally symmetric solution for the propagator only depends on the chordal distance u(x, y), outside of coincidence one obtains the equation A one-parameter family of Hadamard states is given by where the correctly normalised Hadamard parametrix (singular as u → 0) is given by (the logarithmic term that appears for even dimensions is absent for conformal coupling), while the state-dependent regular term reads for α ∈ R. The usual choice of state is given by requiring fast fall-off towards the conformal boundary as u → ∞ [79,80], which is obtained for (4.12) For concreteness, in the following I work in d + 1 = 5 dimensions, and have thus consistent with the dimension of the scalar [φ] = 3 2 .

Two-point OPE coefficients in the free theory
The free OPE coefficients C B,(0) A 1 ···A k can be obtained directly from the defining equation (3.14) and one first needs to determine the renormalised products. Again for simplicity, I choose the expansion point to be y = x k , which gives short expressions for the OPE coefficients. For one argument, the renormalised product is just the normal-ordered expression (2.23) For more arguments, one has to proceed inductively both in the number of arguments and the dimension of the composite operators. For non-coinciding points x = y, by the factorisation condition (2.26) one obtains and since [φ]+[φ] = 3 < 5 = dim M , the renormalisation (extension to coinciding points x = y) is unique. Using the formula (2.11) for the product of two normal-ordered expressions, one obtains φ(y)φ(x) = :φ(y)φ(x): H + H(y, x)1 (4.17) and the local Wick expansion Admittedly, this is quite trivial, but serves to illustrate all steps (which I will abbreviate for the second, more complicated example).
Using the definition of D B x (2.19), I then compute D 1 and all other C B φφ (y, x; x) vanish, where I recall the notation ∇ µ 1 ···µ k ≡ ∇ (µ 1 · · · ∇ µ k ) . Instead of explicitly using D B x , one can also perform a Taylor expansion of the Hadamard normalordered product to obtain and comparing with the OPE the same OPE coefficients are obtained. By the generalised Leibniz rule, we have ...,i ) , (4.25) where the explicit symmetrisation of covariant derivatives is unnecessary because of the factor ξ ⊗k in (4.24) which is already symmetric. By definition, the factor P (i 1 ,...,i ) (2.16) counts how many permutations of the i j give the same composite operator, such that and the non-vanishing OPE coefficients are given by 1φ (y, x; x) = 1 . I furthermore need the coefficients of the OPE of two composite operators which include each more than one basic field φ. In this case, the renormalised products must be non-trivially extended to coinciding points, which I recall corresponds to renormalisation. Concretely, I look at the OPE of φ 4 with φ, φ 2 and φ∇ ρ φ, and thus first determine the corresponding renormalised products. For distinct points x = y, we use the factorisation condition for the renormalised products and the formula (2.11) for the product of two normal-ordered expressions, resulting in Comparing with the local Wick expansion of the renormalised products 29) and imposing that the r commute with covariant derivatives, we obtain in addition to (4.19) r φ (x) = 0 for = 2, 3, 4 , (4.30a) r φ (y) ⊗ φ(x) = 0 for = 2, 3, 4 , (4.30b) r φ (y) ⊗ (φ∇ k φ)(x) = 0 for = 1, 3, 4 , The extension of the second term to coinciding points x = y is unique whenever the degree of divergence is not integer, which is the case for even . For = 3 the degree of divergence of r φ 3 (y) ⊗ φ(x) is given by 4[φ] = 6 > 5, which is integer and greater than the space dimension. Therefore, a renormalisation ambiguity exists, which is a local term (δ and its derivatives), multiplied by the metric or curvature tensors of the correct dimension. However, the only possible term in this case would be ∇ µ δ(x, y), which does not have the correct tensor structure. Also for the third and last term, the degree of divergence is not integer and there is a unique (vanishing) extension to coinciding points. For the fourth term, we first have to define some extension, since the degree of divergence is 4[φ]+k = 6+k > 5. Using the explicit form this can be easily done by computing the naive product for x = y, and then extracting derivatives until the remainder has degree of divergence < 5. Using the relations (4.4), (4.7) and we obtain for x = y (4.33d) Since u(x, y) −2 has degree of divergence 4 < 5, it has a unique extension to coinciding points x = y as a well-defined distribution. Since derivatives acting on a distribution are also always defined, we have found an extension (renormalisation) of all the above expressions. We now have to determine renormalisation ambiguities. For H(y, x) 2 , which has degree of divergence 6, the only possible ambiguity is given by ∇ µ δ(x, y), which however has not the correct tensor structure, such that the renormalisation of H(y, x) 2 is unique.
At this point, it was important that we consider locally covariant renormalised products, where the ambiguities are strongly restricted to only contain polynomials of curvature tensors and their covariant derivatives. Otherwise, a possible ambiguity term for H(y, x) 2 would be −1 δ(x, y), but this is non-polynomial in the Ricci scalar. Expanding the normalordered products in (4.28) around y = x, we can then read off the needed OPE coefficients For the OPE of φ 2 with itself, we have already computed the required extensions. With the local Wick expansion (4.29), we obtain (4.36) using the results for the r distributions (4.19) and (4.30). Expanding the Hadamard normalordered products around y = x, we directly read off the OPE coefficients

Three-point OPE coefficients in the free theory
To determine the first-order corrections, we also need the three-point OPE coefficients where we used the result (4.28) for the renormalised product of two composite operators.
Using the local Wick expansion, we obtain in addition to the results (4.19) and (4.30). Since the degree of divergence of the product H(z, y)H(z, x) is 6 < 2·5, it is already well-defined as a distribution in 3 points. Expanding the normal-ordered products around x, we thus obtain directly the OPE coefficients

The OPE at first order
I now would like to compute the OPE of φ(x)φ(y) up to terms vanishing as d(x, y) 2 , where d(x, y) is the geodesic distance between x and y, to first order in the interaction gφ 4 . Since the theory is massive, one would like to take the adiabatic limit (of constant coupling g), but we will see that this is actually not possible. Nevertheless, as long as x and y are close together and within a region where g is constant, we will see that the OPE coefficients do not depend on the concrete form of the cutoff. In 5 dimensions, the field φ has dimension 3 2 and the interaction has therefore dimension [L] = 4[φ] − 5 = 1. Using the result (3) for the scaling behaviour of the remainder, one needs to take into account the terms up to ∆ = 2 + 2[φ] + [L] = 6. To shorten notation, let me denote in this section the interacting composite fields by capital letters, such that in any Hadamard state ω, where only those OPE coefficients are displayed that do not vanish because of the Z 2 symmetry φ → −φ (of both the free theory and the interaction).
Integrating the recursive formula (3.38) over z with a cutoff function f (z), one obtains at first order (4.42) where the order in g is denoted by the index in parentheses. For the OPE of φ with itself, this reduces to (4.43) where all other contributions vanished either because C φ,(0) φ 4 φ = 0 or by the Z 2 symmetry. We see that for these coefficients, only IR subtractions need to be made, and no UV subtractions. For the coefficients in the OPE (4.41), one then obtains concretely For the first non-vanishing correction C φ 2 ,(1) φφ one obtains using (4.40) and (4.13) (4. 45) In principle, one would like to take the adiabatic limit f → 1, but this would result in an IR divergence close to the boundary ζ → 0 since the integrand diverges there ∼ ζ −2 . One possibility would be to consider normal ordering not with respect to the Hadamard parametrix, but the full propagator; call the corresponding OPE coefficientĈ . For f = 1, the integral is then clearly EAdS-invariant and therefore only depends on the chordal distance u(x, y), such that one can set x = y = 0 to simplify the actual evaluation. With the propagator (4.13), I computê using spherical coordinates d 4 z = 2π 2 |z| 3 d|z|. The integral over ζ is done by splitting the integration range into varions regions depending on whether ζ is larger or smaller than z x and z y , and which of z x and z y is larger. This results in with z + = max(z x , z y ) and z − = min(z x , z y ). Expressing those in terms of u = u(x, y) according to (recall that x = y = 0), one finally obtainŝ where it was used that u(x, y) ∼ d(x, y) 2 as x → y. In terms of the geodesic distance d(x, y) itself, related to u by (4.2), we havê (4.50) For a covariant definition, one would expect a covariant expansion whose coefficients are proportional to the curvature, hence quadratic in the hyperbolic radius . However, since normal-ordering was performed with respect to the full propagator (which is not covariant), this does not hold. This is easily seen by performing an expansion of the propagator (4.13) itself in terms of the geodesic distance: where also odd powers of appear. However, since local covariance is important to restrict the renormalisation ambiguity, this is not the most useful solution. I thus return to the original coefficient (4.45), and observe the following: since the propagator G(x, y) fulfills the equation of motion (4.8), the parametrix fulfills for x and y in the region where f = 1. Note that the second contribution to the integral (4.61) is completely irrelevant for this result and thus the singular terms of the OPE coefficient C φ∇ µ φ, (1) φφ . Instead, it ensures that the integral is convergent in the IR if one considers normal-ordering with respect to the full propagator and takes the adiabatic limit f → 1 (without this term, a logarithmic divergence would appear). Making the ansatz and using the explicit form (4.5) of the tangent vector ξ(x, y) µ , I obtain the differential equation While the first term on the right-hand side is singular as u → 0, the second one is analytic at coincidence. This is seen easiest by using the geodesic distance d(x, y), in terms of which the tangent vector is given by (4.5), and computing , y)) . (4.70) Since the right-hand side of equation ( up to terms analytic at coincidence, from equation (4.68) the differential equations (with the same function K(u) (4.65) as before) are obtained. The coefficient C scales as d(x, y) ∼ u(x, y), such that C(u) scales as √ u andC(u) scales as u − 1 2 . The only solution of the equation for C(u) with this bound on the scaling is then analytic at coincidence, such that the singular part of C(u) vanishes. On the other hand, the equation for C(u) cannot be solved anymore in closed form, but with the known scaling an asymptotic solution to any required order is easily found: Inserting the results in equation (4.41), the OPE with φ with itself reads (4.74) or expressing the chordal distance u in terms of the geodesic distance d according to (4.2) (4.75)

An OPE coefficient at second order
In this subsection, I only want to compute the OPE coefficient C Expanding the normal-ordered products around y = x, the required OPE coefficients are obtained as vanishing. Note that H(x, y) 4 etc. are not well-defined as a distribution because they are too singular at coincidence, but using the same trick as before this will actually be unproblematic. Moreover, the coefficient of φ in the OPE of φ∇ k φ with φ for k = 0, 1, 2, 3 is needed, which is obtained from The required three-point coefficients come from the OPE of φ∇ k φ with φ and φ for k = 0, 1, 2, 3, as well as the OPE of φ 4 with φ 4 and φ. For the first, only the coefficient of the unit operator is required, which in the normal-ordered case is just the expectation value (with the propagator G replaced by the parametrix H) such that expanding around z = x = y one obtains Lastly a four-point coefficient is needed, which is again just the expectation value Putting all together and using the OPE coefficients computed above, it follows that To determine the unknown constants c i in this solution, consider the integral formula for this coefficient, which is given by equation (4.87) with the Hadamard parametrix H replaced by the propagator G, and taking the adiabatic limit f → 1 (which does exist for the propagator but not the parametrix). Looking at the IR behaviour u → ∞ of the integral, I note that the IR fall-off of the propagator G as u → ∞ is ∼ u − 5 2 , and the formula has the fifth power of the propagator and two integrations, such that I expect an IR fall-off of the OPE coefficientĈ , i.e., c multiplies a subleading singularity. This is in fact to be expected, since the value of c is scheme-dependent, and can be changed by another redefinition at first order: Φ 2 → Φ 2 +cg1, which results in c → c +c. Only the most singular part of the OPE coefficient C 1,(2) φφ (4.89) is scheme-independent und uniquely determined.

Outlook
While the formula (3.38) gives in principle a straightforward way to compute OPE coefficients to arbitrary order in perturbation theory, its practical application is hampered by the difficulties of evaluating integrals in curved space. Even if one considers normal ordering with respect to the full propagator, where the adiabatic limit of constant couplings can be taken, a direct evaluation of the integals is seldomly possible. This is not surprising, since in that case (for example) the coefficient of the unit operator is nothing else but the expectation value of the fields appearing in the OPE, and so their evaluation cannot be any simpler than the computation of correlation functions. Even in the simplest case of a maximally symmetric space (Euclidean AdS, or hyperbolic space), computations are very difficult (see e.g. the recent works [82][83][84][85][86][87][88][89] for loop calculations in AdS). On the other hand, for a renormalisable interaction the UV counterterms in this formula are sufficient to remove all non-integrable divergences, such that it is possible to evaluate the OPE coefficients numerically. If a further redefinition of composite operators can be determined such that all non-integrable divergences are removed also for non-renormalisable interactions is currently under investigation.
Nevertheless, the most interesting application would be in the context of the AdS/CFT correspondence with all external points on the AdS boundary. In this case, one can use the Symanzik formula [90] for the AdS integrations, which is a huge simplification compared to the general case with external points in the AdS bulk where no such formula exists. Since quite generally a local bulk field theory with small coupling corresponds to a large-N dual CFT with a gap on the boundary [91,92], one thus obtains a way to systematically compute 1/N corrections to the planar limit N → ∞ of such theories. The flat-space analogue of formula (3.38) has already been applied directly to a conformal field theory [93] (see also [94] and [95]), and allows one to compute corrections in conformal perturbation theory where an existing CFT is perturbed by a strictly marginal operator. Using conformal symmetry the integrals can be performed explicitly, and one obtains a coupled system of ODEs describing the change of the conformal data with the coupling, which can in principle even be solved numerically to obtain non-perturbative results. However, the obvious application where the unperturbed CFT is free (Gaussian) is hampered by the fact that the spectrum of such a CFT is highly degenerate, which violates a non-degeneracy condition that is necessary for the derivation of the ODEs [93]. It is quite probable that this non-degeneracy condition can be more easily fulfilled in the planar limit N → ∞, which is already a non-trivial interacting theory, and that a formula for 1/N corrections could be solved numerically. One issue that needs to be adressed to fulfill this program is that even if the external points are on the AdS boundary, the integration that one needs to obtain corrections to the OPE coefficients from (3.38) still involves lower-order OPE coefficients with points in the bulk. Since the form of the bulk OPE coefficients is quite involved, it seems difficult to directly obtain a formula for the conformal data. I believe that a further redefinition of composite operators is needed to change the cutoff function in this integral in such a way as to involve only boundary points, similarly to the redefinition needed in the flat-space formula of Holland and Hollands [12] to pass from the massive to the massless theory. This is currently under investigation, as well as the extension of formula (3.38) to gauge theories and fermions.
O A 1 = 1. Assume that it has been proven for all k i=1 [O A i ] < D, and consider the difference ∆(x 1 , . . . , where all the appearing r's are known by induction since k i=1 [∂ B i O A i ] < D because the sum is restricted to the case where at least one of the O B i is not the unit operator. By field independence of the renormalised products, one obtains δ δφ(y) ∆(x 1 , . . . , where the restriction on the sum over the B j could be removed because O B i = 1 (otherwise the functional derivative vanishes). Since the renormalised product that appears on the right-hand side has a lower total dimension by induction we know that it has a local Wick expansion. Using that it follows that and ∆ i (x 1 , . . . , x k , y) In the second sum, instead of the sum over all B k , the distribution r depending on ∂ B i O A i and then considering ∂/∂ ∇ j φ O B i (x i ) in the Hadamard-normal ordered product, one can also sum over all B k , the distribution r depending on ∂ B i ∂ ∂∇ j φ O A i and consider O B i (x i ) in the Hadamard-normal ordered product. Therefore, all terms cancel and it follows that δ δφ(y) ∆(x 1 , . . . , x k ) = 0 . (A.8) Hence, ∆ must be proportional to the unit operator, and one can simply define the distribution r to be the proportionality coefficient: