The four-point correlation function of the energy-momentum tensor in the free conformal field theory of a scalar field

We present an explicit momentum space computation of the four-point function of the energy-momentum tensor in 4 spacetime dimensions for the free and conformally invariant theory of a scalar field. The result is obtained by explicit evaluation of the Feynman diagrams by tensor reduction. We work by embedding the scalar field theory in a gravitational background consistently with conformal invariance in order to derive all the terms the correlator consists of and all the Ward identities implied by the requirements of general covariance and anomalous Weyl symmetry. We test all these identities numerically in several kinematic configurations. Mathematica notebooks detailing the step-by-step computation are made publicly available through a GitHub repository. To the best of our knowledge, this is the first explicit result for the four-point correlation function of the energy-momentum tensor in a conformal and non supersymmetric field theory which is readily numerically evaluable in any kinematic configuration.


Introduction
The energy-momentum tensor is the most universal operator for conformal field theories (CFTs). Therefore, its correlation functions are natural objects to study in any CFT.
It has been known for a long time that conformal invariance completely fixes the two-point correlation function of the energy-momentum tensor, whereas the three-point function in 4d is fixed up to three scalar coefficients, which can be inferred by computing the correlator for the field theories of a scalar, a fermion and an abelian gauge field. This was first dealt with in coordinate space for three-point correlators of operators up to spin 2 in [1,2], by solving all the constraints implied by the full conformal group. An analogous program in momentum space, based on the solution of all the conformal group Ward identities to constrain the 3 point functions of scalars, vector currents and energy-momentum tensors modulo a few constants was carried out much more recently [3][4][5] and found to be perfectly consistent with the other approach, although very mathematically very different. A different approach in momentum space to the three-point function of the energy-momentum tensor has been pursued in [6], where an explicit diagrammatic computation was performed (and published in a partially onshell kinematics). The consistency of the two momentum space approaches was later verified by [7].
Much more involved is the case for four-point functions, because the conformal group does not uniquely constrain them. For the case of scalars, Polyakov first showed that the four-point function is fixed up to a function of the coordinates cross ratios [8], (1.1) As for operators with spin, very little is known about four-point correlation functions. One very powerful and popular approach to CFT dynamics is the so called conformal bootstrap [9][10][11], which allows to constrain four-point functions through consistency constraints connecting different conformal blocks representations of the same four-point function. So far, it has been successfully applied mainly to the case of scalar operators. One exception is a very recently proposed approach to conformal blocks for any spin [12,13], which was applied to two, three and four-point fucntions [14][15][16], though not to to the four-point function of the energy-momentum tensor due to its sheer complexity. The number of unrestricted degrees of freedom for the four-point function of the energy-momentum tensor in general dimensions was already determined in [17], where the number of effectively independent constraints provided by the Ward identities stemming from the requirements of both special conformal invariance and energy-momentum conservation was nailed down. There have been works -all of them in the context of the AdS/CFT correspondence -which have investigated the four-point function of the energy-momentum tensor. The first one has explored the structure of some contributions to the four-point function through an OPE analysis [18] . Later, the most general structure in coordinate space of the correlator was pinned down in [19,20]. By different means, analogous structures were conjectured in [21] and later proven to be correct in [22]. 1 So far, however, no explicit result has been presented in momentum space though, except for the case of superconformal N = 4 [23], let alone any which lends itself to numerical evaluation.
With these recent advances in mind, in the present paper we present an explicit perturbative computation of the four-point function of the energy-momentum tensor in the simplest possible conformal field theory in 4d which admits a Lagrangian realization and, thus, a perturbative approach, a free scalar field. It is paramount to notice that the fact that the theory is free does not just make calculations simpler: it is such that the 1 loop results for our four-point function (or lower and higher points as well), which can be computed using perturbation theory, are also exact, because no higher order loops are admitted. In this way, a perturbative calculation produces the same result -just perhaps less compact -that would be produced by non perturbative methods exploiting only the symmetries of the theory.
We provide the results in a set of ancillary Mathematica files stored in the public GitHub repository mentioned in the abstract. In these files, the reader is guided step by step through the computation of the two, three and four-point correlators of the energy-momentum tensor. for the theory at hand, as well as through the tests of the Ward identities stemming from the requirements of energy-momentum conservation and anomalous Weyl symmetry. The Mathematica notebooks code features plenty of comments which make them readily usable and (hopefully) understandable. We exploit the recently developed Package-X [24,25] together with its CollierLink extension for fast numerical evaluation of tensor integrals through the COLLIER library [26][27][28][29].
We did not attempt to track down the few [17] independent tensor structures our correlator should ultimately consist of, by itself a very demanding task given the sheer dimensions of the four-point function. We leave an attempt in this direction for possible future work.
The main reason why we undertook this calculation, beside the fact that we could do it, is the hope that our result could serve as a benchmark check (most probably numerical, though the result itself is fully analytical) for computations carried out with non perturbative methods, e.g. the conformal bootstrap for fields with spin. This paper is organized as follows. In section 2 we describe the setup of our calculation and the diagrammatic expansion of all the energy-momentum tensor correlators through four point, in order to make our discussion as self-contained as possible. In section 3 we derive all the transverse and Ward identities stemming from the requirements of general covariance (named transverse Ward identities hereafter) and Weyl invariance (named trace Ward identities hereafter), the latter of which feature a well known trace anomaly. In section 4 we describe the computation of the 1 loop countertems and the anomaly functional for the energy-momentum tensor correlators through four points and how we can check them against each other independently of the Feynman diagram computations; the countertems are then used for a preliminary check of the diagrammatic expansion, by comparing the UV pole of the correlators to the counterterms themselves. In section 5 we illustrate in detail the Mathematica implementation of our calculation and the procedure followed to analytically test the Ward identities for the two and three-point correlators and to test them numerically for the four-point correlator. Indeed, the last task is too massive to be dealt with analytically. In this section we also give plenty of detail on the implementation of all our computations in Mathematica through Package-X and CollierLink . Section 6 presents our conclusions and perspectives for further work.

