Three twist-two, two spins, two loops

I consider three-point functions of twist-two operators in N=4 SYM, two of which endowed with spin. I supply perturbative data up to twelve units of spins and second perturbative order at weak coupling.


Introduction
As a conformal field theory, N = 4 SYM is defined by the spectrum of its operators and their structure constants. Both problems can be attacked, in principle, with integrability [1,2]. This note focuses on structure constants. In particular, it addresses the computation of three-point functions involving more than one operator with spin. From the perspective of the operator-product expansion (OPE) [3,4], such coefficients emerge from a four-point function involving one un-protected operator, or alternatively from a multiple OPE of a higher-point function of protected ones, such as a five-point correlator. Beyond one-loop order, the computation of these correlators is not as well developed as for four operators, and mostly limited to the work [5]. Recently, the large spin of such higher-point correlators has been analyzed [6,7]. This establishes the exact behaviour of structure constants with various spinning operators at all loops, in the large spins limit, and a duality with null Wilson loops.
The computation presented here lies in the opposite regime, that is small spins. The calculation is performed perturbatively at weak coupling. In this setting, high values of the operators spins are a nuisance rather than a blessing (ten units of spin already means high in this article). In fact, they are the main source of complexity and computational bottlenecks, which this work aims to attack. Being the simplest, twist-two operators are considered. Their spectrum is well-known from both explicit computation [8,9] and integrability [10][11][12][13]. Their three-point functions with two protected operators are also known to vertiginous perturbative precision [14][15][16][17][18][19]. Such results have been derived mostly from the OPE expansion of four-point correlation functions of protected operators.
Less is known for structure constants involving more than one operator with spin. This note intends to address such a paucity of data. The problem has already been considered in [20]. The present work is a non-trivial extension of that. The idea of [20,21] consists in tackling the computation employing an integration trick, first proposed in [22]. This can be interpreted as a soft limit in momentum space. The bottom-line (more details follow in the next section) is that it allows to trade the full three-point function for a much simpler two-point-function-like object. Such a simpler computation is sufficient for isolating the structure constant of the three-point function. The most striking aspect of the simplification lies in the computation of integrals. This is performed through reduction via integration-by-parts (IBP) identities. The latter is a well-developed topic in literature, especially so for this class of propagator integrals. Conformal invariance is instrumental in the extraction of the structure constants. It completely fixes the space-time structure of the correlator. The only dynamical information lies in its coefficient, the structure constant. This is one number at each perturbative order. Subtleties arise, due to regularization, which will be addressed momentarily.
The punchline of [20] is as follows. At one loop the computation is feasible up to high values of the spins, but still in a case-by-case analysis. A solution of the IBP reduction of the relevant integrals for generic spins would grant a complete answer, but was not attempted in [20]. Instead, the amount of data allowed for conjecturing a heuristic formula establishing the dependence of the structure constants on spins and polarization, though in terms of some undetermined coefficients. I was not able to progress further and fix such remaining coefficients. More insightful people managed to do it [7]. Therefore, the problem of computing such structure constants at one loop has been completely solved by [7].
At the next perturbative order a few results were computed in [20], demonstrating that they can be attacked with the method proposed. However, the computation stopped short at rather disappointingly low values of the spins, namely six. This was due to computational complexity, the main bottleneck being integral reduction. Yet, the limitation stemmed mostly from the naivety of the approach taken for completing such a reduction: easy to implement, but too inefficient for attacking a complicated problem. The aim of this paper is to advance further in the two-loop computation and work out a reasonable amount of perturbative data with modest computational resources and time. Thanks to various technical optimizations with respect to [20], a few-weeks computation on an ordinary laptop (plus some aid from a cluster) produced results up to twelve units of operators spins, which are presented here. Such data are too scarce to conjecture an analogous formula to that encountered at one loop. I was hoping to find one at low values of polarizations, where the one-loop formula simplifies considerably, but had no success. Still, this note provides some new solid results, which can be compared with alternative computation methods, such as OPE of higher-point functions and integrability.

