New Modular Invariants in N = 4 Super-Yang-Mills Theory

We study modular invariants arising in the four-point functions of the stress tensor multiplet operators of the N = 4 SU(N) super-Yang-Mills theory, in the limit where N is taken to be large while the complexified Yang-Mills coupling τ is held fixed. The specific four-point functions we consider are integrated correlators obtained by taking various combinations of four derivatives of the squashed sphere partition function of the N = 2∗ theory with respect to the squashing parameter b and mass parameter m, evaluated at the values b = 1 and m = 0 that correspond to the N = 4 theory on a round sphere. At each order in the 1/N expansion, these fourth derivatives are modular invariant functions of (τ, τ̄). We present evidence that at half-integer orders in 1/N , these modular invariants are linear combinations of non-holomorphic Eisenstein series, while at integer orders in 1/N , they are certain “generalized Eisenstein series” which satisfy inhomogeneous Laplace eigenvalue equations on the hyperbolic plane. These results reproduce known features of the low-energy expansion of the four-graviton amplitude in type IIB superstring theory in ten-dimensional flat space and have interesting implications for the structure of the analogous expansion in AdS5×S. August 2020 ar X iv :2 00 8. 02 71 3v 1 [ he pth ] 6 A ug 2 02 0


Introduction
One of the most intriguing features of four-dimensional gauge theories is the possibility of a mysterious duality that exchanges elementary quarks and magnetic monopoles, while relating physics at strong and weak gauge couplings. First conjectured by Montonen and Olive in [1] following the work of Goddard, Nuyts, and Olive (GNO) [2] as a direct generalization of the electric-magnetic duality in Maxwell theories, it is commonly known as the Montonen-Olive duality or S-duality. It was soon realized that such a duality is more likely to hold in a supersymmetric gauge theory rather than in QCD, because supersymmetry provides more control over the spectrum of solitons [3]. In the case of the maximally supersymmetric N = 4 super-Yang-Mills (SYM) theory [4], an overwhelming amount of compelling evidence for S-duality has been provided by analyzing the dyon-monopole bound states [5], from certain topologically twisted partition functions on four manifolds [6], and from embedding this model into type IIB string theory.
Under S-duality, the N = 4 SYM theory with gauge group SU (N ) and complexified gauge coupling is equivalent to the SYM theory with gauge group SU (N )/Z N and gauge coupling τ ∨ = − 1 τ . While the distinction between the gauge group being SU (N ) and SU (N )/Z N is important for studying non-local operators, it does not affect local operators, which are the subject of this work. The S-duality transformation combined with the T-transformation τ → τ + 1 from the periodic identification of the θ-angle gives rise to an SL(2, Z) duality that acts on the complexified coupling as with a, b, c, d ∈ Z and ad − bc = 1. 1 The Coulomb branch of the SYM theory supports an infinite tower of massive BPS particles from W-bosons, monopoles, and their bound states that transform nontrivially under the duality group, while the mass spectrum stays invariant [5]. Correlation functions of local operators have definite SL(2, Z) transformation properties, and correlation functions of half-BPS operators are invariant [8,9]. 1 In general, under S-duality, N = 4 SYM with the gauge group G is mapped to N = 4 SYM with gauge group G ∨ given by the GNO dual of G, and the gauge coupling τ is mapped to τ ∨ = −1/(n G τ ), where n G = 1 when G is simply-laced, n G = 2 for B r , C r , F 4 , and n G = 3 for G 2 . Consequently, for non-simplylaced G, the combination of S-duality and T-duality yields dualities by (extensions of) congruence subgroups of SL(2, Z) [7].
Embedding the N = 4 SYM theory into string theory provides an elegant picture of this non-perturbative duality. The SU (N ) N = 4 SYM theory is the low energy theory on a stack of N coincident D3-branes immersed in asymptotically flat ten-dimensional spacetime.
The gauge coupling τ is identified with the axion-dilaton background τ s = χ + ie −φ (χ being the axion and φ being the dilaton), and the SL(2, Z) duality of the gauge theory is a direct consequence of the SL(2, Z) duality in type IIB string theory [10].
While the string theory perspective is conceptually useful, it is more satisfying to directly investigate the S-duality properties of N = 4 SYM using field theory methods, and this will be our approach here. In fact, one may argue that the field theory methods provide nontrivial support for the duality structures in the quantum gravity theory. Over the past twenty or so years, there have been steady developments on investigating the S-duality properties of N = 4 SYM using field theory methods, including numerous sophisticated checks based on supersymmetric partition functions [6,[11][12][13][14], extensions that incorporate supersymmetric defects [15][16][17][18][19][20], as well as refinements of the duality by keeping track of global structures of the gauge group and topological couplings in the theory [21][22][23].
The goal of this paper is to continue the study began in [24]  is that Z(b, m, τ,τ ) itself can be computed exactly at any N and any coupling (τ,τ ) using supersymmetric localization [25,26] (see also [27][28][29][30][31][32]). Each of the following combinations of derivatives, evaluated at (b, m) = (1, 0), provides, in principle, a different SL(2, Z)-invariant integrated four-point function in N = 4 SYM: 2 Because both m and δb ≡ b − 1 couple in the action to integrated operators that belong to the N = 4 stress tensor multiplet, it should be possible to express all quantities in (1.3) in terms of integrated four-point functions of stress tensor multiplet operators. Of course, it is plausible that not all such integrated correlators are independent, because there may be relations between them that are implied by the N = 4 superconformal symmetry. In fact, as we will discuss in Section 3, one of our main results is a derivation 3 of three linear relations between the quantities in (1.3), as well as the conformal anomaly c, based on the supersymmetric localization results of [25,26].
Taking into account the three linear relations mentioned above, one can take the independent quantities in (1.3) to be . (1.4) The precise relation of these two quantities to integrated correlation functions was explained in [33] and [24], respectively. In slightly more detail, due to the fact that the stress tensor multiplet of N = 4 SYM is a 1/2-BPS multiplet, it can be shown that supersymmetry requires the correlators of any four operators from this multiplet to be algebraically related to a single function T (U, V ) of the conformally-invariant cross ratios U and V [34]. Thus, the two independent quantities in (1.4) should be expressible in terms of integrals T (U, V ) with potentially different integration measures. It is these explicit expressions in terms of integrals of T (U, V ) that were given in [33] and in [24], respectively.
The main question we ask in this work is what modular invariants 4  , this question was answered in [35], building on the work of [33,36] where only the perturbative terms in g YM were studied. Ref. [35] found that this quantity has an expansion in halfthe combinations of derivatives in (1.3) are scheme-independent. In particular, the subtraction of 15∂ 2 b Z in the third quantity is needed for removing such an ambiguity. We will discuss these scheme-dependent ambiguities in Appendix A. 3 For one of these relations, we do not have a full proof, but amass significant evidence in the case where the gauge group is SU (N ). In Appendix B, we will make comments about these relations in N = 4 SYM with a general gauge group. 4 We emphasize here that while all correlators of half-BPS operators are SL(2, Z) invariant, the correlators that involve their superconformal descendants may not be. In particular they would violate the U (1) Y bonus symmetry for five-and higher-point functions [8,9]. For four-point functions, it was conjectured in [8,9] that the U (1) Y bonus symmetry and consequently SL(2, Z) invariance hold for half-BPS operators and their descendants. For four-point functions of stress-tensor multiplet operators, which are of interest here, the U (1) Y invariance follows from the fact that the superconformal Ward identities impose couplingindependent algebraic relations between any four-point function of stress tensor multiplet operators and the U (1) Y -invariant four-point function of the half-BPS superconformal primary of this multiplet. Therefore, while the fourth derivatives in (1.3) and (1.4) are expected to be SL(2, Z) invariant, the modular properties of the higher derivatives will be more complicated. for various half-integer values of s. The perturbative terms in the second quantity in (1.4) were studied in [24], and the non-perturbative contributions will be studied here. As we will show, we find strong evidence that the 1/N expansion of this quantity involves not only the non-holomorphic Eisenstein series but also another class of modular-invariant functions that generalize the non-holomorphic Eisenstein series in the following sense. The Eisenstein series (1.5) satisfies the homogeneous Laplace eigenvalue equation The new modular functions we encounter are solutions to similar Laplace eigenvalue equations but with a source term given by a product of two Eisenstein series: 6 4τ 2 2 ∂ τ ∂τ − r(r + 1) E(r, s 1 , s 2 , τ,τ ) = −E(s 1 , τ,τ )E(s 2 , τ,τ ) . (1.7) In particular, we find that at half-integer orders in 1/N , the second quantity in (1.4) is still written in terms of the Eisenstein series (1.5), while at integer orders in 1/N the expansion is in terms of E(r, s 1 , s 2 , τ,τ ) for various values of r, s 1 , and s 2 . 7 At low orders in the 1/N expansion, these findings are perhaps not entirely surprising, because, as we will explain in Section 4, at these orders one can establish a precise connection between the integrated correlators (1.4), expanded in 1/N , and 10d type IIB superstring scattering amplitudes of gravitons and their superpartners, expanded at low momentum, as a consequence of the AdS/CFT correspondence [37][38][39]. At leading orders in the momentum expansion, the latter quantity contains certain supersymmetric terms that are purely analytic in momentum and whose coefficients are modular functions such as the ones encountered 5 By a half-integer we mean a number in the set Z + 1 2 . 6 The Laplace equation (1.7) and SL(2, Z) invariance do not completely fix E(r, s 1 , s 2 , τ,τ ). In particular, the solution to (1.7) is ambiguous up to a shift by the Eisenstein series E(r + 1, τ,τ ). Later we will fix this ambiguity by specifying the cusp behavior E(r, s 1 , s 2 , τ,τ ) as τ 2 → ∞. See also Appendix C for more details. 7 In particular, we will see that the order 1 N p−2 contributions to the SYM free energy F = − log Z with p ∈ Z ≥0 + 1 2 are given by the Eisenstein series E(s, τ,τ ) with s = p, p − 2, . . . , 3 2 . On the other hand, the order 1 N q contributions with q ∈ Z >0 are controlled by the general modular functions E(r, s 1 , s 2 , τ,τ ) for s 1 , s 2 ∈ Z >0 + 1 2 with s 1 + s 2 = q + 2, q, . . . .
Conversely, if the flat space amplitude is known, it can be used to determine the leading term in the Mellin amplitude. The second fact is that order by order in the 1/N expansion, analytic bootstrap conditions (meaning analyticity, crossing symmetry, and supersymmetry) can be used to write the separated point correlation functions of the stress tensor multiplet operators, encoded in the function T (U, V ) mentioned above, as a finite sum of specific functions of (U, V ) with a priori arbitrary coefficients . The number of coefficients that are not determined by the bootstrap approach grows with the order in the expansion. In particular, the Mellin amplitude corresponding to T (U, V ) is where s, t are the Mellin space variables with u = 4 − s − t, and where we have re-expressed the 1 N expansion as a 1 c expansion in terms of the conformal anomaly coefficient c = N 2 −1 4 of the N = 4 SYM theory, which is more natural from the CFT perspective. In (1.8), the coefficients α, β i , γ i depend on (τ,τ ). 10 To determine these coefficients, the approach proposed in [24,33,35] 11 was to use the integrated correlators that are calculable from supersymmetric localization as well as the flat space limit of the Mellin amplitude. In this case, as we explain 8 See [45] for another perspective on the differential equation (1.7) as coming from constraints of IIB supersymmetry. Analogous arguments have also been applied to higher-point interactions which violate the U (1) symmetry [46]. 9 Functions satisfying (1.7) with various values of s 1 and s 2 were discussed in [47], where they arose in the context of higher-order terms in the low energy expansion of flat-space type II superstring amplitudes. 10 Here M 1−loop is a meromorphic term that corresponds to the regularized supergravity one-loop amplitude in the holographic dual, and it also contains a coefficient that is not determined from the bootstrap approach. We will not discuss this term, but we note that the coefficient mentioned above was determined in [36] from supersymmetric localization. 11 See also [75,76] for similar computations in 3d CFTs. The general method of computing higher derivative in Section 4, the constraints coming from the integrated correlators (1.4), expanded in 1/N , are sufficient to determine all the constants α, β i , γ i in (1.8). 12 As a preview, one finds (1.9) We should note that developing the 1/N expansion of the quantities (1.4) is an onerous task. The supersymmetric localization results of [25,26]  where each factor takes the form of a Nekrasov instanton partition function [25,80,81]. These factors are the ones responsible for the non-perturbative effects we study in this work that are crucial for obtaining the modular functions mentioned above.
The rest of the paper is organized as follows. In Section 2, we start with a general discussion of the large N expansion of the integrated correlators and introduce the new modular functions E(r, s 1 , s 2 , τ,τ ) that generalize the non-holomorphic Eisenstein series. In Section 3, we study in detail the localization constraints coming from (1.4) keeping track of the instanton effects. We apply these results in Section 4 to correlation functions at separated points. In Section 5, we end with a brief summary and discuss future directions. Various technical details are contained in the Appendices. In particular, Appendix C contains some details concerning solutions of (1.7) that are important in guiding the analysis in Section 3.
2 Overview of modular invariants and integrated correlators at large N As mentioned in the Introduction, our interest is in the integrated four-point functions (1.3) of the N = 4 SYM expanded at large N . Before delving into detailed calculations, let us provide an overview of these expansions, review previous results, and state our main results.
corrections using non-trivial CFT was initiated in [77] in 3d, where a certain OPE coefficient computed using localization was used to the fix the R 4 correction. A similar approach was also taken to compute R 4 in 6d in [78], where the nontrivial OPE coefficient was now computed using the 2d chiral algebra subsector of [79]. 12 It is possible that not all relations that reduce the integrated correlators from the list in (1.3) to that

Sphere partition function from supersymmetric localization
The partition function of the mass-deformed SU (N ) N = 4 SYM theory placed on a squashed four-sphere parameterized by squashing parameter b is, up to an overall normalization constant that is independent of b and the mass m, given by [25,26]: 13 The integration is over N real variables a i , i = 1, . . . , N , subject to the constraint i a i = 0. 14 In the integrand of (2.3), |Z inst | 2 captures the contribution from instantons localized at opposite poles of S 4 that we will come to shortly [25,80,81], and the one-loop determinants of SYM fields are written in terms of Υ b (x) which has the following convenient integral definition [84] log .

(2.4)
The bare sphere free energy, F bare = − log Z bare , of the deformed N = 4 SYM theory has logarithmic divergences in addition to power-law divergences, as is the case for all (deformed) sphere free energies in even dimensions. Consequently, the regularized expression in (2.3) has an ambiguity of the form [32,85] log

5)
in (1.4) follow from N = 4 supeconformal symmetry, and if this is the case, then one would be able to determine α, β i , γ i purely from the integrated correlators without the need for the flat space limit. 13 Note that the factors Υ b (0) were missing in the localization formula of [26]. They come from the (regularized) one-loop determinant of the N = 2 vector multiplet associated to each Cartan generator, namely in terms of the Barnes double Gamma function Γ b (x) (see e.g. [82,83] for properties of these special functions). We emphasize that these squashing-dependent factors are crucial for producing the correct CFT correlators by taking derivatives of the SYM partition function with respect to b. 14 The constrained integral over the a i can be implemented, for instance, by an integral over N unconstrained a i 's with a δ( i a i ) insertion.
where the coefficients κ i (τ,τ ) satisfying ∂ τ ∂τ κ i = 0 depend on the regularization scheme (see Appendix A for details). We must therefore be careful to only consider derivatives of log Z such that regularization-dependent terms cancel. For instance, the only well-defined two-derivative term with respect to (b, m, τ,τ ) is ∂ τ ∂τ log Z m=0,b=1 , and it can be seen from the supersymmetric localization result [25] quoted above that it equals is the conformal anomaly for SU (N ) SYM. For three derivatives, we have the non-ambiguous terms which vanish identically. The vanishing of the latter term follows immediately from the fact that Z is an even function of m, while the vanishing of the former comes from the invariance of Z under b ↔ 1/b (see Appendix B for details). At four derivatives, we consider the nontrivial quantities listed in (1.3), which we argue should satisfy the three relations (2.8) We will prove the first two statements, and amass significant evidence for the third. 15 These relations imply that the five quantities in (1.3) can all be written in terms of the two quantities in (1.4). 15 While we focus on the SU (N ) case in the main text, we expect these relations to hold for N = 4 SYM with general gauge groups. See Appendix B for related discussions.

