Far Beyond the Planar Limit in Strongly-Coupled $\mathcal{N}=4$ SYM

When the $SU(N)$ ${\cal N} = 4$ super-Yang-Mills (SYM) theory with complexified gauge coupling $\tau$ is placed on a round four-sphere and deformed by an ${\cal N} = 2$-preserving mass parameter $m$, its free energy $F(m, \tau, \bar \tau)$ can be computed exactly using supersymmetric localization. In this work, we derive a new exact relation between the fourth derivative $\partial_m^4 F(m, \tau, \bar \tau) \big|_{m=0}$ of the sphere free energy and the integrated stress-tensor multiplet four-point function in the $\mathcal{N}=4$ SYM theory. We then apply this exact relation, along with various other constraints derived in previous work (coming from analytic bootstrap, the mixed derivative $\partial_\tau \partial_{\bar \tau} \partial_m^2 F(m, \tau, \bar \tau) \big|_{m=0}$, and type IIB superstring theory scattering amplitudes) to determine various perturbative terms in the large $N$ and large 't Hooft coupling $\lambda$ expansion of the ${\cal N} = 4$ SYM correlator at separated points. In particular, we determine the leading large-$\lambda$ term in the ${\cal N} = 4$ SYM correlation function at order $1/N^8$. This is three orders beyond the planar limit.

In this work, we will derive a new relation between the integrated stress-tensor multiplet correlator and four mass derivatives that was previously studied in [40]. Since F (m, τ,τ ) can be computed using supersymmetric localization [39], both the relation derived here and that of [40] impose non-perturbative constraints on the stress-tensor multiplet correlator for any N and (τ,τ ). As an application, we will use these constraints to derive new terms in the perturbative 1/N and 1/λ expansion of the stress tensor correlator, as we will describe shortly.
In more detail, the stress tensor multiplet of the N = 4 SYM theory contains 42 real scalar operators: 20 of them, which we collectively denote by S, have scaling dimension 2, 1 More precisely, at subleading orders in 1/N there are both contributions from loop diagrams and from tree-level diagrams. Some of the tree-level contributions can be separated out from the loop contributions because they have a different scaling in 1/λ. We will provide examples in the next section.
transform as the 20 of the SU (4) R R-symmetry, and in the Lagrangian description are single trace scalar bilinears; another 20 operators, grouped in the complex combinations P and P , have scaling dimension 3, transform in the 10 and 10, respectively, of SU (4) R , and in the Lagrangian description are single trace fermion bilinears; and lastly 2 operators have dimension 4 that will not be important in this paper. Four-point correlation functions of all these operators, as well as of other operators belonging to the N = 4 stress tensor multiplet, are related to one another by Ward identities, and can all be expressed in terms of a single function T (U, V ) of the conformally-invariant cross-ratios U and V [62].
Coming back to the derivatives F 4 (τ,τ ) and F 2 (τ,τ ) of the S 4 free energy, these quantities can be related to T (U, V ) because each derivative w.r.t. m corresponds to the insertion of a specific linear combination of S, P , and P integrated over the four-sphere, while a derivative w.r.t. τ (orτ ) corresponds [63][64][65] to an insertion of a specific component of S at the north (or south) pole of the sphere. Thus, F 2 (τ,τ ) can be written in terms of SSSS and SSP P , where two of the S operators in the first correlator as well as the P and P operators in the second correlator are integrated. Writing these integrated correlators in terms of integrals of T (U, V ) was achieved in [40], following a similar calculation in the 3d Aharony-Bergman-Jafferis-Maldacena (ABJM) theory [66] described in [32]. In this paper, we will perform the same task for F 4 (τ,τ ), which can be written as a linear combination of SSSS , SSP P , and P P P P , where now all four operators are integrated over the sphere. This calculation has two challenges. The first is to write SSP P and P P P P in terms of T (U, V ) using the Ward identities, which we do following the component method also used in [40,62]. 3 The second challenge is to perform the integrals over the sphere, where, unlike in the case of F 2 (τ,τ ), one now encounters additional divergences that need to be regularized while preserving supersymmetry.
To use the relation between F 4 and T (U, V ) in the holographic regime, we should derive an expansion of F 4 (τ,τ ) at large N and large λ. Using Ref. [39] as a starting point, one can write down F 4 (τ,τ ) as an expectation value of an operator in the free Gaussian matrix model at m = 0. As in [41], this expectation value can then be computed to any order in 1/N at finite λ using topological recursion [68,69], and also at finite N and λ (if we ignore non-perturbative instantons in the Nekrasov partition function) using orthogonal polynomials [70].
The expansion of F 4 in 1/N and 1/λ, combined with various other constraints studied in previous work, can be used to fully determine the function T (U, V ) to higher orders in the double expansion in 1/N and 1/λ than was previously possible. In particular, we determine that the Mellin transform 4 [47,71] of the function T (U, V ), which we denote by M(s, t), takes the form where u ≡ 4−s−t, and where c = (N 2 −1)/4 is the c anomaly coefficient, which is the natural expansion for holographic correlators since it is simply related to the effective 5d Newton's constant. In string theory language, the terms at order 1/c g+1 correspond to genus ≤ g string worldsheets, so the expansion (1.3) contains contributions up to genus three.
The expression (1.3) was determined as follows: • Crossing symmetry and the analytic structure of Witten diagrams in Mellin space [46,47,71,72] determine the s, t dependence of each term in the 1/c and 1/λ expansion in Eq. (1.3) up to undetermined coefficients. In particular, the polynomial terms in s, t, u correspond to contact Witten diagrams, where for a polynomial of degree n, the interaction vertex is schematically of the form D 2n R 4 ; the first term at order 1/c corresponds to the tree-level supergravity amplitude; and the M SG|SG term corresponds to the one-loop supergravity amplitude, which is a non-analytic in s, t, u, and was determined in [51,52,56] using unitarity, up to an additive constant.
• The coefficient of the supergravity term is fixed by the requirement that, when expanding the full correlator in conformal blocks, there are no operators of dimension precisely two [21].
• At each order in the 1/c and 1/λ expansion, one can determine the coefficient of the leading term at large s, t, u from knowledge of the flat space scattering amplitude in 4 The precise definition of the Mellin transform is given in Eq. (2.3) below.
type IIB superstring theory. This was originally done in [49] to fully determine the term of order c −1 λ −3/2 (i.e. the genus zero R 4 term).
• At each order in 1/c and 1/λ, one can determine two coefficients, namely one from F 4 and one from F 2 , when these quantities are also expanded in 1/c and 1/λ. In particular, in [40] the F 2 constraint was used to fully determine the term of order c −1 λ −3/2 (genus zero R 4 ), and to also determine the remaining coefficient in the c −1 λ −5/2 term (genus zero D 4 R 4 ) that remained undetermined after using the flat space limit. In [41], the quantity F 2 was computed to any order in 1/N and 1/λ, and used to also fix the Note that the coefficients corresponding to R 4 and D 4 R 4 can be fully fixed using only the two supersymmetric localization constraints, and they do agree, in the flat space limit, with the scattering amplitude in type IIB superstring theory. They appear at genus zero and genus one in the case of R 4 , and at genus zero and genus two for D 4 R 4 . The match between supersymmetric localization and type IIB scattering amplitudes represents a nontrivial precision test of AdS/CFT at these orders. 5 The terms of order 1/c 4 in (1.3) were obtained by combining the supersymmetric localization and flat space limit constraints, and they represent, to our knowledge, the first known contributions to a holographic correlator at genus three.
The rest of this paper is organized as follows. In Section 2, we discuss the stress tensor multiplet four-point function in the strong coupling limit, and fix the higher order in 1/N and 1/λ terms using the flat space limit, the old F 2 constraint, and the new F 4 constraint. In Section 3, we derive this new integrated constraint. In the Appendices we include many details of the calculation, including the localization calculation of F 4 from topological recursion or orthogonal polynomials. We end with a discussion of our results and future directions in Section 4. Several complicated explicit results are given in an attached Mathematica notebook.