Setting up the computation
This introductory section closely follows Chapter 2 of [30]. For details about the general relation between conformal invariance in flat space and Weyl invariance in curved space see [31]; here we assume that Weyl invariance in curved space is interchangeable with conformal invariance in flat space, as it is actually the case for the theory we are dealing with.
The standard definition of the energy-momentum tensor (EMT in the following) in a classical field theory described by an action S is given in terms of a functional derivative w.r.t. the metric tensor g µν (x) once the theory has been embedded in curved space, i.e.
Introducing the generating functional of the theory in Euclidean conventions, which we call W, where N is a normalization constant, Φ generally indicates the quantum fields of the theory and the second term represents the coupling of matter currents V a1,µ to background gauge fields. Given (2.1), the vacuum expectation value (vev) of the EMT is with the subscript s (which stands for "source") meaning that the background fields are not switched off. Dependence on coordinates will be occasionally dropped if it is not strictly necessary to the understanding of formulas.
In a conformal field theory the trace of the EMT vanishes at the classical level upon using the equations of motion, T µ µ = 0. As for the vev of the EMT, this relationship is modified by the well known trace anomaly [32,33] (2.4) where g and A a are short-hand notations respectively for the background metric and gauge fields. The coefficients κ, β a , β b , β c and β d depend on the field content of the Lagrangian theory and we have a multiplicity factor n I for each particle species. Eq. (2.4) is a reorganization in terms of the squared Weyl tensor and the Euler characteristic of 4d spacetime (see Appendix A.2) of the most general linear combination of the squares of the Riemann tensor and its contractions. The coefficient of R 2 must actually vanish identically since it does not satisfy the Wess-Zumino consistency condition for conformal anomalies [34,35]. In addition, the value of β c is regularization dependent, corresponding to the fact that it can be changed by the addition of an arbitrary local term in the effective action proportional to the integral of R 2 [32]. In particular, the values for the anomaly coefficients which we will use here for one single scalar field, i.e. β a = 3 5760 π 2 , β b = − 1 5760 π 2 (2.6) hold in dimensional regularization, for which one finds the constraint [32,33] For the purpose of this paper, the gauge field sector is not concerned, so from now on we drop it and the trace anomaly functional only depends on the metric tensor, The conformally invariant action for a scalar field coupled to gravity is given by Here χ is the parameter corresponding to the term of improvement obtained by coupling φ 2 to the scalar curvature R. In general dimensions d, χ = 1/4 (d − 2)/(d − 1) gives a CFT, so that χ = 1/6 gives a classically conformal invariant theory in d = 4, whose EMT is The explicit expressions for the vertices involving one or more metric, which can be computed by functional differentiating the action, have been already given in [36] and are also here collected in Appendix B for completeness, beside being explicitly computed in the attached Mathematica file functional derivatives.nb. For further details about our conventions on covariant derivatives, Christoffel symbols and the Riemann tensor, see Appendix A.