The perturbative computation
The setting of the calculation is the same as in [20] and I refer to that for further details. The basic definitions of operators are reported in the appendix. This note focuses on three-point functions with two spinning operators. In a conformal field theory their form is fixed [23] where j 1 and j 2 are the operators spins. The quantities Y and I appearing on the right-hand side are certain invariants whose precise form is spelled out in the appendix. The powers of squared distances are combinations of the conformal dimensions of the operators. The hats symbols denote contractions of the tensor structures with two distinct sets of light-cone vectors. The index l labels the polarization of the three-point functions and to each value 0 ≤ l ≤ min(j 1 , j 2 ) an independent structure constant is attached. The objective is the two-loop computation of such structure constants. They are extracted integrating both sides of (2.1) on one of the integration points of the operators. The left-hand-side is expanded perturbatively in Feynman diagrams. At that level, the integration translates into a soft limit in momentum space. This effectively collapses the three-point function onto a two-point function. The latter, being much simpler than the original three-point problem, can be calculated efficiently, leveraging vast literature and techniques for dealing with two-point integrals. In practice, the soft limit introduces additional powers of some propagators. Conversely, the spins of the operators translate into powers of momenta. Both occurrences are dealt with through integration-by-parts identities (IBP) [24,25], reducing all integrals to a finite set of known master integrals. This is the crucial technical step, allowing for carrying out the computation. More details on the procedure are articulated in the next section.
On the right-hand side of (2.1), the integration over the insertion point of an operator was evaluated in general, as reported in the appendix (A.6). After doing that, (2.1) allows for the extraction of the desired coefficient C, by comparison of both sides. This is the punchline of the calculation method.
There are a few subtleties attached to this seemingly straightforward process. All arise due to the necessity of regularizing divergent intermediate quantities. There are various sources of singularities. First, UV divergence appear. These stem physically from the two-point functions of operators and ought to be renormalized away. In fact, this provides a powerful check for the self-consistency of the computation. Secondly, the soft limit enforced to simplify the computation introduces IR divergences. These are independent of the former and combine with them. Third, individual Feynman diagrams could possess spurious divergences which cancel away from the final result. A regulator is needed for carrying out the computation in practice. Dimensional regularization is used. The advantage of dimensional regularization consists mostly in being compatible with the IBP reduction of integrals, which are carried out in arbitrary dimensions. Moreover this is the regularization method which is usually deployed for renormalization of UV divergences. A different regulator might be employed for IR divergences. A mass regulator for instance could regulate such singularities and separate them clearly from the former. This would introduce a mass in the relevant integrals, hampering the computation considerably. Overall, I reckon dimensional regularization is the only viable option.
In dimensional regularization, IR and UV poles mix and multiply. After renormalizing operators, only IR divergences can be present, arising from the soft limit. Such divergences can behave rather heterogeneously depending on the conformal weights of the operators. Two qualitatively different scenarios may occur. The result may display a divergent behaviour with pole powers increasing with perturbative order; or it can present a fixed lower bound for powers at each perturbative order. The same pattern must appear integrating the generic, conformal form of the three-point function on the right-hand-side of (2.1). Since that integral might in general be divergent, it should be regularized with the same method as in the perturbative expansion of the left-handside, i.e. dimensional regularization, for consistency. If the order of poles is fixed for all perturbative orders, the coefficient of such poles emerges from the d = 4 limit of the integrated expression of the three-point function. This in turn means that only the known expression of the-point function in d = 4 is required. Any subleading in corrections can be safely discarded. The extraction of the structure constant is then correct. On the contrary, if the order of the divergences is not fixed, there is no sensible limit to d = 4 which can be taken. The only possible comparison between the two sides of the equation would involve considering the three-point function in d = 4 − 2 dimensions, whose form is not known in general. This impedes a sensible extraction of the structure constant.
In the problem of this paper, integrating on the insertion point of a protected operator or of a spinning one illustrates this distinction. The former yields a result with a constant 1/ pole, from which the structure constant can be extracted. The latter produces an output with increasing powers of poles, except for the null polarization case, where the order is constant. Indeed for null polarizations a structure constant can be extracted and is equal to the other integration result. For other polarizations it is just impossible to disentangle the perturbative corrections of the structure constant and spurious terms due to undetermined subleading in corrections to the three-point function hitting a higher order IR pole. In fact, a naive comparison fails. In conclusion there is only one sensible choice of the integration point for these structure constants, namely that of the protected operator.
More in general, regulating both sides of (2.1) introduces an order-of-limits issue. On the right-hand-side, using the conformal expression of the three-point function entails taking the d → 4 limit first. Then the soft limit is taken and dimensional regularization is again used to regulate the integration over an operator insertion point. On the left-hand-side, the soft limit is taken before on the Feynman diagrams, where dimensional regularization is applied. Only at the end of the perturbative computation the limit d → 4 is enforced. I could not find any argument why the two limits should commute. In fact, it was shown in [26] that in general they do not, by an explicit example. Luckily, it so happens that the present computation seems not to be affected by this issue. This is easy to verify at one-loop order from direct inspection of the integrals. A posteriori, this is confirmed by the fact that the one-loop structure constants computed in [20] coincide with those determined in [6,7] using an independent method. At two loops the same analysis is difficult due to the complexity of the calculation. Again a posteriori, the method reproduces correctly the known results for three-point functions with one spinning operator [27]. Since there are no qualitatively new diagrams and integrals associated to the computation with two spinning operators (only more momentum powers), I assume the order-of-limits issue to be absent and proceed.