N = 4 stress-tensor four-point function
The main object of study in this work is the stress tensor multiplet four-point function. We begin by discussing general constraints on these correlators coming from invariance under 5 The R 4 and D 4 R 4 coefficients were also fixed in in [32] for the ABJM holographic correlator, which is the N = 4 superconformal algebra. We then discuss the large N strong coupling expansion in Mellin space for the N = 4 SYM theory. Finally, we discuss how to constrain the terms in this expansion from the known Type IIB S-matrix in the flat space limit, as well as using the F 4 (τ,τ ) and F 2 (τ,τ ) introduced in Eqs. (1.1) and (1.2) in the Introduction.

Setup
As mentioned in the Introduction, we denote the bottom component of the stress tensor multiplet by S. This operator is a dimension 2 scalar in the 20 of the SU (4) R ∼ = SO(6) R , and can thus be represented as a rank-two traceless symmetric tensor S IJ ( x), with indices I, J = 1, . . . , 6. However, in order to avoid a proliferation of indices, it is customary to contract them with null polarization vectors Y I , with Y · Y = 0. Superconformal symmetry [62] implies that the four-point function of Here, as before, c is the conformal anomaly coefficient, which for an SU (N ) gauge group equals c = (N 2 − 1)/4; the quantities U ≡ are the usual conformal invariant cross-ratios; and Y ij ≡ Y i · Y j are SO(6) R invariants. Importantly, the only nontrivial information in the correlator (2.1) is encoded in the single function T (U, V ).

Strong coupling expansion
We now restrict our discussion to the case of the SU (N ) N = 4 SYM theory, and discuss the strong coupling 't Hooft limit, where we take N → ∞ (or c → ∞) with λ ≡ g 2 YM N fixed. If we further take λ → ∞, the holographic correlator can be computed from Witten dual to M-theory on AdS 4 × S 7 , using similar localization constraints. In that case, however, while the R 4 coefficient is non-zero, the D 4 R 4 coefficient vanishes. 6 The Mellin transform can also be defined away from the strong coupling limit. For recent work on this topic, see [73]. diagrams in an expansion around AdS 5 × S 5 supergravity. In the strong coupling limit, it is convenient to work with the Mellin transform M of T via 6 Crossing symmetry M(s, t) = M(t, s) = M(s, u) and the analytic properties of the Mellin amplitude (for a detailed description, see [41]) then restrict M(s, t) to have a 1/c and 1/λ expansion of the form which can be transformed to position space using (2.3) to get (2.5) Here, the B's are numerical coefficients that cannot be fixed from symmetry alone. As mentioned in the Introduction, terms at order 1/c g+1 correspond in the flat space limit to genus-g corrections to the Type IIB S-matrix. On AdS 5 × S 5 , these terms receive contributions from l-loop Witten diagrams with l ≤ g. The leading order term is tree-level supergravity, whose expression in Mellin and position space is [49,74] , where the position space expression is written in terms of the functionsD r 1 ,r 2 ,r 3 ,r 3 (U, V ) defined in [8]. The coefficient of M SG is fixed by requiring that the unprotected R-symmetry singlet of dimension two that appears in the conformal block decomposition of the free part S free is not present in the full correlator [21]. In our conventions [40], this amounts to setting the coefficient of M SG to 8/c. To the order considered in (2.5), the only loop term is T SG|SG , which arises from a loop Witten diagram with two supergravity vertices and so scales like 1/c 2 . This term was determined in [52,56,58] using unitarity methods up to a contact term ambiguity, which was further fixed in [41]. Our convention for M SG|SG follows [41], which is the Mellin transform of T SG|SG in [52], although we will not make use of the explicit forms of these quantities.
The remaining terms in (2.5) arise from contact Witten diagrams whose vertices are higher derivative corrections to tree-level supergravity. In particular, the functions M n and T n correspond to vertices of the form D 2n R 4 , and their expressions in Mellin and position space are [28] M 0 = 1 , (2.7) We will now fix the various B's in (2.5) using type IIB string theory and/or the localization constraints.

