Integrated Correlators in $\mathcal{N}=4$ SYM via $SL(2,\mathbb{Z})$ Spectral Theory

We perform a systematic study of integrated four-point functions of half-BPS operators in four-dimensional $\mathcal{N}=4$ super Yang-Mills theory with gauge group $SU(N)$. These observables, defined by a certain spacetime integral of $\langle\mathcal{O}_2\mathcal{O}_2\mathcal{O}_p\mathcal{O}_p\rangle$ where $\mathcal{O}_p$ is a superconformal primary of charge $p$, are known to be computable by supersymmetric localization, yet are non-trivial functions of the complexified gauge coupling $\tau$. We find explicit and remarkably simple results for several classes of these observables, exactly as a function of $N$ and $\tau$. Their physical and formal properties are greatly illuminated upon employing the $SL(2,\mathbb{Z})$ spectral decomposition: in this S-duality-invariant eigenbasis, the integrated correlators are fixed simply by polynomials in the spectral parameter. These polynomials are determined recursively by linear algebraic equations relating different $N$ and $p$, such that all integrated correlators are ultimately fixed in terms of the integrated stress tensor multiplets in the $SU(2)$ theory. Our computations include the full matrix of integrated correlators at low values of $p$, and a certain infinite class involving operators of arbitrary $p$. The latter satisfy an open lattice chain equation for all $N$, reminiscent of the Toda equation obeyed by extremal correlators in $\mathcal{N}=2$ superconformal theories. We compute ensemble averages of these observables and analyze our solutions at large $N$, confirming and predicting features of semiclassical AdS$_5\, \times$ S$^5$ supergravity amplitudes.