Large N expansion of integrated correlators
As mentioned in the Introduction, the first quantity in (1.4) was previously expanded in 1/N in [35]. The expansion took the form: where the non-holomorphic Eisenstein series E(s, τ,τ ) were defined in our normalization in (1.5). In particular, E(s, τ,τ ) is the unique SL(2, Z)-invariant solution of moderate growth (i.e. behaving as τ a 2 , a ∈ R, as τ 2 → ∞) of the Laplace eigenvalue equation (1.6). The differential operator ∆ τ ≡ 4τ 2 2 ∂ τ ∂τ = τ 2 2 (∂ 2 τ 1 + ∂ 2 τ 2 ) appearing in this equation is the hyperbolic Laplacian. As a periodic function of τ 1 with unit period, the Eisenstein series has a Fourier expansion in the form where the divisor sum σ p (k) is defined by σ p (k) = d>0,d|k d p , and K s− . From (2.10), we see that for τ 2 1 the kth and (−k)th modes (k > 0) contribute a term proportional to e 2πikτ + e −2πikτ , which is the sum of a charge-k D-instanton and a charge-(−k) anti-Dinstanton contributions. While the expression (2.9) was derived in [35] by computing the coefficients of e 2πinτ in a power series in 1/τ 2 , in Appendix D.2 we will provide an alternative derivation in which it is not necessary to expand these coefficients in 1/τ 2 .
In the next section we will determine the first few terms of the large-N expansion of the second quantity in (1.4), ∂ 4 m log Z m=0,b=1 . As we will see, this quantity has an expansion in integer powers of N − 1 2 of the general form where a q and b q are rational numbers multiplying powers of π. The functions G(p, τ,τ ) and H(q, τ,τ ) are modular functions (i.e. functions that are invariant under the SL(2, Z) transformations (1.2)). Whereas the large-N expansion of the correlation function (2.9) was an expansion in half-integer powers of 1/N , the expansion (2.11) contains both even and odd powers of N − 1 2 , as was demonstrated in the analysis of the terms that are perturbative in τ −1 2 in [24].
According to the AdS/CFT correspondence, when (τ,τ ) is fixed, the expansion in powers of N − 1 2 corresponds to an expansion in powers of α /L 2 and α D 2 , where L is the curvature radius of anti de-Sitter space, and D denotes, schematically, a space-time derivative. The large-N expansion (2.11) is therefore interpreted as a small curvature and low momentum expansion in the bulk superstring theory. The leading term in (2.11) is proportional to N 2 and corresponds to the contribution of classical IIB supergravity, which is of order (α ) −4 in our conventions.
Half-integer powers of 1/N As we will see, the coefficients of the half-integer powers of N -the functions G(p, τ,τ )bear a strong similarity to the analogous coefficients in (2.9). They are rational sums of non-holomorphic Eisenstein series E(s, τ,τ ) with half-integer indices 3 2 ≤ s ≤ p + 3 2 . As in (2.9), the first term with a half-integer power of N that appears in the large-N expansion is proportional to N 1 2 E( 3 2 , τ,τ ). Setting s = 3 2 in (2.10) we see that this has perturbative contributions proportional to N N is the 't Hooft coupling. This is the order corresponding to the (α ) −1 R 4 interaction in the flat-space type IIB superstring effective action. The next term with a half-integer power of N is proportional to N − 1 2 E( 5 2 , τ,τ ) and corresponds to an α D 4 R 4 interaction in the flat-space type IIB superstring effective action.