The structure of the correlators: topologies and diagrams
Our first step is to exploit the background field method to build the correlation function of n energy-momentum tensors and evaluate it diagrammatically at 1 loop. Since the theory is non interacting, as the scalar field is not coupled to any other quantum field, there are no higher order perturbative contributions, as one can easily verify. Therefore, conformal invariance is such that a 1 loop calculation can tell us all that there is to know about EMT correlators. We work in dimensional regularization. Our definition of the n−EMT correlators (also nT in the following) is that of a symmetric n − th functional derivative of W w.r.t. the metric tensor, Symmetry comes from leaving factors 2/ √ g outside of the derivatives. Given that field theory in Minkowski space is simply given by an analytic continuation from 4d Euclidean theory we work with through δ µν → η µν , we will also refer to this vertex as to the "n-graviton " vertex. We choose to denote such correlators, which contain contact terms, with small angular brackets (< >). Contact terms are characterized in coordinate space by the presence of at least two gravitons at the same spacetime point. Such contact terms are absent by definition in the expression of correlation functions given by the expectation value of the product of n EMT's, which are denoted with large angular brackets ( ) as in This distinction will not apply only to the vev of n insertions of the EMT, but also to contact terms and, in general, to all the correlation functions appearing in this paper. It will also be useful to introduce the following notation to represent the functional derivative with respect to the background metric evaluated in flat space, Since with our conventions all the functional derivatives will be taken w.r.t. the metric tensor with lower indices, which thus produces upper indices, possible functional derivatives with lower indices will mean just (2.14) In order to make tensorial expressions more compact, if we happen to contract two indices with the metric, we will stack the two contracted indices on top of each other, as in With these definitions, a single functional derivative of the action in a correlation function is always equivalent, modulo a factor, to an insertion of a T µν in the flat limit, since To begin with the easiest correlation functions, the definition of Eq. (2.11) implies that the two-point function is The last term on the right hand side of the equation above, which is a massless tadpole, can be set to zero in dimensional regularization, so that the 2T correlation function, obtained by differentiation of the generating functional, coincides with the quantum average of two energy-momentum tensors (2.18) This will not be true for higher order correlation functions, where non vanishing contact terms also appear. For the sake of completeness, the 1 loop contribution to the two-point correlation function is illustrated in Fig. 2.
The 3T correlator functional expansion is given instead by whose right hand side is expressed in terms of ordinary three-point correlators plus extra contact terms. The additional terms obtained by permutation are such to render symmetric the right hand side of the previous equation.
The first term on the right hand side of Eq. (2.19) is an ordinary three-point function, whose connected component is given, at 1 loop, by the triangle diagram of Fig. 3, while the last term is a massless tadpole (see Fig. 1) 2 , which can be set to zero (2.20) In the 3T case, contact terms have the topology of a bubble and are generated by correlators containing insertions of the second functional derivative of the action w.r.t. to the metric. Their general structure is shown in Fig. 3 and each one of them is simply obtained by assigning the momenta to the diagram according to one of the the three possible groupings.
Moving finally to the 4T case, a similar expansion holds and is given by being a massless tadpole contribution. Notice that the left hand side and the right hand side are both symmetric, as they should. In this case the perturbative expansion generates three diagrams of box type, represented by the Green function with 4 EMT insertions on the right hand side of (2.21), plus triangle, bubble and tadpole diagrams generated by the contact terms and graphically represented in Figs. 4 and 5. The analysis of these contributions is more involved compared to the 3T case. In Figs. Fig. 4 and 5 we illustrate the general structures of the four kinds of diagrams involved. The momenta running through them in addition to the loop momentum l must be replaced by the specific assignments detailed below. Later, when we illustrate the implementation of the computation in Mathematica , we will refer to these basic diagrams corresponding to inequivalent topologies as blueprint diagrams. Figure 4. The square and triangle topologies contributing to the < T T T T > correlation function. For the case of fermions, also loop momentum flow in the opposite sense must be considered. All the external momenta are incoming. Figure 5. The two kinds of bubble topology contributing to the < T T T T > correlation function in a scaleless theory. All the external momenta are incoming.
For each topology we have a different different number of contributions, corresponding to cyclically inequivalent orderings of the external momenta.
For the square topology, the 3 distinct contributions can be parametrized, for example, by the following three assignments of momenta (compare to the first diagram in Fig. 4) For the triangle topology, there are the following 6 distinct contributions (2.24) As for the two bubble topologies, we choose to distinguish them through the nomenclature of 22-bubble and 31-bubble, the figures standing for the numbers of gravitons meeting in each of their respective vertices.
The 3 contributions for the 22-bubble are (2.25) The 4 inequivalent contributions for the 31-bubbles are instead The task of computing the three-point function in a completely off-shell configuration and checking at the same time the transverse and trace Ward identities was already performed in [6], but only the result in a partially on-shell configuration was explicitly given. Now we face a more demanding task, which is the computation of the four-point function in a completely off-shell kinematics. This was made possible by the development of Package-X [25], a Mathematica package for tensor algebra and automatic reduction of 1 loop tensor integrals of any rank in arbitrary even dimensions. We detail the computation in section 5. For now, we just mention that also the step-by-step computations of the two and (much more significanly) the three-point functions are made available, so as to give the reader a full overview of the automated procedure in simpler setups, beside the target case of the 4T .
The building blocks of our computation are the scalar-graviton interaction vertices, which we report in Appendix B. Although we were completely confident about their calculation and we tested Package-X extensively by reproducing all of our previous results of [6] and a few other field theory results, the reduction of rank-8 tensor integrals still was an uncharted territory. Just as it was the case for the 3T , we checked all of our computations by testing the transverse and trace Ward identities our four-point correlator was supposed to satisfy, which we derive in the next section.