Introduction
It is reasonable to hope that the four-dimensional maximally-supersymmetric Yang-Mills theory may be exactly solved. In the large N 't Hooft limit, the theory becomes famously integrable, and Yangian symmetries emerge, which has led to the exact determination of various planar observables. At finite N , these special properties are obscured; however, the theory enjoys a non-perturbative S-duality symmetry [1][2][3][4], which for simply-laced gauge groups is an invariance (up to global identifications [5]) under SL(2, Z) transformations of the complexified gauge coupling, (1.1) It seems likely that a solution of this theory (henceforth N = 4 SYM) will be possible only with judicious treatment of S-duality. Even so, determining the exact τ -dependence of N = 4 SYM observables should be very hard in general, precisely because the theory is strongly interacting over most of its parameter space, even after modding out by SL(2, Z).
There is a special family of observables in N = 4 SYM, introduced in [6], which sit askew of this expectation: they vary as functions of τ , and yet are determinable by supersymmetric localization. They are integrated four-point functions of half-BPS operators, whose definition will be recalled further below. In this paper, following [6][7][8], we find strikingly compact solutions for an infinite class of these integrated correlators in the SU (N ) theory, exact in N and τ , and study their general systematics. Tantalizingly, they are optimally understood in a complete basis of SL(2, Z)-invariant functions [8]: in this basis, the integrated correlators are essentially polynomials in the relevant spectral parameter.
To set the context for these results, let us recall the status of (unintegrated) N = 4 SYM four-point functions. We henceforth specialize to the SU (N ) theory. Half-BPS superconformal primary operators in N = 4 SYM, call them O (i) p with p ≥ 2, are spacetime scalars in the [0 p 0] representation of the SU (4) R-symmetry with dimension ∆ p = p. The case p = 2 is the 20 of the stress tensor multiplet. At p > 3 there is degeneracy, which is captured by the label i. Consider the family of four-point functions p . At finite N , this is known only in weak-coupling perturbation theory (g 2 YM 1) through two loops in general [9], and three loops for p = 2 [10] (though integrands are known to four loops [11,12]). In the planar limit, it is known through three loops at weak coupling [13] (with their integrands being known to five-loop order [14] and even through ten loops for p = 2 [15]), while at strong coupling holographic and bootstrap approaches determine the leading strong coupling limit [16][17][18][19][20][21][22][23], as well as some sub-leading 1/λ [6,[24][25][26][27][28][29] and 1/N [30][31][32][33][34][35][36][37][38] corrections. Knowing any of the exactly in the planar limit would be tantamount to knowing the holographic dual of the AdS 5 × S 5 Virasoro-Shapiro amplitudethat is, the four-point scattering of AdS 5 × S 5 supergravitons in classical string theory -and, upon developing its operator product expansion, would furnish a substantial part of the solution of planar N = 4 SYM. Unfortunately, due to the inherent difficulty of intermediate coupling, what is known about this object is essentially perturbative around either weak or strong coupling.
And yet, rather astoundingly, if one integrates these four-point functions against a fairly vanilla spacetime measure, these correlators can be determined exactly and, as we will see, simply. In particular, the integrated correlators of interest here are defined as follows [6]: where H (N |i,j) p (u, v; τ ) is the "dynamical" part of the correlator not fixed by superconformal Ward identities (recalled explicitly in Section 2.1 below). The magic is that G (N |i,j) p (τ ) are computable from derivatives of a free energy on S 4 deformed by sources, which in turn is determined by supersymmetric localization. The aim is then to find simpler expressions upon starting from the localization integrals. In [6], these were computed exactly in the 't Hooft limit for all p, yielding an elegant integral formula for the integrated AdS-Virasoro-Shapiro amplitude. Our goal here is to study these integrated correlators at finite N , for all p.
For p = 2, this was done in [7], who found a conjecturally exact, recursive solution at finite N , an elegant result subjected to many checks. (See [39][40][41][42] for subsequent related work.) That solution was elucidated in [8] from another point of view, one which will be crucial for this paper: namely, that of the SL(2, Z) spectral decomposition. To introduce this briefly, any SL(2, Z)invariant observable of N = 4 SYM admits, by way of finiteness of free N = 4 SYM and the spectral theorem for SL(2, Z), a decomposition into a complete, SL(2, Z)-invariant eigenbasis of the Laplacian on the upper-half τ -plane. There are three branches of eigenfunctions: a continuous branch (the non-holomorphic Eisenstein series, E * s (τ ), with Re s = 1 2 ), a discrete branch (the Maass cusp forms, φ n (τ )), and a constant. The overlap with the constant term is equal to the ensemble average, with respect to the Zamolodchikov metric, over the space of N = 4 theories parameterized by τ . This is a powerful tool in the N = 4 SYM context because the complete eigenbasis "decouples" the τ -dependence from the core information, not made redundant by S-duality, that characterizes a given observable. The entire content of an observable is characterized by its overlaps with these basis elements: τ -space is traded for "spectral space." As in any other physical context where a complete basis is identified, it gives a systematic way to understand the possible τ -dependence of SL(2, Z)-invariant N = 4 SYM observables, and, therefore, to classify these dependences according to the functional complexity of the overlaps. 1 Applied to the p = 2 integrated correlator, G (N ) 2 (τ ), [8] found that the entire integrated correlator is determined by a single polynomial. First, the cusp form overlap vanishes. Second, the Eisenstein overlap is given by a universal gauge-theoretic prefactor times an N -dependent polynomial in the spectral parameter s; this polynomial may be determined solely from the perturbative expansion of G (N |i,j) p (τ ) to finite order. These polynomials solve a short recursion relation in N , equivalent to the recursion originally found by [7]. Finally, the ensemble average, which is simply (one half times) the aforementioned polynomial evaluated at s = 1, is also just a polynomial in N . If one prefers to restore the τ -dependence, the result is where the full information (including the additive constant) is contained in the polynomial f (N ) 2 (s). A consequence of this structure is that the entire G (N ) 2 (τ ) is fixed by the sector of zero total instanton number. The origin of the recursion remains mysterious. However, what the SL(2, Z) decomposition makes clear is that the integrated correlator is about as simple of a non-trivial τ -dependent observable as one could hope for.
In this paper, we carry forward this analysis to the general class of integrated correlators G (N |i,j) p (τ ). The results are highly uniform in p. We will compute the full matrix of correlators for p ≤ 5 -recall that for p ≥ 4 there are multiple trace structures -as well as a certain infinite family that exists for all p. The results will again be very simple in the SL(2, Z) spectral decomposition, with all integrated correlators computed herein determined simply by polynomials in the spectral parameter s. These polynomials obey fascinating recursion relations relating different values of both N and p. They are determined by a finite number of orders in weak coupling perturbation theory. The structural rigidity of these results raises obvious questions about why, exactly, they are so, and begs for a fundamental derivation of these recursion relations. From the AdS 5 × S 5 point of view, the fact that these are known for all N and τ means that the quantum type IIB string theory scattering amplitudes, upon suitable integration over AdS 5 boundary points, are known. To what extent that recasting can be made substantive, by telling us something non-trivial about type IIB string theory, is an intriguing open question for the future.
In Section 2 we introduce the necessary tools and background to study the integrated correlators -half-BPS operators, their correlators both integrated and unintegrated, and the SL(2, Z) spectral decomposition -and review the p = 2 case.
In Section 3 we present our first results, namely, the integrated correlators for p ≤ 5. At both p = 4 and p = 5 there are two half-BPS operators, one single-trace and one double-trace, so each case furnishes a 2 × 2 matrix. The technique is straightforward: using localization, one develops the perturbative expansion in g 2 YM ; infers from this the Eisenstein overlap by matching to the perturbative expansion of the spectral decomposition; and then justifies that the cusp form overlap vanishes. The result, as advertised above, is that the integrated correlators are determined in terms of a single polynomial, f 3 (s), respectively, such that the latter are ultimately fixed by the former in the SU (2) theory (plus a trivial initial condition at N = 1). This phenomenon persists for all cases that we study.
As for the cusp forms, in Subsection 3.3 we provide strong analytical support for their vanishing overlap from first principles, by performing explicit instanton calculations on the localization side 2 In fact, we find that f (N |i,j) p (s) ∝ (2s − 1) 2 , which further reduces the degree of the undetermined polynomial. and checking that they match the spectral result without cusp forms. These calculations use the instanton partition function in the presence of sources for half-BPS operators of higher charge [46], generalizing the original Nekrasov partition function. All told, this means that G (N |i,j) p (τ ) is determined completely by the first N + p 2 − 2 orders in perturbation theory. We do not know the fundamental reason for this remarkable fact. The recursion relations may be immediately uplifted to Laplace difference equations, which we show in Subsection 3.4.
In Sections 4 and 5, we further analyze the p ≤ 5 results. Section 4 collects their ensemble averages. Section 5 analyzes the integrated correlators of single-particle operators 3 at large N , developing the 1/N expansion of the spectral overlaps and the ensemble averages. This may be done algorithmically to arbitrary genus. We observe a match between the large N ensemble averages and the strongly coupled planar results, verifying the general relation established in [8]. Using the same correspondence between ensemble averages and bulk quantities, we also make a prediction for the one-loop AdS 5 × S 5 supergravity result for the integrated correlator at generic p -see (5.15). We explicitly develop the 't Hooft and very strongly coupled large N limits from the genus expansion of the spectral overlaps, noting some interesting functional representations along the way.
In Section 6, we present our second slate of results. We study the integrated correlator for which both operators O p are taken to be composites of O 2 with p/2 constituents, for arbitrary p ∈ 2Z + : We call this the maximal-trace family of integrated correlators, G (N |max) p (τ ). These are a type of generalization of the extremal two-point functions studied in the context of N = 2 SCFTs in four dimensions (e.g. [47][48][49][50][51][52]). We again find an explicit recursive solution for all N and p, with the same features as described earlier, but with one key difference: the recursion relation operates at fixed N , only shifting p. It takes the form of a semi-infinite lattice chain equation where G (N |max) p (τ ) refers to an appropriately normalized G (N |max) p (τ ) defined in (6.5), c = 1 4 (N 2 − 1) is the central charge and p = 2n. This is evocative of the Toda chain equation obeyed by the extremal correlators. The ensemble average admits a closed-form expression in terms of harmonic numbers -see (6.7). The large p limit of this sequence is a large charge limit, the details of which will appear in future work. We also develop the 1/N expansion of these correlators, with more nice polynomials emerging, see e.g. (6.10) and (6.20).
In Section 7, we make an ansatz (7.1) for the general integrated correlator G (N |i,j) p (τ ) with arbitrary trace structures i, j.
In Section 8, we conclude with some open problems and future directions.
Several appendices round out the text with computational details, consistency checks and collected formulas. 3 These operators, dual to single-particle states in AdS5× S 5 , are linear combinations of O (i) p defined to be orthogonal to all multi-trace operators. See (2.5).

Setup
In this section, we introduce the half-BPS operators and their four-point functions; define the integrated correlators and how one computes them from localisation; recall the rudiments of the SL(2, Z) spectral decomposition, a central tool in our study of the integrated correlators; and finally, as a segue to the next sections, review the p = 2 integrated correlator studied in [7,8].

Half-BPS operators and their four-point functions
We will study four-point correlation functions of half-BPS superconformal primary operators in N = 4 SYM theory with gauge group SU (N ). These operators have protected scaling dimension ∆ = p and transform in the [0, p, 0] representation of the R-symmetry group SU (4). The simplest realisation of such operators is in terms of single-trace operators T p , which we define by where Φ I (x), with I = 1, . . . 6, are the six real scalars of the theory and we introduced auxiliary SO(6) null-vectors y I (obeying y 2 ≡ y I y I = 0) to project onto the symmetric traceless representation. In addition to these, we will also consider multi-trace operators, which can be obtained from products of the form with p 1 + · · · + p n = p. Due to this degeneracy in the space of half-BPS operators, we introduce the notation O (i) p to distinguish the different single-and multi-trace operators of a given dimension p, assuming an ordering with increasing number of traces. 4 The species of operators appearing at weight p are in one-to-one correspondence with restricted integer partitions of p that do not include one. In what follows we will sometimes study operators with p ≤ 5, where we encounter at most double-trace operators: Since the operators O (i) p are half-BPS, their two-point functions are protected (meaning they are independent of the coupling τ ) and hence can be computed within the free theory. They are of the general form where g ij = y 2 ij /x 2 ij is the superpropagator and the colour-factors R p (N ) encode the non-trivial N -dependence. For the operators given in (2.3), we record the corresponding colourfactors in Appendix A. Note that in the trace-basis the two-point functions are not diagonal.
In the context of the AdS/CFT correspondence it is also useful to consider so-called singleparticle operators (SPO's), which we denote by O p without the superscript. 5 These operators are 4 For example, since the single-trace operator always comes first, we simply have O (1) p = Tp. Within the same number of traces we adopt lexicographic ordering, such that e.g. T2,4 comes before T3,3, etc. 5 In the following, any quantity without a superscript (i) is evaluated in the SPO basis. dual to single-particle states in AdS 5 × S 5 supergravity, namely the AdS 5 supergraviton multiplet (p = 2) and higher Kaluza-Klein modes (p ≥ 3) on S 5 . They are given by particular linear combinations of the O (i) p and can be elegantly defined for all N by their orthogonality to all multitrace operators [34,53]: At leading order in large N , O p coincides with the single-trace operator T p but receives 1/N suppressed multi-trace admixtures. 6 The first few SPO's are given by While we do not have an explicit formula for the generic colour-factors R (2.4), the N -dependence of the two-point function in the SPO basis is simply given by 7 [54] For the first few SPO's given above, this yields  .3) at fixed p. For instance, 2 T 4 = T 2,2 for N = 2, 3 and 6 T 5 = 5 T 2,3 for N = 3, 4. At large N , R p scales as R p (N ) ∼ N p /p.
Next, let us introduce four-point correlation functions of the half-BPS operators defined above. In particular, we will be interested in the class of correlators given by p . Superconformal symmetry constrains these correlators to take the form [55,56] The factor H (i,j) p is the only part of the correlator which depends on τ , and hence contains all of the non-trivial dynamics. It is multiplied by the factor I which is fixed by superconformal Ward identities to take the factorised form (2.11) 6 The coefficient of an n-trace admixture is suppressed by 1/N n−1 with respect to the single-trace contribution Tp. 7 Due to our choice of normalisation for the external operators (2.1), the colour-factor Rp used here has a factor of p 2 difference with respect to [54], i.e. Rp| here = Rp| there /p 2 . 8 Our conventions for the conformal and SU (4) R-symmetry cross-ratios read We denote the second R-symmetry cross ratio as µ instead of the canonical τ , to avoid confusion with the coupling.
On the other hand, the first term in (2.10) is the correlator in the free theory which can be computed by Wick contractions. It takes a particularly nice form in the SPO basis, which in our conventions reads where c is the central charge c = (N 2 − 1)/4 of the SU (N ) theory. While the structure of G (i,j) p,free for generic (i, j) takes a similar form in terms of the cross-ratios (u, v) and (σ, µ), the N dependence of the different (σ, µ) monomials is in general very non-trivial and the simplicity of the above result is a consequence of the orthogonality property of the SPO's O p [54].

Integrated correlators
Our main focus henceforth will be on a family of observables defined by a Euclidean integration of the dynamical part of It was first shown in [6] that G (N |i,j) p is given by derivatives of S 4 partition function of N = 2 * gauge theory deformed by couplings τ A of higher-weight chiral operators (here and below τ A collectively denotes all irrelevant couplings τ 3 , τ 4 , ... which couple to chiral primary operators of dimension ∆ > 2, whereas τ denotes the exactly marginal coupling). 9 In our conventions, this relation reads (2.14) The S 4 partition function 10 Z N (τ, τ A , m) takes the following form [46,50,57,58] .
Let us unpack these expressions. The integration variables a i are over the Cartan subalgebra of the gauge group SU (N ) -i.e. a i ∈ R and are constrained by the condition i a i = 0 -and a ij := a i − a j . m is the N = 2-preserving mass of the adjoint hypermultiplet. The function H(z) is given by product of two Barnes G-functions (2.16) The instanton factor Z inst (τ, τ A , m, a) in (2.15) is the localized contribution to the partition function from non-trivial instanton configurations of SYM gauge fields. It admits the decomposition where Z (k) inst (τ A , m, a) denotes the contribution from the k th instanton sector. The expansion of Z inst 2 appearing in (2.15), in powers of the instanton counting factor q := e 2πiτ , is The explicit form of Z (k) inst (τ A , m, a) is known for generic k. For our purposes we will only need the expression for k = 1 [46]: We note that this expression from [46] is for U (N ) gauge group. 11 However, it is applicable for computing correlators in N = 4 SYM with SU (N ) gauge group by simply imposing the constraint i a i = 0 on the eigenvalues in the U (N ) matrix model. 12 At the N = 4 conformal point, where m = τ A = 0, Z inst (τ, 0, 0, a) = H(0) = 1, and Z N (τ, 0, 0) is just a Gaussian matrix model in the (special) unitary ensemble that can be explicitly evaluated: An important detail in (2.14) is the operator mixing problem on S 4 . As explained in [50], in curved space, the source for a chiral operator O (i) p can have non-minimal couplings to lowerdimensional chiral operators O (i) p−2n with n ∈ Z + , due to non-vanishing background fields such as the spatial curvature. Therefore, to compute correlators on R 4 from the deformed S 4 partition function, one must solve this operator mixing problem. This is achieved through a Gram-Schmidt orthogonalization process which allows us to compute the complex vectors v i,µ p from the matrix of S 4 two-point functions. These vectors encode the mixing of weight-p chiral operators of species i with lower-weight chiral operators. Insertion of a chiral operator O (i) p in a correlation function on R 4 is then obtained through the combination of derivatives v i,µ p ∂ τ µ , where we sum over the index µ which indexes operators in the following set: p−4 } ∪ · · · (2.21) 11 Our notation is related to those of eq. (2.24) in [46] as follows: 12 This is unlike generic N = 2 theories. We are grateful to Silviu Pufu for a discussion.
where i is a fixed index, whereas the j n indices run over all species of operators of indicated weight, in accord with the discussion in Subsection 2.1. In a slight (but hopefully intuitive) abuse of notation, we write µ = p (j) to refer to the indicated operator O is to differentiate with respect to that source. The derivatives are defined as follows: for a given multi-trace operator O (j) p ≡ T p 1 ,...,pn such that p 1 + · · · + p n = p, with the understanding that τ 2 ≡ τ . Its action on Z N defines differentiation with respect to the source for the composite operators.

The SL(2, Z) spectral decomposition
Many observables of N = 4 SYM are SL(2, Z)-invariant. Such observables, call them O(τ ), admit a spectral decomposition into a complete SL(2, Z)-invariant eigenbasis of the Laplacian on the upper half plane. Let us introduce the main elements of this decomposition. More complete explanations may be found in Sections 2-3 of [8] in the N = 4 SYM context, and in [61,62] in great mathematical detail; see also [63] for a pertinent review. The fundamental domain for SL(2, Z) may be defined as This space is endowed with the hyperbolic metric which defines a Laplacian We define inner products of square-integrable SL(2, Z)-invariant functions via the Petersson inner product, The observable O(τ ), which is square-integrable [8], obeys the invariance equation Restricting to real-valued O(τ ), we develop a Fourier decomposition with respect to x, . (2.28) In N = 4 SYM, where x = θ/2π, the mode number k is the total instanton number. Instantonanti-instanton pairs contribute integer powers of qq = e −4πy , where q := e 2πiτ , with total instanton number zero.
The Roelcke-Selberg spectral decomposition takes the following form: This defines a convergent decomposition for any τ ∈ F. There are three branches here.
The first is the constant function, The bracket notation indicates that this is equivalent to the normalized ensemble average of O(τ ), that is, the average value of O(τ ) over the N = 4 supersymmetry-preserving conformal manifold with respect to the Zamolodchikov measure. This identification between modular and ensemble averages is special to N = 4 SYM, in which the Zamolodchikov metric is exactly hyperbolic due to maximal supersymmetry [64], unlike other sub-maximal theories whose conformal manifolds admit an SL(2, Z) action (e.g. N = 2 SQCD in four dimensions [47,65]).
The second is the continuous branch of non-holomorphic Eisenstein series, The star denotes that we use the "completed" Eisenstein series, whose Fourier expansion with respect to x is where Λ(s) := π −s Γ(s)ζ(2s) is the completed Riemann zeta function and σ 2s−1 (k) is the divisor function. This normalization is convenient because it manifests the reflection symmetry The overlap with O(τ ) is then a rescaled Petersson inner product, which may be simplified to   which follows from E 0 (τ ) = 1.
The third is an infinite discrete branch of Maass cusp forms in the so-called L 2 norm, (φ n , φ n ) = 1. These functions, infinite in number, obey the eigenvalue equation Reality of O(τ ) restricts the basis to the even cusp forms, invariant under x → −x. The cusp forms φ n (τ ) have no zero mode, instead decaying exponentially at the cusp y → ∞, hence contributing to neither the perturbative nor 't Hooft regimes. The cusp forms are rich mathematical objects exhibiting various signs of arithmetic chaos [43][44][45]; however, quite intriguingly, the integrated correlators G (N ) p (τ ) will be seen to have vanishing cusp form overlap, so we will not say anything further about the cusp forms here.
In summary, the entire content of O(τ ) is stored in the spectral overlaps {O, E s } and (O, φ n ).
Looking ahead slightly, our computations will begin from the zero mode, O 0 (y): As shown in [8], N = 4 SYM observables with a consistent perturbative expansion in g 2 YM admit a rather constrained analytic form for their Eisenstein overlaps: The function f p (s) encodes the complete perturbative part of O(τ ), its residues giving the weak coupling (y 1) data. The function f np (s) encodes non-perturbative, instanton-anti-instanton corrections. Both functions must be invariant under s → 1 − s, and are regular for s ∈ C away from s = 0 and its reflection. For the observables we consider, it will turn out that f np (s) = 0, a fact whose origin lies in the Borel summability of their perturbative expansions; this further implies [8] that f p (s) must be regular at s = 0 and is thus entire.

Prelude: Review of p = 2 result
To illustrate the spectral decomposition (2.29) at work -and to set the stage for our computations ahead -let us consider the example of the p = 2 integrated correlator, G That is, the ensemble average, and Eisenstein and cusp form overlaps, are For N > 2, it turns out that the spectral overlaps obey a powerful three-term recursion relation [7] This fully determines G (N ) 2 (τ ) for any SU (N ) in terms of the SU (2) result alone (and its triviality at N = 1).
To reiterate, the complete content of G This observable also admits a (formally equivalent) lattice-integral representation, as originally seen in [7], which may be derived by performing an SL(2, Z) Poincaré sum of the zero mode of (a Borel resummation of) G (N ) 2 (τ ). The same will be true of the p > 2 integrated correlators. 13 As discussed in the Introduction and hopefully made clear by our presentation, instead writing G • The non-perturbative part of the Eisenstein overlap vanishes for all N .
• The perturbative part of the Eisenstein overlap, f The passage between bases also helps to clarify the physical meaning of some of the mathematical aspects of the lattice-integral representation itself, such as the identification of the lattice-integral kernel as an SL(2, Z) Borel transform and the observed properties of G Let us make a few further comments. First, the SU (2) result is, rather literally, the simplest result consistent with the analyticity properties of the overlap: in particular, f (N ) p (s) must be an entire, reflection-symmetric (i.e. even under s → 1 − s), non-constant function of s, the simplest example of which is the degree-two monomial shown above. Second, since the cusp form overlap vanishes, the above recursion relation may be uplifted to a 'Laplace difference equation' for the integrated correlator [7], which reads where we have used (2.31). This differential relation is equivalent to the algebraic relation (2.42). Finally, the large N expansion of G (N ) 2 (τ ) was analyzed in [7,8], in various limits and including non-perturbative corrections in both 1/λ [7] and 1/N [8], where λ := g 2 YM N is the 't Hooft coupling. We refer the reader to those papers for details. Here we just highlight one result, to be generalized to p > 2 later in this work: at λ 1, 13 It is straightforward to prove [8]  The leading term, which gives the value of the integrated correlator in tree-level AdS 5 × S 5 supergravity in our normalization, manifestly equals the ensemble average in (2.40) at large N [8].
With these observations in place, we now proceed to study the more general class of integrated correlators G (N |i,j) p (τ ). The SL(2, Z) spectral decomposition will greatly simplify the physical and technical analysis.

Results I: Integrated correlators for p ≤ 5
We now derive the integrated correlators G (N |i,j) p (τ ) for p ≤ 5 in the SL(2, Z) spectral decomposition. As we will see, their general structure is uniform in p, essentially identical to that of the p = 2 case presented in (2.40)-(2.42). The derivation is a three-step process:

1)
In Subsection 3.1, we reconstruct the full zero(-instanton) mode of the integrated correlator from its weak-coupling expansion. The latter is computed from localization as outlined in Section 2.2 and Appendix B.