Integer powers of 1/N
The coefficients of the terms in (2.11) that have integer powers of 1/N -the functions H(q, τ,τ )-are modular functions that are linear combinations of the generalized Eisenstein series E(r, s 1 , s 2 , τ,τ ) obeying the inhomogeneous Laplace eigenvalue equations (1.7). These equations generalize the equation satisfied by the coefficient of the D 6 R 4 interaction in the flat-space type IIB superstring S-matrix, which has r = 3 and s 1 = s 2 = 3/2 [44,45]. As can be seen from the equation (1.7), the argument r is an integer that labels the eigenvalue r(r + 1), while s 1 , s 2 are indices of the Eisenstein series in the source term subject to the condition s 1 + s 2 ≤ q + 2 and s 1 , s 2 ≥ 3 2 . 16 As previously mentioned, the solution to Laplace equation (1.7) is ambiguous. Here we fix the ambiguity by requiring E(r, s 1 , s 2 , τ,τ ) ∼ τ s 1 +s 2 2 in the limit τ 2 → ∞. In the case that corresponds to the D 6 R 4 interaction discussed above, this condition is required by consistency of the string perturbation expansion [86]. More generally, we demand H(q, τ,τ ) ∼ τ q+2 2 as τ 2 → ∞, as expected for the genus zero string amplitude corresponding to the D 4q+2 R 4 interaction (after transformation to the Einstein frame).
General equations of the form (1.7) were considered in [47] where the method for extracting the perturbative terms (power-behaved in τ 2 ) in the zero Fourier mode (in τ 1 ) of the functions E(r, s 1 , s 2 , τ,τ ) was presented. These terms are summarized as: where the coefficients a r are easy to determine since the corresponding powers of τ 2 arise from the zero Fourier modes of the source, E(s 1 , τ,τ )E(s 2 , τ,τ ). Equating the power-behaved terms on the left-hand and right-hand sides of (1.7) determines the values of a r . The term proportional to τ −r 2 does not arise in the source term, and furthermore it satisfies the homogeneous equation since it is in the kernel of (∆ τ − r(r + 1)), so its contribution is zero on both the left-hand and right-hand sides of (1.7). 17 The value of its coefficient β r may be determined by the procedure in [44,47], which involves multiplying (1.7) by E(r +1, τ,τ ) and integrating over a cut-off fundamental domain of SL(2, Z). The integral over the left-hand side reduces to a boundary term evaluated at the cutoff, while the integral over the product of three Eisenstein series on the right-hand side can be evaluated by the Rankin-Selberg method. In Appendix C we review the details of this procedure and present explicit results for the perturbative terms appearing in (2.12).

2
. For the cases we consider here r ≥ s 1 + s 2 (see (2.13)) and thus such term is forbidden by the boundary condition that E(r, s 1 , s 2 , τ,τ ) ∼ τ s1+s2 2 as τ 2 → ∞. like for the ordinary Eisenstein series E(s, τ,τ ). Indeed, the zero mode consists of the sum of the power-behaved terms in (2.12) together with an infinite series of D-instanton/anti D-instanton pairs with zero net instanton charge. This non-perturbative structure will also be studied in Appendix C. The cases that will be considered in this paper are the following: • The 1/N contribution, which is the order where the flat-space interaction (α ) 2 D 6 R 4 appears. This is a 1/8-BPS interaction of dimension 14 with coefficient the generalized Eisenstein series E(3, 3 2 , 3 2 , τ,τ ) that was determined in detail in [86]. Although we will not perform a complete analysis here, we will find strong evidence that the fourth mass derivative ∂ 4 m log Z m=0,b=1 as computed from supersymmetric localization is also proportional to this generalized Eisenstein series.
• The 1/N 2 contribution, which is the order where the flat-space interaction is (α ) 4 D 10 R 4 .
In this case we will see the function H(2, τ,τ ) contains a rational sum of two new modular functions, E(4, 3 2 , 5 2 , τ,τ ) and E(6, 3 2 , 5 2 , τ,τ ). The power-behaved terms correspond in the flat-space limit to perturbative string theory contributions ranging from genus zero to genus five. Although details of the perturbative sector of these functions are well-understood, we will not discuss the complete expressions for the instanton contributions in the Fourier expansion. Many details of these perturbative terms as well as the k-instanton/anti-instanton pairs have been determined and are presented in Appendix C. This data (and that from other instanton sectors) will be compared with terms arising in the analysis of the localization formula for the SYM free energy, and provides compelling evidence that the 1/N 2 coefficient is proportional to a particular rational linear combination of the two modular invariants mentioned above.
• The 1/N 3 contribution, which is the order where the flat-space interaction is (α ) 6 D 14 R 4 .
In this case we will see that the function H(3, τ,τ ) contains a rational sum of ten modular functions, These consist of the nine functions E(r, 3 2 , 3 2 , τ,τ ), E(r, 5 2 , 5 2 , τ,τ ) and E(r, 3 2 , 7 2 , τ,τ ), where r = 5, 7, 9, together with the function E(3, 3 2 , 3 2 , τ,τ ). The sum of these ten terms contains power-behaved (perturbative) contributions ranging from tree-level to genus-seven. Once again, these perturbative contributions of these functions are well understood but we have not completely analyzed the D-instanton contributions. However, we have obtained sufficient information of particular single Dinstanton contributions, as well as the contributions from instanton/anti D-instanton pairs to compare with the corresponding terms that are obtained from supersymmetric localization. This again provides compelling evidence that the 1/N 3 coefficient is proportional to a particular rational linear combination of the ten modular-invariant functions mentioned above.

Explicit formula
To be concrete, the explicit formula that we find with qualitative features described above where C 0 , C 1 are numerical constants that we will not determine here, while α r , β r , γ r are (2.14) We emphasize that we do not have a complete proof of (2.13). Instead, in the next section, we will provide abundant evidence for this expression by considering various limits of the terms in the 1/N expansion. In particular, we will study perturbative contributions in the zeroinstanton sector, as well as perturbative expansions around certain (n 1 , n 2 ) (anti)instantonpair backgrounds for |n 1 |, |n 2 | ≤ 3. We include a summary of the terms we have computed in Table 1.

Derivatives of deformed S 4 partition function
In this section, we will give evidence for the relations (2.8) and the large N expansion (2.13) by explicit evaluation of the various derivatives of (2.3), which take the form of expectation values in the familiar Hermitian Gaussian matrix model that describes the m = 0, b = 1 order perturbative (n 1 , n 2 ) instanton terms ∼ e 2πi(n 1 τ −n 2τ ) terms n 1 = 1, n 2 = 0 n 1 = 2, n 2 = 0 n 1 = 1, Table 1: A schematic summary of the evidence we have collected for (2.13). An entry g a YM e 2πi(n 1 τ −n 2τ ) in this table means we have computed the coefficient of a term proportional to g a YM e 2πi(n 1 τ −n 2τ ) in the expansion of (2.13). (The perturbative terms correspond to n 1 = n 2 = 0.) A dash means such terms are absent in the expansion. The cases where we have computed all non-zero terms are marked with "all." partition function We will compute these expectation values using topological recursion [87,88]. Since the application of this method to the N = 4 partition function was already explained in detail in [24,36], we will relegate the explicit calculations to Appendix D and only present the results below.