Derivation of the transverse and trace Ward identities
In this section, we derive the Ward identities stemming form the requirements of invariance of the generating functional under general diffeomorphisms and Weyl transformations. We call them transverse and trace Ward identities respectively. We proceed with a derivation of the relevant Ward identities satisfied by the two, three and four-point functions of the EMT.
Invariance under diffeomorphisms is defined by the condition of general covariance on the generating functional W[g], which translates into where in the last step, crucially for the symmetry of the correlators in the derived Ward identities below, we have exploited the cancellation of the last term on the rhs of the second line with the derivative of 1/ √ g x 1 .
The Ward identities for symmetric correlators we are after are obtained by functional differentiation of Eq. (3.1) as many times as it takes to get the correlator we are interested in, followed by taking the flat limit.
For the 2T we have the well known transverse condition, For the 3T we see an already non trivial rhs showing up Much more involved is the case for the 4T correlator, whose transverse Ward identity in coordinate space is given by for the three and four-point functions. The functional derivatives of the Christoffel symbol, obtained from the expansion of the covariant derivative which appear in previous equations, are explicitly given in Appendix B.
Before moving to momentum space, we point out that the Fourier transform formula is defined by (A.1), where all the momenta are incoming. We keep this convention throughout the paper and in all the momentum space computations in our code. By applying(A.1) and some integrations by parts to our Ward identities, we get the momentum space transverse Ward identities, which are a set of algebraic contraints on our correlators.
The momentum space 2T correlator satisfies For the three and four-point functions, the transverse Ward identities are quite cumbersome to write down in a fully expanded way.
We present an expanded version of the first transverse Ward identity for the theepoint function for illustration, whereas we keep the full sets given below implicit and refer the reader to Appendix A for a list of the momentum space explicit forms of the constituting elements For every nT , there is of course one such Ward identity for each EMT. We present the full set here below for the three-point function in compact form, Finally, for the four-point correlator we have 4 transverse Ward identities (this is the last set for which we report all channels; in what follows only one channel for every Ward identity will be explicitly written, though all of them were tested), Deriving the trace Ward identities is simpler. Rewriting (2.4) with n I = 1, κ = 0 as we can functionally differentiate up to three more times and obtain, after Fouriertransforming to momentum space with (A.1), The structure of the anomalies for our correlators is implicitly given by whereas the explicit construction is available in the notebook functional derivatives.nb.
In the next section we discuss in detail their connection with 1 loop counterterms and how we can use the two for a preliminary test of our calculation of the four-point correlator with Package-X . One comment about the anomaly is in order: since the term ∝ R can be removed by either adding an integral propto d d xR 2 [38] or by promoting the numerical coefficients in the Weyl counterterm to be functions of the spacetime dimension d [1], many an author choose to do so. Elsewhere we have discussed these issues in detail too [6], but we will not linger over it here. In this paper, we just perform the calculation in dimensional regularization with no supplemental local counterterm, so that (2.8) is the correct generating functional for the trace anomaly of our correlators.
Finally, there are the Ward identities implied by special conformal transformations, but we will not deal with them in the present work. For a detailed discussion of the derivation and the solution of these identities for two and three-point functions, see [4].

Counterterms, anomalies and a preliminary test of the 4T correlator
In this section we review the derivation of the 1 loop couterterms for EMT correlators and we derive the Ward identities which they must satisfy and which connect them with the respective trace anomalies. We explicitly compute these counterterms and check all of the mentioned Ward identities. Since both the counterterms and the anomalies for the 4T correlator are highly non trivial structures, the successful check of these tests is a very strong indicator that their calculation is likely to be correct. Therefore, if it is possible to evaluate only the ultraviolet pole of our four-point function, after putting together the tensor integrals corresponding to the diagrams of section 2, we can compare it to the independently tested counterterm. If the two expressions match, we have a strong preliminary hint that the diagrams have been assembled correctly. We actually performed this test, as the LoopIntegrate routine of Package-X , which substitutes explicit expressions of scalar coefficients into tensor integrals, has an option for extracting only the ultraviolet pole of each scalar coefficient. Summing up the poles of all diagrams, we did match the rank-8 counterterm for the scalar four-point function.
In order to be as self-contained as possible, we provide here the general formulas for the Passarino-Veltman coefficient functions we employ [39], so as to clarify the meaning of the symbols P V B, P V C and P V D that the reader will encounter in the snippets of code to follow. In the following formula, the symbol r stands for the number of times the metric tensor enters in the symmetric tensor multiplying the coefficient function, whereas n 1 , n 2 and n 3 stand for the number each of the external momenta enters in it, ..µ 2r+n 1 +n 2 +n 3 . The syntax of the routines employed in the code in Fig. 6 is the following • LoopIntegrate is the central routine of the package and performs the reduction of the tensor integrals. Its first argument is the numerator of the tensor integral, the second is the loop momentum, the following pair given as arguments have the structure (mom, mass) where mom is the momentum and m is the mass of the propagating particle for each propagator in the loops. The mass is always 0 for us. Various options can be given, for which we refer the reader to the Package-X guide, included in the Mathematica interactive documentation upon successful installation of the package.
• LoopRefine converts the Passarino-Veltman coefficient functions to analytic expressions which can be readily evaluated numerically, extracting the ultraviolet pole from the scalar two-point function. Various options can be given, among which is Part → UVdivergent, which makes the routine extract only the ultraviolet pole of the expression; for other options we refer the reader to the Package-X guide.  Figure 6. Loading Package-X in Mathematica and extracting the UV pole of the scalar two-point function In the following, we discuss in detail the 1 loop counterterms and the Ward identities they are subjected to. We refer the interested reader to our public repository for details about the calculation not included in the paper. This is actually not the first time that the four-point function counterterms were computed, as they were already worked out in [40], where their relation to the trace anomaly -to be discussed shortly -was exploited to explore dilaton interactions in the conformal limit of the Standard Model. Another application was explored in [41], where a recursive relation allowing to compute fully traced correlators of the energy-momentum tensor was discovered.
The term we must add to the action in our generating functional W in order to renormalize our (unrenormalizable) theory in (2.2) is the following, with d = 4 − 2 , containing the squared Weyl tensor F and the Euler density, defined in Appendix A.
For the general nT correlator, the counterterm action (4.2) generates the n-point vertex where The momentum space couterterms are defined via the same Fourier transform defining the momentum space correlation functions (A.1).
We have derived in detail the first functional variation of the integral of a general expression which is quadratic in the Riemann tensor in a shared appendix of [6,30], to which we refer the interested reader. Here we provide only the final result.
Let K = a R αβγδ R αβγδ + b R αβ R αβ + c R 2 , with a, b and carbitrary real numbers; then From the equation above, the counterterm action (4.2) and taking traces in d dimensions, since we are working in dimensional regularization, one can derive the trace anomaly observing that where we have distinguished the d-dimensional trace with a superscript on the metric tensor.
Using these expressions, the renormalized two, three and four-point correlators in momentum space can be written as From these relations and from (3.5), (3.7) and (3.8) it is apparent that the counterterms are related to each other by the same transverse Ward identities relating the EMT correlators. One can also separately check these identites for Weyl and Euler counterterms just by writing them down and equating the coefficients of β a and β b . This is done in the second part of the file functional derivatives.nb.
Anomalous Ward identities for our counterterms can be derived through up to three more functional derivatives of Eqs. (4.7), so that identities which are completely analogous to (3.10)-(3.12) emerge, the only difference being that now traces are taken in d dimensions. We report them separately for the Weyl and Euler counterterms in the first channel: of course, exactly analogous identities -obtained by proper permutations of indices and momenta-hold in all channels.
It is apparent that all these constraints are not trivial to satisfy. We checked all of them for all anomalies and counterterms and, thus, ensured that these were correctly computed.
As mentioned before, the ultraviolet pole of our four-point function is found to coincide with the overall counterterm for the scalar field 4T . The reason why we perform this preliminary check is twofold.
• The ultraviolet pole of our four-point function is highly non trivial, since there are 16 diagrams contributing to it, so that getting it right is a solid first check of the diagrammatic expansion. Beside, being a polynomial in the external momenta, the ultraviolet pole of our four-point function it is much more handable than the whole correlator, which requires much more memory and computing time. This makes it an ideal tool to quickly scan for potential mistakes.
• The trace Ward identities connect the counterterms to the trace anomalies and, so, are a test for the latter as well. The trace anomalies are present also in turn in the trace Ward identities for the correlation functions (3.10)-(3.12), which we must test in order to ensure the correctness of the four-point function. Thus, having them tested in advance reassures us about the correctness of the trace identities, so that any mismatch emerging from them can be quite surely traced back to the four-point function.
One comment is in order here: every even dimension space with a local metric tensor and an associated Riemann tensor is characterised by a topological quantity, called the Euler characteristic, which is given by the integral of the Euler density over all space. This integral is a pure number and, therefore, its derivatives all vanish. It may surprise, then, that we are checking that its functional derivatives contribute to our anomaly functionals and our counterterms, but we are working in dimensional regularization, so the counterterms inhabit a d = 4 − 2 dimensional space, making it non trivial and non vanishing. Beside, the functional derivatives of the Euler characteristic are zero in d = 4, but "in disguise" [17], meaning that the tensor does not vanish identically, but explicit evaluation of all of its components would give zero. This ensures that symbolic programs cannot tell the difference between such a tensor and any other, so checking the consistency of our computation by looking for such structures on both sides of the identities featuring them is still a further non trivial check of our calculations. To the best of our knowledge, the difference between the Weyl and Euler anomaly due to the topological nature of the Euler density was first pointed out in [42]. Once all these preliminary tests have been successful, we can be reasonably confident about the diagrammatic expansion. What is left it the computation of the whole four-point function and the test of its the transverse and trace Ward identities. We give plenty of details about the procedure, reproducing it in full also for the two and three-point functions in the files correlators calculation.nb and ward identities.nb. so as to present a self-contained overview and let the reader familiarize with the tools we employ in simpler setups.
The purpose of the following section is to go into the necessary technical details about the calculation.