2)
In Subsection 3.2, we leverage the zero mode to deduce the Eisenstein overlap (2.38) by matching to the spectral decomposition. Just as for p = 2, these will have vanishing nonperturbative part, f np (s) = 0, and the remaining piece will be fixed by polynomials obeying powerful recursion relations in N .

3)
Finally, in Subsection 3.3, we upgrade this to the result for the full integrated correlator by arguing that the cusp form overlap vanishes. We do so with several explicit computations. This then implies that the aforementioned recursion relations are relations for the full integrated correlator. We note that they may be recast as Laplace difference equations, an exercise we carry out in Subsection 3.4.

Perturbative expansion
The weak-coupling expansion of G (N |i,j) p (around a zero-instanton background), may be extracted from the localization expression (2.14)-(2.15) by expanding the latter at y 1 (recall that y := Im τ ). Details of this computation may be found in Appendix B for all p ≤ 5. The final result may be written as are given in (B.14) and (B.8) respectively. Let us copy the result from Appendix B for p = 2, 3 below: The results for p = 4, 5 may similarly be assembled from Appendix B.
Before extracting the spectral overlaps from these, we note that these results are in agreement with previous results from two-loop perturbation theory: in particular, following the strategy of [7], in Appendix E we take the two-loop results for the unintegrated correlators from [9], integrate them against the measure (2.13), and compare to the localization result shown above. The results match, as seen in Appendix E. This provides a non-trivial check of the method used here. It also probes a novel aspect of the p = 4, 5 cases, in which there is a matrix of integrated correlators involving double-trace composite operators; the agreement described above verifies that the localization method works even when the external operators are composites. 14