Perturbative sector
We begin by considering the perturbative part Z pert of (2.3) obtained by setting Z inst = 1.
Taking derivatives in m, b and using (2.4), we find where, as in [24], we define the two-and four-body expectation values which are taken with respect to the (b, m) = (1, 0) partition function (3.1). The terms that do not involve I and J come from the factor 3) that does not depend on the integration variables a i . The J terms in (3.2) take the same form for each expression, up to a factor of 3, so the nontrivial difference between the expressions in (3.2) comes from the I terms, which are also easier to evaluate.
From these expressions we see that the first and second relations in (2.8) are identically satisfied even before computing the expectation values, while the third relation follows from integration by parts after taking the expectation value. To justify the latter statement, note that from (3.2) we can write the perturbative contributions to the RHS of the third relation in (2.8) as Now that we have established the identities (2.8) at the perturbative level, let us examine the perturbative contributions to ∂ 4 m log Z pert m=0,b=1 , which include both a two-and four-body term (proportional to I and J , respectively). We can evaluate the two-body term in a large N and large λ expansion as in [36], which we then translate to large N and finite g YM by simply setting λ = N g 2 YM . It is harder to determine the large λ expansion of the four-body term, since as explained in [24], the dependence on the Fourier variables (w, ω) does not factorize and thus the analytic method in Appendix D of [33] cannot be easily applied. Instead, these terms were evaluated numerically to high precision in [24], and from a numerical fit it was possible to extract the first few terms in the 1/λ expansion. In Appendix D.1, we furthermore show that all terms of the form N 2 λ −integer can in fact be computed analytically, since the Fourier variables factorize for these terms. After combining these results and converting λ = g 2 YM N , we get where the f i (g YM ) denote functions that we have not yet been able to compute due to the aforementioned technical difficulties, while the leading small g YM term at any N −integer order can be computed analytically. The terms shown in (3.5) match with the expectation (2.13), after expanding the Eisenstein series using (2.10) and extracting the perturbative parts of the other modular functions shown in (2.12). Note that we have not computed enough perturbative terms to unambiguously fix the modular functions in (2.13), so we need to look at the instanton sector, which we will discuss next.

Instanton sector
We now consider the other parts of log Z in (2.3) that involve Z inst . In particular, we define the non-perturbative contributions to log Z by log Z NP ≡ log(Z/Z pert ). Taking derivatives contains expectation values of n > 2-body terms. As was the case with the J contribution in (3.2), Z appears in each expression with a coefficient proportional to the coefficient of the J term in (3.2). In deriving (3.7), we have used the fact that ∂ 2 b Z inst m=0,b=1 = ∂ 2 m Z inst m=0,b=1 (which we will explain momentarily from the structure of Z inst in (3.10)) to write Z purely in terms of derivatives with respect to m. The second line of (3.7) comes from a mixed term involving derivatives of both Z inst and the Υ factors in (2.3).
The explicit Nekrasov partition function Z inst (m, τ, b, a ij ) was computed in [80,81]. It 18 We emphasize here that Z inst is the U (N ) instanton partition function, which differs from the SU (N ) inst [11,12,81,89], inst can be completely absorbed by the counter-term ambiguities in (2.5), this difference does not affect the physical observables we compute.
can be written in terms of a sum over k-instanton sectors Z and consequently the first two relations in (2.8) hold identically before taking expectation values. Combined with the derivation of the perturbative terms in the previous subsection, this completes the proof of these first two relations for SU (N ) SYM at finite N and gauge coupling τ . 20 For the third relation in (2.8), as well as for computing ∂ 4 m log Z pert m=0,b=1 , we require the explicit expressions for G (k) (m, b, a). These are in general very complicated, so in the next few subsections we will consider just the lowest few values of k. We will compute the expectation values arising from these terms using both the small g YM expansion introduced in [35], as well as a more powerful finite g YM method that uses the full power of topological recursion. 19 The instanton partition function Z (k) inst receives contributions from each N -vector of Young diagrams Y whose total size is k, in the form of a rational expression in a ij , m and squashing parameters 1 = b, 2 = 1 b that is determined by the shapes of the Young diagrams in Y . The factor (m 2 + 1 4 (b − 1/b) 2 ) in (3.10) comes from a universal part of this rational expression for each Y that is due to the outer-corner entries of the non-empty Young diagrams. See Appendix B of [35] for details. 20 For SYM with general gauge groups, see Appendix B for relevant discussions.

One-instanton sector
We begin by considering the one-instanton sector of (3.6), which consists of terms proportional to e 2πiτ , and denote the overall coefficient by log Z NP, (1) with the corresponding derivatives in (b, m). 21 For the two-body (i.e. non-Z) terms in (3.6), we can simply set since the other terms in (3.7) can only contribute to higher instanton sectors. We use the explicit expression for Z (1) inst (m, b, a ij ) in [80,81]: , (3.13) which satisfies the general form in (3.10).
In order to check the third relation in (2.8), we consider the two-body term ( , for which contributions from Z cancel. We can evaluate this expectation value in a large N expansion at finite g YM using topological recursion as shown in Appendix D, and find which satisfies the third relation in (2.8) after comparing to τ 2 2 ∂ τ ∂τ ∂ 2 m log Z m=0,b=1 in (2.9) and extracting the 1-instanton term in the Eisenstein series (2.10).
We can then include Z by considering ∂ 4 m Z NP m=0,b=1 , which includes both the two-body term ∂ 4 m Z inst m=0,b=1 and the n > 2-body term Z (1) m=0,b=1 . It is harder to evaluate Z (1) m=0,b=1 at finite g YM , so instead we will compute it in a small g YM expansion by expanding for small eigenvalue a ij , computing the resulting n-body expectation values using topological recursion, and then performing the ω integral using a large λ expansion, as detailed in Appendix D. Finally, we add the contribution of the small g YM expansion of ∂ 4 m Z inst , which we also compute at finite g YM in the Appendix, to get the full result which matches the 1-instanton sector of the expected combinations of modular invariants in (2.13), which for the Eisenstein series is given in (2.10) and for the other modular functions is given in Appendix C.

Two-instanton sector
The two-instanton sector of (3.6) consists of terms proportional to e 4πiτ whose coefficients we denote by log Z NP, (2) with the corresponding (b, m) derivatives. For the two-body terms, inst , while for Z in (3.7) we have where the first term is new relative to the one-instanton case (3.12). We use the explicit expression for Z (2) inst (m, b, a ij ) in [80,81]: where the three terms correspond to the three vectors of Young diagrams Y given by the superscript (we have only listed the non-empty Young diagrams) and have the following explicit forms and .
which again satisfies the general form in (3.10).
Since the two-instanton expression is much more complicated than the one-instanton expression, we will only perform perturbative in g YM calculations. As shown in Appendix D, we find which satisfies the two-instanton sector of the second relation in (2.8), and (7)) 8388608 which matches the two-instanton sector of the expected modular function (2.13). Here, k denotes the term at order g 2 YM /N that we have not yet computed. 22 Note that the N −1 , N −2 , and N −3 contributions at zeroth order in g YM are new types of terms that did not appear in the one-instanton expression (3.15). The only class of terms that we have not checked so far are the instanton/anti-instanton pairs, which we will consider in the next subsection.

Instanton/anti-instanton sector
The mixed p instanton and q anti-instanton sector of (3.6) consists of terms proportional to e 2πi(pτ −qτ ) with coefficients denoted by the corresponding derivatives of log Z NP,(p,−q) , and only receives contributions from the first two terms in Z in Eq. (3.7): These terms therefore take the same form for all the (m, b) derivatives we consider, so they trivially satisfy the relations in (2.8), but they can be used to nontrivially check the formula for ∂ 4 m log Z m=0,b=1 in (2.13). Conveniently, we can use the expressions for ∂ 2 m Z (p) inst m=0,b=1 22 As explained in Appendix D, the g n YM for even/odd n are computed from different terms in the local-that were already computed for any p in [35]. For the (1, −1) sector, we compute the resulting expectation values in Appendix D in a large N expansion at finite g YM . The answer takes the form of a complicated integral that we write explicitly in Appendix D. We can evaluate this integral for any value of g YM , and find that it matches the relevant term in (2.13). We have also computed contributions from the (p, −q) sectors for p, q ≤ 3 in a small g YM expansion as shown in Appendix D. For instance, the (2, −2) term is while the other terms take a similar form and are given in Appendix D. All these terms agree perfectly with (2.13), which completes the check of that formula.

Four-point function in N = SYM
We will now apply the localization results of the previous section to constrain the four-point function SSSS of the stress tensor superconformal primary, which can also be constrained from its relation to the 10d IIB flat space graviton S-matrix. The superprimary S transforms in the 20 of the SO(6) R R-symmetry, and it can be represented as a traceless symmetric tensor S IJ ( x) with I, J = 1, . . . , 6 as SO(6) R fundamental indices. In order to avoid a proliferation of indices, it is customary to contract them with null polarization vectors Y I ization expression. In particular, the even n terms come from the two body terms ∂ 2 m Z (1) inst 2 m=0,b=1 and ∂ 4 m Z (2) inst m=0,b=1 , while the odd n terms come from Z (2) , and it turns out the latter is easier to compute to higher order in g YM .
Here, U ≡ We would like to study SSSS in the large c expansion at finite τ , which is related to the small momentum expansion of the IIB S-matrix at finite complexified string coupling τ s = χ s + i/g s . In this limit, it is convenient to use the Mellin transform [92,93] M of T , which is defined as [58]: where the coefficients α, β i , γ i , etc. are potentially non-trivial functions of (τ,τ ). The first term corresponds to tree-level supergravity, while M 1-loop is the regularized supergravity one-loop amplitude that can be found in [36] and will not be discussed here. We will instead focus mostly on the 1/c 7/4 , 1/c 9/4 , and 1/c 5/2 terms, which correspond to the R 4 , D 4 R 4 , and D 6 R 4 interaction vertices in type IIB string theory, respectively. At each order in 1/c, one can impose constraints on the coefficients α, β i , γ i , etc. by either comparing with the (super)graviton four-point scattering amplitude in type IIB string theory in the flat space limit or using the quantities (2.9) and (2.13) (or other similar quantities) derived from supersymmetric localization. Let us first discuss the constraints from the flat space scattering amplitude, and then those from supersymmetric localization.

