Non-planar data of $\mathcal N=4$ SYM

The four-point function of length-two half-BPS operators in $\mathcal{N}=4$ SYM receives non-planar corrections starting at four loops. Previous work relied on the analysis of symmetries and logarithmic divergences to fix the integrand up to four constants. In this work, we compute those undetermined coefficients and fix the integrand completely by using the reformulation of $\mathcal{N}=4$ SYM in twistor space. The final integrand can be written as a combination of finite conformal integrals and we have used the method of asymptotic expansions to extract non-planar anomalous dimensions and structure constants for twist-two operators up to spin eight. Some of the results were already know in the literature and we have found agreement with them.


Introduction
The correlation functions of length-two half-BPS operators (also known as 20 operators) in N = 4 SYM have been studied extensively in the literature both at weak [1][2][3][4][5] and strong coupling [6][7][8][9][10][11][12][13][14][15][16][17][18][19]. It is well known that the two-and three-point functions of those operators are protected by supersymmetry [7,20] while higher-point functions receive non-trivial corrections. Each loop order of the correlation function can be computed through the Lagrangian insertion procedure, and because the 20 operator and the Lagrangian belong to the same supermultiplet, the integrand exhibits a hidden permutation symmetry [3]. Imposing also conformal symmetry and restrictions from logarithmic divergences of the correlator in diverse OPE limits, the planar integrand of the four-point correlator has been fixed up to ten loops [4,5]. It is however important to note that knowledge of the integrated correlator is still quite incomplete since conformal integrals are generally not known starting from four loops. Meanwhile, at the non-planar level much less is known. It is easy to see that at one and two loops the non-planar corrections are zero since it is not possible to draw any non-planar Feynman graph at these loop orders. However, the vanishing of the three-loop correction is non-trivial as it results from cancellations of different terms. Finally, at four loops the integrand is known to be a linear combination of four polynomials with constant coefficients. One of the main results of this work is the determination of those four undetermined coefficients. We expect that the methods described below can be adapted to the calculation of other correlators as well.
One of the motivations for computing the four-loop non-planar integrand is that it allows to extract non-planar OPE data by taking a double coincidence limit. In this work we have computed both anomalous dimensions and structure constants for twist-two operators up to spin eight. Even at the planar level, where the tools of integrability are more developed, this method is the best way we currently have for obtaining structure constants of unprotected operators at a high loop order [21][22][23][24]. In contrast, it was only recently that a direct two-loop perturbative computation of planar structure constants was performed [25,26]. The double OPE limit of each conformal integral can be taken by using the method of asymptotic expansions described for example in [21,27]. With this method the integration domains are split into distinct regions, which correspond to the different scales of the problem. Effectively, the conformal four-point integrals can be rewritten in terms of two-point integrals, which are much more tractable. Note that the pseudo-conformal integrals arising in high-loop integrands can also be approximated with this method, as long as one does not assume dependence on cross-ratios. This is however not relevant to this work, as we can write the four-point function in terms of convergent integrals only.
The non-planar anomalous dimensions of twist-two operators up to spin six were also computed in the series of papers [28][29][30] by a direct Feynman diagrammatic calculation and we have found agreement with them 1 . It would be very interesting to find a closed expression for the non-planar anomalous dimensions of twist-two operators for any spin. We expect the result to be given in terms of harmonic sums and Riemann zeta values and to obey the principle of uniform transcendentality. Unfortunately the data obtained in this work is not enough to fix the expression for general spin even if one restricts to a smaller basis consisting only of binomial harmonic sums [31]. The knowledge of the non-planar anomalous dimension for generic spin would allow to compute the non-planar cusp anomalous dimension analytically and also take the BFKL limit. These results are important for understanding non-planar integrability or possible formulations of a non-planar quantum spectral curve [32], see [33] for progress in this direction. One should also stress that the non-planar cusp anomalous dimension was computed numerically by studying Sudakov form-factors with a suitable rewriting in terms of uniformly transcendental integrals [34][35][36].
Another motivation for computing non-planar structure constants is to further the understanding of non-planar integrability. Three-point functions can be computed at the planar level as a product of two integrable hexagon form-factors [37], and it was later understood that higher-point functions can also be decomposed into a weighted product of hexagon form-factors, both at the planar and non-planar level [38][39][40][41][42][43]. The integrability setup was tested at two loops for long operators and at one loop for short operators such as the 20 . In the case of four-point functions of length-two operators, all the non-planar corrections in the integrability setup come from the stractification procedure described in [41]. More specifically, one embeds the tree-level planar graphs in higher genus surfaces and properly subtracts possible boundary terms. It would be interesting to test the stratification procedure at higher loops (at the moment four loops seems a difficult task) or at least understand why non-trivial non-planar corrections to this correlator first show up at four loops. This paper is organized as follows. In the remaining part of the introduction, we review what is known in the literature about four-point functions of 20 operators up to four loops. In section 2 we introduce the reformulation of N = 4 SYM in twistor space and explain how to compute correlation functions in that framework. We then describe the strategy used to fix the four-loop non-planar integrand. The OPE analysis of the correlator is performed in section 3, where we present the non-planar OPE data of twist-two low-spin operators. We conclude in section 4 and refer the reader to the appendices for conventions and examples of analytic computations using twistors.