Into the full calculation
It is time to illustrate in detail how we employed our tools to nail down the whole four-point function. In the first part we review the Feynman diagram computation, in the next two we give a survey of the checks of the Ward identities. Snippets of code are provided throughout this section, but we encourage the interested reader to explore the notebooks, which are documented in great and hopefully sufficient detail.

The calculation of the four-point correlator
The calculation is performed by the file correlators calculation.nb. The code employs Package-X both to automate the Passarino-Veltman reduction of the many tensor integrals one encounters.
The notebook starts by loading the vertices and counterterms stored in the all functional derivatives file, which must be generated beforehand by running functional derivatives.nb. After that, the notebook is divided into three parts, one for each of the three computed correlators. The structure of three sections is similar and unfolds as follows.
• The numerators of each contributing diagrams are constructed: see Fig. 7. As detailed in section 2, there is only 1 for the 2T , while there are 4 for the 3T and 16 for the 4T . An image of each computed Feynman diagram is loaded in the notebook right before the line computing its numerator. In the case of the twopoint function, the known result for general values of the improvement parameter (see Eq. 2.9) is presented and rederived for the reader's convenience. The case with χ = 1/6 corresponds to the conformally invariant case and is eventually selected.
• The LoopIntegrate routine is invoked to perform the tensor reduction of the diagrams, see Fig. 8. For the three bubbles in the three-point function bubbles and for all the diagrams contributing to the four-point function, only blueprint diagrams with generic momenta are computed. The actual diagrams are obtained by replacing the generic momenta with the sets classified in section 2. This is topical for the four-point function, due to the memory it requires (see next point).

In[!]:=
Unfortunately, this is not the case for the tests of the Ward identities for the 4T correlator, as we will explain in a while. The next two sections deal with the content of the file ward identities.nb.