Constraints from the flat space limit
The IIB four-point scattering amplitude of 10d gravitons and superpartners are restricted by supersymmetry to be proportional to a single function f (s, t) where A SG tree is the tree-level four-point supergravity amplitude, 23 s and t are the Mandelstam invariants. We will also define u ≡ −s − t. In turn, this function has an expansion at small momentum (more correctly, the expansion is for small values of the dimensionless product between momentum and the string length s ) of the form where the coefficient function that appears at each order in the expansion may be a nontrivial function of the complexified string coupling τ s . The functions f R 4 , f D 4 R 4 , and f D 6 R 4 can be written in terms of the modular functions introduced in Section 2 as [41][42][43][44]86] where the non-holomorphic Eisenstein series was defined in (2.10) and the other modular function was defined as the SL(2, Z)-invariant solution of the inhomogeneous equation (1.7).
The relation between the function f (s, t) in (4.5) and the Mellin amplitude (4.4) is given by the flat space limit formula [33] f (s, t) = stu 2 11 π 2 g 2 where κ > 0. 24 This relation, as well as the AdS/CFT dictionary allow us to fix the leading s, t terms in (4.4), such as (4.10)

Constraints from supersymmetric localization
As explained in [24,33], the localization quantities τ 2 2 ∂ τ ∂τ ∂ 2 m log Z m=0,b=1 and ∂ 4 m log Z m=0,b=1 impose constraints on SSSS integrated over S 4 . In the large c expansion, these constraints take the form where C 1-loop , C 1-loop are constants that depends on the precise form of the M 1-loop amplitude that we will not study here, and ∂ τ ∂τ log Z m=0,b=1 was given in (2.6). The right-hand sides of (4.11) are obtained by integrating the Mellin amplitude (4.4) with certain integration measures that produce the integrated correlators which are accessible by localization. For ∂ τ ∂τ ∂ 2 m log Z, the integration measure was first obtained in [33], and for ∂ 4 m log Z m=0,b=1 the measure was derived in [24]. Explicitly, the two integration measures are given in equations (2.15) and (2.16), respectively, of [24]. As for the left-hand sides of (4.11), we use the explicit localization results in (2.9) and (2.13). After converting the expressions in (2.9) and (2.13) into the 1/c expansion using c = (N 2 − 1)/4, the constraints (4.11) fix the coefficients to be (4.12) 24 When evaluating this integral, it is useful to note that Note that the values of α and β 2 match those computed from the flat space limit in (4.10), which is a non-perturbative in τ check of AdS/CFT to this order in 1/c. We can then combine the flat space limit and localization constraints to fix all the coefficients shown in (4.4) and obtain which is one of our main results. Note that we could not yet make use of the localization quantities involving derivatives of squashed parameter b, since the integrated constraints for those have not yet been derived.

Conclusion
Let us start with a summary of our results, and afterwards discuss several future directions.
In this paper, we have studied integrated correlators of four operators from the stress tensor Our first result was that among five possible (nontrivial) combinations of derivatives of Z with respect to (m, b, τ,τ ), only two are independent and can be taken to be τ 2 2 ∂ τ ∂τ ∂ 2 m log Z m=0,b=1 and ∂ 4 m log Z m=0,b=1 . The former was studied in [35] where strong evidence was presented that, beyond the leading term that scales as N 2 , this quantity has an expansion only in half-integer powers of 1/N whose coefficients are linear combinations of non-holomorphic Eisenstein series. (In Appendix D.2, we presented an alternative method that improved on the one in [35].). In this paper, we focussed on the other integrated correlator, ∂ 4 m log Z m=0,b=1 , for which the 1/N expansion contains both half-integral and integral powers of 1/N . Based on our computations, we conjectured that the coefficients of the half-integer powers of 1/N are again linear combinations of non-holomorphic Eisenstein series, while the coefficients of integer powers of 1/N are generalized Eisenstein series, which obey inhomogeneous Laplace eigenvalue equations of the form (1.7). In particular, the modular invariant coefficient at order 1/N is proportional to E(3, 3 2 , 3 2 , τ,τ ), the well-known coefficient of D 6 R 4 in the lowenergy expansion of the flat-space type IIB superstring amplitude. The terms at order 1/N 2 and 1/N 3 have coefficients that are linear combinations of generalized Eisenstein series with rational coefficients. See Table 1 where we summarized our evidence for our conjectures for the specific modular functions appearing in the 1/N expansion.
Lastly, in Section 4 we discussed the relation between the integrated correlators we computed and the type IIB superstring low-energy effective action on AdS 5 × S 5 which encodes scattering amplitudes of bulk (super)gravitons. As shown in [33], the separated point correlation functions of the same operators can be determined by general consistency conditions up to one, two, and three undetermined coefficients at orders N 1 2 , N − 1 2 and N −1 , respectively. At these orders, the correlators are determined from contact R 4 (at order N 1 2 ), 1 L 4 R 4 and D 4 R 4 (at order N − 1 2 ), and 1 L 6 R 4 , 1 L 2 D 4 R 4 , and D 6 R 4 (at order 1/N ) interaction vertices, and in Mellin space they asymptote to the flat space scattering amplitudes corresponding, respectively, to the R 4 , D 4 R 4 , and D 6 R 4 contact interactions in type IIB string theory. This information from the flat space limit combined with the integrated correlators ∂ τ ∂τ ∂ 2 m log Z m=0,b=1 or ∂ 4 m log Z m=0,b=1 allows us to uniquely determine the separated point correlators at orders N 1 2 , N − 1 2 , and N −1 . As we discussed, it is plausible that one of the relations (2.8) does not follow from superconformal symmetry and that it thus imposes an additional non-trivial constraint on the separated-point correlation function. If this is the case, then one would be able to determine the separated-point correlator at orders up to order 1/N , and from it derive the flat space scattering amplitude corresponding to the R 4 , D 4 R 4 , and D 6 R 4 contact interactions. These are precisely the terms that are also determined by supersymmetry in flat space.
The structure of the integrated correlators beyond order 1/N is worth highlighting. At these orders, the integrated correlators that can be computed using supersymmetric localization do not provide enough constraints to determine the separated-point correlators, so the Eisenstein and generalized Eisenstein series that we find do not completely characterize superstring scattering in AdS or in flat space. Nevertheless, they do represent supersymmetryprotected interactions in AdS 5 × S 5 . It is interesting to analyze their perturbative structure by examining the powers of g YM that appear in the zero-instanton expansion of the Eisenstein series in Eqs. (2.9) and (2.13). For terms of order N (1−2n)/2 for n = 0, 1, 2, . . . , which correspond to D 4n R 4 vertices in the bulk, the Eisenstein series have perturbative terms that can arise from up to genus-(n + 1) string worldsheets. For terms of order N −n for n = 1, 2, 3, . . . , which correspond to D 2(2n+1) R 4 vertices in the bulk, we find contributions of up to genus-(2n + 1). Interestingly, these features match previous observations about the type IIB S-matrix on flat space in [96,97]. Indeed, in [97] it was argued that if the duality between M-theory and string theory is naively assumed to hold exactly for all terms in the effective action, then one would conclude that the D 2k R 4 interaction vertices in type IIB string theory receive contributions from up to genus k worldsheets, in agreement with the observation we made above about the integrated four-point functions. It is believed, however, that the relations implied by the M-theory/string theory duality hold only for supersymmetry-protected interactions and are violated otherwise. 25 Thus, one expects that for the non-supersymmetric D 2k R 4 interactions with k > 3, both in flat space as well as their corresponding AdS Mellin amplitudes, there should be no restriction on the genus of the string worldsheets that contribute. However, the integrated correlators we study here are much simpler quantities than the full interaction vertex or the full Mellin amplitude, and, as mentioned above, these integrated correlators are supersymmetric. It is thus not entirely surprising that the arguments based on the M-theory/string theory duality seem to apply to them and restrict the genus of the worldsheet contributions in a manner consistent with our explicit computations in N = 4 SYM theory. Whether this observation is a coincidence or whether it can be made more precise are questions that we leave for future work.
Another future direction is the study of modular functions that appear in higher-point correlators. 26 Whereas the four-point correlators studied in this paper conserve the bonus U (1) Y symmetry of [8,9], n-point correlators may violate U (1) Y by a maximum of 2(n − 4) units. Maximal U (1) Y -violating n-point correlators of operators in the stress tensor multi- 25 We emphasize that this instance of M-theory/string theory duality (see [97] for details) assumes that the 11d supergravity description continues to be valid from large radius in the 11th dimension to small radius (to make contact with type IIA string theory) or from a large two-torus to small two-torus (to make contact with type IIB string theory). While this turns out to be true for BPS interactions protected by supersymmetry, in general non-perturbative M-theory effects (e.g. from M2 and M5 branes) that become large in these continuations cannot be ignored. 26 The five-point function of the stress tensor multiplet superconformal primary was considered in [98] in the supergravity approximation.
plet are holographically dual to type IIB superstring n-particle amplitudes that violate the U (1) R-symmetry maximally. The coefficients of terms in the large-N expansion of these correlators transform as SL(2, Z) modular forms with modular weights related to their U (1) charges. In a forthcoming paper [99] these correlators are determined up to order 1/N by a recursion relation analogous to the soft dilaton relations of flat-space superstring amplitudes. These relate the higher-point correlators to the four-point correlators determined in this paper and in [35], and make contact with the results in [46] concerning flat-space type IIB maximal U (1)-violating superstring amplitudes.
It would be interesting to construct n-point correlators that violate the bonus U (1) Y symmetry directly from the localization procedure by generalizing the analysis of this paper to cases where one takes more than four derivatives of the N = 2 * partition function. In this manner we would hope to determine expressions for the modular form coefficients to any order in the large-N expansion of the integrated n-point correlators.
An important loose end of our work is the proof of the last relation in (2.8) as well as determining whether or not any of these relations are consequences of supersymmetry (see Appendix B for evidence for these relations for general gauge groups). The Ward identities relating the four-point functions of various operators in the stress tensor multiplet were solved in [34]. 27 The derivatives in the relations (2.8) involving squashing are directly related to correlators of operators with spin in the stress tensor multiplet (namely the stress tensor, the R-symmetry current, as well as a rank-two anti-symmetric tensor operator), and one would have to use the Ward identity solution in [34] to relate such correlators to those of the stress tensor multiplet superconformal primary. It would be very valuable to perform this analysis, because it could have applications beyond 1/N perturbation theory, for instance in numerical bootstrap studies.
Lastly, let us point out that the large N expansion in this paper is asymptotic, as can be seen already from the all orders in 1/N expressions for ∂ τ ∂τ ∂ 2 m Z m=0,b=1 in [36], and so is expected to have exponentially small in N corrections. In the bulk, we speculate that these exponentially-suppressed corrections can be interpreted as boundary-anchored strings and branes in AdS 5 ×S 5 . It would be interesting to understand these contributions as well as their dependence on (τ,τ ). More generally, one might hope that our 1/N expansion supplemented by these exponential corrections can be resummed into a finite-N modular function. For the perturbative terms in g YM , the finite-N integrated correlators were computed using the method of orthogonal polynomials for both ∂ τ ∂τ ∂ 2 m Z m=0,b=1 [36] and ∂ 4 m Z m=0,b=1 [24], but such an analysis would be more challenging for the instanton terms. We nevertheless hope to come back to these issues in the near future.