Spectral overlaps and recursion relations
To deduce the Eisenstein overlap from the perturbative expansion, one must develop the latter from the spectral decomposition. This is straightforwardly done by closing the integration contour in (2.37) to the left, yielding 15 where we recall the definition of the completed Riemann zeta function Λ(n + 1 2 ) below (2.32). This is then compared to the weak-coupling expansion of G Using the perturbative expansions derived in the previous subsection and Appendix B, we find that in all cases we considered, where g (N |i,j) p (s) are reflection-symmetric polynomials of degree 2N + 2 p 2 − 6, for all labels (i, j). The perturbative expansions are Borel summable. This typically (but not always [66,67]) implies, in the context of resurgence, that non-perturbative corrections in g 2 YM are not present. This was confirmed in [7] for p = 2, implying that the non-perturbative part of the Eisenstein overlap (2.38) vanishes [8]. One can confirm this property for the p > 2 cases considered here as well: 16 f 14 This is also a further verification that the Konishi operator does not contribute to G (N |i,j) p (τ ), as was discussed in [6] in the 't Hooft limit and [7] for finite N . Therefore, (3.5) gives the full Eisenstein overlap. Summarizing so far, the Eisenstein overlaps of G (N |i,j) p for p = 3, 4, 5 take the form , E s } may be completely fixed by the first N + p 2 − 2 orders in weakcoupling perturbation theory.
Given this simplicity, and inspired by the p = 2 case (2.42), we are encouraged to seek recursive relations for p > 2. To do so, we use the above observations on the general structure of the spectral overlaps to construct the most restrictive polynomial ansatz, allowing us to compute the overlaps f (N |i,j) p (s) for p ≤ 5 and N ≤ 20 from the perturbative expansions to order y −20 . As we will discuss in detail in the following, we are then able to find many interesting relations generalising the p = 2 recursion to relations not only among theories with different values of N , but also among integrated correlators of different values of p.

p = 3
In the p = 3 case, the perturbative expansion vanishes for N = 1, 2 (due to vanishing of the operator O 3 itself) and hence we have f 3 (s) = 0. For concreteness, we display the expressions for the next few values of N which read 3 (s) = 1 120 (2s − 1) 2 (s 6 − 3s 5 + 95s 4 − 185s 3 + 1824s 2 − 1732s + 6240). (3.8) Considering many more cases up to N = 20, we find that the above overlaps satisfy the following striking recursion relation: This is a remarkable formula: it allows one to recursively compute the p = 3 overlaps for generic N from knowledge of the p = 2 overlaps f (N ) 2 (s), which in turn are fully determined by the SU (2) theory alone, together with the trivial initial condition f (1) 3 (s) = 0! Let us note that this recursion is consistent starting from the value N = 1. In particular, the fact that f 2 (s) = 0 in the first term on the RHS thanks to the presence of an (N − 1) factor. An interesting feature of the above relation is that it has no explicit dependence on s, in contrast with the p = 2 recursion relation (2.42). This is concordant with the fact that the spectral overlaps f (3.10) Note that due to the factors of (N − 2)(N − 3) on the RHS, the recursion starts at N = 4, and as a consequence of being a four-term recursion relation one needs to supply f (N ) 3 (s) with N = 2, 3, 4 as initial data, making it a somewhat less powerful statement than the previous relation (3.9).
Interestingly, while for N = 1 the perturbative expansions vanish, we find that for N = 2, 3 the three expansions degenerate and are simply proportional to each other termwise in their 1/y expansions. In terms of spectral overlaps, this yields This phenomenon is a non-trivial consequence of considering the SU (N ) theory and is explained by the fact that for low values of N the single-trace operator T 4 is no longer linearly independent of the double-trace operator T 2,2 : for N = 2, 3 one has T 4 ∝ T 2,2 and therefore also their correlators become proportional to each other. In view of the different intricate N -dependence of the weakcoupling expansions (B.14c)-(B.14e), the degenerations for N = 2, 3 seem highly non-trivial (and further support our claim that the localisation integrals correctly compute multi-trace insertions).
For N ≥ 4, the spectral overlaps f (N |i,j) 4 (s) start to differ for different (i, j). For example, the case (i, j) = (2, 2), corresponding to the correlator with two double-trace operators, is the one with simplest s-dependence: (3.12) We have obtained analogous expressions for the (i, j) = (1, 1) and (1, 2) cases. However, we refrain from presenting more explicit data since, by considering many cases of different N , we discover that the p = 4 spectral overlaps likewise obey interesting recursion relations, which we present now.
(i, j) = (2, 2): Let us start with the simplest case given by the correlator involving two doubletrace operators. In this case we find the very direct relation 2 (s) only. This will be elaborated upon in Section 6.
(i, j) = (1, 2): For the mixed correlator with one single-trace and one double-trace operator, This, too, determines f Recursion relations also exist for the overlaps of the correlator involving two singletrace operators albeit in somewhat more complicated form. For this case, we find The new feature relative to the previous relations is the presence of f . Note that this particular linear combination is precisely such that f  17 Here and in what follows, we leave implicit the trivial vanishing conditions of the spectral overlaps at N = 0, 1. 18 We have also found a recursion relation which does not involve the p = 3 overlaps, but its structure as well as the N -dependence of the coefficients is more complicated and we refrain from giving it here. Notably, the RHS of that relation has a pole at N = 3 and hence the recursion can be used only for N ≥ 4. 19 Because the overlaps f (N |i,j) 4 (s) become proportional to each other for N = 2, 3 (cf. (3.11)), the mixing coefficient β must take the same value for N = 2 and N = 3 in order to be consistent with f  p (s), whose precise form we record in Appendix C.
As for p = 4, we observe that the perturbative expansions degenerate for low values of N : as a consequence of the two half-BPS operators becoming proportional to each other, i.e. T 5 ∝ T 2,3 for N ≤ 4, we have that the weak-coupling expansions vanish for N = 1, 2 (due to O 3 which ceases to exist for N = 1, 2), and that for N = 3, 4 the expansions degenerate. As mentioned for p = 4, the latter property is not at all obvious from the different N -dependent coefficients in the perturbative expansions, see equations (B.14f)-(B.14h). For the spectral overlaps we therefore have and only from N ≥ 5 onwards their s-dependence differs. Comparing again cases with varying N , we are able to find recursion relations for the p = 5 spectral overlaps, which we list below: where the above recursion relations start from N = 1 with initial values f (N |i,j) p (s) = 0 for N ≤ 1. As noted earlier in the p = 3 case, the factors of (N − 1) in front of the f (N +1) 2 (s) terms on the RHS are required by consistency of these formulae for N = 1.
(SPO, SPO): Lastly, the p = 5 spectral overlaps of the correlator of SPO's are given by an analogous equation as given in (3.21), where the mixing coefficient now reads . This is the right linear combination which yields f

Including instantons: completion to the full correlator
Having computed the zero mode of G (N ) p (τ ), and therefore the Eisenstein spectral overlap, what remains is to determine the cusp form overlap. We claim that it vanishes: That is, the full integrated correlator is where recall that the average is G . Note how strongly constraining this is: because g (N |i,j) p (s) is an even polynomial of degree given in (3.7), the integrated correlator may be completely fixed by the first N + p 2 − 2 orders in weak-coupling perturbation theory.
We will provide strong support for (3.22) with explicit instanton computations below. The basic idea behind the check is as follows. First, one can compute the perturbative expansion of the integrated correlators in nonzero instanton sectors as predicted by (3.23). This result is then to be compared with a direct computation from localisation formula (2.14) using the relevant instanton partition function. For p = 2, the absence of cusp forms was extensively checked in this manner in [7,8]. In this subsection, we carry out the check for several p > 2 integrated correlators using the instanton partition function (2.19). 21 Happily, the two results agree.

Calculation
The first step of extracting the perturbative expansion around k instantons from (3.23) is straightforward. (The anti-instanton result, proportional to e −2πikτ , is identical.) One just inserts the k'th Fourier mode of E * s (τ ) from (2.32) and the y 1 expansion Borel summability allows us to swap the s integral in (3.23) with the sum over n [8]. Closing the contour towards the left, and taking k = 1, we find where we have used the symmetry property a n (−m) = a n (1 + m). For fixed n, the sum in parentheses is divergent, since both f (N |i,j) p (s) and a n (s) are polynomials in s. The sum can be regulated in a standard manner using an exponential or zeta function regulator.
Let us first apply this strategy to the p = 3 integrated correlator, G (N ) 3 (τ ). Plugging its spectral overlap into (3.25) gives the prediction where h 3,n (N ) are polynomials in N of degree n + 2. These may be obtained algorithmically for any n via recursion. For first few values of n, We now compare (3.26) to the exact localization calculation of the same quantity via the definition (2.14). To proceed, we use the expression (2.19) for the one-instanton partition function and set the higher-weight couplings τ p>3 = 0, for which we have Plugging this expression into (2.15), we develop a weak-coupling (large y) expansion for the p = 3 integrated correlator in the one-instanton sector (exactly in the same manner as we did for the zero-instanton computations) using the definition (2.14). A finite generic N computation is technically involved due to the form of the instanton partition function. But for low values of N , the computation is straightforward. Performing this calculation for N = 3, 4, 5 through the first several orders in 1/y, we find precise agreement with (3.26) as obtained from the spectral decomposition. 22 The above strategy applies for any integrated correlators G (N |i,j) p (τ ). Analogous results for the full matrix of p = 4, 5 correlators are deferred to Appendix D for which we have also checked in this way that the cusp form overlaps vanish at low values of N . 23 There is yet another class of integrated correlators -to be introduced more thoroughly in Section 6 -which is especially simple. This is what we call the maximal-trace family of integrated correlators. For p ∈ 2Z + , these are defined by taking both operators O (s) was computed in (3.13). Inserting this into (3.25), we find where h (3.31) We now compare (3.30) to the exact localization calculation of the same quantity using the instanton partition function with vanishing higher-weight sources (τ p≥3 = 0), The authors of [68] have, in work to appear, independently studied the p = 3 integrated correlator at low values of N and performed numerical checks of the vanishing cusp form overlap at finite y, using the convergent expansion of the deformed instanton partition function. This further complements our conclusion. We thank Shai Chester for discussions. 23 It would be interesting to reproduce and independently confirm these perturbative results by explicitly integrating the one-instanton correction to the unintegrated O2O2O correlators against the measure (2.13). The latter may be computed in principle using methods based on the ADHM construction of instantons [69][70][71][72][73], though very few explicit results for instanton corrections to the unintegrated correlator exist in the literature [69,74,75].
It is convenient to first compute ∂ 2 m Z (1) inst (0, a) in a weak-coupling expansion (the expectation value is in the ensemble (2.20)). This result was obtained for finite N in [7,76], with the first few orders being Plugging the above into (2.14) (using (B.7) for solving the operator mixing problem) we precisely find (3.30). The same calculation may be repeated for higher-p maximal-trace integrated correlators, studied further in Section 6; having done so for many p > 4 again yields perfect agreement.
Altogether, these results provide substantial evidence of vanishing cusp form overlap for integrated correlators G (N |i,j) p (τ ). We leave a deeper physical understanding of this intriguing mathematical property for the future.