Constraints from flat space type IIB string theory
Following the general approach of [47], one can relate the holographic correlator SSSS on AdS 5 × S 5 as written in terms of the Mellin amplitude M(s, t) in (2.3) to the type IIB S-matrix. The scattering amplitude of four gravitons (or superpartners) in type IIB string theory takes the form where A SG is the tree-level supergravity amplitude, and s, t, u = −s − t are the Mandelstam invariants. The full amplitude as well as the tree-level supergravity amplitude in (2.8) depend on the momenta and polarizations of the scattered particles, which is information that we suppress in writing down (2.8). The function f (s, t) has been computed in a small g 2 s expansion to genus-two for finite s [75,76], and to genus-three [77] to the lowest few orders in s . We will consider the following terms in the small g s and s expansion: .
Higher orders in s can come from contact terms of higher derivative correction to supergravity, which are analytic in s, t, u and have an expansion in g s , as well as loops, which are non-analytic in s, t, u. The first few higher derivative terms are R 4 , D 4 R 4 , and D 6 R 4 . These are the only protected terms. They receive corrections at genus-zero for R 4 , genus two for D 4 R 4 , and up to genus three for D 6 R 4 . These take the form (2.10) The only loop term shown in (2.9) is the one-loop term with two supergravity vertices, which can be computed from the genus-zero supergravity term using unitarity cuts [78].
The Mellin amplitude M(s, t) is then related to the function f (s, t) according to the flat space limit formula [30,32,40,46,47]: where the momenta of the flat space S-matrix are restricted to lie within five of the ten dimensions. (When taking this limit, one uses the AdS/CFT dictionary 7 to first write the correlation function in terms of g s , s , and L, and then one takes L/ s to infinity as in (2.11).) We can then use the known terms in (2.10) for the type IIB S-matrix 7 In the strong coupling limit we consider in this paper, the θ angle does not appear.
to fix the leading s, t terms B m m in the AdS 5 × S 5 correlator: Constraints from flat space limit: , where the constraints on the R 4 and D 4 R 4 coefficients were already derived in this way in [49]. Note that the R 4 term is thus entirely fixed from the flat space limit alone.

Constraints from supersymmetric localization
As mentioned in the Introduction, we can also constrain SSSS just from the mass-deformed sphere free energy F (m, τ,τ ), which Ref. [39] expressed as an N -dimensional matrix model integral using supersymmetric localization. We have two such constraints, one coming from The first one was shown in [40] to take the form 8 (2.14) As we will show in the next section, we can simplify this expression using crossing symmetry, The second constraint, whose detailed derivation we postpone until Section 3, takes the form The right-hand sides of Eqs.
respectively, which, when evaluated on the functions of position defined in (2.6)-(2.7) are (2.17) The , and I 2 [T 2 ] integrals were first computed numerically in [40] using the position space expressions (2.7). In Appendix B, we confirm these results analytically using the simplified version (2.15) of the constraint (2.14), and we also compute analytically the The LHS of (2.15) was computed to leading order in the 't Hooft limit in [79], and used along with the flat space limit in [40] to fix the coefficients of both R 4 and D 4 R 4 at genus zero. In [41], this computation was extended to O(N −6 ) and finite λ using topological recursion. As reviewed in Appendix A, the finite λ result takes the form of a single Fourier integral, which can then be expanded analytically to any order in 1/λ following Appendix D of [40]. This can be used to fix the remaining non-zero genus terms in R 4 and D 4 R 4 , as well as the one-loop ambiguity term M 0 , giving Constraints from F 2 : In Appendix A, we similarly use topological recursion to compute the LHS of (2. 16) to O(N −6 ) for finite λ. 9 The result now involves two Fourier integrals, which cannot be analytically expanded in 1/λ as in Appendix D of [40] unless the Fourier integrals factorize.
Instead, we had to resort to a numerical large λ expansion, which we show in (A.33). Without using the flat space limit, the two integrated constraints can then be used to fix more coefficients beyond those in (2.18), namely we completely determine the D 4 R 4 coefficients and determine relations between the D 6 R 4 ones: Constraints from F 4 and F 2 : .
In addition, we can check that the constraint coming from F 4 by itself is sufficient to determine the R 4 coefficients constrained using F 2 in ( to get which is one of our main results, and which we also quoted in the Introduction. In particular, we have fixed the genus-zero and genus-three terms in D 6 R 4 , where the latter scales as λ 3 /N 8 and is the first result in N = 4 SYM at three orders beyond the planar limit! Note that we cannot yet fix the genus-one and genus-two terms in D 6 R 4 , since we have been unable to accurately expand the localization constraint at large λ to the required orders yet. Now that SSSS has been fixed to the order shown in (2.20), we can use it to extract any CFT data to this order that we like. For instance, we find the anomalous dimensions γ j of the unique lowest twist even spin j double trace operators [S∂ µ 1 . . . ∂ µ j S] to be where the three O(c −1 ) terms were computed in [80,81], [49], and [40], respectively, while the one-loop supergravity term was computed in [41,52]. Contact terms with n-derivatives only contribute to operators up to spin n/2 − 4, as explained in [20]. For higher twist there are many degenerate double trace operators, so one would need to compute many different half-BPS correlators to determine their anomalous dimensions [51,52].

Constraints from the sphere partition function
In this section, let us complete our discussion by relating T (U, V ) defined in (2.1) to the derivatives of the S 4 partition function of the mass-deformed N = 4 SYM theory that were introduced in Eqs (1.1)-(1.2). Both quantities involve special types of supersymmetric operators on S 4 that were considered in [40]. The first type are specific components of S IJ that are Coulomb branch operators from an N = 2 point of view, and that are placed at the poles of S 4 . Adding them to the action is equivalent to changing the gauge coupling [63][64][65].
The second type of supersymmetric operators considered in [40] are integrated operators that couple to an N = 2-preserving real-mass deformation. As mentioned in the Introduction, Ref. [40] considered only the quantity 2)) that equals the four-point function of two operators of the first type mentioned above and two operators of the second type. This quantity gave the constraint on T given in (2.14). In this section we consider the constraint coming from the quantity F 4 (τ,τ ) ≡ ∂ 4 m F m=0 (defined in Eq. (1.1)). This quantity corresponds to four integrated insertions, and, for reasons that will become clear, this integrated correlator is more difficult to analyze than the correlator that led to (2.14).