Analytical checks of the Ward identities for the 2T and 3T
Checking the transverse and trace Ward identities for the two-point function, Eqs. (3.5) and (3.10) is trivial. The first two lines in fig. 11 should be self explanatory, once the reader has understood the essential purpose of the LoopRefine routine explained in section 4 and that the A2 object is just the anomaly on the rhs of (3.10).
It is now time to explain how to check the naive scaling dimension of the correlators, the underpinning idea being very simple. This dimensions is called naive because proper scaling accounts for the effect of logarithms depedent on the ultraviolet scale ∝ log(p 2 /µ 2 ), which come from the two-point scalar integrals. Since every term in our EMT correlators must be of dimension (momentum) 4 , we can replace every term in our correlator with its dimension, which we choose to denote with λ in the code. So, for instance, p µ → λ, g µν → λ 0 . Once all the replacements have gone through, we expect to get just a real number times λ 4 , which is what is shown in the snippet below and what is found for the two higher-point functions as well. = p2 ν2 T2[μ1, ν1, μ2, ν2, p2 Figure 11. Check of the transverse and trace Ward identities as well as of the naive scaling dimensione for the two-point function.

In[!]:
One might ask why we do not undertake a check of the full dilatation identity, which would be, in d = 4, something like The reason is that checking the trace identities (3.10)-(3.12) already implies satisfying the scale invariance requirement. Indeed, the trace anomaly equation (2.4) is obtained by studying the transformation of the generating functional of a field theory embedded in curved space under a local rescaling of the metric tensor, which is called a Weyl transformation. This is a more general transformation than a global rescaling of coordinates and, for Lagrangian field theories which are at most quadratic in their fields derivatives, is effectively equivalent to full-fledged conformal invariance [31]. This means in turn that, since this symmetry is broken only by the trace anomaly, which can be seen as a by-product of the ultraviolet renormalization of the EMT, then it follows that, if an EMT correlator satisfies its anomalous trace identity, the correlator must also be scale invariant in the sense of (5.1). In particular this implies that, if its momenta are uniformly rescaled by a global factor p → λp, it must change by an overall factor λ 4 plus additional contributions due to the logarithms associated to the regularization scale µ 2 , introduced to consistently regularize ultraviolet divergencies ∝ log(p 2 /µ 2 ) → log(p 2 /µ 2 ) + log λ 2 . These are the contributions that, acted upon by the differential operators in (5.1), would render the anomalous term on the rhs (5.1). The fact that in the last line of Fig 11 we are just cheking the naive scale dimension is due to the replacement of the two-point scalar integrals with λ 0 , which does not take into account the logarithmic contribution. Now, the fact is that the check of the naive scale dimension can be done analytically also for our most complicated correlator, the four-point function and, as such, will be useful in the last part of this section, when numerical stability of our results will be discussed.
Next we come to the tests of the Ward identities for the three-point function, which is more involved but was already done in [6] in pretty much the same way, which we reproduce here. The idea behind the procedure for both transverse and trace identities is the same and quite simple.
• Load the tensor basis which spans the space of rank-7 or rank-6 tensors dependent on 2 independent momenta to which belongs the correlator contracted with the momentum in (3.7) or with the metric tensor in (4.10); assign to them the correct indices (those which survive after contraction).
• Build the lhs and rhs for each Ward identity separately.

Numerical checks of the Ward identities for the 4T correlator
This is the most slippery and time consuming part of our project, because the whole 4T correlator is just too big to be dealt with analytically within reasonble time with an ordinary computer, even one with a 64GB RAM like the one at our disposal. As we mentioned above, the blueprint of the square diagrams contributing to the 4T are the culprits, occupying alone over 7GB of memory space, which more than triple when three of them are summed after the proper momenta assignments. This is because the latter imply that the second argument of each diagram is the sum of two external momenta (see the notebook and the diagrams of Section 2) which, considering the amount of terms it consists of, significantly increases the necessary memory. The first feature which is required from a computer is quite some RAM: at least 40 available GB, in order for the computation not to crush. The second crucial feature is to have more than one available CPU core to perform calculations, in order to reduce the required time to a reasonable window through parallelization. In our case, we had 8 kernels and a 64GB RAM at our disposal on our machine, which allowed us to complete the test of all the 9 Ward identities for the 4T correlator in slightly less than 1 hour. 4 In order to speed up our numerical effort, we resorted to the CollierLink package, which connects Package-X with the COLLIER library (written in C language) in order to provide fast numerical evaluation of the scalar integrals. A fundamental component of the numerical evaluation procedure is to isolate the ultraviolet poles from the rest of the correlators.
The steps of the testing procedure explained for the 3T Ward identities are perfectly valid here as well, with two caveats. First, the tensor bases to be used are now much bigger due to the increased number of indices from 5 to 7 for transverse Ward identities and from 4 to 6 for trace Ward identities, beside the fact that we have one more independent momentum in our correlator w.r.t. the three-point case; second, the matching of the coefficients must now be checked numerically in order to guarantee sufficient speed and, for most common computers, not to crush.
On the ground of what was explained so far, one might think that subtracting the counterterm for the 4T after applying the LoopRefine routine to the computed correlator -even parallelizing it-would do the job: it would indeed, but that would be much too slow.
What one can do, instead, is to resort to CollierLink , which can perform direct numerical evaluation of the Passarino-Veltman coefficient functions bypassing the analytic expansion. In Fig. 13 is a snippet showing the 16 diagrams arranged in a list (notice the different momenta routings and compare to the diagrams in section 2). Here, the ParallelMap routine maps the pure function which is its first argument in parallel to all diagrams, each one of which is handled by one subkernel. The parallelized task is in turn a sequence of two steps.
• The SeparateUV routine of CollierLink isolates the UV pole in each coefficient function.
• The numrep set of rules, defined above in the notebook, replaces all the scalar products of momenta with the corresponding numerical values obtained after choosing a completely off-shell configuration which respects p 1 + p 2 + p 3 + p 4 = 0.
Once the numerical evaluation of the scalar coefficients has been performed, what is left in order to test the Ward identities is to contract the listed diagrams with either one of the momenta or the metric tensor, perform the few more numerical replacements on the scalar products this produces, perform the same operations on the correspondong rhs (see Eqs. 3.8 and 3.12) and compare the scalar coefficients one by one. Then, if all the differences are acceptably small, we can declare the test successful. The corresponding snippet of code is given in Fig. 14, which we comment below • lhscov1 is the lhs of (3.8) and the set of rules in the curly brackets are the (faster) equivalent of the contraction with the p 1 momentum. When momentum conservation is employed at the end, it can affect only tensor structures, as it must in order not to miss any terms, since the tensor basis tens371 depends only on 3 independent momenta • rhscov1 is the rhs of Eq. (3.8) , the Γ 1 and Γ 2 symbols being defined in Appendix B (or, equivalently, in the notebook functional derivatives.nb), whereas the T 2 and the T 3 functions are obviously the two and three-point correlators.
• checkcov1 is a table, evaluated in parallel, made up by the differences of all scalar coefficients of the lhs and the rhs over the basis of rank 7 tensors with 3 independent momenta. The additional function Labeled just marks every difference with the number of the tensor in the list, which could be needed after an unsuccessful test to track down the tensors whose coefficients on the lhs and rhs of the identity do not match.
• The final line simply removes from the list all those numbers which were set equal to 0 by the Chop routine, because both their real and imaginary parts were under the prec threshold, which we set to 10 −8 for the first three momentum configurations listed in (5.2). A dedicated discussion of the fourth momentum configuration, for which the precision threshold we manage to reach is not even close, being just 10 −1 , aims at proving that this is definitely an expected issue of numerical instability. To this end, the fact that we can easily check analytically naive scaling dimension is useful.
when specific momenta are plugged into them. This is done in the very last few lines of the notebook.

Discussion about numerical stability
The last point we wish to discuss before coming to our conclusions is the issue of numerical stability, which is raised when one tries to check the four-point functions Ward identities for momenta configurations like the last one in (5.2), which was purportedly chosen to be a rescaled version of the back-to-back third configuration, for which the test is precise through 8 decimal figures, just as for the former 3 as well. This is when the tests of the naive scale invariance of our correlators comes in handy, as we finally explain, reassuring us that the problem is only due to numerical instability. Now, the naive scaling dimension of every EMT correlator is (momentum) 4 , as one can make sure analytically for all correlators (see the code). This means that both the lhs and the rhs of every identity in (3.8) and (3.12) should naively scale by the same factor λ. To be completely precise, one should not forget to take into account the logarithmic contributions associated with scalar two-point functions which, as mentioned in section 3, behave under rescaling like log(p 2 /µ 2 ) → log(λ 2 p 2 /µ 2 ) = log(p 2 /µ 2 ) + log λ 2 .
If it comes to numerically checking the rescaled identities, we can take care of this problem in two ways: • by rescaling the regularization scale as well by µ → λµ, using the dedicated initialization option of the CollierLink package. The snippet of code in Fig.  15 illustrates the invariance of the two-point function after rescaling both the momentum and the regularization scale.
• by leaving everything as it is, since the same two-point functions are present on both the lhs and rhs of the Ward identity, so that extra terms should coincide as well.
In  So this proves that, if the identities are satisfied for a given momentum configuration, they should be identically satisfied for the rescaled configuration as well. Testing both possibilities is useful to realize that the extra logs coming from the rescaling of the momenta are not the source of the numerical instability, or at least not the only one.
Indeed, what we find for our fourth momentum configuration, whether we take care of µ in the first way suggested above or not, is that the precision up to which the Ward identity is numerically satisfied does not exceed 10 −1 for some scalar coefficients in the tensor bases. Nevertheless, since we proved that this result must vanish exactly too, the fact that the numerical agreement is excellent for such diverse momenta configurations as the first three leads us to the conclusion that our calculation is correct, although its numerical evaluation suffers from numerical instability, which is nevertheless to be expected for such huge expressions.
Of course, one can think of ways to improve our numerical tests. The naive way could be to pursue a brute force approach: one would feed the numerical routines higher and higher precision inputs (Mathematica can numerically evaluatae exact numbers to arbitrary numerical precision), but this would take an even heavier toll on the memory of the machine, which is already quite stretched with the required accuracy standards. A batch of several machines on which we could parallelize our numerical tests could do the job, of course. We do not have access to an integrated Mathematica deployment of this kind, but we firmly believe that the discussion above has clarified that this would add neither any further understanding nor improvement to our result, not least because, once one has a much higher available computing power, the identities could be tested analytically right off the bat.
A second way to go could be to compile our output into optimized C++ or Fortran code. This is presumably feasible, but it goes behind the scope of the present work.
These problems will stay, of course, only until a more clever way to compute the four-point function correlator is devised, which can yield a more compact result cutting off all the redundancies our method necessarily implies, but for the time being we consider ourselves satisfied with our result.

Conclusions and perspectives
We have computed the four-point correlation function of the energy-momentum tensor and provided a set of tools which allow the interested reader to reproduce the perturbative computation of the two, three and four-point correlation functions, together with the explicit and detailed construction of counterterms and anomalies and the test of the transverse and trace Ward identities. The results for the correlators are fully analytical, expressed in terms of Passarino-Veltman coefficient functions which can be easily either fully reduced to scalar integrals or evaluated numerically to arbitrary precision.
Considering the sheer dimension of our results, the main purpose of this work should be understood as the delivery of a benchmark for numerical checks of more compact results to be obtained in the future with different strategies, most probably the conformal bootstrap for fields with spin.
It would also be interesting to perform the same calculation for other Lagrangian conformal field theories, such as a fermion or a gauge field: actually, we already did this, but decided not to publish these results because they are even more massive than the scalar case and checking the Ward identities for the corresponding four-point correlation functions is computationally even more demanding. Furthermore, unlikely the case of the three-point function, the results for the scalar, fermion and gauge fields do not allow to account for the full set of constants which would allow one to reconstruct any four-point correlator, so one does not gain much further insight into the structure of the CFT by computing them. Anyway, should we make progress on speeding up such checks, we will certainly update our repository.
Given the result by Dymarsky [17], it would be also very useful to identify the 22 independent structures making up the 4T correlation function in 4d in momentum space, perhaps after successfully accomplishing the task in 3d, where the number of independent structures is just 5. Once and if this task is accomplished, it would be very useful to check the expected one-to-one correspondence between the momentum space result and the coordinate space results of [18][19][20][21][22]. The mapping between correlators in coordinate and momentum space has been investigated in some detail in [6,30] and shown to lead to some non trivial consistency conditions, which make this effort worth a separate work.
Also checking conformal Ward identities for our four-point correlator in momentum space would be interesting, but the task requires the implementation of second order differential operators and their application to a very complicated object. Since the Ward identities for the four-point correlator we did test were manageable only numerically, which was already a highly non trivial task, we did not try to figure out a way to perform the same numerical task as efficiently when working with differential operators, not least because our perturbative results for the lower point functions were shown in [7] to coincide with the non perturbative results of [4,5] The latter were obtained by solving the full set of conformal group Ward identities, so the work of implementing differential operators in our codes for the three-point function without being able to extend the procedure to the four-point function would have added no original output to our effort. Again, should we make progress on this point, we will update our repository. the development of this project. I am very grateful to Hiren Patel for sharing insightful suggestions about Package-X , to Andreas van Hameren for pointing it to my attention and for making several useful comments about the numerical tests. Special thanks go to Misha Lublinsky for allowing me to deploy this computation on a powerful computer of the Ben Gurion University of the Negev. Useful discussions with Claudio Corianò and Matteo Maria Maglio are also gratefully acknowledged.

A.1 Conventions for signs and momentum space correlators definition
The definition of the Fourier transform of the correlation function of n EMT's, which holds for any other n-point coordinate space function as well, is given by . . . T µnνn (p n (A.1) where all the momenta are conventionally taken to be incoming.
The covariant derivatives of a contravariant vector A µ and of a covariant one B µ are respectively with the Christoffel symbols defined as Γ α βγ = 1 2 g ακ [−∂ κ g βγ + ∂ β g κγ + ∂ γ g κβ ] . (A.4) Our definition of the Riemann tensor is The Ricci tensor is defined by the contraction R µν = R λ µλν and the scalar curvature by R = g µν R µν .

A.2 Weyl invariant and Euler density in 4d
It is well known that the object one has to deal with in order to construct Weyl-invariant objects for general dimensions d is the traceless part of the Riemann tensor, called the Weyl tensor, defined by (g αγ g δβ +g αδ g γβ )R . In 4 dimensions, the only quantity which is Weyl invariant when multiplied by √ g, is the Weyl tensor squared, which is The other quantity with which one can construct an integral which is Weyl invariant (in fact, a constant) for general even dimensions d = 2 k is the Euler density, defined as E 2k = 1 2 k δ µ 1 a 1 ν 1 b 1 ...µ k a k ν k b k R µ 1 ν 1 λ 1 κ 1 . . . R µ k ν k a k b k . (A.12) The antisymmetric Kronecker symbol is given by δ ν 1 a 1 ν 2 a 2 ...νnan = n! P(a 1 ,...,an) (−1) T P g ν 1 P(a 1 ) . . . g νnP(an) , (A. 13) where T P is the number of inversions in the permutation P of the n numbers a 1 , . . . a n . By applying the general definition A.12, we find that in 4 dimensions we have

B Functional derivatives B.1 Basic functional derivatives
Here we give the basic momentum space functional derivatives, which are building blocks needed to for all of the counterterms, vertices and anomalies. Once these few formulas are understood, the construction of all the rest is a matter of (very) careful bookkeeping and a lot of applications of the chain rule for functional derivatives: please notice that in the following formulas terms which vanish in the flat spacetime limit are already dropped, but one must take them into account for higher order derivatives.

B.2 Interaction vertices of the scalar field with gravitons
We provide here the explicit forms of the three vertices used in the Feynman diagrams. The computation of the vertices can be done by taking at most three functional derivatives of the scalar action in curved space with respect to the metric tensor and, of course, the scalar field. This is because the vev's of the fourth order derivatives correspond to massless tadpoles, which are set to zero in Dimensional Regularization, so no fourth derivative is needed.