Laplace difference equations
Having given evidence for vanishing cusp form overlap and thereby completion to the full correlator, the various recursion relations from Section 3.2 should be regarded as non-trivial identities between the full integrated correlators. Said another way, relations for f Prime among these relations is the p = 3 recursion relation (3.9), which uplifts to the following powerful difference equation obeyed by the integrated correlator G (N ) 3 (τ ): The absence of an inhomogeneous term coming from the ensemble averages G (N |i,j) p follows from consistency of the recursion relations evaluated at s = 1.
One may, if desired, recast all of the previously derived recursion formulas as Laplace difference equations. Let us just give one more example explicitly, namely, the very simple relation (3.13) between G (N |2,2) 4 (τ ) and G We leave other applications of (3.34) as an exercise for the enthusiastic reader.

Ensemble averages
We now consider G (N |i,j) p , the ensemble averages of the integrated correlators with respect to the Zamolodchikov measure, cf. (2.30). As can be seen from equation (2.35), these are obtained by simply evaluating the Eisenstein spectral overlaps f (N |i,j) p (s) at s = 1.
For p = 2, this yields the result already given in [8], which for convenience we reproduce here: We find that for higher p the ensemble averages continue to be given by rational functions of N . Explicitly, we find (in slight abuse of matrix notation) We observe that at leading order in large N the diagonal entries scale as N p , while the offdiagonal ones scale as N p−1 : Furthermore, note that the ensemble averages of correlators with two single-trace operators (i.e. the (i, j) = (1, 1) components) at large N exhibit p-dependence given by For the correlators with one single-trace and one double-trace operator, the available data for the ensemble averages at large N is consistent with the formula where the label j denotes a specific double-trace operator T p 1 ,p 2 with p 1 + p 2 = p.

Ensemble averages for SPO's
To make contact with correlators describing scattering amplitudes in AdS 5 ×S 5 , let us also consider the ensemble averages of integrated correlators involving SPO's. For p = 2, 3 these are equal to the averages for single-trace operators already given in equations (4.1) and (4.2) above. On the other hand, for p = 4, 5 we find (4.8) Due to the gamma factors, these have the expected property that they vanish for N = 1, 2, . . . , p−1.
While at first sight these expressions look rather messy, a surprising simplification occurs upon pulling out a factor of the SPO two-point function normalisation R p , given in (2.7), and performing a decomposition into partial fractions: Recalling that at large N we have R p ∼ N p /p, one recognises that the first term in the brackets is again consistent with the leading contribution (4.6), which is then followed by a string of 1/N suppressed terms. Moreover, a clear pattern is visible in the above expressions, namely with the variable x p taking the values x p = 0, 0, 1, 5 for p = 2, 3, 4, 5, respectively. Let us emphasise that the above pattern is based on data up to p = 5 only. It would be interesting to consider higher values of p to either confirm that this nice structure persists or to find a more general structure of 1/N corrections.
We will further analyze the large N structure of G (N |i,j) p in Section 5.2.

Large N
In this section we study the large N limit of the integrated correlators G (N ) p (τ ). In order to make contact with the dual AdS 5 × S 5 supergravity, we restrict ourselves in this Section to correlators of single-particle operators (SPOs) O p . As such, we drop the (i, j) labels.
All of the N -dependence of a given observable O(τ ) lies purely in its spectral overlaps, so the large N expansion of the former follows from that of the latter. The large N expansion of G Anticipating the genus expansion of the 't Hooft limit in (5.16), we call f (g) p (s) the genus-g spectral overlap.
Developing the genus expansion of f p (s). The latter may be computed in various ways: for example, by studying the leading-order N -dependence order-by-order in the 1/y weak-coupling expansion and reconstructing the s-dependence; or, by using the previously known results for p = 2 to seed the recursion.
In Subsection 5.1 we present results for f (g) p (s) through g = 2, though it is straightforward to compute to arbitrary genus using the recursion formulas. Taking their s → 1 limit, it is then trivial to compute ensemble averages at large N (cf. (2.35)) and make contact with AdS 5 × S 5 supergravity. We do this in Subsection 5.2. Finally, plugging in the large N expansion (5.1) of the spectral overlaps into the spectral decomposition (2.29) allows us to make manifest two different large N limits of G

Genus expansion of spectral overlaps
At genus zero, we find the following result for p = 2 which reproduces the result of [7,8]. The arbitrary p generalization of this result is given by For a derivation of the equation we refer the reader to Appendix G. The presence of simple poles at s = 1 2 − m for m ∈ Z + (in addition to zeros at s ∈ Z − ) is consistent with the allowed general properties of spectral overlaps in the genus expansion. 25 Let us make a curious observation. If we write (5.3) in the form 24 The canonical development of the genus expansion of a general observable O(τ ) starts at order N 2 . While a given observable may have vanishing contributions at low genus -for example, a conformal dimension which starts at order N 0 -it is nevertheless convenient to refer to its leading contribution as genus zero, instead adjusting the leading power of N accordingly. Our conventions for G (N ) p (τ ) are, for various reasons of clarity, unnormalized, which gives them a leading power N p . Upon normalizing by the two-point functions of O2 and Op using the colour factors R2 and Rp given in (2.7), their leading power would be 1/N 2 , which is the correct scaling of a normalized connected four-point function in a CFT obeying large N factorization. 25 We recall from [8] that whereas the finite N overlap f we find that the numerator polynomials n (0) p (s), of degree 2 p−2 2 , obey a reflection symmetry for all p ≥ 2: that is, This symmetry is not required by SL(2, Z)-invariance of the spectral decomposition.

(5.8)
Let us point out the occurrence of certain zeroes in these polynomials, which have non-trivial consequences for the large N expansion of the ensemble averages discussed in the next section. In particular, note that for p = 2 the spectral overlap f (g) p (s) vanishes at s = 0, 1 for g ≥ 1, while for p = 3 the same is the case for g ≥ 2. On the other hand, the p = 4, 5 overlaps do not vanish at s = 0, 1, and the reason for this will be explained next.

Ensemble averages at large N and supergravity
We now make contact with semiclassical AdS 5 × S 5 supergravity.
It was shown in [8] that for any SL(2, Z)-invariant observable O(τ ) that admits a genus expansion in the 't Hooft limit, the ensemble average O and the large 't Hooft coupling limit of O(τ ) are equal at leading order in large N : This correspondence extends to all genera: defining the leading-order term at genus-g as there is an equivalence Translated into bulk terms, the genus-zero relation implies an equivalence between the large N ensemble average of O(τ ) and its value in tree-level AdS 5 × S 5 supergravity, while the g > 0 relation relates the higher-genera averages to loop-level AdS 5 × S 5 supergravity after string theory regularization of divergences. 26 (The g > 0 statement of (5.11) will be elaborated upon below.) At g = 0, our results confirm this equivalence. As foreshadowed in (4.6), we find After accounting for differences in normalisations, this agrees with the λ → ∞ limit at g = 0 [6].
At g = 1, the results (4.9) for p ≤ 5 are consistent with the following pattern: (5.13) Via (5.11), this makes a prediction for the finite term in the λ → ∞ limit at genus one. This amounts to a holographic computation of the 1-loop supergravity computation of G (N ) p (τ ). To put this in conventional supergravity language, we note that the bulk loop expansion proceeds not in powers of 1/N but in powers of 1/c = 4/(N 2 − 1). In addition, let us divide by the product of colour factors R 2 R p to normalize the correlator, 27 Putting things together leads to This sets a benchmark for a direct bulk computation of the integrated correlators in regularized 1-loop AdS 5 × S 5 supergravity. 28 We reiterate that this is the finite result that remains after the 26 We note that one can perform a mutual consistency check of this dictionary and our results: namely, matching the large N expansion of the finite N ensemble averages G (N ) p to the genus expansion of the spectral overlaps fp(s). We do so in Appendix F. 27 In this convention, the tree-level supergravity result is G (unambiguous) string theory regularization of 1-loop supergravity divergences. It is notable that the ensemble average is sensitive to the UV details of this regularization; in particular, string theory regularizes divergences unambiguously, choosing a specific renormalization scheme, and the SL(2, Z) average identifies precisely this scheme.
Finally, let us make the side remark that (4.7) gives a prediction for tree-level supergravity, this time for the integrated correlator involving one single-trace operator T p and one double-trace operator T p 1 ,p 2 with p 1 + p 2 = p.

Integrated correlators at strong coupling
As emphasized throughout this section, the spectral overlaps in the 1/N expansion contain complete information about the large N integrated correlators themselves: simply plug the formulas for the overlaps f (g) p (s) into the spectral representation (2.29), and expand in the desired limit. Having said that, for convenience of future study, we assemble these ingredients into explicit expressions for G (N ) p (τ ) in two strongly coupled limits of interest: the 't Hooft limit, and the very strongly coupled (VSC) limit. We will briefly review salient aspects of these limits in the spectral language; a systematic treatment of these limits for general SL(2, Z)-invariant observables can be found in [8].