Four-point function of 20 operators
We will now review some well-known results about correlation functions in N = 4 SYM. We will focus our attention on the four-point function of length-two half-BPS operators, which can be written as with Φ i the six real scalars of the theory and y i a null polarization vector which projects the operator into the symmetric traceless representation. The correlator admits a double expansion in the effective coupling constant a = g 2 YM N c /(4π 2 ) and in the number of colours N c Starting at one loop, N = 4 superconformal symmetry [44] implies that at any order in N c the correlator G (g, ) 4 factorizes in the following way where all depence on the polarization vectors of the external operators is encoded in the R factor Notice that all dynamical information is contained in the functions F (g, ) , which crucially multiply all six R-symmetry structures. We can therefore work with a particular choice of polarizations where only y 12 and y 34 are non-zero, so that a single term in the R factor survives. This will significantly reduce the number of graphs to be computed, as explained in more detail later in the next section, but one can unambiguously reconstruct the correlator for generic polarizations.
Loop corrections for the correlator can be obtained by the Lagrangian insertion procedure, where the integrand of the -loop four-point function is viewed as a Born-level (4 + )-point function. We can then rewrite the dynamical function as where the integrand f (g, ) carries conformal weight +4 in all external and internal points. Furthermore, an analysis of the possible OPE singularities indicates that the integrand is a rational function which diverges at most as a simple pole in the coincidence limit x ij → 0, which implies that we can rewrite it as Finally, P (g,l) is a linear combination of polynomials which have the following properties: 1. It is homogeneous in x 2 ij .
3. It is invariant under the permutation of all its arguments, i.e. under the group S 4+ . Property 2 follows from conformal symmetry, while property 3 reflects a hidden permutation symmetry, which follows from the fact that the Lagrangian operator is in the same supermultiplet of the external 20 operators.
At each loop order there is a finite number of polynomials P ( ) i that satisfy the properties listed above and each P (g,l) is a linear combination of those with constant coefficients. Note that the properties above are independent of g, so the basis P ( ) i which solves these constraints can be used to construct the numerator of the integrand at any order of the genus expansion. At one-, two-, three-and four-loops there are 1, 1, 4 and 32 independent polynomials respectively. We will not write them explicitly in this work so we refer the reader to reference [4]. In order to fix the integrand completely, one only needs to find the coefficient multiplying each polynomial. A powerful method to fix these coefficients is to study the asymptotic behaviour of the correlator either in the Euclidean double shortdistance limit, where both x 1 → x 2 and x 3 → x 4 , or the Minkowski light-cone limit where x 2 12 , x 2 23 , x 2 34 , x 2 41 → 0. In these limits the logarithm of the correlator must develop soft logarithmic singularities, which imposes strong constraints on the coefficients. These constraints, together with the conformal Gram determinant relations 2 were powerful enough to fix the planar four-loop result and to reduce the non-planar corrections at four loops to only four unknown coefficients.
Before writing down the form of the non-planar integrand, let us clarify the classification of the polynomials regarding their planarity. For each homogeneous polynomial P that of its associated f -graph and it can only contribute to the integrand P (g, ) if the genus obtained does not exceed g.
The analysis of the non-planar integrand in [4] showed that corrections to G 4 first appear at four loops, but the constraints were not sufficient to fix it uniquely. At genus one, the integrand is given up to four undetermined coefficients Each term is a linear combination of the 32 four-loop polynomials P (4) i (see equations (5.9) and (C.1) of [4] for definitions) with where the short-hand notation 0 n corresponds to a list of n zeros.
One of the results of this work is the determination of the coefficients c i . We obtained c 1 = c 2 = c 3 = 0 and c 4 = −6, and in that way we fixed the non-planar integrand at four loops completely. The method we have used relies on the reformulation of N = 4 SYM in twistor space, which is the subject of the next section.

Twistors
In this section, we first review how to compute correlation functions of the stress-tensor supermultiplet in N = 4 SYM using twistor space, see [45] for further details. One of the advantages of this formalism is that each Feynman diagram in twistor space has manifest N = 4 superconformal symmetry apart from some reference twistor. Then we explain how the four-loop non-planar calculation was performed for a particular polarization of the external operators. The necessary graphs were generated with the open source program Sage [46].