Integrals treatment
The problem described above boils down to the reduction of various integral topologies with powers of momenta in their numerators. The latter are contracted with two different sets of null vectors z 1 and z 2 . This gives rise to polynomials in the scalar product z 1 · z 2 , with degree up to the minimum between the powers of z 1 and z 2 contractions with loop momenta. At two-loop order three-loop momentum integrals are needed.
IBP reduction. The reduction can be performed in various ways. The integrals can be IBP reduced including directly z 1 and z 2 contractions. This introduces new external null momenta z 1 and z 2 , with their respective IBP identities. The advantage of such an approach is that it is directly implementable on IBP reducers on the market. In particular, I used FIRE6 [28,29] with LiteRed solved rules [30,31] (which are more effective than Laporta reduction). The reduction works fine and fast for two-loop integrals and for three-loop integrals, up to certain powers. The disadvantage is that FIRE6 tends to crash (at least on my computer) for high enough powers of numerator momenta of some three-loop integrals. Not surprisingly, the onset of this behaviour seems to happen for lower powers of momenta, when including z 1 and z 2 dependence and their additional IBP rules, than the case of integrals with no such external momenta contractions. In particular the reduction fails for some integrals involved in the threepoint function with two spin-6 operators. That is the main reason why the earlier computation of [20] stopped short around this complexity level.
Tensor reduction An alternative approach consists in performing a reduction of tensor integrals to scalar ones, later projecting the numerator onto the z 1 and z 2 contractions. This process generates scalar integrals involving only products with the external momentum p. These are more easily and rapidly reduced via IBP reductions, than the aforementioned case with additional external momenta. In this case another bottleneck is the tensor reduction itself. The computational load can be alleviated considering the symmetry properties of both numerator momenta and of the resulting contraction with null vectors. This is however a case-by-case analysis. In order to speed up such a reduction I used the following method. Each numerator structure of momenta is defined by the set of indices to be contracted with z 1 or z 2 factors. Let's say that there are n 1 contractions with z 1 and n 2 contractions with z 2 . The indices to be contracted with z 1 and z 2 are simmetrized (and traceless) by construction. Therefore the reduction corresponds to a case with symmetric numerator powers where the indices α and β identify loop momenta in the integral. A generic ansatz for the right-hand-side of such a tensor integral is of the form . . . g µ i+k ν j+k p µ i+k+1 . . . p µn 1 p ν j+k+1 . . . p νn 2 + all perms among µ and ν indices with the same coefficient for various tensor structures, as indicated, thanks to symmetry. Eventually the permutations all evaluate to the same result after contractions with the null vectors, giving rise to a combinatorial factor. The important information resides in the coefficients c, which depend on the particular integral. They are indexed according to the number of different possible contractions among z 1 -and z 2to-be-contracted indices and mixed ones. These are labeled by i, j and k in the above formula. Only the coefficients c(n 1 , n 2 ) 0,0,k survive the contraction with null vectors. The system can be inverted, in principle, after multiplying both sides with suitable tensor structures. The inversion coefficients are by construction polynomials in z 12 and rational functions of d, which multiply powers of external momenta and metric tensors, of a similar structure, eventually to be contracted with the original numerator momenta of the integral. Hence, each coefficient c reads . . p σn 2 + all perms among ρ and σ indices where the permutations take care of the symmetry of numerator indices. The inversion coefficients d(n 1 , n 2 ) l,m,n i,j,k multiply tensor structures contracting indices in the numerator of the to-be-reduced integral. Such structures are generic for given n 1 and n 2 and do not depend on the specific numerator or integral topology. Therefore they can be computed once and recycled into other integrals with the same n 1 and n 2 . The inversion process to determine such coefficients can turn costly for high powers of numerators. Moreover, all solutions for the coefficients c may have in principle to be derived, though only the limited subset c(n 1 , n 2 ) 0,0,k is relevant for the computation. The index k is eventually attached to a term with the corresponding power of z 12 in the result.
In practice, I took an alternative route. I fixed the relevant inversion coefficients d(n 1 , n 2 ) l,m,n 0,0,k heuristically, by comparing a sufficient number of independent reduced tensor integrals with a direct IBP reduction including z 1 and z 2 scalar products. Since the inversion coefficients must not depend on the particular integral, just on n 1 and n 2 , I evaluated a test reduction on the simplest three-loop integral topology, the triple bubble, with a bunch of different numerator combinations. Then I compared the result with a FIRE6 reduction including z 1 and z 2 contractions. This is feasible and fast also for high momentum powers, thanks to the simplicity of the integral topology. The two sets of results must coincide, independently for each z 12 power. This in turn produces a linear system of equations (potentially with a high number of redundant constraints), which can be solved to compute the inversion coefficients. Moreover, for each z 12 power, these are expected to be rational functions of the dimension d. With some experimentation, a rough upper bound estimate of the maximal power of d in the numerator and denominator of such rational functions can be established. Let's say that these powers combined have upper bound n for given n 1 and n 2 . Then, one just needs to evaluate the numeric coefficients of such powers of d. For this, the relevant systems of equations can be evaluated at n independent numerical values of d (or a bit more for a check), where the system does not loose rank. This could be integers, or rationals, so that the final systems possess only numeric (rational) coefficients and can be inverted more rapidly. If this is the case, the IBP reduction part can also be performed directly with rational values of the dimension, instead of generic d, if that produces a faster process. One trial is sufficient for restricting a sufficient set of integrals to be IBP reduced, that supplies enough independent equations to invert the system. This eliminates the redundancy in the system, mentioned above for the subsequent reductions at fixed values of d.
This is the method that I applied and that I found the quickest for computations involving operators with spin 6 and higher. For larger values of spins, it becomes more efficient to evaluate the inversion coefficients at the relevant integer dimension d = 4 and then solve for corrections perturbatively (fixing the exact in d expressions takes too much time and is practically useless). The needed order at two loops for such en expansion is 2 , given the maximal order of the possible poles of the relevant integrals emerging from Feynman diagrams.
In a few situations, one can determine analytically the inversion coefficients, for generic values of some of their parameters. For instance, this is the case when only a single null momenta is present, i.e. n 2 = 0. The inversion coefficients only depend on the power of numerator momenta n 1 . A heuristic formula for the inversion coefficients can be derived for generic n 1 , in terms of combinatorial factors.
As said, some parts of this process are a bit heuristic. For instance, the choice of numerator momenta for deriving inversion coefficients, was highly redundant. I have not studied what a minimal selection of numerators giving rise to independent constraints would be. The same goes for the exact upper bound of d powers in the coefficients. In practice, this entails some efficiency loss, and there is room for possible optimizations. However, for the computation I am presenting, the method worked sufficiently fast and there were other, more pressing bottlenecks.
Whence all the relevant inversion coefficients are known for a numerator with powers n 1 and n 2 , the explicit tensor reduction of a generic integral can be carried out. This step generates a larger and larger number of scalar integrals, for higher and higher powers of numerator momenta. A first optimization consists in leveraging the further symmetries of the reduction owing to the particular numerators. In general, numerators contain repeated powers of the loop momenta k 1 , k 2 and k 3 (for the three-loop case) which introduce subsets of symmetric indices. This symmetry can be used, which reduces the number of independent contractions in (2.4). The output (containing scalar integrals yet to be reduced by IBP identities) can still be very bulky. At that stage the remaining scalar integrals are reduced, expanded in up to the required order (2 at two loops) and substituted in the expressions for tensor integrals.
An initial implementation of this algorithm in Mathematica turned out to be too slow starting at spin 10. Hence I performed most of the simplifications in FORM [32], which efficiently deals with large expressions. For reducing scalar integrals via IBP reductions I mostly used FIRE6. However, some reductions failed around the level of complexity of a three-point function with two spin-10 operators, in the non-planar topology. I used Mincer [33], for dealing with such cases. The process described above is rather roundabout, but it works far more efficiently than the approach in [20] and allows to push the computation to higher values of the spins. At spin 12 a few hundred thousands integrals were needed to be reduced. That required deployment on a small cluster of around 300 cores for completing the process in a reasonable amount of time.