't Hooft limit
In the 't Hooft limit, where g YM → 0 and N → ∞ with λ := g 2 YM N fixed, the 1/N expansion organises into a genus expansion, up to non-perturbative corrections in 1/N . Because instantons and anti-instantons are non-perturbatively suppressed in 1/N , the genus expansion is an expansion of the zero-instanton mode: p,0 (λ) + (non-perturbative) .

(5.16)
We henceforth drop the "0" subscript with the understanding that we are working perturbatively in 1/N . The parameter g differs from the previous definition of genus in (5.1) because of a nontrivial repackaging of the 1/N expansion by the SL(2, Z) spectral decomposition, a fact we return to momentarily.
We want to perform a spectral decomposition of this quantity. The zero mode of completed Eisenstein series may be written in terms of a 't Hooft coupling as As the constant term in the spectral decomposition, the ensemble average G (N ) p does not depend on λ. There are two terms in brackets, each of a different nature. Large N forces us to close the contour of the second term in brackets, carrying an N 2s−1 , to the left; this generates terms with positive powers of λ, regardless of the value of λ. These terms were denoted 'renormalization terms' in [8], due to their strong coupling interpretation as bulk UV divergences regularized by the string scale cutoff. 29 The integration contour of the first term in brackets may be closed either to the left to develop the strong coupling expansion (in negative powers ofλ); to the right to develop the weak-coupling expansion (in positive powers ofλ); or not deformed at all, thus giving an expression which is exact in λ.
For example, at g = 0, 1 one finds p (s) were given in (5.3) and (5.6), respectively. The first term in (5.21) is the 'renormalization term' discussed above, and descends from the second term in brackets in (5.19) for g = 0. At risk of repetition, we stress anew that this spectral representation is exact in λ: no expansion has been made. G (N ) p (τ ) contains non-perturbative effects at both large λ and large N . This follows from the analysis in [8], where they were determined for p = 2, and the fact that the p-dependence of f respectively, where λ S := 16π 2 N 2 /λ is an S-dual 't Hooft coupling. Both corrections are present in the 't Hooft limit. The p-dependence of the non-perturbative series, entering only in the expansion coefficients, would be nice to examine in future work.
It is formally interesting to derive two alternative presentations of these same objects G (g) p (λ). First, it is possible to show (see Appendix G.3) that the genus-g correlator can be written entirely in terms of genus-g spectral overlaps as follows: This rewriting involves shifting the contour of (5.19) and making use of identities relating different genera -namely, eq. (7.9) of [8] -that are required by consistency of the weak-coupling expansion in integer powers of λ. The novelty of this expression is that the integrand involves only the genus-g overlap, unlike what (5.19) seems to give. Note the absence of the second line for g < 2.
Secondly, in [6] the genus-zero result G (g=0) p (λ) was derived for all λ as an integral expression involving squared Bessel functions. In Appendix G, we prove its equivalence to (5.20). Indeed, a "Bessel representation" at all genera maybe derived by deforming the contour in (5.23) to the right, yielding This is convergent for |λ| ≤ π 2 for all finite g and p. The expansion may be resummed using the integral identity The resulting expression is a one-dimensional integral of products of Bessel functions. For example, we present such a result for G (g=1) 3 (λ): For the sake of completeness, we provide the results at genus one for p = 2, 4, 5 in Appendix G. Similar expressions for G (g) p (λ) for various g and p are easily obtained from (5.24).

Very strongly coupled limit
The VSC limit is the limit of N → ∞ with τ fixed. In contrast to the 't Hooft limit, instanton contributions are not suppressed. This limit is trivially extracted from the spectral decomposition (2.29) upon inserting the large N expansion of the overlaps (5.1) into (2.29), and deforming the integration contour to develop the 1/N expansion. Summing over residues in this way gives [8] G (N ) The τ -dependence appears only via Eisenstein series, and only down by odd-half-integer powers of 1/N . 30 Both of these features, like many others in this work, follow from the relations f For example, applying the above formulae to the p = 3 case using the genus expansion (5.6) of its spectral overlap through g = 2 yields the following VSC expansion:

Results II: Integrated maximal-trace correlators for all p
So far, we have mostly restricted ourselves to integrated correlators with p ≤ 5. Recall that our approach was to obtain the weak-coupling expansion of G (N |i,j) p (τ ) using the localisation relation (2.14), from which we then inferred the spectral overlaps. When trying to generalise the localisation computation to higher values of p for most choices of (i, j), one is met with some computational expense, to be elaborated upon in Section 7.
However, there exists a special class of integrated correlators, defined for arbitrary p, which we now determine exactly: the so-called maximal-trace correlators.
The key point is that, as shown in [50], one can recursively construct a basis of operators on the sphere such that certain families of operators become orthogonal to each other. The simplest such family of operators with this property is given by maximal-trace operators, introduced earlier in (3.29), which for even p are given by They mix only among themselves, leading to an effective reduction in the size of mixing matrices. 31 The solution to the sphere mixing is given in terms of the vectors v j p , recall the discussion around 30 There can be terms down by powers of 1/N 2 as well, due to the ensemble average terms G (g) p , but these are τ -independent. 31 The same reasoning applies also to the odd p family of maximal-trace operators , {T3, T2,3, T2,2,3, T2,2,2,3, . . .}, but in what follows we restrict ourselves to the even p case.

(2.21). For the first few cases of p we find
(6.2) The above vectors admit the closed form expression is built from products of T 2 only, one only needs to act with repeated ∂ τ and ∂τ derivatives in order to bring down insertions of these operators in the localisation computation. As such, no new matrix integrals have to be performed when considering correlators of these maximal-trace operators.
In what follows, we will denote the integrated correlators following the notation of (2.10) and (2.13). 32 The corresponding colour-factor appearing in the localisation relation (2.14) is recorded in (A.2). We also introduce a hatted notation to denote normalization by R (max) p (N ), e.g.
and likewise for other quantities.
Let us present, then elaborate upon, the result: The ensemble average is where H n denotes the harmonic number, and the overlaps obey the polynomial recursion We first address the spectral overlaps. These are deduced as usual from the weak-coupling expansion obtained from localisation. 33 The usual function f (N |max) p (s) parameterizing the Eisenstein overlap à la (2.38) is found to be of the now-familiar form where the polynomials g (N |max) p (s) are even and of indicated degree. Considering many cases of different N and p leads to (6.8). 34 The central property here is that there are no shifts in N : in contrast to the previously-studied cases, this is a recursion relation in p alone, with g (N ) 2 (s) as the initial condition (together with the trivial condition g (N |max) 0 (s) = 0). Indeed, the solution of the recursion relation for general p and finite N is simply proportional to the p = 2 solution: (6.10) The functions F p (N, s) have the structure where the complete s-dependence is captured by the h   s) . (6.13) 33 As mentioned earlier, one can verify that the weak-coupling expansion agrees with the explicit results from perturbation theory upon integration. For the maximal-trace correlators, we have performed this check to two-loop order and up to p = 6, see Appendix E. 34 As we have observed in previous cases, for low values of p the above formula neatly degenerates in a consistent way. This can be most easily seen when considering the recursion relation for the unnormalised overlaps g (s) = 0: p = 0 corresponds to an external identity operator, for which the connected correlator vanishes. Moreover, the other coefficients combine to correctly recover relation (3.13). Finally, for p = 6, the two terms in the second line of (6.8) degenerate and in particular their coefficients combine to give the correct recursion relation consistent with the data.
As is evident from (6.6), the cusp form overlap vanishes: (6.14) This assertion follows from the computations discussed in Section 3.3: for many values of p, we have directly checked to the first few perturbative orders in 1/y around one instanton that (6.6) agrees with the first-principles localisation computation.
Finally, we turn to the ensemble averages G . For general p we have found an explicit formula (6.7), not just a recursion relation. Since p ∈ 2Z + , an equivalent representation of (6.7) for the ensemble averages is given by 35 At leading order in large N (with finite p), For p = 2 and p = 4, we recover equations (4.6) and the lower-right entry of (4.5), respectively, after multiplying by R (max) p (N ) at large N . 36

Comments and connections to N = 2 extremal correlators
The maximal-trace family of integrated correlators, G (N |max) p (τ ), has been determined above for all N and arbitrary p. Compared to the general trace structures studied earlier in this work, its solution is especially clean. It also has interesting connections to various other quantities studied in the literature.
First, note that G (N |max) p (τ ) may be thought of as a generalization of the so-called "extremal correlators" studied in 4d N = 2 SCFTs. Let us summarize the latter. Extremal correlators involve several chiral operators O m and a single anti-chiral operator O n . The most well-studied extremal correlators are two-point functions, These depend non-holomorphically on the complexified gauge coupling τ , and may be determined by localisation. In an N = 2 SCFT of arbitrary rank, one may take O n to be n-fold composites of the weight-two chiral primary, what we called O 2 . Such an O n is precisely our family O (max) p defined in (6.1), with the identification 2n = p. The correlators G 2n (τ ) solve differential recursion relations in n [47,64]. In rank-one N = 2 SCFTs, such as SU (2) SQCD with four flavors, the 35 It may be worth noting that the unnormalised average, G , is actually a polynomial in N of total degree p. This follows from the explicit formula (A.2) for R At higher rank, differential relations are believed to exist, but to take a more involved form that couples different families of operators besides merely composites of O 2 [50].
We can now state the connections between G 2n (τ ) and G (N |max) p (τ ) with 2n = p. Both are correlators involving the composite operators O (max) p . While the former are two-point functions, the latter are four-point functions, but integrated over positions, leaving only τ -dependence. Strikingly, both obey a similar form of differential recursion: as a consequence of (6.14), the recursion (6.8) uplifts to a differential relation for G This was written in terms of the central charge and n = p/2 in (1.6). Despite obvious differences, this equation is intriguingly similar to the decoupled semi-infinite Toda chain equations obeyed by the extremal two-point functions G 2n (τ ) in the SU (2) theory [47][48][49][50]. 38 Note, though, that (6.19) holds for all N ! Evidently, the extra power of N = 4 supersymmetry makes the integrated fourpoint function even more soluble than its N = 2 extremal counterparts. Altogether, this suggests that one may think of the Laplace difference equations for integrated N = 4 SYM correlators in general, and (6.19) in particular, in the same spirit as the Toda equations (6.18), and that they may be derivable by some generalization of the tt * method that led to (6.18).
The other point we wish to make is the following: the large p limit of G (N |max) p (τ ) is a large charge limit [51,80]. Since G (N |max) p (τ ) is determined explicitly for all p, this regime is easy to access, given the simplicity of the recursion (6.8). For example, together with the spectral decomposition (6.6), it is quick to see that G (N |max) p 1 (τ ) limit admits a 't Hooft-like expansion with g 2 YM p fixed, for any finite N . This generalizes a main result of [52], which established the 't Hooft-like limit of the extremal correlators (6.17) in rank-one SCFTs at large charge (also previously studied in [81][82][83]). Non-perturbative effects in large p or large λ are also straightforward to determine along the lines of [8], again using the explicit recursion above. Moreover, there are several interesting regimes of p and N scaling together. Work in these directions is ongoing [84].