Acknowledgments
We As shown in [100], with the chiral superfields W αβ , Φ, and A, there are three classes of composite fields F of chiral weight −2 and Weyl weight 2, corresponding to three counter-terms in (A.1). Here Φ is an auxiliary chiral superfield that can be identified with the compensating vector multiplet in N = 2 supergravity [101], and T ∝D 4 is an antisymmetric combination of the four anti-chiral superspace covariant derivatives [100]. Importantly the corresponding counter-term does not depend on the value of Φ when A take constant values [100] which is the case for the deformations considered here.
When evaluated on the supersymmetric mass deformed background [101], the first term

B Comments on relations between fourth derivatives of log Z(b, m) for general gauge groups
Here we provide some further evidence for the three relations (2.8) between various fourthderivatives of the SYM free energy with respect to the mass and squashing deformations as well as the complexified gauge coupling (τ,τ ) for N = 4 SYM with general gauge groups.
We will also argue for the first relation in (2.8) based on superconformal Ward identities.
The N = 2 * partition function of SYM with a general gauge group G on a squashed sphere is given by 28 28 We emphasize again that the factors of Υ b (0) that have been missing the previous works (e.g. [26]) carry nontrivial b dependence, where γ is the Euler's constant, and are thus crucial to produce the correct CFT free energy in the presence of squashing deformations.
where r denotes the rank of G, W is the Weyl group, ∆ is the set of roots and (·, ·) defines the standard Killing form on the Lie algebra g. The instanton contributions at the two poles of S 4 are captured by the factor Z inst (m, τ, b, a) = 1 + k≥1 q k Z (k) inst (m, b, a) and its conjugate respectively, whose explicit forms are available for G = SU (N ) and are used extensively in the following sections. For more general classical Lie groups of BCD types, (m, b, a) admits a contour integral expression at each instanton number k thanks to the ADHM construction of the instanton moduli space and an equivariant localization procedure thereof [80,81,[103][104][105]. The instanton contribution to Z for exceptional Lie groups is still an open question, though by the AGT correspondence they are related to torus one-point blocks of the corresponding W-algebras (with twist for the non-simply laced cases) [11,106].
We note the following simple properties of Z The first equality is due to m being the mass parameter for an SU (2) 29 To verify this for the perturbative contributions, the following identities of the Upsilon function is useful, Let us define the one-loop contribution in (B.2) from a pair of root vector α and its Weyl reflection −α as

This implies that
, (B.6) with z = α(a) ∈ R. Then by using the integral expression of the Upsilon function (2.4), one Note that the Upsilon function Υ b (iz) has a simple zero at z = 0. Consequently, the relations (B.7) continue to hold with H(b, m, z) replaced by .

(B.8)
Putting them together, we conclude and thus we verify the first two relations in (2.8) for the perturbative contributions in (B.2) (before [d r a] integral). 32 In the main text, we have further proved the first two relations of (2.8) non-perturbatively for G = SU (N ) at finite N by using the explicit form of the instanton partition function which can be found in Appendix B of [35].
Concerning the last relation of (2.8), we have verified it perturbatively for G = SU (N ) at finite N using integration by parts as in (3.4). The derivation extends trivially to general gauge group G after replacing I(w) by  H(b, m, z). 32 Said differently, in a weak coupling expansion, Z pert captures the perturbative contributions to the full Furthermore, we have provided evidence for this relation of (2.8) in the main text at the non-perturbative level for G = SU (N ).
Let us now comment on (2.8)  harmonic ambiguities (2.5) in (τ,τ ) that can be removed by taking ∂ τ ∂τ derivatives. Since the mass and squashing couplings are introduced in a theory-independent way in the localization setup [26], we conclude that ∂ τ ∂τ ∂ 2 m log Z m=0,b=1 and ∂ τ ∂τ ∂ 2 b log Z m=0,b=1 must be proportional up to a theory independent constant. A quick calculation in the abelian SYM theory which has partition function 33 confirms that this proportionality constant is one and thus the desired relation follows.

C Solutions of inhomogeneous Laplace equations
In this appendix we will describe some properties of the generalised Eisenstein series that satisfy equations of the form (1.7) that arise in the coefficients H(q, τ,τ ) of even terms in the 1/N expansion (2.11) up to order 1/N 3 . The function E(3, 3 2 , 3 2 , τ,τ ) (the coefficient of SYM partition function Z. Thus we have verified the first two relations in (2.8) up to instanton effects. 33 Here we have included the abelian instanton contributions for completeness [89], though they do not affect the physical observables that come from fourth derivatives of the SYM free energy. the D 6 R 4 interaction in flat-space type IIB superstring theory) was completely determined in [86], Certain properties of more general functions satisfying (1.7) were presented in [47] but these were mainly restricted to the perturbative terms, whereas we are here also interested in detailed properties of the D-instanton terms for the specific functions appearing in the 1/N expansion.
We will first review the structure of this modular invariant based on the solution of the Laplace equation [44] 34 Following [86], this equation may be solved in terms of its Fourier modes defined by It is important to understand the boundary conditions imposed on the Fourier modes that are necessary in order for the complete function to be SL(2, Z) invariant. According to a theorem proved in [86]: SL(2, Z) invariance of a function that grows as τ a 2 as τ 2 → ∞ implies that its Fourier modes are bounded by τ 1−a 2 in the τ 2 → 0 limit.
Since E(r, s 1 , s 2 , τ,τ ) = Writing the source term on the right-hand side of (C.1) as the kth mode satisfies the equation The general solution to the above differential equation can be found in [86], where it is expressed as the sum of a particular solution and a solution of the homogeneous equa- 34 The overall normalisation of the source term on the right-hand side of this equation has been arbitrarily set equal to −1. Since E(3, 3 2 , 3 2 , τ,τ ) is the coefficient of the D 6 R 4 interaction in the low energy expansion of tion. F k (τ 2 ) = F P k (τ 2 ) + F H k (τ 2 ). The coefficient of the homogeneous solution is uniquely determined by imposing the τ 2 → 0 boundary condition described above.
The particular solution that was determined in [86], was written in the form F P k (τ 2 ) = k 1 It is useful to consider the solutions in several sectors: (a) k 1 = −k 2 (so k = 0); (b) These terms contribute to the k = 0 mode. The term with k 1 = k 2 = 0 is a sum of powers of τ 2 that is given by The first three terms in this expression originate from f P 0,0 (τ 2 ) and are easily obtained by equating the coefficients of the powers of τ 3 2 , τ 2 and (τ 2 ) −1 on both sides of (C.5). The term βτ −3 2 is a solution of the homogenous equation, and its coefficient β was determined in [44] by multiplying both sides of (C.5) by E(4, τ,τ ) and integrating over a fundamental domain of SL(2, Z). A detailed analysis can be found in [44], which leads to β = 4ζ(6)/27 so that In addition to the (0, 0) term, the k = 0 mode receives contributions from a sum over an infinite number of terms with k 1 = −k 2 = 0, which represent D-instanton/anti D-instanton pairs. These terms are bilinear in K-Bessel functions and are given by the flat-space type IIB superstring action its normalisation is simple to fix from the string theory scattering amplitude.
where the coefficients q i,j 3 are given by Making use of the weak coupling (τ 2 → ∞) expansion of the K-Bessel functions, we see that the expression (C.8) is suppressed by a factor proportional to e −4π|k 1 |τ 2 , which is characteristic of an instanton/anti-instanton pair.
The complete zero mode is given by . In order to check the small-τ 2 boundary condition we note that in the small-τ 2 limit where we have used the Ramanujan identity Using (C.6), (C.7) and (C.11) we see that F 0 (τ 2 ) = F P 0 (τ 2 ) + F H 0 (τ 2 ) = O(τ −2 2 ), which is the required boundary condition.
This depends only on the sum of the source mode numbers, k 1 + k 2 , and α k is determined by imposing the boundary condition at τ 2 = 0. It turns out that α k = 0 for k = 0 (although this is a property of the solution that was not noticed in [86]). Since F H k (τ 2 ) = 0 for k = 0 the solution for mode k is identified with the particular solution, F k (τ 2 ) = F P k (τ 2 ), so we can drop the superscript P in the following.
(b) k 1 = k = 0, k 2 = 0 and k 2 = k = 0, k 1 = 0 These modes are solutions of (C.5) given by f k,0 (τ 2 ) = f 0,k (τ 2 ) = 8 σ 2 (|k|) 9 π |k| 3 × q 0 k,0 (π|k|τ 2 )K 0 (2π|k|τ 2 ) + q 1 k,0 (π|k|τ 2 )K 1 (2π|k|τ 2 ) , where the coefficients are given by (C.14) Using the weak coupling expression (C. 10) we see that f k,0 (τ 2 ) has the form of a charge-k D-instanton contribution with a characteristic e −2π|k|τ 2 suppression factor together with an unlimited number of perturbative corrections (powers of τ −1 2 ). In order to find the complete expression for the kth mode we need to sum an infinite number of terms of the form F k (τ 2 ) = In these cases the solution of (C.5) was found in [86] to have the form where q i,j k 1 ,k 2 (z) are specific polynomials with powers of z ranging from z to z −2 which we will not display here. There are two distinct cases to consider: The solution contains D-instantons of charge k = k 1 + k 2 > 0, or anti D-instantons when k < 0. These are characterized by an exponentially suppressed behavior of the form e −2π|k|τ 2 in the τ 2 → ∞ limit.
This is again a contribution that has total D-instanton charge equal to k, but the solution describes a D-instanton of charge k 1 together with an anti D-instanton of charge k −k 1 (when we assume that k 1 > 0). The large-τ 2 behavior is characterized by an exponential suppression factor of e −2π(|k 1 +|k 2 |)τ 2 . Since |k 1 | + |k 2 | > |k| these terms do not contribute to the leading exponential behavior in the weak-coupling limit.
Given the complete solution it is straightforward to compare with our analysis of the 1/N terms in the large-N expansion of the localized integrated correlator that are determined in Section 3 and Appendix D. For example, among these terms there are certain perturbative contributions to the leading exponential dependence in the k = 1 D-instanton term, which is the sum of the (1, 0) and (0, 1) components (the k = 1 components (2, −1), (3, −2), . . . are exponentially suppressed relative to the leading term). In particular, the first few terms in the perturbative expansion around the k = 1 D-instanton contribution matches terms in the expansion of F 1 (τ 2 ) = f 1,0 (τ 2 ) + f 0,1 (τ 2 ), which is given by where τ −1 2 = g 2 YM /4π. Similarly, the leading exponential contribution to the k = 2 mode gets contributions from the sum of the (1, 1), (2, 0), and (0, 2) components. The perturbative expansion around the leading exponential dependence of the k = 2 contribution should therefore match the expansion of F 2 (τ 2 ) = f 2,0 (τ 2 ) + f 0,2 (τ 2 ) + f 1,1 (τ 2 ), which has the form (C.17) An intriguing aspect of this expansion is that it contains a sum of odd powers of 1/τ 2 ∼ g YM that come from the expansion of f 2,0 (τ 2 ) + f 0,2 (τ 2 ) and even powers of 1/ √ τ 2 that come from the expansion of f 1,1 (τ 2 ).
A note on the rôle of the τ 2 = 0 boundary condition.
An expression such as (C.17) is uniquely determined by a large-τ 2 expansions of the exact solution. The uniqueness is associated with the fact that the solution is valid for all τ 2 and builds in the τ 2 = 0 boundary condition. It is important to stress that simply solving (C.5) for the (k 1 , k 2 ) component (with k 1 + k 2 = k) by means of a perturbation expansion around the large-τ 2 limit does not determine the expansion coefficients uniquely since such an expansion does not address the boundary condition at τ 2 = 0. This is reflected in the arbitrary coefficients of the solutions of the homogeneous equations for the (k 1 , k−k 1 ) sectors, . It is enlightening to illustrate this by considering the perturbative solution around the k = 2 D-instanton that is defined by the solution of (C.5) with k 1 + k 2 = 2. We have seen that e 4πτ 2 (f 2,0 (τ 2 ) + f 0,2 (τ 2 )) has an expansion in half-integer powers of τ −1 2 . However, the perturbative expansion of e 4πτ 2 f H k 1 ,k 2 (τ 2 ) (where f H k 1 ,k 2 (τ 2 ) satisfies (C.5) with the source term set to zero), is in integer powers of τ −1 2 . Therefore, there is no ambiguity in the perturbative solution in the (2, 0) + (0, 2) sector. However, the perturbative expansion of e 4πτ 2 f 1,1 (τ 2 ) (the (1, 1) sector) is in integer powers of τ −1 2 , and so f 1,1 (τ 2 ) mixes with the expansion of the solution of the homogeneous equation, f H 1,1 (τ 2 ). More explicitly, we can extract the k = 2 instanton power behavior by setting k 1 = k 2 = 1 in (C.5) and writing The function C(τ 2 ) satisfies a differential equation that is easily solved to any given order in a perturbation expansion in τ −1 2 , but the solution has one arbitrary constant that cannot be determined without additional information. Since (C.17) is an expansion of the exact solution, it builds in the τ 2 = 0 boundary condition, so there is no ambiguity involved in comparing (C.17) with the results of the localization calculation in the main text.
Other terms of order 1/N were also determined from the localized correlator in section 3 and appendix D. In particular, the analysis of in appendix D determines the exact expression for the (1, −1) component of the k = 0 mode rather than simply its perturbative expansion. C.2 E(4, 3 2 , 5 2 , τ,τ ) and E(6, 3 2 , 5 2 , τ,τ ) The modular functions E(4, 3 2 , 5 2 , τ,τ ) and E(6, 3 2 , 5 2 , τ,τ ) satisfy the inhomogeneous Laplace equations, From here on we will write the Fourier expansion of E(r, s 1 , s 2 , τ,τ ) using the notation E(r, s 1 , s 2 , τ,τ ) = k F r,s 1 ,s 2 k (τ 2 )e 2πikτ 1 = k 1 k 2 f r,s 1 ,s 2 k 1 ,k 2 (τ 2 )e 2πi(k 1 +k 2 )τ 1 , (C.20) 35 In this notation the modes f k (τ ) in the previous subsection would be denoted f where we have introduced the superscripts r, s 1 , s 2 to indicate the eigenvalue and the source term. 35 We have not determined the complete solution in these cases, but we have determined many features that can be correlated with the results of the 1/N 2 contribution to the large-N expansion of the localized integrated correlator. Some of these properties of the solutions are summarised as follows.
We have determined the complete zero-instanton sector, which again includes terms that are power-behaved in τ 2 , and (k 1 , −k 1 ) D-instanton/anti D-instanton terms.
The power-behaved terms (k 1 = k 2 = 0) with eigenvalues r = 4, 6 are given by where the coefficients q i,j r are symmetric, q i,j r = q j,i r . We find that the expressions for q i,j 4 (z) for r = 4 are given by whereas for the case r = 6, they are given by q 2,2 6 (z) =

(C.24)
Once again it is important that the τ 2 → 0 boundary condition is satisfied.
(ii) k = 0 We have determined perturbative expansions around the leading exponential behavior in various instanton sectors that are needed in order to compare with the results obtained from the 1/N 2 contribution to the localized integrated correlator discussed in Section 3. We will discuss the explicit perturbative expansions around the leading exponential behavior of the charge k = 1 and k = 2 D-instantons to the first few orders in powers of 1/τ 2 .
For the one-instanton contributions, i.e. the (1, 0) + (0, 1) components, the expansions of (C.25) For the case of k = 2, the leading exponential contributions to F r, 3 2   C.3.1 E(r, 3 2 , 3 2 , τ,τ ) We described this function with r = 3 in detail in Section C.1 in order to compare with the coefficient of the 1/N term in the large-N expansion of the localized correlator. In the r = 5, 7, 9 cases the perturbative terms (the terms power behaved in τ 2 ) are given by The coefficients are given by the following polynomials: • r = 5 11π . (C.29) • r = 7 . (C.30) • r = 9 .
The k = 0 terms with r = 5, 7, 9 that are power-behaved in τ 2 are given by (C.32) The non-perturbative k = 0 terms are given by the sum of (k 1 , −k 1 ) D-instanton/anti Dinstanton contributions that takes the following form The coefficients in this equation are given by, the following polynomials.

D Topological recursion
In this appendix we will show the details of the localization calculations whose results were discussed in the main text. All of these calculations involve computing expectation values with respect to the m = 0, b = 1 free gaussian matrix model in (3.1). In fact, as explained in [36], if an expectation value only depends on the difference of eigenvalues a ij , as all the ones we consider do, then we can equivalently take the expectation value with respect to the where we now integrate over N eigenvalues with no constraint, unlike the SU (N ) matrix model in (3.1). In the following we will for simplicity take all expectation values with respect to (D.1). We will then compute these expectation values using topological recursion, which we will briefly review following [36].
Let us begin by defining the n-point operator The expectation value of this operator with respect to (D.1) can be shown to obey recursion relations in n and 1/N , which are called topological recursion. It is customary to write down these recursion relations in terms of the connected correlators W n (y 1 , . . . , y n ) ≡ N n−2 R n (y 1 , . . . , y n ) conn = N n−2 which in a slight abuse of notation we will refer to as resolvents. These resolvents can then be expanded in 1/N 2 as W n (y 1 , . . . , y n ) ≡ and each genus-m term W n m can be computed for finite λ using a recursion formula in n, m [87,88] starting with the base case W 1 0 , as described e.g. in [36]. We use resolvents up to n + m ≤ 5, which were given in Mathematica files attached to [24,36], except one should set y 2 i − λ/(4π 2 ) → y i 1 − λ/(4π 2 y 2 i ) in all expressions given there, so that the resolvents have the correct properties as y → ∞. In the following subsections, we will relate the expectation values we are interested in to these resolvents, which allows us to compute their 1/N expansion.

D.1 Details of perturbative calculation
The goal of this subsection is to compute (3.5) starting from the expectation values in (3.2).
We start by reviewing the calculation of [24,36], where the former computed the two-body operator I(ω), and the latter computed the four-body operator J (ω, w).
where we define We then take the inverse Laplace transform in (D.7) of the explicit resolvents to get the 1/N 2 expansion at finite λ for I(ω) and J (ω, w) in terms of integrals over the Fourier variables w, ω shown in (3.2). For instance, at leading order in 1/N 2 we need only consider the genus-zero resolvents in J 0 (ω, w) and I 0 (ω), which give (D.10) We can then plug these expressions into (3.2) to get the leading order in N 2 result at finite λ: and the higher order in 1/N 2 terms take a similar form of integrals of two Bessel functions for the 2-body terms, and four Bessel functions for the 4-body terms. We need to take the large λ expansion of these results, which will correspond to the large N expansion after we set λ = g 2 YM N . As described in Appendix D of [33], the first step is to express products of Bessel functions in their Mellin-Barnes form c+∞i c−∞i ds Γ(−s)Γ(2s + µ + ν + 1) 1 2 x µ+ν+2s Γ(s + µ + 1)Γ(s + ν + 1)Γ(s + µ + ν + 1) . (D.12) For the two-body terms, we can them perform the resulting integrals over ω in (3.2) using the identity ∞ 0 dω ω a sinh 2 ω = 1 2 a−1 Γ(a + 1)ζ(a) . (D.13) After doing these ω integrals, we can then do the s integral in (D.12) by closing the contour to the left, which gives an expansion in 1/λ.
For the 4-body term in (D.11), we can now apply (D.12) twice to get 14) where note that the w, ω dependence does not factorize due to the w 2 −ω 2 in the denominator.
While in general it is difficult to compute the s, t integrals by closing the contour to the left, since there are likely poles that can only be seen after doing the w, ω integrals, for the poles at s, t = − 3 2 , − 5 2 , . . . we find that the residues at each order in λ factorize in w, ω. The w, ω integrals can then be computed with (D.13) analytically continued to negative even integers (recall that this quantity is only divergent for a = 1, −1, −3, −5, . . . ). These poles correspond to the N 0 λ −integer terms discussed in the main text, which is why we can compute all of them analytically. Unfortunately, this factorization after taking poles does not apply to all the expected large λ terms, such as the N 0 λ − integer 2 terms that we know to exist from the numerical results of [24], nor does it apply to any terms at higher orders in 1/N 2 . D.2 τ 2 2 ∂ τ ∂τ ∂ 2 m log Z m=0,b=1 at large N and finite g YM Before we discuss the instanton sector contribution to the relations in (2.8) and ∂ 4 m log Z m=0,b=1 in (2.13), we first introduce a new large N and finite g YM method that we will use for these calculations, by demonstrating it in the simpler case of τ 2 2 ∂ τ ∂τ ∂ 2 m log Z m=0,b=1 in (2.9). This result was previously computed in [35] to the first couple orders in 1/N at finite g YM , and at subsequent orders in 1/N in a small g YM expansion. Here, we complete this derivation by computing all orders in 1/N at finite g YM .