N =4 SYM in twistor space
The supertwistor space [47,48] is the complex projective superspace CP 3|4 . An element Z A of this space has four bosonic and four fermionic coordinates and it is defined up to the equivalence relation Z A ∼ cZ A , with c ∈ C * . These variables are parametrized in the following way with α,α = 1, 2 and χ a , a = 1, . . . , 4, the fermionic coordinates. A nice property of these variables is that they transform linearly under the action of all generators of the complexified super conformal group SL(4|4; C), see for example [49,50] for the explicit form of the generators. In addition, these variables can be related to the usual superspace variables (x αα , θ aα ,θα a ). We are interested here only in the chiral superspace, i. e. we are going to set allθα a to zero, and in this case These are called incidence relations and they map a point in chiral superspace to a line in supertwistor space.
The first relation in (11) can be understood as follows. For simplicity, let us consider the bosonic components of the supertwistors, in which case the complexified conformal group is SL(4; C). The supersymmetric case is a simple generalization. The twistors Z I = (λ α , µα), with I = 1, . . . , 4, transform in the fundamental representation of this bosonic group. We can define a null antisymmetric tensor X IJ as (see the Appendix A for conventions) These tensors are also homogeneous, with X IJ ∼ cX IJ , and the set of null rays is in correspondence with the original 4d spacetime coordinates x αβ . This identification is known as the embedding formalism. Because of the null condition given in (12), the matrix X IJ has rank two and it can be written in terms of two twistors as As mentioned before, this implies that a spacetime point is mapped to a line in twistor space. The line connects the two twistors {Z I 1 , Z J 2 } which have linearly independent values for λ α and the component µα given by by (11).
It is possible to reformulate N = 4 SYM in supertwistor space and in that way we gain an alternative method for computing correlation functions of the stress-tensor multiplet. The fields of N = 4 SYM sit inside a one-form superfield A living in supertwistor space. Accordingly, the action can be written as a function of this superfield in the following way (see [51,52] for details) In order to perform calculations, it is convenient to choose a gauge in which a component of the superfield A vanishes in the direction of a reference twistor Z , so that the kinetic term becomes quadratic and the interaction term simplifies to where the superfields are integrated along a line in twistor space with Z α satisfying the incidence relations (11) for different λ α . The measure is Dσ = σ dσ and the bracket notation stands for σ i σ j = αβ σ α i σ β j . While non-trivial, it has been shown that any physical quantity is independent of the reference twistor Z .
The spacetime equations of motion can be obtained from the action above by expanding the superfield A in components In the formula above the fields on the right-hand side depend on the bosonic twistors Z,Z parametrizing a line. Moreover {a, a } are the two gluon helicity states, {ψ a , ψ d } are the gluinos and φ ab are the six scalars. It is very important to notice that the action above is chiral and it contains the topological term iFF , where F is the field strength. While this term is not important in perturbation theory because it is a total derivative, it is going to contribute to the integrand we want to compute by introducing a term proportional to the spacetime tensor µνρλ . These terms have the wrong parity and they integrate to zero [53].