Large N expansion
Let us proceed to develop the large N , fixed p expansion of the (normalised) maximal-trace overlaps g (N |max) p (s). 37 In the special case of N = 4 SYM, this may be solved for G2n(τ ) [50], and matches the color factor R (max) p (N ) defined in (A.2), up to a p 2 factor due to normalization differences (see footnote 7). In this case, G2n(τ ) ∝ y −2n , so the RHS is manifestly τ -independent.
In light of (6.10), the large N expansion follows from that of the functions F p (N, s) and g (N ) 2 . The large N expansion of the latter 39 has been discussed in [7,8], with the genus-g overlaps f (g) 2 (s) being of the general form depicted in (5.6).
On the other hand, the genuinely new information about the p-dependence of the maximal-trace overlaps lies entirely in the functions F p (N, s). For some fixed value of p, we note that both the sum in the numerator as well as the Pochhammer factor in the denominator of (6.11) give rise to even polynomials in N of degree p − 2. Therefore, the large N expansion of F p (N, s) starts at order N 0 and proceeds in even powers of 1/N . Explicitly, for the first few orders we find the expansion (6.20) At order N −2m (with m ≥ 1) we find an overall factor (s+1)(s−2) times a degree 2m−2 polynomial in s, whose coefficients depend only polynomially on p. We have verified that this simple structure persists to high order in m. Moreover, as a consequence of (6.13), every term in this expansion is invariant under s → 1 − s.
With this expansion of F p (N, s) at hand, it is clear that the normalised overlaps g (N |max) p (s) for general p admit the large N expansion Let us reiterate that the genus-g overlaps g (g|max) p (s) can be computed using (6.10) from the expansion (6.20) of F p (N, s) together with the known genus-g overlaps f (g) 2 (s) studied previously in [7,8]. Importantly, since the F p (N, s) are polynomial in s, the pole structure of g (g|max) p (s) is identical to their p = 2 counterparts, which we discussed in Section 5.1. As such, the general structure of the expansion in the 't Hooft or the very strongly coupled limit remains unchanged, the only effect being that the coefficients acquire some p-dependence. Let us illustrate this by explicitly developing the corresponding large N expansions.

't Hooft and very strongly coupled limits
We start by considering the 't Hooft limit of the normalised correlators G As mentioned previously, (anti-)instanton contributions are non-perturbatively suppressed in 1/N and hence this is an expansion in the zero-instanton sector only. 40 Note that since we consider the normalised correlator, there is a difference of N p in the overall power of N compared to (5.16).
Since the genus-zero overlaps are directly related to the p = 2 overlaps by an overall factor of p, 2 (s), it follows that the g = 0 term is simply given by where G (g=0) 2 (λ) was given in Section 5.3.1. At higher orders in 1/N such a factorisation no longer happens and some polynomial dependence on p will appear. For example at the next order, for g = 1, one has

(6.24)
Like the SPO integrated correlators considered earlier, this admits a Bessel representation, which we provide in Appendix G.2. Its strong coupling expansion is Lastly, let us also consider the large N expansion in the VSC limit. Following the same steps as described in Section 5.3.2, we arrive at This takes a remarkably simple form considering that this describes an infinite family of integrated correlators. Compared to the corresponding expansion of the (unnormalised) p = 2 correlator, as given in e.g. eq. (5.67) of [7], one notes that also higher-order even powers of 1/N are present, which comes from the non-trivial large N expansion of the ensemble averages G (N |max) p , cf. (6.7). 40 Again, we will drop the '0' subscript and it is understood that we work perturbatively in 1/N .

General p ansatz
the ensemble averages G

Future directions
The integrated correlators G (N |i,j) p (τ ) are remarkable observables: they are non-trivial functions of the complexified gauge coupling τ , and yet tractable enough not only to be exactly determined, but to be determined simply as polynomials in an appropriate functional basis. This basis is, moreover, one in which S-duality invariance is manifest. This link suggests that there may be lessons in the offing for how to better understand τ -dependence of N = 4 SYM observables at large.
These properties raise many conceptual and fundamental questions. We close by discussing some of these, but we first describe some more concrete open problems. There are several computations that are ripe for completion, which include the following: • Our formulas allow the explicit determination of non-perturbative corrections, both in 1/λ and 1/N , to G (N |i,j) p (τ ) in the large N limit. This was analyzed previously for the p = 2 case, where corrections in powers of the two scales (5.22) -identified with fundamental and D-string worldsheet instantons in AdS 5 × S 5 , respectively -are easily identified from the SL(2, Z) spectral integral, the two being related by S-duality. It would be interesting to understand the p-dependence of these corrections. We also note that the analysis of [8] did not compute all non-perturbative corrections in N , in particular leaving the sum over genera to future study; this is an interesting problem because the genus sum can in principle generate higher-order (e.g. "doubly") non-perturbative effects.  [34,35], another worthwhile computation.
• The ensemble averages of the (SPO|SPO) correlators through p = 5 were seen in (4.9)-(4. 10) to obey an intriguing pattern. Can a closed-form expression be derived, or perhaps inferred from more data points, for generic p?
• The SL(2, Z) spectral decomposition faciliates the study of the ensemble statistics of observables over the N = 4 conformal manifold M. The variance of the integrated correlators over M follows immediately from their explicit spectral overlaps and a formula derived in [8]: This was studied at both finite and large N in [8] for the p = 2 case. It would be interesting to ask how the general p integrated correlator statistics behave as well.
There are various generalizations of our work -e.g. extension to N = 4 SYM with other gauge groups, analysis of other families of integrated correlators G If so, new ideas are needed to determine the supersymmetric integration measure and to relate it to some as-yetunknown generalization of the localization methods based on N = 2 * partition functions. It should also be possible to extend the integrated correlator construction to SCFTs in two spacetime dimensions using localisation techniques.
Most fascinating are various foundational questions about these observables, centered around the question of why they exhibit such simplicity, and the relation to manifest S-duality invariance.
Our results for G (N |i,j) p (τ ) are characterized by strong uniformity in the half-BPS charge p. This calls for a deeper explanation. At large N , an obvious possible touchstone is the 10-dimensional "hidden conformal symmetry" that is known to govern at least some aspects of planar correlators at weak and strong coupling [23,86]. The origin and regime of applicability of the hidden conformal symmetry itself are still unknown. It is tempting to ask whether, and how, the integrated correlators G (N |i,j) p (τ ) are organized in some way by this symmetry, perhaps at large N only. If so, one may wonder whether the integrated correlators are subject to stronger constraints than the unintegrated ones.
Is there a natural bulk description of these objects besides as integrals over boundary points? That point of view is rather clumsy from the bulk spacetime perspective. The integrated correlators have an alternative formulation as derivatives of the free energy deformed by sources, but that quantity breaks the N = 4 superconformal symmetry. Perhaps there is a nicer picture which preserves the full symmetries without using the unintegrated correlator or free energy as an intermediary.
Of course, since the G (N |i,j) p (τ ) are soluble for all N and τ , we are really after their bulk string theory description. The formulas herein for G (N |i,j) p (τ ) may be regarded as the exact AdS 5 × S 5 quantum string theory answers. As emphasized in the introduction, the integrated correlators seem to be picking out some privileged subsector of the AdS 5 × S 5 four-point string amplitude. Whether this is a useful perspective that can be leveraged to teach us something more general about string theory, either in AdS 5 × S 5 or in 10-dimensional flat space, is an open question. It would obviously be very interesting to make this more explicit from the worldsheet point of view, perhaps first in various simplifying limits.
The central role of SL(2, Z) in our results underscores not only the power, but also the practical utility, of S-duality. That the integrated correlators are so streamlined in the SL(2, Z) eigenbasis suggests that the SL(2, Z) spectral decomposition may be generally useful in studying unprotected observables of N = 4 SYM such as operator dimensions or OPE coefficents. This deserves further investigation in the context of the superconformal bootstrap and integrability (and their combination [29,87,88]). In the present case, that the G (N |i,j) p (τ ) can be written down in such simple form for finite N and τ is partly 44 due to the absence of Maass cusp forms in their SL(2, Z) spectral decomposition. Why do the cusp form overlaps vanish? What is the physical principle that tells us when an SL(2, Z)-invariant observable in N = 4 SYM has vanishing cusp form overlap? Observables with this property are "non-chaotic" in a precise mathematical, if not (yet) physical, sense. It seems important to bridge this gap. This may help understand whether the mathematical notion of arithmeticity is a useful physical criterion in characterizing CFT observables. 45 Indeed, the radical simplicity of the integrated correlators -exact for all N and τ , with ingredients no more complicated than rational functions of N and p -deserves to be understood in the deepest terms possible. The fundamental origin of the recursion relations derived herein, generalizing the G (N ) 2 (τ ) recursion relation of [7], is a major open question. We made some suggestions about their provenance in Section 6.1. Can they be derived from a 6d (2,0) theory point of view? Can we bootstrap them directly without the input of localization? And finally, what are these equations telling us about unintegrated CFT data at finite N and τ , and can we somehow extract these data by cleverly combining the whole family of integrated correlators indexed by p?
The integrated correlators seem to be beautiful and rich observables, standing out even in the wonderful world of N = 4 SYM. They are, in our view, well-deserving of intensive further study.