D.2.1 One-instanton sector
We start by considering the one-instanton contribution to τ 2 2 ∂ τ ∂τ ∂ 2 m log Z m=0,b=1 , and for simplicity we will consider just ∂ 2 m log Z m=0,b=1 , since the τ derivatives can be trivially applied to the result.
For this calculation, it is useful to express Z (1) inst (m, b, a ij ) in (3.13) as a contour integral where the integration contour is counter-clockwise around the poles at z = a j + i, and the subtraction of 1 from the integrand does not contribute to the final result, but makes the integrand decay as 1/z 2 at |z| → ∞, so that the contour can be taken to be the real line.
We can then take the m derivatives to get where we define 17) in terms of the resolvent operator R given in (D.2). We now take the expectation value, and use the cumulant expansion A m conn m! (D. 18) to get These integrals can then be performed as described in [35] to get (D. 24) We can then take the τ derivatives and compare to the one-instanton term in (2.9).

D.2.2 Higher instanton sector
We can similarly compute the k > 1 instanton terms. As described in [35], these instantons are described by rectangular Young diagrams of height p and length q, which will correspond to the partition of unity in the divisor sum that defines the Eisenstein series. Following [35], we thus define for integers p, q such that k = pq. This I p×q was given in [35] as 1)) .

(D.26)
We can then expand at large N and perform the integrals similarly to the one-instanton case to get We can then take the τ derivatives, take the p, q sum in (D. 25), and compare to the relevant instanton term in (2.9), which is the complete finite g YM derivation of this result to any order in 1/N .