Correlators in N = 4 SYM
Each integrated insertion is a linear combination of specific components of S IJ and specific components of the ∆ = 3 operators P AB and their conjugatesP AB transforming in the 10 and 10 of SU (4) R , respectively. Here, the indices A, B = 1, . . . , 4 are placed in a lower (upper) position when they correspond to the 4 (4) of SU (4) R . Since the 10 and 10 are rank-two symmetric products of 4 and 4, respectively, the operators P AB andP AB are symmetric tensors. For concreteness, we normalize S and P such that vectorsX and X, respectively. The dot product in (3.1) and in subsequent equations stands for contraction using the Kronecker delta symbol. Before discussing in detail which components of S, P , andP participate in the mass deformation, let us point out that the four-point function PP PP is restricted by conformal symmetry and R-symmetry to take the form where a basis for the three distinct SU (4) R invariants can be taken to be Thus, the PP PP correlator involves three functions P i (U, V ), with i = 1, 2, 3. We will also need the mixed correlator SSP P , which also involves three functions that we denote by R i (U, V ), with i = 1, 2, 3: where, in this case, the basis of SU (4) R invariants can be taken to be Here {X, X, Y 1 , Y 2 } is an SU (4) R invariant that can be formed in the product 10 ⊗ 10 ⊗ 20 ⊗ 20 defined in Eq. (A.8) of [40].
The coefficients P i (U, V ) and R i (U, V ) are related by the supersymmetric Ward identity to the function T (U, V ) appearing in the SSSS correlator. These relations are tedious to derive and quite complicated, so they are relegated to Appendix C. Separating out the contribution from the free theory as in Eq. (2.1), these relations take the form where the free theory contributions are and the differential operators R and P are given in (C.1) and (C.5)-(C.7), respectively.

N = 2-preserving mass deformation
Let us now discuss the mass deformation on S 4 in more detail. It is customary to describe the N = 4 SYM theory as an N = 1 gauge theory with three adjoint chiral multiplets with scalar and fermionic components that we will denote by Z i and χ i , respectively, where i = 1, 2, 3, and with canonical kinetic terms. In the decomposition of the N = 4 vector multiplet into an N = 2 vector multiplet and hypermultiplet, we can consider the Z i and χ i with i = 1, 2 as forming the hypermultiplet, and Z 3 and χ 3 as part of the N = 2 vector multiplet. The mass deformation we consider corresponds to giving a mass to the hypermultiplet fields. On S 4 , this mass deformation takes the form where we assumed that the radius of S 4 is set to one, and where the operators J, K, and L are given by (3.9) While the operators K and L multiplying m and m 2 , respectively, are familiar from a flatspace mass deformation, the operator J (with its coefficient being inversely proportional to the radius of S 4 ) is present here in order to preserve N = 2 supersymmetry on S 4 .
The operators appearing in (3.9) can be written in terms of specific components of operators with well-defined transformation properties under the SU (4) R symmetry. In particular, J can be written in terms of S IJ and K can be written in terms of P AB andP AB : where, as before, the anomaly coefficient equals c = (N 2 − 1)/4.