A Colour factors
Here we list the colour factors R (i,j) p (N ) used in the main text. For the cases p = 2, 3, these are just given by the general formula (2.7). For p = 4, 5 we have The colour-factor of the two-point function of maximal-trace operators O (max) p = T 2,...,2 has been derived in [50]. In our conventions, for even p it reads At large N , the counting simplifies. In particular, computation of colour factors may be computed by large N factorization. Single-trace colour factors become slightly different context, see [92]. This is especially useful for computation of colour factors involving multi-trace operators T p 1 ,··· ,pn , defined in (2.2). The large N colour factor of T p 1 ,...,pn with n distinct indices is where p = p 1 + · · · + p n = p. For example, taking p 1 = 2 and p 2 = 3 reproduces the large N asymptotic of R

B Localization details B.1 Solving operator mixing on S 4 and integrated two-point functions
In this Section we focus on computing the denominator of (2.14) for p = 2, 3, 4, 5. Let us define where we recall that µ belongs to the index set (2.21). One of the main problems is to determine the vectors v i,µ p which encode the mixing of weight-p chiral operators of species i, denoted by O Next, we construct the matrix of connected two-point functions on S 4 for a given chiral operator of weight p with every other chiral operator of weight p − 2, p − 4, . . . with which it mixes. We will list below all such matrices for operators up to p = 5. First, we introduce some short-hand notation which will be convenient. In what follows we restrict our treatment and notation to singleand double-trace operators only, but everything can be straightforwardly generalized to arbitrary multi-trace insertions following the prescription of [50]. Let us introduce the notation O p and O p,q for single-and double-trace chiral operators, respectively, which couple to respective sources τ p and τ p,q . Two-point functions on S 4 are constructed via .
(B.2) 46 On S 4 , operators of even conformal dimension also mix with the identity operator. However, there is no need for us to consider this in the unmixing problem because the final results differ only by one-point functions on R 4 , which vanish.

Note that O p = O
(1) p is the weight-p single-trace chiral operator introduced in the main text, while O p,q is a single representative of the set of weight-(p + q) double-trace chiral operators O (j) p+q indexed by j. Let us also introduce a standard notation for the connected two-point functions on S 4 (one-point functions do not vanish on S 4 ), obtained by differentiation with respect to the free energy: and likewise for the other two-point functions in (B.2).
Through p = 5, we thus have the following mixing matrices: • p = 2: There is just one single-trace operator which does not mix with any other operator.
So there is nothing to unmix here and we have M 2 = O 2Ō2 c .
• p = 3: Again there is nothing to unmix, and we have M 3 = O 3Ō3 c .
• p = 4: There are two mixing matrices to consider here since operators of weight four have degeneracy two. For the single-trace weight-four operator O 5 as a function of N and τ by evaluating the matrix integrals using methods outlined in Section B.3. The next step is Gram-Schmidt orthogonalization of the matrices M (i) p . As mentioned above there is nothing to unmix in the p = 2, 3 cases. In particular this gives ∂ τ 3 = ∂ τ 3 . For p = 4, 5, the Gram-Schmidt process gives the new orthogonal vectors ∂/∂τ 4 (j) , ∂/∂τ 5 (j) (which brings insertions of chiral operators O where j = 1, 2 corresponds to single-and double-trace operators, respectively. We find the following results for the v-vectors: This data is sufficient to solve the sphere mixing problem for p = 2, 3, 4, 5 and therefore allows to compute both the numerator and denominator of (2.14).
As for the denominator, assembling everything, we find the following results for the matrix of integrated two-point functions D (N |i,j) p on R 4 as defined in (B.1): in a τ,τ independent single-particle operator basis (discussed in Section 2.1) but for our current purpose we do not need such an expression. It is straightforward to check that the N -dependence of these results are precisely those given by computing the flat space two-point function using Wick contractions in the free N = 4 SYM with gauge group SU (N ).

B.2 Integrated four-point functions
In this Section, we compute the numerator of (2.14) for p = 2, 3, 4, 5 in the zero instanton sector.
where we recall that µ belongs to the index set (2.21). We first note that ∂ τ µ ∂τ ν ∂ 2 m log Z N gives the following five terms that compute the connected integrated correlator Single mass derivatives of the N = 2 * partition function do not appear above since those terms vanish for vanishing mass. At vanishing sources and with µ, ν refering to single-trace species p 1 , q 1 respectively, the above derivatives give the following insertions in the matrix model integral: For µ, ν refering to multi-trace species, as explained in Section 2.2, the precise insertion depends upon the species under consideration. The first line in (B.11) is easy to see from (2.15). On the other hand, insertions from mass derivatives resulting from (2.15) are complicated. A pragmatic approach to evaluate them is in the large y expansion. One way to do this is to expand the quantity in small a i . This expansion is represented by the sum in (B.11). Performing the matrix integrals term-by-term in the sum gives an asymptotic expansion for the integrated correlators at large y. So computing ∂ τ µ ∂τ ν ∂ 2 m log Z boils down to computing the matrix integral (2.20) with insertions dictated by the derivatives in (B.11) and their multi-trace generalization. For p ≤ 5, the most general such insertion involves computing the following 47 where the expectation value is computed in the Gaussian ensemble (2.20) whose details we present in Section B.3. Putting together various ingredients we find the following results for C (N |i,j) p (in the expressions below the matrix elements D (N |1,1) p were computed in (B.8)) • p = 2: 64y 5 π 5 − 117117N 2 11N 4 + 40N 2 + 9 ζ(13) 64y 6 π 6 + · · · (B.14b) 47 For correlators involving just single-trace operators, (B.13) is still the most general insertion for generic p, but when multi-trace operators are involved, there can be several more insertions of the type N i=1 a q i .
• p = 4: Having computed these results, we remark that while for p = 2 and 3, where the non-planar corrections start at four-loop order, for higher p, non-planarity sets in at lower loop orders. Specifically, for the matrix of p = 4, 5 correlators, we see that the non-planar corrections already appear at two loops which is in agreement with known results in the literature [9,12].
With these results, we now have all the ingredients to assemble the full matrix of integrated correlators G (N |i,j) p (τ ) defined in (2.14) for p = 2, 3, 4, 5 at finite N in the zero-instanton sector in a weak-coupling expansion.

B.3 Expectation values in the Gaussian (special) unitary ensemble
Here we demonstrate how we evaluate expectation values of the type (B.13) in the Gaussian special unitary ensemble (2.20). We found [93][94][95][96] as useful references for the material presented in this section.
One of the technical complications in evaluating such expectation values is to properly deal with the delta function constraint on the eigenvalues in (2.20). When the insertions are just differences of eigenvalues then such a constraint simply gives an N -dependent prefactor with the expectation values to be evaluated in the usual Gaussian unitary ensemble. However, when the insertions are not differences of eigenvalues (like in (B.13)), a systematic way to handle the constraint is to use the integral representation of the delta function, δ( a i ) = dµ 2π e iµ i a i . 48 After doing a linear change of variables and completing the squares in the exponent, we get that (B.13) evaluates to (where for convenience we defined α := 2πy) The it is easiest to handle these expectation values by first providing sources to them and performing the µ integral, i.e., we first do the following replacement above A practical way to evaluate these expectation values is as follows. Using the orthogonal Hermite polynomials [97], it is straightforward to show that (B.21) can be written as the following sum over a certain determinant A convenient way to represent the sum over the determinant in (B.22) is to think of it as sum over products of traces of products of matrices Q a i ,a j . Working out the traces one finds that the result takes the form of an exponential times an even polynomial in γ i : where the coefficients C (k) i 1 ,i 2 ···i k (N ) are functions of N and completely symmetric in the indices i 1 , i 2 · · · i k . We do not know them in closed form as a function of N for generic indices but nevertheless can determine for a specific set of indices, which is all that we really need for calculations. (C.4)

D One-instanton results
Here we provide results for the one-instanton contribution to the integrated correlators G (i,j) p (N ; τ ) for p = 2, 4, 5 in a weak-coupling expansion. The result for p = 3 was already presented in the main text in Section 3.3.
where Φ (l) (u, v) = − 1 x−x φ (l) (x ,x ) with φ (l) (x,x) being the well-known l-loop ladder integrals. The C (i,j) n,p (N ) with n = 1, 2, 3 are certain colour-factors encoding the N -dependence of the one-and 49 See also the recent work [41] for a comparison up to four-loop order and interesting connections to periods of Feynman graphs. p and is given on a case-by-case basis in [9] (for p = 3 see their eqs. (44)- (46), for p = 4 see Table 2, and for p = 5, 6 see Tables 6 and 7, respectively). 50 Furthermore, this matching also indicates that we have correctly solved the operator mixing problem on S 4 , which we performed explicitly up to p = 5 in Appendix B.1 and for the maximaltrace correlators for all even p in Section 6, thus extending the results of [6] beyond the single-trace sector and to finite N .

G.1 Genus-zero spectral overlap for general p
First we derive the genus-zero, generic p spectral overlap (5.3). The idea is to make use of the known genus-zero correlator for generic p, recast it into the spectral form and read off the overlap. The genus-zero correlator is given by