Correlation functions of the stress-tensor multiplet
Our aim is to compute the correlation function of four 20 operators defined in (1). This operator is the lowest component of the stress-tensor supermultiplet T , whose top component is the Lagrangian. The fact that this is a short multiplet implies that T depend only in half of the odd variables θ aα ,θα a .
In order to define the relevant fermionic degrees of freedom, it is convenient to introduce the auxiliary harmonic variables u b a ≡ (u +b a , u −b a ) which parametrize the coset SU (4)/(SU (2) × SU (2) × U (1)). The indices a, b are fundamental indices of SU (4), while b, b are fundamental indices of SU (2) and SU (2) respectively, with the signs indicating the U (1) charge. These variables and their complex conjugates satisfy several unitary and completeness conditions which follow because u b a is in SU (4). The harmonic variables allow us to write manifestly SU (4) invariant expressions. If one defines then one can decompose every θ a α as whereū are the complex conjugates of the harmonic variables u. The stress-tensor superfield T (x, θ + ,θ − , u) depends on half of the odd variables, both chiral and anti-chiral, but it is useful to focus on the chiral part of the multiplet, where allθ − vanish, so that we have and we have omitted the other powers of θ + in the expansion. The field L(x) is the chiral Lagrangian and O ++++ (x) = Tr(φ ++ φ ++ ) is a representation of the half-BPS operator defined in (1) with φ ++ = φ ab u +b a u +c b bc . The connexion between u a b and y i is made with the following particular parametrisation of the harmonic variables with y 2 = −y b a y a b /2 and the indices a, a are raised and lowered as usual with the epsilon tensors.
The correlation functions G n = T (1) . . . T (n) have a series expansion in θ + i as a consequence of (20), and we can extract the correlation function of four 20 operators by computing G 4 and reading its lowest component, or equivalently, by sending all θ + i to zero. For four-and higher-point functions there is a dependence on the coupling a, which can be made precise through the Lagrangian insertion method or, more generally, On the other hand, from equation (14) we can also derive the following insertion formula which hints at the following representation of the stress-tensor superfield in twistor space In order to extract the -loop four-point function, we will then have to compute the following (4 + )-point correlator in twistor space Looking at equation (14) we see that each L int comes with at least two superfield insertions along the twistor line. Moreover, according to equation ( The computation of G n at tree level is made by summing all relevant diagrams and for each graph in twistor space we have to use the Feynman rules which were derived in [45]: 1. A propagator connecting vertices i and j provides a factor d ij = y 2 ij /x 2 ij and a colour delta function δ a i a j , 2. A bivalent vertex contributes a colour factor Tr(T a 1 T a 2 ) = δ a 1 a 2 , 3. Higher-valence vertices are associated with the factor R i with T a the generators of the gauge group and the R factor defined by The delta functions are fermionic and therefore, by construction, R i j 1 j 2 ···j k has Grassmann degree 2k − 4. The σ ij originate from the integrations in equation (15), which are localised by the twistor propagators and become where Z j,1 and Z j,2 are the bosonic components of the twistors parametrizing a line which corresponds to the spacetime point x µ j . Finally, by setting the fermionic components of the auxiliary supertwistor to zero we have The R factors defined in (27) have several important properties and satisfy some identities which can be found in [45]. In what follows we will only need two of these identities. First, since the numerator of (27) does not depend on the ordering of {j 1 , . . . , j k } then the effect of a permutation ρ on the indices is simply (1) .
In this way we can rewrite R factors in a canonical way and reduce the number of fermionic integrations we need to perform. Second, a multi-index R i j 1 j 2 ···j k can always be factorized as follows which implies that each diagram can be rewritten in terms of a fundamental building block R i j 1 j 2 j 3 , which takes the following form after the integration over the θ − Once we sum all relevant diagrams, we obtain the component of G 4+ with fermionic degree 4 . And since we want to perform the d 4 θ + integrations at the Lagrangian insertions to obtain the loop-level four-point function, effectively we need to send the θ + i at the external points to zero. This implies that the A a ij defined in (29) can only give a non-zero contribution Figure 1: The skeleton graphs used for the one-loop computation of a four-point function.
The ribbon graphs are obtained from these by adding colours traces. Only the first graph contributes to the particular choice of polarization, and the only configuration which does not vanish is when the middle vertex corresponds to the integrated Lagrangian insertion. One can easily compute it analytically with equation (65).
if at least one of the indices corresponds to an internal point, which in turn means that we only have dependence on y −1 ij if at least one of the indices is from an integrated point. Despite the obvious simplicity of this statement, it does imply that if there is a propagator between external points k and l, the resulting y 2 kl factor can never be cancelled as an effect of the fermionic integrations. Consequently, if we choose a particular polarization where y kl vanishes, then we can neglect all diagrams which contain a propagator between those two points. This is a great simplification because the factorized form of the correlator in equation (3) allows us to select external polarizations such that only y 12 and y 34 are different from zero, thus greatly reducing the number of twistor space diagrams we need to evaluate. In the next subsection, we explain how to compute a four-point function with this assumption. It is possible to perform several intermediate analytical computations as well and we give examples in the Appendix B.

Four-loop four-point function
The first step in the calculation of the four-point function is to generate all relevant graphs. As discussed previously, at loops we need graphs with 4 + vertices and 4 + 2 propagators. These graphs can be easily constructed at lower loops but the number of graphs increases very fast with the loop order, so we used the program Sage [46] to generate them. At one and two loops, the skeleton graphs are shown in figures 1 and 2 respectively, while the number of skeleton graphs up to four loops is shown in Table 1. At one loop one can draw graphs loops one two three four # of skeleton graphs 3 11 63 513  Figure 2: The skeleton graphs at two loops generated by Sage. The final set of ribbon graphs is obtained from these by adding the colour factors and summing over all inequivalent assignments of external and internal points to the vertices.
The next step in the computation is to generate the ribbon graphs from the skeleton graphs obtained with Sage. In other words, each vertex of valence v is assigned both a colour trace Tr(T j 1 . . . T jv ) and an R i j 1 j 2 ···jv factor as defined in (27), and each propagator supplies an additional factor of d ij . In general, each skeleton graph obtained with Sage leads to a large number of ribbon graphs. This happens because any inequivalent permutation of the indices {j 1 , . . . , j v } appearing in the colour traces and R factors gives rise to a different ribbon graph. More precisely, a skeleton graph with n vertices of valences {v 1 , . . . , v n } produces n i=1 (v i − 1)! ribbon graphs, corresponding to the non-cyclic permutations at each vertex.
The graph propagators also provide colour delta functions so that the indices in the colour traces of the vertices are fully contracted. Effectively the color structure of the ribbon graph simplifies to a polynomial in N c through successive application of the fission and fusion rules Tr (T a BT a C) = Tr(B) Tr(C) − γ N c Tr(BC) , where B and C are arbitrary matrices and the cases γ = 0, 1 correspond to U (N c ) and SU (N c ) respectively. We have checked that our four-loop four-point correlator is independent of γ. This seems a bit surprising at the non-planar loop level, but it is true for lenght-two operators. The U (1) part of the gauge group is free and the external operators only give a trivial colour factor Tr(T a 1 T a 2 ) = δ a 1 a 2 in this case. 3 For correlators involving half-BPS operators of higher weight, we expect that generically the result will depend on γ, see for example Appendix A of [41]. The evaluation of ribbon graphs can be further simplified by using the property of the R factor given in (30), which allows us to rewrite all R factors in a chosen canonical order, up to factors of σ ij σ ik . This implies that all ribbon graphs which originate from the same skeleton graph will produce the same canonical R factors and differ only by a simple prefactor. In conclusion, to each skeleton graph we associate a canonical ribbon graph. With this knowledge, and without performing any fermionic integration, we could already rederive the well-known result that the correlator of four 20 operators can only receive a non-planar correction at four loops. This happens because for lower loop orders all skeleton graphs produce the simple colour factor of (N 2 c − 1). Finally, from the canonical ribbon graphs we can generate the final set of diagrams by associating vertices with external and internal positions in all inequivalent ways. Naively this would generate (4 + )! graphs from each canonical ribbon graph contributing to the -loop correlator. However, most graphs possess a non-trivial automorphism group, so that the number of inequivalent permutations is effectively much smaller. Moreover, once we have assigned positions to the vertices of the graph we can check if it is possible to produce an (θ + ) 4 factor for each of the internal points, and it turns out that this test considerably decreases the number of allowed graphs. The number of final graphs for a generic polarization of the external operators is shown in the first line of table 2. Luckily, as explained in subsection 2.2, we can perform the computation for a particular choice of polarizations and reconstruct the final result unambiguously. If polarizations are such that only y 12 and y 34 are non-vanishing, then any graph with external-to-external propagators other than d 12 or d 34 is identically zero. This drastically reduces the number of diagrams to compute, as can be seen in table  2. loops one two three four # of graphs generic polarization 45 1417 75141 6019618 # of graphs particular polarization 1 73 7939 715350 Table 2: The final number of diagrams which is obtained from canonical ribbon graphs by associating vertices to internal and external positions. In the generic polarization all y ij are non-vanishing, while in the particular polarization any scalar product other than y 12 and y 34 is zero. The four-loop correlator in this paper was evaluated in the particular configuration, for which the number of graphs is greatly reduced.
After generating all the graphs and their prefactors, the final step in the calculation is to replace the R factors by their expression (27) and perform the fermionic integrations in θ + at the Lagrangian insertions. At four loops, all graphs become a product of eight R i j 1 j 2 j 3 factors, thus producing many terms with the correct number of θ + i . An analytical computation is then very hard, especially because the dependence on the auxiliary twistor Z only disappears after summing many graphs, see [45] for details. Therefore, in this work we have performed the computation numerically by giving integer values to all components of the polarization and position vectors. While it is not necessary to restrict to integer numbers, we found that was a practical way of avoiding numerical fluctuations and errors.
The planar integrand had already been fixed in [4], and we successfully reproduced their result up to four loops with our diagrammatic expansion in twistor space. This important comparison provides a good cross-check of our implementation. At the non-planar level the integrand was written as a linear combination of four polynomials with undetermined coefficients. By computing each twistor diagram for twelve different sets of numerical values we produced an overcomplete system of equations and were able to fix those coefficients. Each numerical evaluation of the graphs took approximately 3 days in a single computer with 20 cores. The result is given in the following section, where we cross-check it further against some available non-planar data.
If we consider configurations of points living in four dimensions, then there is a technical issue that arises. As discussed at the end of subsection 2.1, the twistor action has the topological term iFF and so the evaluation of the graphs generates terms involving the tensor µνρλ . These are spurious contributions at the level of the integrand and they must be absent in the final result, i.e. they necessary multiply functions that integrate to zero as iFF is a total derivative and it cannot give any perturbative contribution. Note that in Lorentzian signature the terms with an odd number of epsilon tensors will be imaginary, which means that we can single out their contribution in our calculations. In order to fix them we found the complete basis of pseudoscalar conformal integrands with extended permutation symmetry. At first there seem to be many structures one can form, but when we implement permutation symmetry we must include minus signs to compensate the antisymmetry of the epsilon tensor, and at the end there are only 4 such polynomials. Furthermore, two of those correspond to pseudoscalar conformal Gram polynomials, leaving only two degrees of freedom. We obtained an imaginary non-planar component for five of our numerical data points, thus producing an overcomplete system of equations. Effectively, the result is a shift to the polynomial P (1,4) bỹ where we defined an c tensor by using six-dimensional embedding vectors (1, By constructing the pseudoscalar polynomials with this six-dimensional epsilon tensor we guarantee that the polynomial in (35) is conformal. In addition, recall that a product of ordinary epsilon tensors can always be reduced to a sum of products of Kronecker deltas Thus, if spurious terms are generate in this way, the conformal and permutation symmetries ensure they can be written as a linear combination of the 32 four-loop polynomials from [4]. We have to be confident that our numerical result is not polluted by this kind of spurious terms. As mentioned before, the terms generated with tensors have to integrate to zero, but we believe there is no such combination apart from the conformal Gram polynomials. Thus we are confident that our twelve data points still give an overcomplete set of equations if we consider an extended basis which have these additional polynomials that integrate to zero. Another check of our result can be made by noticing that the bosonic twistors, which are written in terms of λ α and xα β due to the incidence relations (11), can only contribute to the R-factors via the σ α ij defined in (28). Thus, all the matrices σ µ αα are always contracted with a vector x µ i . This implies that if a graph generates µνρλ through simplifications of products of σ µ αα , the resulting epsilon tensor will always be fully contracted with a set of x µ i 's and therefore vanish identically if all points are in a three-dimensional subspace. Thus these possible spurious terms belong to the space spanned by the 6 exclusively three-dimensional conformal Gram polynomials. We have fit our numerical expressions to an extended basis consisting of the four non-planar polynomials defined in (9) together with the 3d conformal Gram polynomials, but our result remained the same. Finally, we performed a few additional numerical tests by checking invariance under a change of polarizations for the internal vertices and independence of the result on the auxiliary twistor.

OPE analysis
The diagrammatic computation of the previous section fixed the non-planar polynomial to be with Q 4 (x i ) defined in (8). However, some of the terms in the polynomial lead to pseudoconformal integrals. The weight at each integrated point is +4, but those integrals are divergent in four dimensions and once we introduce dimensional regularization they lose their conformal properties. Since it is more cumbersome to deal with pseudo-conformal integrals, we chose to add a conformal Gram polynomial to the integrand and rewrite it as with P The conformal Gram polynomials parametrize a three-dimensional subspace of the allowed polynomials, and we used two of its degrees of freedom to ensure that none of the pseudoconformal integrals contributes in (40). The last parameter is chosen in a way that eliminates some of the more difficult conformal integrals.
While the twistor space method allowed us to obtain the integrand of the correlation function, we are still faced with the integration of the Lagrangian insertions. Conformal symmetry reduces the complexity of the problem by restricting to functions of two crossratios only The four-loop ladder diagram is known exactly [54], and recently the method of differential equations was used to fix a different topology [55], but in general the evaluation of four-loop four-point integrals is a daunting task.
Therefore, in this work we will focus on the Euclidean coincidence limit, which is tractable. We use the freedom of conformal symmetry to send x 4 to infinity, so that effectively we deal only with three-point integrals, and then we let x 1 approach x 2 . In that case the cross-ratios become Each of the conformal integrals has now two distinct scales |x 12 | |x 13 |, which means that we can approximate the integrals with the method of asymptotic expansions. The idea is that for each integration variable x i we can now divide the integration domain into two regions, one where x 1i is of the order of x 12 , and another where it is of the order of x 13 . In that way, an -loop conformal integral will split into 2 terms, corresponding to different distributions of the integration variables in the two regions. For example, if x 1i is of the order of x 12 and x 1j of the order of x 13 , we can Taylor expand the propagator This Taylor expansion is only convergent for |x 1i | ≤ |x 1j |, but we extend the integration domain in this region to the whole space. That makes the integrals strictly divergent, which can be resolved by introducing dimensional regularization and requiring that scaleless integrals vanish [56]. In order to find the lowest order in the small u expansion we can ignore the x 2 1i terms from the numerator of (43), and higher powers in the Taylor expansion will contribute to higher orders in the small Y expansion.
Looking at equation (43) we realize that the two regions are effectively disentangled. This indicates that a term with a {k, − k} split of the integration variables into the {x 12 , x 13 } regions will lead to a product of k-and ( − k)-loop two-point integrals, with external points x 1 and x 2 , or x 1 and x 3 , respectively. This is a considerable simplification because we are able to approximate the conformal integral with single-scale integrals. The two-point integrals will generically have numerators with tensor structures, but we can follow the strategy of [21] to reduce them to scalar integrals. Once that is accomplished we use LiteRed [57] and FIRE [58] to implement IBP identities and reduce to a basis of simpler master integrals, which were obtained in [59].
For simplicity we focused on a parametrization of the integrand which involved solely convergent integrals, but that was not strictly necessary. The method of asymptotic expansions can also be applied to pseudo-conformal integrals and the only caveat is that one must disregard conformal symmetry. More specifically, different terms in the integrand lead to integrals which are the same up to a permutation of the external points. When they are convergent one can rely on the invariance of the cross-ratios under the group of double-transpositions {id, (12)(34), (13)(24), (14)(23)} (44) to relate them. However, the pseudo-integrals are divergent and in dimensional regularization they are not simply a function of cross ratios, which means that we must perform the asymptotic expansion independently for each permutation of the external points.
where σ and τ are the R-symmetry cross-ratios, Y n,m is the R-symmetry block for the exchange of an operator in the SU (4) representation [n − m, 2m, n − m], and g ∆,l is the conformal block for an operator of dimension ∆ and spin l. In the Euclidean OPE limit that was described above the conformal block simplifies to which means that only the lowest-twist non-protected operators contribute to the leading u behaviour of the four-loop correlator. Furthermore, a suitable choice of the external polarization vectors y i allows us to single out the [0, 2, 0] representation. Since there is a single twist-two operator in the 20 representation for each spin, we are able to extract all their OPE coefficients and anomalous dimensions.
The OPE data is written as a double expansion on the genus and coupling constant Note that with the twistor space calculation we reconstruct the full N c dependence of the four-loop four-point function, which shows that the genus expansion truncates at the first non-planar order. The non-planar anomalous dimensions are therefore These results match the perturbative computation of Velizhanin [28][29][30] for spin 2, 4 and 6. 4 Meanwhile the expression for spin 8 is new and its leading transcendental piece matches Velizhanin's conjecture for general spin with S 1 (l) the harmonic sum. This seems to imply the behaviour log 2 (l) for large spin, which is in contradiction wih the expected one. However, this is only part of the result and cancellations can occur at large spin. We also extracted the non-planar correction to the OPE coefficients, all of which are novel results, The transcendental structure of the spin 2 structure constant is quite interesting due to the absence of rational and ζ 3 terms, but we do not currently understand why that happens. It is possible to derive a closed expression for the ζ 7 part of the OPE coefficient, which is given by 5 α with the tree-level planar coefficient given by α Finally, we were able to perform the asymptotic expansion of the integrand up to spin 12, but the IBP reduction of the resulting two-point integrals was not possible. The data we have available is not sufficient to reconstruct an expression for generic spin in terms of harmonic sums, which means that we cannot access the large spin limit of the anomalous dimensions. Perhaps the method of differential equations can be used to find a higher-order expansion of the conformal integrals in a more efficient way.

Conclusion
In this work, we have fixed the four-loop non-planar integrand of the four-point function of length-two half-BPS operators by a direct computation in twistor space. We performed the Grassmann integrations for numeric vales of the position and polarization vectors and fit the results obtained against a polynomial ansatz with four unknown coefficients. In this formalism each individual graph preserves N = 4 superconformal symmetry, apart from the reference twistor. The calculation was done with a particular choice of external polarizations in order to reduce the number of graphs, but the general result can be unambiguously reconstructed due to the factorized dependence on the polarizations in equation (4). In principle, one can also use our method to compute higher-point and higher-loop integrands, both at the planar and non-planar level. However, the number of diagrams can grow a lot in those cases and the complexity of the fermionic integrations can also increase considerably. In our case, the integration took an average of three days on 20 cores for each set of numerical values. In order to compute more complicated correlators it is perhaps necessary to implement a method which mixes numerical and analytical methods. More specifically, one can first perform some of the fermionic integrations analytically, as in the examples of Appendix B, and then complete the calculation numerically.
Let us now stress a technical detail of our method that could also show up for other correlators. The twistor action from equation (14) is chiral and contains the topological term iFF when expanded in components. This implies that terms with µνρλ can be generated in twistor space calculations. In this work we were able to write down an ansatz for these terms and found the contribution given in (35). This terms must integrate to zero and they do not appear in the final integrand, but in principle it can be difficult to isolate the epsilon terms from numerical calculations. However, when working in Lorentzian signature, the terms with an odd number of epsilon tensors have an imaginary contribution and they can be easily isolated. On the other hand, terms with an even number of epsilons give a real contribution which can be rewritten in the usual polynomial basis. Notice that these terms are identically zero when we restrict to a three-dimensial subspace, which implies that this possible spurious contribution can be written in terms of three-dimensional conformal Gram polynomials. A careful analysis using an extended basis showed that our result is not contaminated by this type of terms.
In this work we focused on the non-planar corrections to the four-point function of 20 operators, which start at four loops. There are also results in the literature for correlation functions of length-k half-BPS operators and for generic k the non-trivial non-planar corrections usually start at lower loops. At the planar level, all four-point functions of such operators are know up to five loops [60,61] while at the non-planar level correlators of four length-k operators are known up to two loops only [62,63]. It would be very interesting to compute non-planar corrections to more general four-point functions of half-BPS operators at two and higher loops. In that way we would further test the integrability approach to non-planar correlators developed in [40,41], and it could also help understand why nonplanar corrections show up at specific loop orders in this approach. Maybe some of these computations can be done using the twistor reformulation of N = 4 SYM used in this paper together wih the prescription for composite operators in twistor space of [64][65][66].
Finally, it would be extremely nice to obtain a closed-form expression for the non-planar anomalous dimension and structure constants of twist operators. This would allow to extract the non-planar cusp anomalous dimension and it would give valuable data for an integrability approach to the non-planar spectrum. We have made an OPE analysis of the four-point function and we obtained a few data points. However, that data was not sufficient to completely fix the anomalous dimension at generic spin. It seems hard to push the OPE expansion further with the method of asymptotic expansions and so a new strategy is very likely needed. We hope to return to this point in the future.
for valuable comments on the manuscript. T.F would like to thank the warm hospitality of the Arizona State University where this work was finished. This work was supported by the Serrapilheira Institute (grant number Serra-1812-26900). R.P. is supported by SFI grant 15/CDA/3472.

A Conventions
In this work, we have used the following conventions to raise and lower SU (2) indices Similarly, for the y a b variables introduced in (21) The epsilon tensors are defined with 12 = 12 = 1, so that they obey ab ac = δ c b .
For both the spacetime and R-symmetry matrices we use the following short-hand notation Using the properties of the Pauli matrices it is possible to show that xα α y α α = −2 x · y , xα α xβ α = x 2 αβ , We can manipulate the last equation above to see that Analogously, one has Concerning the integrations of the Grassmann variables θ ±a i,α , we use the convention This implies that θ ±1 1 θ ±2 1 θ ±1 2 θ ±2 2 integrates to one.

B Examples of fermionic integrations
In this work we have computed all necessary diagrams numerically, i.e we have set all components of the position and polarization vectors to integer values. Nevertheless, for a given diagram it might be possible to perform some or even all the fermionic integrations analytically. Eventually one can combine the two methods in an efficient way and reduce the complexity of many graphs. Here, we give some examples of integrations that can be performed easily. First, consider the θ + i integration when the internal point i is a bivalent vertex. In that case the Grassmann variables can be found in the product of R factors from the adjacent vertices j and k R j a 1 ...ami R k b 1 ...bni .
From equation (27) we can see that the relevant terms for the integration of θ + i come solely from which means that the only term with four θ + i is Therefore, the fermionic integration gives σ ik σ ij 2 y 2 ki y 2 ji σ jam σ ja 1 σ jam σ ji σ ji σ ja 1 σ kbn σ kb 1 σ kbn σ ki σ ki σ kb 1 .
Another example of a simple θ + i integration consists of the internal point sitting at a trivalent vertex connected to two bivalent vertices j and l and one higher-valence denoted by k. In that case the necessary fermionic variables are provided by the product of two R factors which integrates to d 4 θ + i R i jkl R k a 1 ...ani = R k a 1 ...an y 2 jl σ ij σ ik σ ik σ il y 2 ij y 2 ik y 2 il σ ij σ il σ kan σ ka 1 σ kan σ ki σ ki σ ka 1 .
As the final example, let us now consider the θ + i integration when the internal point i is a quartic vertex connected to bivalent vertices only. In that case the Grassmann variables originate from the R factor at the internal vertex i and we have d 4 θ + i R i 1234 = 1 y 2 i1 y 2 i2 y 2 i3 y 2 i4 y 2 12 y 2 34 σ i1 σ i3 σ i2 σ i4 σ i1 σ i2 σ i3 σ i4 + y 2 13 y 2 24 + y 2 14 y 2