Four mass derivatives
The four-point function contribution to the fourth mass derivative of the free energy is where I d ∆ A ,∆ B [G] denotes the integrated correlator on S d of four operators of dimensions (∆ A , ∆ A , ∆ B , ∆ B ), which was studied in detail in [32]. To write it down explicitly, first note that on R d , such a correlation function takes the form G(U, V )/(| x 12 | 2∆ A | x 34 | 2∆ B ). On 10 They are −12 a round S d of unit radius parameterized in stereographic coordinates such that the line element is ds 2 = Ω( x) 2 d x 2 , the analogous correlator is , (3.14) The integrated correlator on S d is then The quantity (3.15) was evaluated in [32], where, for d = 4, it was found that (3.16) Here, theD r 1 ,r 1 ,r 2 ,r 2 function is related to a contact Witten diagram in AdS 4 of four scalar fields dual to operators of dimensions r 1 , r 1 , r 2 , r 2 . While one can write explicit positionspace expressions for theD functions we need, for our purposes, however, the most useful definition theD function is through the Mellin transform: (3.17) Note that directly plugging (3.16) into (3.13) is problematic, because I 4 ∆ A ,∆ B [G] contains a factor of Γ(6 − ∆ A − ∆ B ), which diverges for the I 4 3,3 terms from the second line of (3.13). We will thus have to find a way to regularize this divergence.
The expression (3.13) can be split into two parts: one that is independent of T (U, V ) corresponding to the free theory, and one that is linear in T (U, V ) and its derivatives: The first term, −(∂ 4 m F ) free m=0 can be calculated by plugging S = S free , R = R free , P = P free defined in (2.2) and (3.7) into (3.13). When performing this calculation, there are various divergences that arise and one has to be careful to regularize them properly. We will not do that here, and instead calculate −(∂ 4 m F ) free m=0 using a different method that avoids some of these complications. This alternative method relies on the observation that where F H (m) is the S 4 free energy of a hypermultiplet of mass m. The relation (3.19) holds true because in the zero coupling limit, the SU (N ) SYM theory of central charge c = (N 2 − 1)/4 has N 2 − 1 = 4c such hypermultiplets of mass m.
The massive hypermultiplet free energy F H (m) can be determined as follows. We start by writing down the theory of a single hypermultiplet with scalar and fermionic components (Z i , χ i ), as in (3.8). The action is  This expression can be simplified and then regularized: where G is the Barnes G-function and γ is the Euler-Mascheroni constant. The normalization of Z H (m) in (3.22) was chosen such that Z H (0) = 1. The function H(m) appeared in the supersymmetric localization computation of [39], and indeed, the result (3.22) can be also read off from [39]. The equation (3.22) is imprecise, however, partly because the regularization of (3.22) possesses ambiguities, and partly because we dropped an unambiguous overall coefficient that depends on the radius of the sphere, as required by the conformal 11 The eigenvalues of the bosonic operator are (n + 1 + im)(n + 2 − im) and (n + 1 − im)(n + 2 + im), with n = 0, 1, 2, . . ., each with degeneracy D n = 1 6 (n + 1)(n + 2)(2n + 3). The eigenvalues of the fermionic operator are n + 2 + im n + 2 − im, each with degeneracyD n = 1 3 (n + 1)(n + 2)(n + 3).

T ] + (2-and 3-pt function contributions)
, where the 2-point and 3-point function contributions here are the subset of the ones from (3.12) that were not accounted for in 4c copies of the free theory. We will not write them in detail because, as we will discuss, we believe that the boundary terms from the integration by parts we will be performing shortly precisely cancels them. Such a phenomenon was observed also in [32] in 3d.
Let us study the first three terms in (3.24) separately, and let us aim to write them in the "canonical form" where u = −2−s−t, where each term will have a different function Q(s, t, u). The expression For the first term in (3.24), we combine (3.16) with (3.17) and shift the integration variables as appropriate to obtain For the second term, we plug (C.2) into (3.16), and then we integrate by parts to put all the derivatives on theD 2,2,1,1 function. We obtain (3.27) Then, using the expression (3.17) for theD 2,2,1,1 function in Mellin space, we can write (3.27) in the canonical form (3.25): For the third term in (3.24), we should follow a similar procedure. Since the prefactor in (3.16) diverges when ∆ A = ∆ B = 3, we should evaluate this quantity in a way that avoids this divergence. This can be done by first considering ∆ A = 3 and ∆ B = 3 − , both in the prefactor and in the expression for theD function in (3.16). 13 We then use the Ward identity (C.8) (which, as will be justified shortly, should also hold for non-zero ), integrate by parts, and then take → 0. This procedure gives The fact that we can use the same Ward identity (C.8) when = 0 can be justified as follows. of SU (4) R , respectively, and have scaling dimensions ∆ p + 1. The form of any of the mixed correlation functions with two operators from the p = 2 multiplet and two operators from 13 We are grateful to Thomas Dumitrescu for extensive discussions about this issue. a p = 2 multiplet is precisely the same as when p = 2, except that there is an additional factor of (Y 3 · Y 4 ) p−2 . Moreover, the Ward identity relations relating PP P pPp to SSS p S p is independent of p. Analytically continuing PP P pPp in p to p = 2 − leads to ∆ B = 3 − as above.
Combining the above results, we have (3.30)

Simplification using crossing symmetry and final formula
, or (V, U ) (along with corresponding changes in (r, η)), and then using the relations in (3.31). Averaging over these six possibilities (the original expression (3.25) as well as the five expressions obtained as above), one obtains a similar expression to (3.25), with the only difference being that the factor U s 2 V t 2 is symmetrized in s, t, and u: where u = −2 − s − t. One can then rename s, t, and u to rewrite (3.25) such that Q(s, t) is replaced by the symmetrized expression Q(s, t) → Q(s, t) + Q(s, u) + Q(u, t) + Q(t, s) + Q(t, u) + Q(s, u) 6 . (3.33) After this symmetrization, Eqs. (3.26), (3.28), and (3.29) imply that we can replace Q SSSS , Q SSP P , and Q P P P P by Q P P P P = 2π 7 9 (s 2 + t 2 + u 2 + 4) .