Results and conclusions
In the following tables the results are resumed for the two-loop corrections to three-point functions of twist-two operators with two spinning ones of spins j 1 and j 2 and polarization l. Such structure constants are normalized by the two-point functions of the corresponding operators. Said another way, the relevant operators form an orthonormal basis with respect to two-point functions. Namely, the coefficients multiplying the space-time conformal structure (A.9) are 1 for two-point functions of the same operator, and 0 otherwise. Moreover, the ratio of each quantum correction C (2) j 1 ,j 2 ,l with the corresponding tree-level value C (0) j 1 ,j 2 ,l is reported. The two-loop corrections include a contribution proportional to ζ(3), which reads 24|S 1 (j 1 ) − S 1 (j 2 )| ζ(3). For brevity, I removed such a part from the results below, which are finally understood to represent the quantitiesC and possibly outperform the results presented here, from that angle. This would also constitute an important consistency check of the data presented here. A first test that I have implemented on my results is the gauge invariance of the underlying Feynman diagram expansion, which is inherited from [20]. Moreover the fact that the threepoint functions renormalize correctly and through the known anomalous dimensions of the relevant operators provides a further consistency check. These checks however limit the tests to the gauge dependent part and divergent terms of the computation, respectively. A further check on the finite part comes from the expected behaviour of the transcendental terms, though their coefficients are tied to the divergent term, hence not independent. In conclusion, there are several indications that the results are correct. Still, an additional independent computation would be desirable. The extension of the present computation to the next perturbative order, perhaps restricted to lower values of the spins, constitutes a possible interesting development, especially so since not much information on higher-point functions of protected operators is presently available at three loops.