D.3 Details of instanton calculation
We now continue with the calculation of the expectation values that show up in the relations (2.8) and ∂ 4 m log Z m=0,b=1 in (2.13), and address the instanton terms. For some of these calculations, we will use the the large N and small g YM method introduced in [35], while for others we will use the new large N and finite g YM method that we demonstrated in the previous section. We follow the main text and discuss the one-instanton sector, then the two-instanton sector, and finally the mixed instanton/anti-instanton sector.

D.3.1 One-instanton sector
We start by detailing the large N and finite g YM calculation of (3.14). Consider the contour integral representation of Z (1) inst (m, b, a ij ) given in (D.15). We can then take derivatives in m, b to get ∂ 4 m Z (1) inst m=0,b=1 = 24 π dz e Q(z) (∂ z R(z) − 1) − 1 , where R is the resolvent operator given in (D.2), and Q was defined in (D.17). We then take the expectation value and use the cumulant expansion (D.18) to get where we introduced the derivatives of s to put all terms in (D.31) into the exponential.
From (D.17), we see that this expression is written in terms of connected correlators of R, i.e. resolvents with the known 1/N 2 expansion described in previous sections. We can then expand (D.32) at large N and perform the integrals, just as in Section D.2, to get which we combine to get (3.14).
Next, we compute ∂ 4 m Z NP m=0,b=1 , which consists of the two-body term ∂ 4 m Z (1) inst m=0,b=1 computed above as well as the higher-body term Z (1) in (3.12). For Z (1) , it is difficult to perform the large N and finite g YM calculation due to the e 2ωa ij terms and the Fourier integral over ω. Instead, we will perform a large N and small g YM expansion by expanding ∂Z (1) inst m=0,b=1 at small eigenvalue, which corresponds to small g 2 YM , to get an infinite series of n-body terms. We will then compute their expectation value with e 2ωa ij in a large N expansion at finite λ using topological recursion, and then do the large λ expansion as we did with the perturbative terms of Section D.1. After setting λ = g 2 YM N , these steps give a consistent large N and small g YM expansion.
We start by expanding ∂ 2 m Z inst m=0,b=1 in (3.13) at small eigenvalue to get 3 (N )C 2 C 4 + f

D.3.2 Two-instanton sector
The calculation in the two-instanton sector is similar to the one-instanton sector, except all the expressions are much more complicated, so we only do calculations in the large N and small g YM expansions. For the two-body terms ∂ 4 m Z inst m=0,b=1 and ∂ 2 m ∂ 2 b Z inst m=0,b=1 , we expand to leading order in small eigenvalue to get N -dependent coefficients that satisfy complicated recursion relations, similar those found at k > 1 instantons in the small eigenvalue expansion of ∂ 2 m Z inst m=0,b=1 in [35]. We can then expand these recursion relations at large N and perform the trivial expectation value (since their is no eigenvalue dependence to leading order) to get , which satisfies the second relation in (2.8) for the two-instanton sector. For the higher body term Z (2) , we note that the first term in (3.16) can be computed to leading order in g 2 YM by simply squaring the leading order expression in (D.35) to get .