(3.34)
Plugging (3.34) into (3.30), we see quite nicely that the dependence on s, t, and u inside the argument of F disappears, and we simply have Recalling that evaluating F at Q = 1 means replacing the second line of (3.25) with with I 4 [T ] defined in (2.16). Combining this expression with (3.23), we immediately obtain (2.16). This is our final formula for F 4 (τ,τ ).
There are two loose ends to be tied up. The first concerns the 2-and 3-point function contributions in (3.24). We note that the mass parameter m has dimension 1, so m must couple to an operator of dimension precisely 3 and m 2 must couple to an operator of dimension precisely 2. Away from zero coupling, the only such operators present are operators in the SU (N ) N = 4 stress tensor multiplet, whose canonically-normalized 2-and 3-point functions are proportional to c. Thus, the entire 2-and 3-point function contribution to (3.12) must be proportional to c as well, so any such contribution that we did not take into account would simply modify the first term in (2.16). However, we did check that with the formula (2.16) as written, there is agreement between the leading large c supersymmetric localization result for F 4 and the explicit evaluation of the RHS of (2.16) using the known supergravity amplitude (2.6). This is a strong check that the coefficient of the term proportional to c in (2.16) is as written.
The second loose end concerns the simplified formula (2.15) for the mixed derivative F 2 (τ,τ ). As mentioned before, Ref. [40] derived the relation (2.14) between the mixed derivatives of the N = 2 * partition function and an integral of T (U, V ). One can simplify this formula using crossing symmetry, namely by taking the average between (2.14) and the expression obtained after replacing (U, V ) in (2.14) by

Discussion
The main result of this work was a new exact relation between four derivatives of the mass deformed sphere free energy, F 4 (τ,τ ) ≡ ∂ 4 m F (m, τ,τ ) m=0 , and an integral of the four point function SSSS of the superconformal primary S of the stress tensor multiplet in N = 4 SYM theory. For gauge group SU (N ), we applied this constraint in the strong coupling 't Hooft limit at large c ∼ N 2 and large λ, where the N = 4 SYM theory is holographically dual to type IIB string theory on AdS 5 × S 5 . In combination with the constraint coming from F 2 (τ,τ ) = ∂ τ ∂τ ∂ 2 m F (m, τ,τ ) m=0 derived in [40], the F 4 (τ,τ ) constraint allowed us to completely fix the D 4 R 4 contributions to the SSSS correlation function. (This contact interaction vertex contributes non-trivially at genus zero and genus two.) In the flat space limit, we matched these contributions to the known D 4 R 4 terms in the Type IIB S-matrix.
Using the constraint on SSSS from the known flat space S-matrix combined with the two constraints from F 2 and F 2 , we were able to further fix the genus-zero and genus-three D 6 R 4 term in SSSS , where the latter scales as λ 3 c −4 and is the first known contribution to SSSS that has been computed at order 1/c 4 .
Looking ahead, it would be useful to develop an analytic large λ expansion of F 4 (τ,τ ), as was achieved for F 2 (τ,τ ) in [41]. In the latter case, the large N and finite λ expressions obtained by topological recursion were given in terms of a single Fourier integral, which could then be analytically expanded to any order in λ using the method described in Appendix D of [40]. For F 4 (τ,τ ), however, the large N and finite λ expressions that we derived in this work are given in terms of two Fourier integrals, which were not amenable to the method of Appendix D of [40] unless the Fourier integrals factorized. Instead, we had to resort to a numerical large λ expansion, which only gave precise results at low orders in λ. In particular, we were unable to compute the N 0 λ 0 , N 0 λ − 1 2 , and N −2 λ contributions to F 4 (τ,τ ). The first term could be used to confirm the derivation of the one-loop constant ambiguity B SG|SG 0 that was previously fixed in [41] using F 2 (τ,τ ). The N 0 λ − 1 2 and N −2 λ terms could be used to derive the genus one and two contributions, respectively, to D 6 R 4 in SSSS , which would complete the derivation of the D 6 R 4 term begun in this work.
While in this work we considered the strong coupling 't Hooft limit, one could also consider the holographic limit where N → ∞ and τ = θ 2π + 4πi g 2 YM is finite. In the flat-space limit, this strong coupling limit of the SSSS correlator should match the small s and finite τ s = χ s + ig −1 s expansion of the type IIB S-matrix, where χ s is the expectation value of the type IIB axion. The coefficients of the various powers of 1/c in SSSS for each expansion must be SL(2, Z) invariants of τ andτ . In the flat space limit, the coefficients of 1/c 7/4 , 1/c 9/4 , and 1/c 5/2 correspond to the protected R 4 , D 4 R 4 , and D 6 R 4 contact amplitudes that were derived in [82][83][84][85][86]. In [42], the mixed mass derivative F 2 (τ,τ ) was studied in this limit, and combining the integrated constraints with the flat space limit, it was possible to completely determine the SSSS correlator at orders 1/c 7/4 and 1/c 9/4 . It would be very interesting, but much harder, to extend this analysis to F 4 (τ,τ ). We leave this topic for future work [87].
In addition to the constraints on the SSSS correlator considered here, one could also consider new integrated constraints that come from derivatives in terms of the squashing parameter b for the free energy F (b, m, τ,τ ) on the squashed sphere, which was also computed in terms of a matrix model using localization in [88]. Of the three possible constraints: we expect that only one of the first two constraints to be linearly independent from the two already considered, which is exactly enough constraints to fix the D 6 R 4 term in SSSS purely from CFT. These three localization constraints could also be combined with the known type IIB S-matrix in the flat space limit to fix the four ambiguities in the one-loop term M SG|R 4 genus-0 with one supergravity vertex and one genus-zero R 4 vertex [55,56], which scales like λ − 3 2 N −4 . One could similarly fix the D 8 R 4 contact term to genus two. Lastly, while the application of integrated constraints and localization to holographic correlators has been to the large N expansion in this paper and previous work [40][41][42], these relations are in fact non-perturbative, and so could be applied to the numerical bootstrap for N = 4 SYM [89,90]. For this purpose, the finite N formula for the perturbative part of the mass deformed free energy, as derived using orthogonal polynomials in Appendix A.3, will be especially useful, especially if one could augment it with a similar formula for the contribution from the Nekrasov partition function. These constraints could allow one to impose the values of τ andτ in the numerical bootstrap for finite N , just as N was imposed in the original studies [89,90] using the conformal anomaly c, and thereby solve N = 4 SYM numerically for all τ ,τ and N .