Acknowledgments
This work was supported by Fondo de Investigación VIDCA 2020. I also acknowledge computing at the Magi cluster at ICFM, Universidad Austral de Chile and thank its administrator Simón Poblete for assistance and maintenance after several crashes. I thank the authors of [6] for a conversation which stimulated this work.

A Setting
I work in N = 4 SYM with SU (N ) gauge group. The results presented in this article, that is up to second order in perturbation theory, do not exhibit any non-planar correction. Therefore I use ubiquitously the 't Hooft coupling λ = g 2 N 16π 2 as the perturbative expansion parameter (g is the Yang-Mills coupling constant). In a conformal field theory such as N = 4 SYM, the form of three-point functions is fixed thanks to conformal symmetry. The focus of this note is on three-point functions with two spinning operators, of spins j 1 and j 2 , and a third scalar operator. Their three-point function is constrained to exhibit the general form where the invariants read and the following notation is used The vectors z 1 and z 2 are two null vectors (z 2 i = 0) onto which the Lorentz tensor structure of the operators is projected, in such a way to automatically form the symmetric traceless combinations required for the correct representations of the operators. Since there are two independent operators, two null vectors are used, which give rise to a non-trivial invariant z 12 ≡ z 1 · z 2 , parameterizing the various polarizations of the three-point function.
The coefficients C l are min(j 1 , j 2 ) + 1 independent structure constants, which are functions of the coupling constant g 2 and the rank of the gauge group N , through the t' Hooft coupling λ since, as mentioned, no color-subleading corrections appear up to two loops. The perturbative expansion of structure constants is thus As mentioned in the main text, a space-time integral is considered of (A.1), which streamlines the extraction of the structure constant. The (dimensionally regulated in d = 4 − 2 ) integral of the right-hand-side of (A.1) can be computed for general spins: The structure constants obtained via this method are not normalized in the standard form yet. The normalization is fixed in such a way that the two-point function of the relevant operators are orthonormal, which is a canonical basis of operators in a conformal field theory. In this note I consider three-point functions of twist-two operators O j of spin j, whose bare expression reads schematically where X is one of the complex scalars of N = 4 SYM and D are covariant derivatives. The correct spin representation can be attained contracting the covariant derivatives on the light-cone via null vectors z i , giving rise tô O j ≡ Tr(z 1 ·D) k X (z 1 ·D) j−k X)+. . . ,Ô k (x) ≡ Tr((z 2 ·D) kX (z 2 ·D) j−kX )+. . . (A.8) The scalar fields are chosen with such flavours that the three-point function is nonvanishing. Such operators mix non-trivially under renormalization but in a closed form, that is they mix with operators of the same sl(2) sector. For twist-two operators mixing only occurs with descendants ∂ k O j−k of lower spin operators with the same total spin j. For each even value of the spin there is a single primary operator (odd spin operators are all descendants).
After re-normalization, an eigenbasis of operators for the dilatation operator is selected, whose operators are finally the conformal ones appearing in (A.1). Their two-point function is fixed thanks to conformal symmetry Ô j (x 1 )Ô k (x 2 ) = C(g 2 , N ) δ jk I j 12 |x 12 | 2∆ (A. 9) where ∆ is the conformal dimension of the operator, including quantum corrections. The tensor structure is encapsulated in the quantity I. The null vectors z 1 and z 2 are in principle distinct for the two operators, but they may be set equal to simplify the computation, without affecting the computation of the coefficient C. The latter is an overall normalization, depending on the perturbative definition of the operators themselves. This is chosen canonically in such a way that the operators form an orthonormal basis. In dimensional regularization this includes subleading in corrections as well, to the required order. The structure constants for such a basis of operators are those appearing in the results.