Acknowledgments
We thank Thomas Dumitrescu and Ofer Aharony for useful discussions, as well as Damon A F 4 (τ,τ ) from supersymmetric localization In this appendix, we show how F 4 (τ,τ ) can be computed from the supersymmetric localization result of [39], following a similar calculation for F 2 (τ,τ ) in [41]. We will start by writing F 4 (τ,τ ) as an expectation value of an operator in a Gaussian matrix model. We then evaluate this expectation value to any order in 1/N 2 using topological recursion [68,69], or for finite N and λ = 4πN Im τ (ignoring instantons) using orthogonal polynomials [70].

A.1 Matrix model expectation value
As shown by Pestun [39], the S 4 partition function Z = exp(−F ) of the SU (N ) N = 2 * is given by where we denoted a ij ≡ a i − a j , and where the delta function enforces the SU (N ) constraint that the eigenvalues sum to zero. The function H(m) appearing in (A.1) was already defined in the main text in Eq. (3.22). The quantity |Z inst (m, τ )| 2 represents the contribution to the localized partition function coming from instantons located at the North and South poles of S 4 [91][92][93][94], and can be ignored in the 't Hooft limit because it is non-perturbative when The LHS of the perturbative part of the integrated constraint (3.36) is then where K(z) ≡ − H (z) H(z) , and where the expectation values are taken in the Gaussian matrix model The function K (z) can be simply expressed using its Fourier transform [79] K (z) = − The 2-body term I(ω) also occurs in the matrix model computation of the LHS of (2.15), whose non-instanton part can be written as [41] c 8 This quantity was actually already computed in [41] to all orders in 1/N 2 in the 't Hooft limit for finite λ (and to any order in 1/λ) using topological recursion [68,69], and for finite N, λ (ignoring instantons) using orthogonal polynomials [70]. We will apply these methods to the 4-body term J (ω, w), and then combine with the known results for I(ω) to compute (A.7).
A.2 1/N 2 expansion from topological recursion with γ i chosen so that the contour lies to the right of all singularities in the integrand, we then have In the formalism involving resolvents, a similar computation shows that the inverse Laplace transforms of an n-point function of resolvents in the SU (N ) and U (N ) matrix models are related by Thus, as long as the arguments of R n sum to zero, as is the case in (A.11), there is no difference between the U (N ) and SU (N ) theories, so we will drop the subscript U (N ) in what follows.
The correlators of resolvents obey various relations similar to Ward identities in QFT. In particular, the change of variables a i → a i + δa i with δa i = 1/(z − a i ), with infinitesimal, leads at first order in to the relation The more complicated change of variables corresponding to δa i = 1 Eqs. (A.14)-(A.15), combined with large N factorization properties, lead to recursion relations that allows one to determine R p recursively in p and in 1/N . It is customary to write down these recursion relations in terms of the connected correlators W n (y 1 , . . . , y n ) ≡ R n (y 1 , . . . , y n ) conn = N n−2 . (A. 16) In terms of the inverse Laplace transforms of W n , we have and where we define The resolvents can then be expanded in 1/N 2 as W n (y 1 , . . . , y n ) ≡ ∞ m=0 1 N 2m W n m (y 1 , . . . , y n ) , (A. 22) and each genus-m term W n m can be computed for finite λ using a recursion formula in n, m [68,69] starting with the base case W 1 0 , as described e.g. in [41]. We use resolvents up to n + m ≤ 5, which we give in an attached Mathematica file. 14 We then take the inverse Laplace transform in (A.19) to get the 1/N 2 expansion at finite λ for J (ω, w) in terms of integrals over the Fourier variables w, ω from (A.4). For instance, at leading order in 1/N 2 we need only consider the genus-zero resolvents in J 0 (ω, w), which give 14 Some of these resolvents were already shown in Appendix B of [41].
We can then plug this expression, along with the leading order term in I(ω) as given in Appendix B of [41], into (A.2) to get the leading order in N 2 result at finite λ: In the attached Mathematica file, we give explicit formulae for J (ω, w) to order O(N −4 ), which similarly take the form of four Bessel functions, while the expressions for I(ω) were already given in Appendix B of [41] and consist of two Bessel functions.
We would also like to take the large λ expansion of these results, so that we can apply them to the strong coupling expansion of the integrated correlator. For the first term in (A.7), which depends on I(ω), the large λ expansion can be performed just as in [41], and Note that none of these terms have the right powers of π compared to the holographic correlator, so all of them must be cancelled against corresponding terms in the second term in (A.7). For these terms, which depend on J (ω, w), every W n m except W 2 0 factorizes in terms of their argument y i , so the inverse Laplace transform can be easily taken and gives products of four Bessel functions of the form for various integers a, b, n i . As described in Appendix D of [40], we can take the large λ limit of these Bessel functions using the Mellin-Barnes form where the integrals over w, ω can be separately done with factors csch 2 w, csch 2 ω from (A.4) by twice using the identity (A. 28) and then the contours can be closed to the left to get an expansion in 1/λ. For these factorizable terms, using the finite λ expressions for J in the attached Mathematica file, we where we have only showed the terms that we will use, and recall that there is no leading order factorizable term.
The only exception to factorizability is the genus-zero 2-body resolvent W 2 0 (y 1 , y 2 ) 15 which takes the form and whose inverse Laplace transform for distinct imaginary arguments is This term shows up in both J 0 (ω, w) and J 1 (ω, w), and so appears at every order in the large N 2 expansion of J (ω, w). We do not know how to take the large λ expansions of such 15 Note that W 2 0 (y 1 , y 2 ) also appears in I(ω), but in this case there is only one Fourier variable w so the integral was always of the form (A.28).
terms, since the w, ω dependence does not factorize due to the (w + ω) in the denominator of (A.31). For these terms, we instead performed the large λ expansion numerically by evaluating the w, ω integrals at many value of λ at high precision and fitting a curve. Using the expressions in the Mathematica file, we get .
. to the holographic correlator. Note that for the leading N 2 term we were able to do the numerical large λ expansion to many orders, but for subleading terms in 1/N 2 we were only able to accurately read off a couple orders in large λ sufficient to get the terms shown here.

A.3 Finite N from orthogonal polynomials
We can also compute J (ω, w) at finite N and λ in terms of four finite sums using the method of orthogonal polynomials [70], as was already done for I(ω) in [41]. We start by writing I(ω) and J (ω, w) as I(ω, w) = N (N − 1) cos(2ω(a 1 − a 2 )) + N , J ijkl ≡ cos(2ω(a i − a j )) cos(2w(a k − a l )) − cos(2ω(a i − a j )) cos(2w(a k − a l )) . (A.34) We then introduce a family of polynomials p n (a) using the Hermite polynomials H n (x): which are orthogonal with respect to the Gaussian measure As shown in [41], these orthogonal polynomials can be used to write the expectation value of an n-body operator as an n-dimensional integral: was computed in this way in [41] and yields where L b a (x) are generalized Laguerre polynomials. The expression for J (ω, w) similarly involves sums over L b a (x), but takes a rather complicated form that we give in the attached Mathematica file. After plugging these terms into (A.7) we can perform the sums and Fourier integrals for any finite N and compare to the topological recursion results (as given in the attached Mathematica file). We find that they match even down to N = 2 for a large range of λ, which is a nontrivial check.

B Evaluating I 2 analytically
In this Appendix, we describe how to compute The integrals over r, θ are standard one-loop integrals in four dimension, which can be done explicitly to get The integrals over s, t can be done for polynomial M(s, t) (or M SG (s, t)) by twice applying the Barnes lemma: which holds for contours for which the poles of each Gamma function lie either to the left or to the right of the contour. Applying this to (2.6) and (2.7) we get the first line in (2.17).
C Ward identities C.1 Ward identity for SSP P As mentioned in the main text, conformal and R-symmetry invariance implies that the four-point function SSP P takes the form given in (3.4). The non-trivial information is encoded in the functions R i (U, V ), which are related by SUSY Ward identities to the functions S i (U, V ) defined in (2.1). The relations are given in Eq. (B.6) of [40].
Since the SSSS correlator is split into a free part and a part depending on a single function T (U, V ), one can also write R i (U, V ) reflecting this split, as we did in (3.6). The non-free part of R i (U, V ) is thus encoded in three differential operators R i (U, V, ∂ U , ∂ V ) that act on the function T (U, V ) from the SSSS correlator. From (B.6) of [40], we deduce that these differential operators are: (C.1) Note that the differential operator R 1 + 3R 3 takes the form C.2 Ward identity for P P P P Let us now move on to discussing PP PP . As mentioned in the main text, conformal symmetry and R-symmetry imply that this correlation function can be written as in (3.2) in terms of three functions P i (U, V ). These functions must be related by Ward identities  to the functions appearing in the SSSS correlator. To derive these relations, we use the component field method of [32,40,62]. This method was already discussed in a closely related context in [40], so we will only present an outline of the derivation here.
In general, we can derive these Ward identities by first determining the most general forms of the four-point functions that are consistent with conformal symmetry and R-symmetry, and then imposing invariance only under the Poincaré supercharges. For P P P P , the relevant Ward identity takes the schematic form 0 =δ P P P χ = P P P P + ∂χP P χ + P P ∂χχ + P P P F + P P P λ , (C.3) whereδ denotes the action of the supercharge, and the other operators in the stress tensor multiplet are summarized in Table 1. This Ward identity will give P P P P in terms of powers and derivatives of U, V of χP P χ , which must be related to other correlators in a chain that will eventually reach SSSS . These variations are 0 =δ SSSχ = χSSχ + SχSχ + SSχχ + SSSj + SSS∂S , 0 =δ SSP χ = χSP χ + SχP χ + SS∂χχ + SSP P + SSP F , 0 =δ SP P χ = χP P χ + SλP χ + SP ∂χχ + SP P j + SP P ∂S , where the first two were already considered in [40]. Finally, we write SSSS in terms of T (U, V ) using (2.1) to obtain the split of the quantities P i (U, V ) into a free part and an interacting part dependent on T (U, V ) as in (3.6). Following the procedure outlined above, we find that the differential operators P i (U, V, ∂ U , ∂ V ) appearing in (3.6) are Note that the differential operators 2P 1 + 2P 2 + P 3 appearing in (3.24) can be simplified to (C.8)