###### Abstract

Integrated correlation functions in ${\cal N}=4$ supersymmetric Yang–Mills
theory with gauge group $SU(N)$ can be expressed in terms of the localised
$S^{4}$ partition function, $Z_{N}$, deformed by a mass $m$. Two such cases
are ${\cal C}_{N}=({\rm
Im\,}\tau)^{2}\partial_{\tau}\partial_{\bar{\tau}}\partial_{m}^{2}\log
Z_{N}|_{m=0}$ and ${\cal H}_{N}=\partial_{m}^{4}\log Z_{N}|_{m=0}$, which are
modular invariant functions of the complex coupling $\tau$. While ${\cal
C}_{N}$ was recently written in terms of a two-dimensional lattice sum for any
$N$ and $\tau$, ${\cal H}_{N}$ has only been evaluated up to order $1/N^{3}$
in a large-$N$ expansion in terms of modular invariant functions with no known
lattice sum realisation. Here we develop methods for evaluating ${\cal H}_{N}$
to any desired order in $1/N$ and finite $\tau$. We use this new data to
constrain higher loop corrections to the stress tensor correlator, and give
evidence for several intriguing relations between ${\cal H}_{N}$ and ${\cal
C}_{N}$ to all orders in $1/N$. We also give evidence that the coefficients of
the $1/N$ expansion of ${\cal H}_{N}$ can be written as lattice sums to all
orders. Lastly, these large $N$ and finite $\tau$ results are used to
accurately estimate the integrated correlators at finite $N$ and finite
$\tau$.

QMUL-PH-23-20

Relations between integrated correlators in

${\cal N}=4$ Supersymmetric Yang–Mills Theory

Luis F. Alday1, Shai M. Chester2,3,4, Daniele Dorigoni5, Michael B. Green6,
Congkao Wen7

1 Mathematical Institute, University of Oxford, Woodstock Road, Oxford OX2
6GG, UK

2 Jefferson Physical Laboratory, Harvard University, Cambridge, MA 02138, USA

3 CMSA, Harvard University, Cambridge, MA 02138, USA

4 Blackett Laboratory, Imperial College, Prince Consort Road, London, SW7 2AZ,
UK

5 Centre for Particle Theory, Durham University, Stockton Road, Durham DH1
3LE, UK

6 DAMTP, Wilberforce Road, Cambridge CB3 0WA, UK

7 Centre for Theoretical Physics, Queen Mary University of London, London, E1
4NS, UK

###### Contents

  1. 1 Introduction
    1. 1.1 Brief review of integrated correlators
    2. 1.2 Summary of results
    3. 1.3 Outline
  2. 2 Free energy at large $N$
  3. 3 Half-integer powers of $1/N$ at finite $\tau$
  4. 4 Integer powers of $1/N$ at finite $\tau$
    1. 4.1 Holomorphic cusp forms and their cancellation
    2. 4.2 Lattice sum representation
  5. 5 Numerical estimate of free energy at finite $N$
  6. 6 Constraints on stress tensor correlator
  7. 7 Discussion
  8. A Matrix model details
    1. A.1 $O(N^{2})$ at large $\lambda$
    2. A.2 $O(N^{0})$ at large $\lambda$
  9. B Generalised Eisenstein series
  10. C Properties of lattice sums
    1. C.1 A family of lattice sum integrals
    2. C.2 Modular local polynomials: definition and properties
  11. D From lattice sums to generalised Eisenstein series

##  1 Introduction

Four-dimensional maximally supersymmetric Yang-Mills (${\cal N}=4$ SYM) theory
with gauge group $SU(N)$ is of great interest for a variety of reasons. It is
a non-trivial four-dimensional superconformal gauge theory parameterised by
the complex coupling
$\tau\equiv\frac{\theta}{2\pi}+i\frac{4\pi}{g_{{}_{YM}}^{2}}$ with
Montonen–Olive duality group $SL(2,{\mathbb{Z}})$ [1];111In this paper we will
focus on the theory with gauge group $SU(N)$. S-duality involving more general
gauge groups is discussed by Goddard, Nuyts and Olive in [2]. it is a toy
model of QCD and an integrable system in the planar limit. Importantly for our
present considerations it provides a non-perturbative formulation of type IIB
string theory in an $AdS_{5}\times S^{5}$ background [3]222Similar holographic
dualities exist for other gauge groups [4].. In this holographic
interpretation, type IIB string theory with coupling $\tau_{s}=\chi+i/g_{s}$
and string length scale $\ell_{s}$ in a background $AdS_{5}\times S^{5}$ of
length scale $L$ is identified with ${\cal N}=4$ SYM with gauge group $SU(N)$,
where $\tau_{s}=\tau$ and $(L/\ell_{s})^{4}=g_{{}_{YM}}^{2}N$. Therefore, the
gravity description at small $\ell_{s}$ is dual to SYM at large $N$ and finite
$\tau$.333At large $N$ and finite $\tau$, the bulk is weakly curved because
$\ell_{s}$ is small, but the string theory is still strongly coupled because
$\tau_{s}$ is finite. One can further take $g_{s}$ to be weakly coupled, which
is then dual to the familiar large-$N$ large-$\lambda\equiv g_{{}_{YM}}^{2}N$
limit of SYM. Super-graviton scattering amplitudes in string theory on
$AdS_{5}\times S^{5}$ are then dual to correlation functions of the stress
tensor multiplet. While this duality in principle gives a non-perturbative
definition of string theory in terms of a well defined CFT, in practice it is
difficult to use the duality to study string theory in the gravity regime
since the CFT is strongly coupled at large $N$. We thus need to study ${\cal
N}=4$ SYM non-perturbatively to make progress in studying quantum gravity via
AdS/CFT.

###  1.1 Brief review of integrated correlators

Certain properties of ${\cal N}=4$ SYM can be obtained analytically for all
values of $N$ and $\tau$ using supersymmetric localisation (see [5] for a
review). This method was originally used to compute non-local quantities such
as supersymmetric partition functions and Wilson loops [6]. More recently, it
was understood [7, 8] that certain integrals of the correlator of four
superconformal primary operators in the ${\cal N}=4$ stress tensor multiplet
are obtained from the following derivatives of the $S^{4}$ partition function
$Z_{N}(m,\tau)$ deformed by the $\mathcal{N}=2$ preserving mass $m$:

$\begin{split}\mathcal{C}_{N}(\tau)\equiv\frac{1}{4}\Delta_{\tau}\partial_{m}^{2}\log
Z_{N}(m,\tau)\big{|}_{m=0}\,,\qquad\mathcal{H}_{N}(\tau)&\equiv\partial_{m}^{4}\log
Z_{N}(m,\tau)\big{|}_{m=0}\,,\end{split}$ (1.1)

where $\Delta_{\tau}\equiv 4\tau_{2}^{2}\partial_{\tau}\partial_{\bar{\tau}}$
is the hyperbolic laplacian and $\tau=\tau_{1}+i\tau_{2}$ with
$\tau_{1}={\theta}/{(2\pi)}$, and $\tau_{2}={4\pi}/{g_{{}_{YM}}^{2}}$. The
quantities ${\cal C}_{N}$ and ${\cal H}_{N}$ are identified with integrals
over the insertion points of the operators in the four-point correlator. The
choice of integration measure distinguishes the choice of integrated
correlators.

The relations (1.1) are useful because $Z_{N}(m,\tau)$ can be computed using
supersymmetric localisation in terms of an $(N-1)$-dimensional matrix model
integral [6]. The integrand is a product of a classical contribution that
depends on $g_{{}_{YM}}$, a one-loop contribution that is $\tau$-independent,
and a sum over an infinite number of instanton contributions that depend on
$\tau$ and $\bar{\tau}$ [9, 10, 11, 12]. In order to use these integrated
correlators to study the gravity regime of string theory via AdS/CFT, we need
to compute these $(N-1)$-dimensional integrals explicitly for large $N$.

The perturbative sector of $\mathcal{C}_{N}(\tau)$ was first computed for both
finite $N$ and in a large-$N$ expansion in [13]. The instanton sectors were
then also computed at large $N$ in [14], and the large-$N$ expansion of
$\mathcal{C}_{N}(\tau)$ was found to take the following form

$\begin{split}\mathcal{C}_{N}(\tau)&=\frac{N^{2}}{4}-\frac{3\sqrt{N}}{2^{4}\
}E({\scriptstyle{3\over
2}};\tau)+\frac{45}{2^{8}\sqrt{N}}E({\scriptstyle{5\over
2}};\tau)+\frac{1}{{N}^{\frac{3}{2}}}\left[-\frac{39}{2^{13}}E({\scriptstyle{3\over
2}};\tau)+\frac{4725}{2^{15}}E({\scriptstyle{7\over 2}};\tau)\right]\\\
&\quad+\frac{1}{{N}^{\frac{5}{2}}}\left[-\frac{1125}{2^{16}}E({\scriptstyle{5\over
2}};\tau)+\frac{99225}{2^{18}}E({\scriptstyle{9\over
2}};\tau)\right]+O(N^{-\frac{7}{2}})\,,\end{split}$ (1.2)

where $E(s;\tau)$ is the non-holomorphic Eisenstein series, that may be
defined by the two-dimensional lattice sum444The normalisation of $E(s;\tau)$
differs from that of [14] in a manner that eliminates factors of $\pi$ in the
coefficients of the series in (1.2) and subsequent expressions.

$\displaystyle
E(s;\tau)=\frac{1}{\pi^{s}}\sum_{(m,n)\,\neq\,(0,0)}\frac{\tau_{2}^{s}}{|m+n\tau|^{2s}}\,,$
(1.3)

which is the unique $SL(2,\mathbb{Z})$ invariant solution of moderate growth
to the homogeneous Laplace eigenvalue equation

$\begin{split}\left(\Delta_{\tau}-s(s-1)\right)E(s;\tau)=0\,.\end{split}$
(1.4)

It was shown in [15, 16] that, for all $N$ and $\tau$, ${\cal C}_{N}(\tau)$
can be written as a two-dimensional lattice sum,

$\displaystyle{\cal
C}_{N}(\tau)=\sum_{(m,n)\in\mathbb{Z}^{2}}\int_{0}^{\infty}dt\,B_{N}(t)\exp\Big{(}-t\pi\frac{|m+n\tau|^{2}}{\tau_{2}}\Big{)}\,,$
(1.5)

where $B_{N}(t)={\cal Q}_{N}(t)/(t+1)^{2N+1}$ and ${\cal Q}_{N}(t)$ is a
polynomial of degree $(2N+1)$. Although the expression for $B_{N}(t)$ was
arrived at in [16] by analysing the perturbative and instanton contributions
for many values of $N$, a derivation of its form was given in [17] based on
the construction of a generating series, ${\cal
C}(z;\tau)\equiv\sum_{N=1}^{\infty}{\cal C}_{N}(\tau)z^{N}$. The expression
(1.5) was shown to satisfy a Laplace difference equation

$\begin{split}\Delta_{\tau}{\cal C}_{N}(\tau)=N^{2}\big{[}{\cal
C}_{N+1}(\tau)-2\,{\cal C}_{N}(\tau)+{\cal C}_{N-1}(\tau)\big{]}-N\big{[}{\cal
C}_{N+1}(\tau)-{\cal C}_{N-1}(\tau)\big{]}+2\,{\cal
C}_{N}(\tau)\,.\end{split}$ (1.6)

Using the initial condition ${\cal C}_{1}(\tau)=0$, this recursion relation
can be used to efficiently compute all ${\cal C}_{N}(\tau)$ in terms of ${\cal
C}_{2}(\tau)$.

The $1/N$ expansion (1.2) is non-convergent, and the explicit modular
invariant trans-series containing the large-$N$ completion was obtained by
analysing the generating series [17]. In the large-$\lambda$ limit, the
large-$N$ completion behaves as $e^{-2\sqrt{\lambda}}$ as obtained first using
resurgence [16, 18]. Many other results concerning the integrated correlator
${\cal C}_{N}(\tau)$ as well as its extensions for other gauge groups and more
general 1/2-BPS operators can be found in the literature, see e.g. [19, 20,
21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32].

Much less is known about the second correlator ${\cal H}_{N}(\tau)$, which is
the focus of this paper. While the perturbative sector was computed at finite
$N$ in [8], the large-$N$ expansion of these perturbative terms is much more
challenging and only the first few orders were computed in that study. In
[33], it was shown that instanton terms could be computed more easily at large
$N$, which were used to compute the first few terms at large $N$ and finite
$\tau$. The resulting series has the following form

$\begin{split}&{\cal H}_{N}(\tau)=6N^{2}+{\cal H}^{h}_{N}(\tau)+{\cal
H}^{i}_{N}(\tau)+\dots\,,\end{split}$ (1.7)

where ellipsis denotes terms that are exponentially suppressed in $N$. The
series ${\cal H}^{h}_{N}$ contains terms with half-integer powers of $1/N$,

$\begin{split}&{\cal H}^{h}_{N}(\tau)=\sqrt{N}E({\scriptstyle{3\over
2}};\tau)-\frac{9}{2\sqrt{N}}E({\scriptstyle{5\over
2}};\tau)+\frac{1}{{N}^{\frac{3}{2}}}\left[\frac{117}{2^{8}}E({\scriptstyle{3\over
2}};\tau)-\frac{3375}{2^{10}}E({\scriptstyle{7\over 2}};\tau)\right]\\\
&\quad+\frac{1}{{N}^{\frac{5}{2}}}\left[\frac{675}{2^{10}}E({\scriptstyle{5\over
2}};\tau)-\frac{33075}{2^{12}}E({\scriptstyle{9\over
2}};\tau)\right]+O(N^{-\frac{7}{2}})\,,\\\ \end{split}$ (1.8)

where the coefficients are rational multiples of non-holomorphic Eisenstein
series, just as in the case of ${\cal C}_{N}$ in (1.2). The series ${\cal
H}^{i}_{N}$ contains integer powers of $1/N$. The terms up to $O(1/N^{3})$
were evaluated in [33] and take the following form

$\displaystyle{\cal H}^{i}_{N}(\tau)=C_{0}+\frac{27}{2^{3}N}{\cal
E}(4,{\scriptstyle{3\over 2}},{\scriptstyle{3\over
2}};\tau)+\frac{1}{N^{2}}\left[C_{1}-\frac{14175}{704}{\cal
E}(7,{\scriptstyle{5\over 2}},{\scriptstyle{3\over
2}};\tau)+\frac{1215}{88}{\cal E}(5,{\scriptstyle{5\over
2}},{\scriptstyle{3\over 2}};\tau)\right]$ (1.9) $\displaystyle+$
$\displaystyle\,\frac{1}{N^{3}}\Big{[}\alpha_{3}{\cal
E}(4,{\scriptstyle{3\over 2}},{\scriptstyle{3\over
2}};\tau)+\sum_{s=6,8,9}[\alpha_{s}{\cal E}(s,{\scriptstyle{3\over
2}},{\scriptstyle{3\over 2}};\tau)+\beta_{s}{\cal E}(s,{\scriptstyle{5\over
2}},{\scriptstyle{5\over 2}};\tau)+\gamma_{s}{\cal E}(s,{\scriptstyle{7\over
2}},{\scriptstyle{3\over 2}};\tau)]\Big{]}\\!+O(N^{-4})\,,$

where $\alpha_{s},\beta_{s},\gamma_{s}$ are known rational constants, while
the values of the constants $C_{0}$ and $C_{1}$ require further knowledge of
the perturbative sector. The function ${\cal E}(s,s_{1},s_{2};\tau)$ is a
generalised Eisenstein series, which is an $SL(2,\mathbb{Z})$ invariant
solution of the inhomogeneous Laplace eigenvalue equation 555Note that in [33]
a different convention was used in which $\hat{{\cal
E}}(s,s_{1},s_{2};\tau)=-\pi^{s_{1}+s_{2}}{\cal E}(s+1,s_{1},s_{2};\tau)$,
where $\hat{{\cal E}}(s,s_{1},s_{2};\tau)$ denotes the generalised Eisenstein
series of [33].

$\begin{split}\left(\Delta_{\tau}-s(s-1)\right){\cal
E}(s,s_{1},s_{2};\tau)=E(s_{1};\tau)E(s_{2};\tau)\,.\end{split}$ (1.10)

An example of such a function, ${\cal
E}(4,{\scriptstyle\frac{3}{2}},{\scriptstyle\frac{3}{2}};\tau)$, entered as
the coefficient of $d^{6}R^{4}$ in the low energy expansion of the type IIB
effective action in [34, 35] where it was expressed as a four-dimensional
lattice sum and its Fourier expansion was studied in [36]. More general
examples of such lattice sums were obtained in [37] 666We will later see that
such lattice sums require more careful definitions than those given in [35]
and [37].. Starting at eigenvalue $s=6$ it was shown, in [38] for cases when
$s_{1}$ and $s_{2}$ are integers and in [39] when $s_{1}$ and $s_{2}$ are
half-integers, that the Fourier-mode decomposition of ${\cal
E}(s,s_{1},s_{2};\tau)$ in general contains terms which are related to
holomorphic cusp forms. This fact will play an important in our analysis.

We will now summarise the main results of this paper, which will shed much
more light on the structure of the large-$N$ expansion of ${\cal H}_{N}$ and
hint at its possible structure for all values of $N$.

###  1.2 Summary of results

Most of the explicit results in this paper apply to the large-$N$ expansion of
${\cal H}_{N}$. This involves developing techniques for determining the
explicit form of the perturbative terms in ${\cal H}_{N}$ to any desired order
in $1/N$. When combined with the large-$N$ instanton expressions obtained
following [33], this determines the full ${\cal H}_{N}$ to any desired order
in $1/N$.

Several interesting relations between ${\cal H}_{N}$ and ${\cal C}_{N}$ will
emerge from these explicit large-$N$ results. For example, we find that the
half-integer power terms, ${\cal H}^{h}_{N}$ in (1.7) satisfy a Laplace
difference equation

$\displaystyle\Delta_{\tau}{{\cal H}_{N}^{h}(\tau)}=$
$\displaystyle\,N^{2}(N^{2}-1)\Big{[}\frac{{\cal
H}_{N+1}^{h}(\tau)}{(N+1)^{2}}-2{{\cal H}_{N}^{h}(\tau)\over
N^{2}}+\frac{{\cal H}_{N-1}^{h}(\tau)}{(N-1)^{2}}\Big{]}$ (1.11)
$\displaystyle+\frac{16N}{N^{2}-1}\Big{[}\frac{3N^{3}}{4}-(N-1)^{2}{\cal
C}_{N+1}(\tau)+3N{\cal C}_{N}(\tau)+(N+1)^{2}{\cal C}_{N-1}(\tau)\Big{]}\,,$

which allows us to efficiently compute ${\cal H}^{h}_{N}$ to all orders in
$1/N$. This recursion relation is similar to that of ${\cal C}_{N}$ in (1.6),
except now both integrated correlators appear on the right hand side, and
furthermore we cannot use (1.11) to compute the full ${\cal H}_{N}$ at finite
$N$ because ${\cal H}_{N}^{i}$ does not satisfy this recursion relation.

We also derived several new results concerning the integer power terms ${\cal
H}_{N}^{i}$. Firstly, the $\tau$-independent terms $C_{0}$, $C_{1}$, $\dots$
that appear at even powers of $1/N$ in (1.9) can now be evaluated. These
coefficients can be used to constrain contact-term ambiguities of loop
corrections to the holographic correlator. In particular, we will use the
$C_{0}$ term to give a consistency check on the unique contact-term ambiguity
of the one-loop graviton exchange term [40, 41], which was originally fixed
using ${\cal C}_{N}$ in [13]; while $C_{1}$ can be used to fix one of the four
contact-term ambiguities of the 2-loop graviton exchange term [42, 43].

Although we suspect that the full integrated correlator, ${\cal H}_{N}$,
satisfies some differential equation analogous to the Laplace difference
equation satisfied by ${\cal C}_{N}$, this remains to be understood. However,
we find that the laplacian acting on ${\cal H}_{N}^{i}$ is related to ${\cal
C}_{N}$ via a nonlinear relation

$\begin{split}\Delta_{\tau}{\cal
H}^{i}_{N}(\tau)\big{|}_{E}=\frac{96}{N^{2}-1}\,\left({\cal
C}_{N}(\tau)-\frac{N^{2}}{4}\right)^{2}\,,\end{split}$ (1.12)

where $\Delta_{\tau}{\cal H}_{N}|_{E}$ means that the expression is restricted
to the terms that are bilinear in non-holomorphic Eisenstein series that
appear when acting with the laplacian on the generalised Eisenstein series
${\cal E}$, as appeared in the right-hand side of (1.10). Although this
relation cannot recursively fix each order of ${\cal H}_{N}^{i}$ as (1.11) did
for ${\cal H}_{N}^{h}$, we will see that it does strongly restrict the space
of ${\cal E}$ that arise in ${\cal H}^{i}_{N}$ at each order in $1/N$
expansion.

As we commented earlier, an individual generalised Eisenstein series ${\cal
E}(s,s_{1},s_{2};\tau)$ with eigenvalue $s\geq 6$ in general contains
contributions which originate from holomorphic cusp forms of modular weight
$2s$. However, we will show that the combinations of ${\cal
E}(s,s_{1},s_{2};\tau)$ that appear at each order of $1/N$ in (1.9) and its
extension to higher orders are very special ones in which the holomorphic cusp
forms contributions cancel. Closely related to this, the combinations
appearing in the $1/N$ expansion are precisely the combinations of generalised
Eisenstein series that admit four-dimensional lattice sum representations.

More precisely we will see that at each order in $1/N$ up to order $O(N^{-7})$
that we have checked explicitly, the combination of generalised Eisenstein
series that appears has the form of an integral of a four-dimensional lattice
sum of the type that has appeared in [37, 44, 36] in the study of the low
energy expansion of type IIB string amplitudes. Such specific sums of
generalised Eisenstein series (modulo linear combinations of single Eisenstein
series) have the form

$\displaystyle{\cal
E}^{w}_{i,j}(\tau)=\\!\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\neq
0}}\int_{0}^{\infty}d^{3}t\,B_{i,j}^{w}(t_{1},t_{2},t_{3})\,\exp\Big{(}\\!-\frac{\pi}{\tau_{2}}\sum_{i=1}^{3}t_{i}|p_{i}|^{2}\Big{)}\,,$
(1.13)

where $p_{i}=m_{i}+n_{i}\tau$ with $m_{i},n_{i}\in{\mathbb{Z}}$. The function
$B_{i,j}^{w}(t_{1},t_{2},t_{3})$ is a symmetric function of
$(t_{1},t_{2},t_{3})$ which satisfies a Laplace eigenvalue equation777Here and
in subsequent expressions we use the shorthand notation $B(t)\equiv
B(t_{1},t_{2},t_{3})$.

$\displaystyle(\Delta_{t}-s(s-1))B_{i,j}^{w}(t)=0\,,$ (1.14)

within the domain $t_{i}>0$, where $s=3i+j+1$ with $i,j\in\mathbb{N}$ and
$\Delta_{t}$ is a laplacian expressed in terms of derivatives with respect to
$t_{i}$, as given in (4.22). The value of the index $w\in{\mathbb{N}}$
determines the specific combinations of generalised Eisenstein series that
arise in ${\cal E}_{i,j}^{w}$. Some details of this expression will be
described later. In previous work, such as [37, 44] the function
$B_{i,j}^{w}(t)$ was expressed in the form $V(t)^{w}A_{i,j}(\rho(t))$, where
$\rho=\rho_{1}+i\rho_{2}$ is a simple complex function of $t_{1},t_{2},t_{3}$,
as will be reviewed in appendix C.

The motivation for combining generalised Eisenstein series into expressions
such as (1.13) can be seen already at $O(N^{-3})$ in (1.9). As we will show,
the coefficients $\alpha_{s},\beta_{s},\gamma_{s}$ (with $s=6,8,10$) of the
generalised Eisenstein series at order $1/N^{3}$ in (1.9) satisfy a constraint
for each value of $s$. We will further see that this pattern continues to all
the orders in $1/N$ that we have evaluated – at each order the coefficients
are in the basis defined by the integrated lattice sums, ${\cal E}_{i,j}^{a}$,
which is a smaller basis than that defined by generalised Eisenstein series.
Since the non-holomorphic Eisenstein series in ${\cal H}_{N}^{i}$ are already
written as lattice sums in (1.3), this shows that the full ${\cal H}_{N}$ can
be written as a lattice sum to all orders in $1/N$, up to at least $O(N^{-7}$)
which we have verified explicitly.

All the explicit results concern properties of ${\cal H}_{N}$ at large $N$ and
finite $\tau$, and we have not yet found an analytic expression at finite $N$
and $\tau$, as was found for ${\cal C}_{N}$ in [16]. However, we will show
that the asymptotic large $N$ expansion for both ${\cal H}_{N}$ and ${\cal
C}_{N}$ can be used to accurately estimate these quantities at finite $N$ for
all $\tau$ in the fundamental domain. We do this by combining the exact $N$
expressions for the perturbative terms from [13, 8] with the instanton sectors
of the $1/N$ expressions in this work truncated to a certain order in $1/N$.
The result is shown to numerically match $N=2,3$ expressions for all $\tau$ to
good precision, which have been evaluated numerically from the matrix model
integrals in [45].

###  1.3 Outline

The rest of this paper is organised as follows. In Section 2 we review the
localised expression for ${\cal H}_{N}$ with emphasis on the zero instanton
sector. We then determine the large-$\lambda$ expansion (where
$\lambda=g_{{}_{YM}}^{2}N$ is the ’t Hooft coupling) of the correlator,
explicitly displaying the first few terms in the $1/\lambda$ expansion at
order $N^{2}$ and $N^{0}$.

In Section 3 we will use the above methods to compute the coefficients of
half-integer powers of $1/N$ that arise in ${\cal H}^{h}_{N}$. We will see
that this expression satisfies a recursion relation, which generates ${\cal
H}^{h}_{N}$ to all orders in $1/N$ and at finite $\tau$, thereby enabling us
to determine the expression for ${\cal H}_{N}^{h}$ to all orders. The
structure of this expression is remarkably similar to that of ${\cal C}_{N}$.

In Section 4 we consider terms with integer powers of $1/N$ in the large-$N$
expansion of ${\cal H}^{i}_{N}$. We will see that the coefficients of each
term in the $1/N$ expansion is a particular combination of generalised
Eisenstein series in which contributions from holomorphic cusp forms cancel.
We further show that ${\cal H}^{i}_{N}$ admits a four-dimensional lattice sum
representation at each order in the large-$N$ expansion.

In Section 5 we show how to use the large-$N$ expansion of ${\cal H}_{N}$ and
${\cal C}_{N}$ to accurately estimate the value of the integrated correlator
at finite values of $N$ and generic values of $\tau$.

In Section 6, we show how our new results for the $\tau$-independent terms in
${\cal H}_{N}$ can be used to constrain higher loop corrections to the stress
tensor correlator at loop level.

We conclude in Section 7 with a discussion of how these results might suggest
further insights. Technical details are relegated to the appendices, and we
also include various results in an attached Mathematica file.

##  2 Free energy at large $N$

In this section, we will review the form of the $S^{4}$ partition function
$Z_{N}(m,\tau)$ of ${\cal N}=2^{*}$ SYM, which reduces to the partition
function of ${\cal N}=4$ SYM with gauge group $SU(N)$ in the $m\to 0$ limit.
We will be particularly interested in computing the large-$N$ expansion of
$\mathcal{H}_{N}$. The emphasis in this section is on the large-$N$ and
large-$\lambda$ expansion, which is not sensitive to instantons. In subsequent
sections we will then combine this new data as well as instanton results
obtained following [33] to fix the large-$N$ finite-$\tau$ expansion to very
high orders.

The ${\cal N}=2^{*}$ SYM $S^{4}$ partition function can be determined using
supersymmetric localisation and can be expressed as a matrix model integral
[6]

$\begin{split}Z_{N}(m,\tau)\\!=\\!\int
d^{N-1}a\,\big{(}\prod_{i<j}a_{ij}^{2}\big{)}e^{-\frac{8\pi^{2}}{g_{\text{YM}}^{2}}\sum_{i}a_{i}^{2}}\frac{\prod_{i\neq
j}H(a_{ij})}{H(m)^{N-1}\\!\prod_{i\neq j}H(a_{ij}+m)}\left\lvert
Z_{\text{inst}}(m,\tau,a_{ij})\right\rvert^{2}\,,\end{split}$ (2.1)

where $a_{ij}\equiv a_{i}-a_{j}$. The integration is over $N$ real variables
$a_{i}$, $i=1,\ldots,N$, subject to the constraint $\sum_{i}a_{i}=0$. The
function $H(z)=e^{-(1+\gamma)z^{2}}G(1+iz)G(1-iz)$ is the product of two
Barnes G-functions, and $Z_{\text{inst}}$ is the contribution from instantons
localised at the poles of $S^{4}$. We can write the instanton partition
function as

$\begin{split}Z_{\text{inst}}(m,\tau,a_{ij})=\sum_{k=0}^{\infty}q^{k}Z_{\text{inst}}^{(k)}(m,a_{ij})\,,\qquad{\rm
with}\qquad q=e^{2\pi i\tau}\,,\end{split}$ (2.2)

where $Z_{\text{inst}}^{(k)}(m,a_{ij})$ represents the contribution of the
$k$-instanton sector and is normalised such that
$Z_{\text{inst}}^{(0)}(m,a_{ij})=1$. Explicit expressions for
$Z_{\text{inst}}^{(k)}(m,a_{ij})$ can be found in [11, 12]. We can then take
mass derivatives to get the perturbative contribution

$\begin{split}\partial_{m}^{4}\log
Z_{N}\big{|}_{m=0}^{pert}=-12\zeta(3)+\sum_{i,j}\langle
K^{\prime\prime\prime}(a_{ij})\rangle+3\sum_{i,j,k,l}\left[\langle
K^{\prime}(a_{ij})K^{\prime}(a_{kl})\rangle-\langle
K^{\prime}(a_{ij})\rangle\langle
K^{\prime}(a_{kl})\rangle\right]\,,\end{split}$ (2.3)

where $K(z)\equiv-\frac{H’(z)}{H(z)}$ and expectation values are taken with
the Gaussian matrix model measure so that

$\begin{split}\langle{\cal O}(a_{ij})\rangle=\frac{1}{\mathcal{N}}\int
d^{N}a\,\delta\big{(}\sum_{i}a_{i}\big{)}\,\big{(}\prod_{i<j}a_{ij}^{2}\big{)}e^{-\frac{8\pi^{2}N}{\lambda}\sum_{i}a_{i}^{2}}\,{\cal
O}(a_{ij})\,,\end{split}$ (2.4)

where the normalisation factor $\mathcal{N}$ is chosen such that $\langle
1\rangle=1$.

The instanton terms in (2.1) are contained in Nekrasov partition function,
$Z_{\text{inst}}(m,\tau,a_{ij})$, whose small-$m$ expansion has been studied
in [14, 33]. In [33], it was already shown how to compute these instanton
terms to any order in $1/N$, so in the following we will focus on the
perturbative terms (2.3) which we denote by $\partial_{m}^{4}\log
Z_{N}|_{m=0}^{pert}$.888Note that the perturbative terms have $k=0$ (zero
instanton number), but there is also an infinite set of $k=0$ instanton/anti-
instanton terms, which do not contribute to the perturbative sector.

To evaluate $\partial_{m}^{4}\log Z_{N}|_{m=0}^{pert}$, we write the function
$K^{\prime}(z)$ in terms of the integral representation [46]

$\begin{split}K^{\prime}(z)=-\int_{0}^{\infty}d\omega\frac{2\omega[\cos(2\omega
z)-1]}{\sinh^{2}\omega}\,,\end{split}$ (2.5)

so that the perturbative sector in (2.3) can be written as

$\begin{split}\partial_{m}^{4}\log
Z_{N}\big{|}_{m=0}^{pert}=&-12\zeta(3)+\int_{0}^{\infty}d\omega\frac{8\omega^{3}\mathcal{I}(\omega)}{\sinh^{2}\omega}+\int_{0}^{\infty}d\omega_{1}d\omega_{2}\frac{12\omega_{1}\omega_{2}\mathcal{J}(\omega_{1},\omega_{2})}{\sinh^{2}\omega_{1}\sinh^{2}\omega_{2}}\,,\end{split}$
(2.6)

where we define the 2-body and 4-body expectation values

$\begin{split}\mathcal{I}(\omega)\equiv\sum_{i,j}\langle e^{2i\omega
a_{ij}}\rangle\,,\qquad\mathcal{J}(\omega_{1},\omega_{2})\equiv\sum_{i,j,k,l}\left[\langle
e^{2i\omega_{1}a_{ij}}e^{2i\omega_{2}a_{kl}}\rangle-\langle
e^{2i\omega_{1}a_{ij}}\rangle\langle
e^{2i\omega_{2}a_{kl}}\rangle\right]\,.\end{split}$ (2.7)

In [13, 8], it was shown how these expectation values can be computed to any
order at large $N$ and finite $\lambda$ using topological recursion, where the
expansion takes the form

$\begin{split}\mathcal{I}(\omega)&=N^{2}\mathcal{I}_{2}(\omega)+N^{0}\mathcal{I}_{0}(\omega)+N^{-2}\mathcal{I}_{-2}(\omega)+\dots\,,\\\
\mathcal{J}(\omega_{1},\omega_{2})&=N^{2}\mathcal{J}_{2}(\omega_{1},\omega_{2})+N^{0}\mathcal{J}_{0}(\omega_{1},\omega_{2})+N^{-2}\mathcal{J}_{-2}(\omega_{1},\omega_{2})+\dots\,.\\\
\end{split}$ (2.8)

The $1/N^{2}$ coefficients $\mathcal{I}_{i}(\omega)$ are expressed as products
of two Bessel functions with the same argument. For instance, the first two
terms are

$\begin{split}\mathcal{I}_{2}(\omega)=\frac{4\pi^{2}J_{1}(\frac{\sqrt{\lambda}\omega}{\pi}){}^{2}}{\lambda\omega^{2}}\,,\qquad\mathcal{I}_{0}(\omega)=\frac{\sqrt{\lambda}\omega
J_{1}(\frac{\sqrt{\lambda}\omega}{\pi}){}J_{2}(\frac{\sqrt{\lambda}\omega}{\pi}){}}{12\pi}\,,\end{split}$
(2.9)

while higher terms are given in [13]. The coefficients
$\mathcal{J}_{i}(\omega_{1},\omega_{2})$ can be expressed as products of four
Bessel functions. For instance,

$\displaystyle\mathcal{J}_{2}(\omega_{1},\omega_{2})$
$\displaystyle=\frac{8\pi
J_{1}(\frac{\sqrt{\lambda}\omega_{1}}{\pi})J_{1}(\frac{\sqrt{\lambda}\omega_{2}}{\pi})}{\sqrt{\lambda}(\omega_{2}^{2}-\omega_{1}^{2})}\left[\textstyle\omega_{1}J_{0}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)-\omega_{2}J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{0}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)\right]\,,$
$\displaystyle\mathcal{J}_{0}(\omega_{1},\omega_{2})$
$\displaystyle=\mathcal{J}^{\text{non-
fac}}_{0}(\omega_{1},\omega_{2})+\mathcal{J}^{\text{fac}}_{0}(\omega_{1},\omega_{2})\,,$
(2.10)

where we defined the auxiliary functions

$\displaystyle\mathcal{J}^{\text{non-fac}}_{0}(\omega_{1},\omega_{2})$
$\displaystyle\equiv\frac{\lambda}{12\pi^{2}(\omega_{2}^{2}{-}\omega_{1}^{2})}\\!\left[\omega_{1}\textstyle
J_{0}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\textstyle
J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)-\omega_{2}\textstyle
J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\textstyle
J_{0}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)\right]$ (2.11)
$\displaystyle\qquad\qquad\qquad\qquad\times\left[\omega_{2}^{3}\textstyle
J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\textstyle
J_{2}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)+\omega_{1}^{3}\textstyle
J_{2}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\textstyle
J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)\right]$
$\displaystyle+\frac{\lambda\omega_{1}\omega_{2}(\omega_{1}^{2}+\omega_{2}^{2})}{2\pi^{2}(\omega_{2}^{2}{-}\omega_{1}^{2})^{2}}\left[\omega_{1}\textstyle
J_{0}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\textstyle
J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)-\omega_{2}\textstyle
J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\textstyle
J_{0}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)\right]$
$\displaystyle\qquad\qquad\qquad\qquad\times\left[\omega_{2}\textstyle
J_{0}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\textstyle
J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)-\omega_{1}\textstyle
J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\textstyle
J_{0}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)\right]\,,$
$\displaystyle\mathcal{J}^{\text{fac}}_{0}(\omega_{1},\omega_{2})$
$\displaystyle\equiv-\frac{7\lambda\omega_{2}^{2}J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right){}^{2}J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right){}^{2}}{12\pi^{2}}-\frac{\lambda\omega_{2}^{2}J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right){}^{2}J_{0}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right){}^{2}}{2\pi^{2}}$
(2.12)
$\displaystyle-\frac{7\lambda\omega_{1}^{2}J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right){}^{2}J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right){}^{2}}{12\pi^{2}}+\frac{\lambda\omega_{1}\omega_{2}J_{0}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{0}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)}{2\pi^{2}}$
$\displaystyle+\frac{\lambda\omega_{1}\omega_{2}J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{2}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{2}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)}{12\pi^{2}}-\frac{\lambda\omega_{1}^{2}J_{0}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right){}^{2}J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right){}^{2}}{2\pi^{2}}\,.$

Note that at order $O(N^{0})$ we have separated the terms into
$\mathcal{J}^{\text{fac}}_{0}$, which factorises in terms of $\omega_{2}$ and
$\omega_{1}$, and $\mathcal{J}^{\text{fac}}_{0}$, which does not factorise.
Higher order terms take a similar form; they can be found in [8] and up to
order $O(N^{-8})$ in the attached Mathematica file.

The large-$\lambda$ expansion of the 2-body term $\mathcal{I}(\omega)$ in
(2.6) is straightforward at all orders in $1/N^{2}$. As shown in [7], we start
by writing bilinears of Bessel functions in the Mellin-Barnes form

$\begin{split}J_{\mu}(x)J_{\nu}(x)=\oint\frac{ds}{2\pi
i}\frac{\Gamma(-s)\Gamma(2s+\mu+\nu+1)\left(\frac{x}{2}\right)^{\mu+\nu+2s}}{\Gamma(s+\mu+1)\Gamma(s+\nu+1)\Gamma(s+\mu+\nu+1)}\,,\end{split}$
(2.13)

where the contour separates the poles of each Gamma function. We then perform
the $\omega$ integrals in (2.6) using the identity

$\begin{split}\int_{0}^{\infty}d\omega\
\frac{\omega^{a}}{\sinh^{2}\omega}=\frac{1}{2^{a-1}}\Gamma(a+1)\zeta(a)\,,\end{split}$
(2.14)

and then close the contour to the left to get an expansion in $1/\lambda$. For
instance, using the explicit expressions in (2.9) we get

$\begin{split}&\int_{0}^{\infty}d\omega\frac{8\omega^{3}\mathcal{I}(\omega)}{\sinh^{2}\omega}=N^{2}\Big{[}\frac{16\pi^{2}}{\lambda}-\frac{32\pi^{2}}{\lambda^{\frac{3}{2}}}+\frac{24\pi^{2}\zeta(3)}{\lambda^{\frac{5}{2}}}+O(\lambda^{-\frac{7}{2}})\Big{]}\\\
&\qquad\qquad\qquad\qquad+N^{0}\Big{[}\frac{4\pi^{2}\sqrt{\lambda}}{15}-\frac{13\pi^{2}}{16\lambda^{\frac{3}{2}}}-\frac{75\pi^{2}\zeta(3)}{32\lambda^{\frac{5}{2}}}+O(\lambda^{-\frac{7}{2}})\Big{]}+O(N^{-2})\,,\\\
\end{split}$ (2.15)

where the expansion at each order in $1/N^{2}$ including a finite number of
positive powers of $\lambda$, plus an infinite number of negative powers.

The large-$\lambda$ expansion of the 4-body terms
$\mathcal{J}(\omega_{1},\omega_{2})$ in (2.6) is much harder to evaluate
because the dependence on $\omega_{1}$ and $\omega_{2}$ does not in general
factorise. As a result, even if we write each pair of Bessel functions in
Mellin-Barnes form, it is still difficult to perform the
$\omega_{1},\omega_{2}$ integrals. In [8], the 4-body terms were instead
computed at large $\lambda$ by numerically evaluating the
$\omega_{1},\omega_{2}$ integrals in (2.6) for many values of finite
$\lambda$, and then using a numerical fit to guess the large $\lambda$ values.
This method only worked for the lowest few orders in $1/\lambda$ at each order
in $1/N^{2}$, where it gave

$\displaystyle\int_{0}^{\infty}$ $\displaystyle
d\omega_{1}d\omega_{2}\frac{12\omega_{1}\omega_{2}\mathcal{J}(\omega_{1},\omega_{2})}{\sinh^{2}\omega_{1}\sinh^{2}\omega_{2}}=$
$\displaystyle
N^{2}\Big{[}6+\frac{96\zeta(3)}{\lambda^{\frac{3}{2}}}-\frac{288\zeta(5)}{\lambda^{\frac{5}{2}}}-\frac{144\zeta(3)^{2}}{\lambda^{3}}-\Big{(}\frac{16\pi^{2}}{\lambda}-\frac{32\pi^{2}}{\lambda^{\frac{3}{2}}}+\frac{24\pi^{2}\zeta(3)}{\lambda^{\frac{5}{2}}}\Big{)}+O(\lambda^{-\frac{7}{2}})\Big{]}$
$\displaystyle+N^{0}\Big{[}\Big{(}2-\frac{4\pi^{2}}{15}\Big{)}\lambda^{\scriptstyle\frac{1}{2}}+O(\lambda^{0})\Big{]}+N^{-2}\Big{[}\Big{(}\frac{\pi^{2}}{504}-\frac{1}{60}\Big{)}\lambda^{\frac{3}{2}}+O(\lambda)\Big{]}$
(2.16)
$\displaystyle+N^{-4}\Big{[}-\frac{\lambda^{3}}{120960}+O(\lambda^{\frac{5}{2}})\Big{]}+O(N^{-6})\,.$

Note that the $\pi^{2}$ terms here exactly cancel those from the 2-body term
in (2.15) up to the order shown.

We now show how to evaluate the 4-body terms at large $\lambda$ analytically
to any order. Let us begin with the leading $N^{2}$ term
$\mathcal{J}_{2}(\omega_{1},\omega_{2})$ in (2):

$\displaystyle I^{(0)}(\lambda)$
$\displaystyle\equiv\int_{0}^{\infty}\\!d\omega_{1}d\omega_{2}\frac{12\omega_{1}\omega_{2}\mathcal{J}_{2}(\omega_{1},\omega_{2})}{\sinh^{2}\omega_{1}\sinh^{2}\omega_{2}}$
(2.17)
$\displaystyle=\int_{0}^{\infty}\\!d\omega_{1}d\omega_{2}\frac{96\pi\omega_{1}\omega_{2}J_{1}\\!\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\\!J_{1}\\!\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)}{\sqrt{\lambda}(\omega_{2}^{2}{-}\omega_{1}^{2})\sinh^{2}\omega_{1}\sinh^{2}\omega_{2}}\left[\textstyle\omega_{1}J_{0}\\!\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\\!J_{1}\\!\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)\\!-\omega_{2}J_{1}\\!\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)\\!J_{0}\\!\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)\right]\,.$

We start by applying the Bessel kernel identity

$\begin{split}\frac{\omega_{1}J_{0}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)-\omega_{2}J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{0}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)}{\omega_{2}^{2}-\omega_{1}^{2}}=\sum_{\ell=1}^{\infty}\frac{4\ell\pi
J_{2\ell}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{2\ell}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)}{\omega_{1}\omega_{2}\sqrt{\lambda}}\,,\end{split}$
(2.18)

to write (2.17) as an infinite sum of terms

$\begin{split}I^{(0)}(\lambda)=\sum_{\ell=1}^{\infty}\frac{384\pi^{2}\ell}{\lambda}\int_{0}^{\infty}d\omega_{1}d\omega_{2}\frac{J_{1}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)J_{2\ell}\left(\frac{\sqrt{\lambda}\omega_{1}}{\pi}\right)}{\sinh^{2}\omega_{1}}\frac{J_{1}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)J_{2\ell}\left(\frac{\sqrt{\lambda}\omega_{2}}{\pi}\right)}{\sinh^{2}\omega_{2}}\,.\end{split}$
(2.19)

Since this expression now factorises in $\omega_{1}$ and $\omega_{2}$, we can
convert each pair of Bessel functions to Mellin-Barnes form using (2.13),
shift the contour integral variables by $1-\ell$, and perform the $\omega_{1}$
and $\omega_{2}$ integrals using (2.14), which gives

$\begin{split}I^{(0)}(\lambda)=\oint\frac{d{s}}{2\pi i}\frac{d{t}}{2\pi
i}p({s},{t})c({s},{t})\,,\qquad
c({s},{t})=\sum^{\infty}_{\ell=1}c_{\ell}({s},{t})\,,\end{split}$ (2.20)

where the $\ell$-independent prefactor is

$\begin{split}p({s},{t})=\frac{3\lambda^{{s}+{t}+2}\sin(\pi{s})\sin(\pi{t})\Gamma(2{s}+4)^{2}\Gamma(2{t}+4)^{2}\zeta(2{s}+3)\zeta(2{t}+3)}{2^{4{s}+4{t}+3}\pi^{2({s}+{t}+2)}}\,,\end{split}$
(2.21)

while the $\ell$-dependent terms are

$\begin{split}c_{\ell}({s},{t})=\frac{\ell(\ell-{s}-3)(\ell-{t}-3)\Gamma(\ell-{s}-3)\Gamma(\ell-{s}-1)\Gamma(\ell-{t}-3)\Gamma(\ell-{t}-1)}{\pi^{2}\Gamma(\ell+{s}+2)\Gamma(\ell+{s}+3)\Gamma(\ell+{t}+2)\Gamma(\ell+{t}+3)}\,.\end{split}$
(2.22)

We can now perform the contour integrals by closing the contours to the left,
which includes poles where ${s}$ and ${t}$ are each negative numbers, as well
as poles where the sum ${s}+{t}$ is a negative number. The details of this
pole counting are summarised in appendix A. The result matches the leading
large-$N$ term in (2), but now we can extend the series to any order in
$1/\lambda$ systematically.

The 4-body term at order $N^{0}$ is more complicated, as
$\mathcal{J}_{0}(\omega_{1},\omega_{2})$ in (2) now includes twice the non-
factorised denominator $\omega_{2}^{2}-\omega_{1}^{2}$.999This additional
complication is only relevant for the $O(N^{0})$ term. To make the
$\omega_{1},\omega_{2}$ factorise, we must therefore apply the Bessel kernel
identity (2.18) twice, so the resulting expression will be a double infinite
sum. In appendix A, we show how the above method is modified to compute the
large-$\lambda$ expansion in this case. After combining these 4-body results
with the 2-body expression (2.15), we find that that the expansion of the full
perturbative sector (2.6) is given by

$\displaystyle\partial_{m}^{4}\log$ $\displaystyle
Z_{N}\big{|}_{m=0}^{pert}=N^{2}\Big{[}6+\frac{96\zeta(3)}{\lambda^{\frac{3}{2}}}-\frac{288\zeta(5)}{\lambda^{\frac{5}{2}}}-\frac{144\zeta(3)^{2}}{\lambda^{3}}-\frac{3375\zeta(7)}{4\lambda^{7/2}}-\frac{1080\zeta(3)\zeta(5)}{\lambda^{4}}+O(\lambda^{-\frac{9}{2}})\Big{]}$
$\displaystyle+\Big{[}2\sqrt{\lambda}-\frac{96\zeta(3)+23}{10}-\frac{3\zeta(3)}{\lambda}+\frac{117\zeta(3)}{16\lambda^{\frac{3}{2}}}-\frac{45\zeta(5)}{8\lambda^{2}}+O(\lambda^{-\frac{5}{2}})\Big{]}+O(N^{-2})\,.$
(2.23)

We note that the $\pi^{2}$ terms from the 2-body contribution (2.15) cancel at
all orders, and the $\lambda^{0}N^{0}$ term also receives a contribution from
the first term in (2.6). All higher $1/N^{2}$ corrections can be evaluated
similarly. In the attached Mathematica file, we show explicit results to
higher orders in both $1/N^{2}$ and $1/\lambda$. They take a similar form to
the above expressions, except the number of positive powers of $\lambda$ grows
with each order in $1/N^{2}$, and all subsequent $\lambda^{0}$ terms are
simple fractions. In the following sections we will use these perturbative
data, combined with the instanton results obtained following [33], to
determine the coefficients of the large-$N$ finite-$\tau$ expansion.

##  3 Half-integer powers of $1/N$ at finite $\tau$

We now consider the large-$N$ and finite-$\tau$ expansion. As discussed around
(1.7), this expansion can be divided into terms with half-integer powers of
$N$ contained in ${\cal H}_{N}^{h}(\tau)$, as well as terms with integer power
of $N$ in ${\cal H}_{N}^{i}(\tau)$. In this section we will consider the form
of ${\cal H}^{h}(\tau)$. We can detemine the coefficients of the powers of
$1/N$ using the perturbative sector of the localisation expression (2.3), as
well as the instanton terms. As discussed earlier, the instanton terms were
obtained in [33], so we will mainly discuss perturbative terms.

The perturbative sector at large $N$ and finite $\tau$ only depends on
$\tau_{2}$, and can be obtained from the large-$N$ and large-$\lambda$
expansion in (2) by simply setting $\lambda={4\pi N}/{\tau_{2}}$ and
reorganising the large-$N$ expansion. In [33] it was observed that ${\cal
H}_{N}^{h}(\tau)$ is entirely written in terms of non-holomorphic Eisenstein
series $E(s;\tau)$, which have the well-known Fourier-mode decomposition

$\begin{split}E(s;\tau)&=\frac{2\zeta(2s)}{\pi^{s}}{\tau_{2}^{s}}+\frac{2\Gamma(s-\frac{1}{2})}{\pi^{s-\frac{1}{2}}\Gamma(s)}\zeta(2s-1)\tau_{1}^{1-s}\\\
&{}+\frac{4\sqrt{\tau_{2}}}{\Gamma(s)}\sum_{k\neq 0}\left\lvert
k\right\rvert^{s-{\scriptstyle\frac{1}{2}}}\sigma_{1-2s}(k)\,K_{s-\frac{1}{2}}(2\pi\left\lvert
k\right\rvert\tau_{2})\,e^{2\pi ik\tau_{1}}\,,\end{split}$ (3.1)

where the divisor sum $\sigma_{p}(k)$ is defined as
$\sigma_{p}(k)=\sum_{d|k}d^{p}$ with $d>0$, and $K_{s-\frac{1}{2}}$ is a
Bessel function of the second kind. The contribution to ${\cal
H}_{N}^{h}(\tau)$ coming from a particular non-holomorphic Eisenstein series,
$E(s;\tau)$, is in fact uniquely fixed by the form of the two $k=0$ terms in
the Fourier series (3.1). We can therefore determine ${\cal H}_{N}^{h}(\tau)$
to any order in $1/N$ and finite $\tau$ by simply comparing to the
perturbative data from (2) (and the higher order terms given in the attached
Mathematica notebook). We can then double check that the $k>0$ terms in (3.1)
are consistent with the $k>0$ terms derived in [33] (as well as new data
obtained following the methods outlined in [33]). We find that ${\cal
H}_{N}^{h}(\tau)$ takes the following form

$\begin{split}{\cal
H}_{N}^{h}(\tau)=\sum_{j,\ell=0}^{\infty}N^{{{\scriptstyle\frac{1}{2}}-\ell-2j}}d_{\ell}^{j}\,E({\scriptstyle\frac{3}{2}}+\ell;\tau)\,,\end{split}$
(3.2)

where the coefficients $d_{\ell}^{j}$ for the lowest few $j$ are found to be

$\displaystyle d_{\ell}^{0}$
$\displaystyle=-\frac{(\ell+1)(\ell+3)\Gamma\left(\ell-\frac{1}{2}\right)\Gamma\left(\ell+\frac{3}{2}\right)^{2}}{4^{\ell-1}\pi^{\ell+\frac{9}{2}}\Gamma(\ell+1)}\,,$
(3.3) $\displaystyle d_{\ell}^{1}$
$\displaystyle=\frac{(\ell+1)(2\ell+13)(\ell(2\ell+7)+9)\Gamma\left(\ell+\frac{1}{2}\right)\Gamma\left(\ell+\frac{3}{2}\right)\Gamma\left(\ell+\frac{5}{2}\right)}{3\cdot
4^{\ell+2}\pi^{\ell+\frac{11}{2}}(\ell+2)\Gamma(\ell+1)}\,,$ (3.4)
$\displaystyle d_{\ell}^{2}$
$\displaystyle=-\frac{(\ell+1)^{2}(\ell(2\ell+7)+15)(4\ell(5\ell+52)+219)\Gamma\left(\ell+\frac{1}{2}\right)\Gamma\left(\ell+\frac{3}{2}\right)\Gamma\left(\ell+\frac{11}{2}\right)}{45\cdot
4^{\ell+5}\pi^{\ell+\frac{13}{2}}\Gamma(\ell+4)}\,,$ (3.5) $\displaystyle
d_{\ell}^{3}$
$\displaystyle=\frac{(\ell+1)^{2}(\ell(2\ell+7)+21)(560\ell^{4}{+}10752\ell^{3}{+}59272\ell^{2}{+}124800\ell{+}75771)\Gamma\\!\left(\ell+\frac{1}{2}\right)\\!\Gamma\\!\left(\ell+\frac{3}{2}\right)\\!\Gamma\\!\left(\ell+\frac{15}{2}\right)}{2835\cdot
4^{\ell+8}\pi^{\ell+\frac{15}{2}}\Gamma(\ell+5)}\,.$

The integrated correlator ${\cal C}_{N}(\tau)$ can be written in a similar
form

$\begin{split}{\cal
C}_{N}(\tau)=\frac{N^{2}}{4}+\sum_{j,\ell=0}^{\infty}N^{{{\scriptstyle\frac{1}{2}}-\ell-2j}}{\tilde{d}}_{\ell}^{j}\,E({\scriptstyle\frac{3}{2}}+\ell;\tau)\,,\end{split}$
(3.6)

where the lowest few coefficients are [16]

$\displaystyle\tilde{d}_{\ell}^{0}$
$\displaystyle=\frac{(\ell+1)\Gamma\left(\ell-\frac{1}{2}\right)\Gamma\left(\ell+\frac{3}{2}\right)\Gamma\left(\ell+\frac{5}{2}\right)}{4^{\ell+1}\pi^{\ell+\frac{9}{2}}\Gamma(\ell+1)}\,,$
(3.7) $\displaystyle\tilde{d}_{\ell}^{1}$
$\displaystyle=-\frac{(\ell+1)^{2}(2\ell+13)\Gamma\left(\ell+\frac{3}{2}\right)\Gamma\left(\ell+\frac{5}{2}\right)^{2}}{3\cdot
2^{2\ell+7}\pi^{\ell+\frac{11}{2}}\Gamma(\ell+3)}\,,$ (3.8)
$\displaystyle\tilde{d}_{\ell}^{2}$
$\displaystyle=\frac{(\ell+1)^{2}\left(20\ell^{2}+208\ell+219\right)\Gamma\left(\ell+\frac{3}{2}\right)\Gamma\left(\ell+\frac{5}{2}\right)\Gamma\left(\ell+\frac{11}{2}\right)}{45\cdot
2^{2\ell+13}\pi^{\ell+\frac{13}{2}}\Gamma(\ell+4)}\,,$ (3.9)
$\displaystyle\tilde{d}_{\ell}^{3}$
$\displaystyle=-\frac{(\ell+1)^{2}\left(560\ell^{4}+10752\ell^{3}+59272\ell^{2}+124800\ell+75771\right)\Gamma\\!\left(\ell+\frac{3}{2}\right)\\!\Gamma\\!\left(\ell+\frac{5}{2}\right)\\!\Gamma\\!\left(\ell+\frac{15}{2}\right)}{2835\cdot
2^{2\ell+19}\pi^{\ell+3}\Gamma(\ell+\frac{15}{2})}\,.$

While the coefficients in both (3.3) and (3.7) may seem quite complicated, it
turns out for all $j$ they are simply related by

$\begin{split}\left(3+8\ell+4\ell^{2}\right){d_{\ell}^{j}}=-32\left(3+6j+7\ell+2\ell^{2}\right){\tilde{d}_{\ell}^{j}}\,.\end{split}$
(3.10)

These simple polynomials in $\ell,j$ can be interpreted as coming from
derivatives with respect to $\tau$ and $N$ acting on (3.2) and (3.6), which
implies the formal relation

$\begin{split}\Delta_{\tau}{\cal
H}_{N}^{h}(\tau)=-8\left(3-3\partial_{N}+2\Delta_{\tau}\right){\cal
C}_{N}(\tau)+6N^{2}-12N\,,\end{split}$ (3.11)

where the polynomial in $N$ cancels the $N^{2}/4$ term in ${\cal
C}_{N}(\tau)$. We can eliminate $\partial_{N}$, whose interpretation is
confusing for integer $N$, by acting with $\partial_{N}$ on the recursion
relation for ${\cal C}_{N}(\tau)$ in (1.6) and combining with shifted versions
of (3.11), which yields the recursion relation given in (1.11). This recursion
relation, or equivalently (3.10), can be used to generate ${\cal
H}_{N}^{h}(\tau)$ to any desired order in $1/N$ from the known ${\cal
C}_{N}(\tau)$ [16]. The new recursion relation essentially puts ${\cal
H}^{h}_{N}(\tau)$ on the same level as $C_{N}(\tau)$, which is known to any
order in $1/N$. However, the full second integrated correlator ${\cal
H}_{N}(\tau)$ does not obey the recursion relation (1.11). In particular, the
integer powers in $1/N$ expansion ${\cal H}^{i}_{N}(\tau)$ behave differently,
as we will discuss in the next section.

##  4 Integer powers of $1/N$ at finite $\tau$

We will now discuss ${\cal H}_{N}^{i}(\tau)$, which contains the integer
powers of $1/N$ in the large-$N$ expansion. The terms up to $1/N^{3}$ were
obtained in [33] and are reproduced in (1.9). The coefficients of these terms
are sums of generalised Eisenstein series ${\cal E}(s,s_{1},s_{2};\tau)$ with
$s_{1},s_{2}\in{\mathbb{N}}+{\scriptstyle\frac{1}{2}}$, as well as
$\tau$-independent pure numbers. We find that the generalisation of (1.9) to
arbitrary orders in $1/N$ takes the following form,

$\begin{split}{\cal
H}_{N}^{i}(\tau)=\sum_{j=0}^{\infty}\frac{C_{j}}{N^{2j}}+\sum_{j=0}^{\infty}\sum_{\ell=0}^{\infty}\sum_{r=0}^{\lfloor\ell/2\rfloor}\sum_{s=4+\ell,6+\ell,\dots}^{3\ell+6j+4}\frac{\alpha_{j,\ell,r,s}}{N^{{1+\ell+2j}}}{\cal
E}(s,{\scriptstyle\frac{3}{2}}+r,{\scriptstyle\frac{3}{2}}+\ell-r;\tau)\,.\end{split}$
(4.1)

The $\tau$-independent terms $C_{j}$ are associated with supergravity
contributions, and terms containing ${\cal E}(s,s_{1},s_{2};\tau)$ are due to
the stringy effects. The choice of the eigenvalues $s$ is to ensure that the
terms at order $1/N^{n}$ obeys a boundary condition such that the
corresponding (integrated) $d^{2(2n+1)}R^{4}$ vertices in the holographic
string amplitudes contribute only up to genus-$(2n+1)$, which was previously
observed in [33]. Whereas the choice of the sources, which are labelled by
$s_{1},s_{2}$, can be understood as a consequence of the relation (1.12).

The $\tau$-independent terms $C_{j}$ in (1.9) were not computed in [33],
because they can only be fixed from the perturbative sector. Using our new
perturbative data, we find that

$\begin{split}C_{0}=-\frac{96\zeta(3)+23}{10}\,,\end{split}$ (4.2)

while for $j>0$ they all appear to be rational numbers. For instance,
$C_{1}=-\frac{1}{10}\,,$ and $C_{2}=\frac{1}{14}$. It is straightforward (but
laborious) to compute higher orders from the perturbative results.

As for the ${\cal E}(s,s_{1},s_{2};\tau)$ terms in ${\cal H}_{N}^{i}(\tau)$,
we first note that the ${\cal E}(s,s_{1},s_{2};\tau)$ have the Fourier
expansion

$\begin{split}{\cal E}(s,s_{1},s_{2};\tau)=\sum_{n,m=0}^{\infty}{\cal
E}^{(n,m)}(s,s_{1},s_{2};\tau_{2})q^{n}\bar{q}^{m}\,.\end{split}$ (4.3)

The coefficients ${\cal E}^{(n,m)}(s,s_{1},s_{2};\tau_{2})$ correspond to
terms in the $n$-instanton/$m$-anti-instanton sector with total instanton
number $k=n-m$. Further discussion of ${\cal
E}^{(n,m)}(s,s_{1},s_{2};\tau_{2})$ can be found in appendix B. The
coefficients of terms in the $1/N$ expansion of ${\cal H}_{N}^{i}(\tau)$ can
in principle be fixed directly by comparing with explicit results obtained
from localisation for certain choices of $(m,n)$. This was the procedure used
in [33] where all the terms up to $O(1/N^{3})$ were determined. This however
becomes more and more challenging as we consider higher-order terms. In order
to simplify the computation, we will assume that ${\cal H}_{N}^{i}(\tau)$
satisfies certain conditions at all orders in $1/N$ that are satisfied by low-
order terms. We will then check that the resulting expression indeed
reproduces the localisation result. The conditions that will be imposed are
the following:

  * •

Up to order $1/N^{3}$ we observed in (1.12) that the action of $\Delta_{\tau}$
on ${\cal H}_{N}^{i}(\tau)$ results in terms bilinear in Eisenstein series
that are proportional to those in the square of ${\cal C}_{N}(\tau)$. Assuming
that this continues to hold at higher orders in $1/N$ provides strong
constraints on the coefficients of the generalised Eisenstein series that
arise in the $1/N$ expansion.

  * •

The second condition is that the coefficients $\alpha_{j,\ell,r,s}$ of the
expansion (4.1) in generalised Eisenstein series are related for each
eigenvalue $s$ in such a manner that contributions coming from holomorphic
cusp forms cancel. As will be discussed in the next subsection individual
${\cal E}(s,s_{1},s_{2};\tau)$ with $s\geq 6$ generically have such cusp
forms, but they are not present in the localisation result of ${\cal H}_{N}$.
Such cancellation was present (but not appreciated) in the order $1/N^{3}$
expression in [33]. Here it will be extended to arbitrary order in $1/N$.
Furthermore, we will show that the combination of ${\cal
E}(s,s_{1},s_{2};\tau)$ for which the cusp forms cancel is one for which there
is lattice sum representation.

These supplementary conditions significantly reduce the number of independent
coefficients $\alpha_{j,\ell,r,s}$ at each order in $1/N$. For example, at
order $1/N^{3}$, although there are apparently $10$ coefficients only $4$ of
these are independent once the above extra conditions are satisfied. At order
$1/N^{4}$ there appear to be $13$ independent coefficients, but this is
reduced to $7$ after imposing the above conditions, and at order $1/N^{5}$ the
$34$ apparently independent coefficients are reduced to $16$. Having
constrained the coefficients in the $1/N$ expansion in this manner, we will
then use new perturbative and non-perturbative data obtained in this paper
from the localisation computation to fix these coefficients. In this way, we
have determined ${\cal H}_{N}^{i}(\tau)$ completely up to order $O(N^{-7})$.
It is important to stress that we have checked that the results agree with
much more localisation data than is needed to fix the coefficients, This
procedure therefore provides non-trivial consistency checks on the properties
observed in the lower-order results.

Some examples of the coefficients of powers of $1/N$ will be displayed later
in this section, but in general the explicit expressions at higher orders in
$1/N$ are rather lengthy, so we relegate them to the attached Mathematica
file.

###  4.1 Holomorphic cusp forms and their cancellation

In [33], then generalised in [47], it was shown how one can construct
particular solutions, ${\cal E}_{{\rm p}}(s,s_{1},s_{2};\tau)$, to the Laplace
equation (1.10). However, in the earlier works [48, 38] it was demonstrated
that when $s_{1},s_{2}\in{\mathbb{N}}$ (relevant for modular graph functions),
such particular solutions do not always correspond to modular invariant
functions. In such cases it was shown that the failure of modularity of ${\cal
E}_{{\rm p}}(s,s_{1},s_{2};\tau)$ requires the addition of a suitable
homogeneous solution.

Using the results of [39], a similar argument can be given for the case
$s_{1},s_{2}\in\mathbb{N}+\frac{1}{2}$ that is of present interest. It can be
shown that the generalised Eisenstein series ${\cal E}$ differs from the
particular solution ${\cal E}_{{\rm p}}$ constructed via [33, 47] in the
purely instantonic or anti-instantonic sector, i.e. for $(n,m)=(n,0)$ or
$(n,m)=(0,m)$, and takes the form

${\cal E}^{(n,0)}(s,s_{1},s_{2};\tau_{2})={\cal E}^{(n,0)}_{{\rm
p}}(s,s_{1},s_{2};\tau_{2})+\sum_{\Delta\in{\cal
S}_{2s}}\lambda_{\Delta}(s,s_{1},s_{2})\frac{a_{\Delta}(n)}{n^{s}}\sqrt{n\tau_{2}}e^{2\pi
n\tau_{2}}K_{s-\frac{1}{2}}(2\pi n\tau_{2})\,.$ (4.4)

The asymptotic expansion at the cusp of ${\cal E}^{(n,0)}_{{\rm
p}}(s,s_{1},s_{2};\tau_{2})$ produces a perturbative expansion in
$\tau_{2}^{-1}$, which are the terms that we matched to the instanton sector
of ${\cal H}_{N}^{i}(\tau)$ here and in [33]. The second term in (4.4) was
missed in the earlier analysis of [33, 47]101010As already mentioned, and will
be shown in detail later, the holomorphic cusp forms cancel out in the final
expression of the integrated correlator, therefore the final result presented
in [33] is not affected by this.  and it involves a sum over the vector space
${\cal S}_{2s}$ of all holomorphic cusp forms with modular weight $2s$, i.e.
we need to sum over all $\Delta\in{\cal S}_{2s}$ of the form

$\Delta(\tau)\equiv\sum_{n=1}^{\infty}a_{\Delta}(n)q^{n}\in{\cal S}_{2s}\,,$
(4.5)

with Hecke normalisation $a_{\Delta}(1)=1$, and again $q=\exp(2\pi i\tau)$.

Following [48, 38] we consider the contribution of a fixed cusp form,
$\Delta\in{\cal S}_{2s}$, to ${\cal E}(s,s_{1},s_{2};\tau)$ coming from the
second term in (4.4) when summed over all instanton and all anti-instanton
sectors, which takes the form

$H_{\Delta}(\tau)\equiv\sum_{n=1}^{\infty}\frac{a_{\Delta}(n)}{n^{s}}\sqrt{n\tau_{2}}K_{s-\frac{1}{2}}(2\pi
n\tau_{2})\big{(}e^{2\pi in\tau_{1}}+e^{-2\pi in\tau_{1}}\big{)}\,.$ (4.6)

This function is not modular invariant but it does provide for a homogeneous
solution to the Laplace equation (1.10). We refer to the combination
$\lambda_{\Delta}(s,s_{1},s_{2})H_{\Delta}(\tau)$ as the cusp form $\Delta$
contribution to ${\cal E}(s,s_{1},s_{2};\tau)$, see appendix B for more
details.

For the case $s_{1},s_{2}\in{\mathbb{N}}$ relevant for modular graph
functions, the coefficient $\lambda_{\Delta}(s,s_{1},s_{2})$ was determined in
[48, 38], while for $s_{1},s_{2}\in\mathbb{N}+\frac{1}{2}$ we use the results
of [39] and find that $\lambda_{\Delta}(s,s_{1},s_{2})$ takes the schematic
form

$\lambda_{\Delta}(s,s_{1},s_{2})=q(s,s_{1},s_{2})(-1)^{s_{1}}\frac{\Lambda(\Delta;s+1-s_{1}-s_{2})\Lambda(\Delta;s+s_{1}-s_{2})}{\pi\langle\Delta,\Delta\rangle}\,,$
(4.7)

where $q(s,s_{1},s_{2})$ is a rational number depending on $s,s_{1},s_{2}$,
and symmetric under the exchange $s_{1}\leftrightarrow s_{2}$. We denote by
$\Lambda(\Delta;t)$ the completed ${\rm L}$-value of $\Delta$,

$\Lambda(\Delta;t)\equiv(2\pi)^{-t}\Gamma(t){\rm L}(\Delta;t)\,,$ (4.8)

which is defined for $t\in\mathbb{C}$ via analytic continuation of the
Dirichlet series

${\rm L}(\Delta;t)\equiv\sum_{n=1}^{\infty}\frac{a_{\Delta}(n)}{n^{t}}\,,$
(4.9)

and the denominator in (4.7), $\langle\Delta,\Delta\rangle$, is the Petersson
norm of $\Delta\in{\cal S}_{2s}$.111111The precise definition of the Petersson
norm is not essential for the following discussion. Although not manifest from
(4.7), for $\Delta\in{\cal S}_{2s}$ the coefficients $\lambda_{\Delta}$ are
symmetric under the exchange $s_{1}\leftrightarrow s_{2}$ as a consequence of
the functional equation

$\Lambda(\Delta,t)=(-1)^{s}\Lambda(\Delta,2s-t)\,.$ (4.10)

Therefore, whenever the eigenvalue $s$ is such that $\dim\mathcal{S}_{2s}\neq
0$, the particular solution to the inhomogeneous Laplace equation (1.10), as
constructed from the particular solution ${\cal E}^{(n,0)}_{{\rm
p}}(s,s_{1},s_{2};\tau_{2})$ [33, 47], in general does not lead to a modular
invariant solution. Only by adding the second term in (4.4), do we then arrive
at the modular invariant solution to (1.10) that we denote by ${\cal
E}(s,s_{1},s_{2};\tau)$.

The dimension of the vector space of holomorphic cusp forms $\mathcal{S}_{2s}$
is

$\dim\mathcal{S}_{2s}=\left\\{\begin{array}[]{lc}\left\lfloor\frac{2s}{12}\right\rfloor-1&2s\equiv
2\,{\rm mod}\,12\,.\\\\[5.69054pt]
\left\lfloor\frac{2s}{12}\right\rfloor\phantom{-1}&\mbox{otherwise}\end{array}\right.$
(4.11)

Since $\dim\mathcal{S}_{12}=1$, we see from our ansatz (4.1) that cusp forms
become relevant starting at order $O(N^{-3})$ which is the first instance
where the eigenvalue $s=6$ appears. In this case, the unique Hecke normalised
cusp form of weight $12$ is given by the Ramanujan cusp form
$\Delta_{12}=\sum_{n=1}^{\infty}\tau(n)q^{n}$ where $\tau(n)$ (not to be
confused with the coupling $\tau$) denotes the Ramanujan tau function.

Let us consider the $s=6$ terms at order $O(N^{-3})$ as an example, which were
denoted in (1.9) as

$\begin{split}{\cal
H}_{N}^{i}(\tau)\big{|}^{s=6}_{N^{-3}}=\frac{1}{N^{3}}\Big{[}\alpha_{6}\,{\cal
E}(6,{\scriptstyle{3\over 2}},{\scriptstyle{3\over 2}};\tau)+\beta_{6}\,{\cal
E}(6,{\scriptstyle{5\over 2}},{\scriptstyle{5\over 2}};\tau)+\gamma_{6}\,{\cal
E}(6,{\scriptstyle{7\over 2}},{\scriptstyle{3\over
2}};\tau)]\Big{]}\,.\end{split}$ (4.12)

From [39], we find that the contributions from the Ramanujan cusp form to the
relevant generalised Eisenstein series are

$\displaystyle{\cal
E}(6,{\scriptstyle\frac{3}{2}},{\scriptstyle\frac{3}{2}};\tau)\big{|}_{\rm
cusp}$
$\displaystyle=\frac{2\Lambda(\Delta_{12};6)^{2}}{\pi\langle\Delta_{12},\Delta_{12}\rangle}\times
H_{\Delta_{12}}(\tau)\,,$ $\displaystyle{\cal
E}(6,{\scriptstyle\frac{5}{2}},{\scriptstyle\frac{5}{2}};\tau)\big{|}_{\rm
cusp}$
$\displaystyle=-\frac{128\Lambda(\Delta_{12};6)^{2}}{225\pi\langle\Delta_{12},\Delta_{12}\rangle}\times
H_{\Delta_{12}}(\tau)\,,$ (4.13) $\displaystyle{\cal
E}(6,{\scriptstyle\frac{3}{2}},{\scriptstyle\frac{7}{2}};\tau)\big{|}_{\rm
cusp}$
$\displaystyle=\frac{32\Lambda(\Delta_{12};6)^{2}}{105\pi\langle\Delta_{12},\Delta_{12}\rangle}\times
H_{\Delta_{12}}(\tau)\,,$

where we have used a classic result due to Manin [49], which implies that

$\Lambda(\Delta_{12};6):\Lambda(\Delta_{12};8):\Lambda(\Delta_{12};10)=1:\frac{5}{4}:\frac{12}{5}\,,$
(4.14)

so that all the numerators are expressed in terms of $\Lambda(\Delta_{12};6)$,
The coefficients in (4.12) were determined in [33], to be

$\quad\alpha_{6}=\frac{135}{52}\,,\qquad\beta_{6}=\frac{30375}{832}\,,\qquad\gamma_{6}=\frac{42525}{832}\,.$
(4.15)

With these particular coefficients we find that the sum of the cuspidal
contributions (4.13) to (4.12) is zero.

The cancellation of the cusp forms contributions is indeed expected, since
they do not appear in the localisation computation. So for each $s$, they must
cancel between the different ${\cal E}(s,s_{1},s_{2};\tau)$ that appear in
${\cal H}_{N}^{i}(\tau)$, as we demonstrated explicitly for the above example.
This property leads to the aforementioned constraints on the coefficients
$\alpha_{j,\ell,r,s}$. We have verified and exploited similar cancellations
for all the other generalised Eisenstein series that appear in the large-$N$
expansion of ${\cal H}_{N}^{i}(\tau)$ up to order $O(N^{-7})$ and explicit
results can be found in the attached Mathematica file.

###  4.2 Lattice sum representation

The linear relations between the coefficients $\alpha_{j,\ell,r,s}$ in (4.1)
that lead to the cancellation of cusp forms suggests that there is a more
refined basis of modular functions than the ${\cal E}(s,s_{1},s_{2};\tau)$.
Following [37], which we summarise in appendix C, we introduce a new class of
modular invariant objects defined via lattice sums over four integers

$\displaystyle{\cal E}^{w}_{i,j}(\tau)\equiv$
$\displaystyle\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\neq
0}}\int_{0}^{\infty}d^{3}t\,B_{i,j}^{w}(t)\,\exp\Big{(}-\frac{\pi}{\tau_{2}}\sum_{i=1}^{3}t_{i}|p_{i}|^{2}\Big{)}\,,$
(4.16)

where $p_{i}=m_{i}+n_{i}\tau$ with $m_{i},n_{i}\in{\mathbb{Z}}$. The function
$B_{i,j}^{w}(t)$ is a symmetric function of $(t_{1},t_{2},t_{3})$ that can be
expressed in the form

$B_{i,j}^{w}(t)\equiv Y(t)^{\frac{w-3}{2}}A_{i,j}(\rho(t))\,,$ (4.17)

where $Y(t)=t_{1}t_{2}+t_{1}t_{3}+t_{2}t_{3}$ and

$\displaystyle\rho_{1}(t)=\frac{t_{1}}{t_{1}+t_{2}}\,,\qquad\rho_{2}(t)=\frac{\sqrt{Y(t)}}{t_{1}+t_{2}}\,,\qquad{\rm
with}\qquad\rho(t)=\rho_{1}(t)+i\rho_{2}(t)\,.$ (4.18)

The domain $t_{i}\geq 0$ translates into the domain $0\leq Y(t)\leq\infty$ and

$\displaystyle 0\leq\rho_{1}\leq 1\,,\ \
\qquad(\rho_{1}-{\scriptstyle\frac{1}{2}})^{2}+\rho_{2}^{2}\geq\frac{1}{4}\,,$
(4.19)

which is the fundamental domain of $\Gamma_{0}(2)$ illustrated in figure 1.

$\rho_{1}$$\rho_{2}$$0$$\frac{1}{2}$$1$$-\frac{1}{2}$$-1$$\frac{1}{2}$$1$$(123)$$(213)$$(132)$$(231)$$(312)$$(321)$$t_{1}=0$$t_{1}=t_{2}$$t_{2}=0$$t_{2}=t_{3}$$t_{1}=t_{3}$$t_{3}=0$
Figure 1: In the $\rho$-plane, the domain $t_{i}>0$ is mapped to the region
shaded in grey, isomorphic to a fundamental domain of $\Gamma_{0}(2)$. Here
the labels $(ijk)$ denote the ordering $t_{i}<t_{j}<t_{k}$ of the Schwinger
parameters.

The function $A_{1,0}(\rho)$ originally arose in the construction of the
$d^{6}R^{4}$ coefficient ${\cal E}(4,\frac{3}{2},\frac{3}{2};\tau)$ in [35],
and a more general discussion in [37] presented the general procedure for
constructing linear combinations of ${\cal E}(s,s_{1},s_{2};\tau)$. The
expressions for $A_{i,j}(\rho)$ are Laurent polynomials in $\rho_{2}$ with
coefficients which are polynomial in $\rho_{1}$. They satisfy homogeneous
Laplace equation with respect to the Laplace-Beltrami operator
$\Delta_{\rho}=\rho_{2}^{2}(\partial_{\rho_{1}}^{2}+\partial_{\rho_{2}}^{2})$

$\left[\Delta_{\rho}-s(s-1)\right]\,A_{i,j}(\rho)=0\,,$ (4.20)

inside the domain121212The upper half $\rho$-plane is denoted by
$\mathcal{H}=\\{\rho\in\mathbb{C}\,|\,\rho_{2}>0\\}$, while the congruence
subgroup $\Gamma_{0}(2)$ is defined as
$\Gamma_{0}(2)=\\{\left(\begin{smallmatrix}a&b\\\
c&d\end{smallmatrix}\right)\in SL(2,\mathbb{Z})\,|\,c\equiv
0\,\rm{mod}\,2\\}$. $\rho\in\mathcal{H}\backslash\Gamma_{0}(2)$ and where
$s=3i+j+1$ (for example, $s=4$ for $A_{1,0}$).

In appendix C.2 we present the systematic construction of $A_{i,j}$ given by
Don Zagier in unpublished notes that are expanded on in section 5.2 of [44],
where they are called ‘modular local polynomials’, see also [50].

When translated into $t$-variables (4.20) implies that the function
$B_{i,j}^{w}(t)$ defined in (4.17) satisfies the homogeneous Laplace equation

$\left[\Delta_{t}-s(s-1)\right]\,B^{w}_{i,j}(t)=0\,,$ (4.21)

where again $s=3i+j+1$ and the laplacian $\Delta_{t}$ is defined by

$\displaystyle\Delta_{t}[F(t)]=-2\sum_{i=1}^{3}\partial_{i}[t_{i}\,F(t)\,]+\sum_{i,j=1}^{3}\partial_{i}\partial_{j}\\{[t_{i}t_{j}+(2\delta_{ij}-1)Y(t)]\,F(t)\\}\,.$
(4.22)

Furthermore it is straightforward to show that

$\Delta_{\tau}\exp\Big{(}-\frac{\pi}{\tau_{2}}\sum_{i=1}^{3}t_{i}|p_{i}|^{2}\Big{)}=\Delta_{t}\exp\Big{(}-\frac{\pi}{\tau_{2}}\sum_{i=1}^{3}t_{i}|p_{i}|^{2}\Big{)}\,,$
(4.23)

so by acting with $\Delta_{\tau}$ on (4.16) and combining (4.23) and (4.21) we
deduce

$\big{[}\Delta_{\tau}-s(s-1)\big{]}{\cal E}^{w}_{i,j}(\tau)={\mbox{boundary
terms}}\,,$ (4.24)

where the boundary terms arise from integrating by parts the Laplacian
$\Delta_{t}$. In appendix D, we show that the boundary terms produce the
“source terms” in the inhomogeneous Laplace eigenvalue equation,

$\displaystyle[\Delta_{\tau}-s(s-1)]{\cal
E}_{i,j}^{w}(\tau)=\sum_{\delta=\frac{2-s}{2}}^{\frac{s-2}{2}}c_{i,j}(w,\delta)E\left(\frac{w}{2}+\delta;\tau\right)E\left(\frac{w}{2}-\delta;\tau\right)+d_{i,j}(w)E(w;\tau)\,,$
(4.25)

for particular values of the coefficients $c_{i,j}(w,\delta),\,d_{i,j}(w)$.

We note a few important facts. Firstly, when the weight $w$ and the eigenvalue
$s=3i+j+1$ have opposite parity (as in our case), the source terms are
bilinears in Eisenstein series with half-integer indices. Secondly, the
modular invariant function ${\cal E}_{i,j}^{w}(\tau)$ satisfies an
inhomogeneous Laplace equation where the source terms have total
“trascendental weight” given by $w=s_{1}+s_{2}$. It is always possible to use
the functional equation

$\Gamma(s)E(s;\tau)=\Gamma(1-s)E(1-s;\tau)\,,$ (4.26)

to rewrite (4.25) so that all Eisenstein series in the source term appear with
positive indices, with the drawback of spoiling uniform trascendentality in
weight.

Lastly, since we are interested in eigenvalues $s\geq w+1$ it is always
possible to use (1.4) and invert the laplacian131313We can check from the
lattice sum representation (4.16) that for the range of eigenvalues of present
interest, the coefficient of the homogeneous solution $E(s;\tau)$ vanishes.
This implies that there is no issue in inverting the Laplace operator in
(4.25). over the single Eisenstein series that appear as source terms to get

${\cal
E}_{i,j}^{w}(\tau)=\sum_{\delta=\frac{2-s}{2}}^{\frac{s-2}{2}}c_{i,j}(w,\delta){\cal
E}\\!\left(s,\frac{w}{2}+\delta,\frac{w}{2}-\delta;\tau\right)+\frac{d_{i,j}(w)}{w(w-1)-s(s-1)}E(w;\tau)\,.$
(4.27)

Note that for eigenvalue $s$ such that ${\rm dim}\,{\cal S}_{2s}=0$ it is
possible to find an isomorphism between the two vector spaces spanned by
either ${\cal E}_{i,j}^{w}(\tau)$ or by ${\cal E}(s,s_{1},s_{2};\tau)$ with
$s=3i+j+1$ and fixed weight $w=s_{1}+s_{2}$ (modulo the addition of single
Eisenstein series). For example, in the first few cases relevant to (1.9) up
to order $O(N^{-2})$ we have

$\displaystyle{\cal E}_{1,0}^{3}(\tau)$ $\displaystyle=-\frac{9}{10}\pi{\cal
E}(4,{\scriptstyle\frac{3}{2}},{\scriptstyle\frac{3}{2}};\tau)-\frac{3}{10}\pi
E(3;\tau)\,,$ $\displaystyle{\cal E}_{1,1}^{4}(\tau)$
$\displaystyle=-\frac{81}{28}\pi{\cal
E}(5,{\scriptstyle\frac{3}{2}},{\scriptstyle\frac{5}{2}};\tau)-\frac{99}{140}\pi
E(4;\tau)\,,$ (4.28) $\displaystyle{\cal E}_{2,0}^{4}(\tau)$
$\displaystyle=-\frac{15}{11}\pi{\cal
E}(7,{\scriptstyle\frac{3}{2}},{\scriptstyle\frac{5}{2}};\tau)-\frac{1}{14}\pi
E(4;\tau)\,.$

We can use these relations and similar ones at higher order in $1/N$ to
express ${\cal H}^{i}_{N}(\tau)$ in terms of ${\cal E}^{w}_{i,j}(\tau)$ (and
$E(s;\tau)$). For example, up to $1/N^{3}$, we have

$\displaystyle{\cal H}^{i}_{N}(\tau)\big{|}_{N^{-1}}$
$\displaystyle=-\frac{15}{4\pi}{\cal
E}^{3}_{1,0}(\tau)-\frac{9}{8}E(3;\tau)\,,$ (4.29) $\displaystyle{\cal
H}^{i}_{N}(\tau)\big{|}_{N^{-2}}$ $\displaystyle=C_{1}+\frac{945}{64\pi}{\cal
E}_{2,0}^{4}(\tau)-\frac{105}{22\pi}{\cal
E}^{4}_{1,1}(\tau)-\frac{297}{128}E(4;\tau)\,,$ (4.30) $\displaystyle{\cal
H}^{i}_{N}(\tau)\big{|}_{N^{-3}}$ $\displaystyle=\frac{645}{572\pi}{\cal
E}^{3}_{1,0}(\tau)-\Big{[}\frac{945}{104\pi}{\cal
E}^{5}_{1,2}(\tau)+\frac{945}{208\pi}{\cal
E}^{3}_{1,2}(\tau)-\frac{63}{104\pi}{\cal E}^{3}_{0,5}(\tau)\Big{]}$ (4.31)
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!+\Big{[}\frac{31185}{544\pi}{\cal
E}^{5}_{2,1}(\tau)+\frac{23625}{544\pi}{\cal
E}^{3}_{2,1}(\tau)-\frac{23625}{4352\pi}{\cal E}^{3}_{1,4}(\tau)\Big{]}$
$\displaystyle\\!\\!\\!\\!\\!\\!\\!\\!-\Big{[}\frac{45045}{512\pi}{\cal
E}^{5}_{3,0}(\tau)+\frac{183645}{2048\pi}{\cal
E}^{3}_{3,0}(\tau)-\frac{10395}{1024\pi}{\cal
E}^{3}_{2,3}(\tau)\Big{]}-\frac{496125}{146432}E(3;\tau)+\frac{5805}{512}E(5;\tau)\,.$

We thus see that ${\cal H}_{N}^{i}(\tau)$ can be written in terms of the more
refined basis of modular functions ${\cal E}^{w}_{i,j}(\tau)$ and $E(w;\tau)$,
which manifestly does not include the unwanted cusp form contributions that
appear in generic sums of ${\cal E}(s,s_{1},s_{2};\tau)$. Furthermore, they
are the objects that are naturally written as lattice sums. This pattern
generalises to all the higher orders, which we have verified explicitly up to
order $O(N^{-7})$. The explicit expression for the lattice sum representation
is described in some detail in appendix C.

##  5 Numerical estimate of free energy at finite $N$

We will now discuss how the large-$N$ finite-$\tau$ results for ${\cal C}_{N}$
and ${\cal H}_{N}$ can be used to accurately estimate these quantities at
finite $N$ and $\tau$. This is especially important for ${\cal H}_{N}$, for
which we have no exact expression. These finite-$N$ numerical results will
provide an important input for the numerical bootstrap study of
$\mathcal{N}=4$ SYM following [45].

Our strategy is to divide ${\cal C}_{N}$ and ${\cal H}_{N}$ into perturbative
and non-perturbative terms as $g_{{}_{YM}}^{2}\to 0$, or equivalently for
$\tau_{2}\gg 1$, arising in the matrix model (2.1) and (2.2)

$\begin{split}{\cal C}_{N}={\cal C}_{N}\big{|}_{pert}+{\cal
C}_{N}\big{|}_{non-pert}\,,\qquad\qquad\qquad{\cal H}_{N}={\cal
H}_{N}\big{|}_{pert}+{\cal H}_{N}\big{|}_{non-pert}\,,\end{split}$ (5.1)

which is a meaningful distinction for all $N$. The perturbative sector only
depend on $\tau_{2}$, and was already computed for finite $N$ and $\tau_{2}$
in [13, 8] using the method of orthogonal polynomials. For instance, for
${\cal C}_{N}$ we have [13]141414A slightly simplified version of this formula
can be found in (A.12) in [17].

$\begin{split}{\cal
C}_{N}\big{|}_{pert}=-\tau_{2}^{2}\partial_{\tau_{2}}^{2}\int_{0}^{\infty}dw\frac{e^{-\frac{w^{2}}{\pi\tau_{2}}}}{2\sinh^{2}w}\Big{[}L_{N-1}^{(1)}({\scriptstyle\frac{w^{2}}{\pi\tau_{2}}})-\sum_{i,j=1}^{N}(-1)^{i-j}L_{i-1}^{(j-i)}({\scriptstyle\frac{w^{2}}{\pi\tau_{2}}})L_{j-1}^{(i-j)}({\scriptstyle\frac{w^{2}}{\pi\tau_{2}}})\Big{]}\,,\end{split}$
(5.2)

where $L_{i}^{(j)}(x)$ are generalised Laguerre polynomials, and note that for
all $N$ we need only perform a single integral and some finite sums. The
analogous expression for ${\cal H}_{N}|_{pert}$ from [8] involves two
integrals and some finite sums for all $N$. Since it is much more complicated,
we do not show it here explicitly, but we include it in the attached
Mathematica file.

For the non-perturbative sector, we will use the large-$N$ finite-$\tau$
expressions from the previous sections, which were written in terms of
Eisenstein series and generalised Eisenstein series. The advantage of only
using the non-perturbative terms of these functions is that we avoid the
divergence at the free theory point $\tau_{2}\to\infty$ that occurs for their
perturbative terms, as seen from the Fourier expansions (3.1) and (B). From
the Fourier decomposition of the non-holomorphic Eisenstein series (3.1) we
see that non-perturbative terms correspond to Fourier modes with $k\neq 0$
while from (4.3) we see that the non-perturbative terms of the generalised
Eisenstein series correspond to terms with $(n,m)\neq(0,0)$. The non-
perturbative sectors converge extremely quickly for $\tau$ in the fundamental
domain

$\begin{split}|\tau|\geq 1\,,\qquad|{\rm
Re\,}(\tau)|\leq\frac{1}{2}\,,\end{split}$ (5.3)

so that only a few values of $k\neq 0$ and $(n,m)\neq(0,0)$ need be included
from the Fourier expansions to achieve good accuracy. Since the large-$N$
expansion is an asymptotic expansion, we expect that the expansion will
actually get worse at some order. We observe that the best approximation is
given by keeping terms up to $O(1/N^{\frac{5}{2}})$ in the non-perturbative
sector for ${\cal C}_{N}$,151515For ${\cal C}_{N}(\tau)$ we may use the exact
result in [16]. However, as we will show, the simple approximation given here
already provides an extremely good estimate. And in this way, we treat ${\cal
C}_{N}(\tau)$ and ${\cal H}_{N}(\tau)$ uniformly.

$\begin{split}{\cal C}_{N}\big{|}_{non-
pert}&\approx\Bigg{[}-\frac{3\sqrt{N}}{2^{4}}E({\scriptstyle{3\over
2}};\tau)+\frac{45}{2^{8}\sqrt{N}}E({\scriptstyle{5\over
2}};\tau)+\frac{1}{{N}^{\frac{3}{2}}}\big{[}-\frac{39}{2^{13}}E({\scriptstyle{3\over
2}};\tau)+\frac{4725}{2^{15}}E({\scriptstyle{7\over 2}};\tau)\big{]}\\\
&\qquad+\frac{1}{{N}^{\frac{5}{2}}}\big{[}-\frac{1125}{2^{16}}E({\scriptstyle{5\over
2}};\tau)+\frac{99225}{2^{18}}E({\scriptstyle{9\over
2}};\tau)\big{]}\Bigg{]}_{k\neq 0}\,,\\\ \end{split}$ (5.4)

and for ${\cal H}_{N}$:

$\begin{split}&{\cal H}_{N}\big{|}_{non-
pert}\approx\Bigg{[}6\sqrt{N}E({\scriptstyle{3\over
2}};\tau)-\frac{9}{2\sqrt{N}}E({\scriptstyle{5\over
2}};\tau)+\frac{27}{2^{3}N}{\cal E}(4,{\scriptstyle{3\over
2}},{\scriptstyle{3\over 2}};\tau)\\\
+&\,\frac{1}{{N}^{\frac{3}{2}}}\big{[}\frac{117}{2^{8}}E({\scriptstyle{3\over
2}};\tau)-\frac{3375}{2^{10}}E({\scriptstyle{7\over
2}};\tau)\big{]}+\frac{1}{N^{2}}\big{[}-\frac{14175}{704}{\cal
E}(7,{\scriptstyle{5\over 2}},{\scriptstyle{3\over
2}};\tau)+\frac{1215}{88}{\cal E}(5,{\scriptstyle{5\over
2}},{\scriptstyle{3\over 2}};\tau)\big{]}\\\
+&\,\frac{1}{{N}^{\frac{5}{2}}}\big{[}\frac{675}{2^{10}}E({\scriptstyle{5\over
2}};\tau)-\frac{33075}{2^{12}}E({\scriptstyle{9\over
2}};\tau)\big{]}\Bigg{]}_{\begin{subarray}{c}\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!\\!k\neq
0\\\ (n,m)\neq(0,0)\end{subarray}}\,,\end{split}$ (5.5)

where on the RHS we only keep the $k\neq 0$ instanton terms in the Fourier
expansions for $E(s;\tau)$ in (3.1) and the instanton/anti-instanton sectors
$(n,m)\neq(0,0)$ for ${\cal E}(s,s_{1},s_{2};\tau)$ given in(4.3). The
explicit expressions are given in the attached Mathematica file.

Figure 2: Comparison of large-$N$ approximation to exact $N$ for ${\cal
C}_{N}$ and ${\cal H}_{N}$ for two paths in the fundamental domain of $\tau$:
from the free theory $\tau\to i\infty$ to the $\mathbb{Z}_{2}$ self-dual point
$\tau=i$ along the line $\tau_{1}=0$, and from $\tau=i$ to the
$\mathbb{Z}_{3}$ self-dual point $\tau=e^{i\pi/3}$ along the arc $|\tau|=1$.

We can check our approximations for ${\cal C}_{N}$ and ${\cal H}_{N}$ by
comparing to low values of $N$, where the exact expression can be computed
from the matrix model integral (2.1) by doing $N-1$ integrals and truncating
the instanton expansion in (2.2), which converges very quickly in the
fundamental domain. This finite-$N$ calculation was performed for $N=2,3$ in
[45]. In Figure 2, we compare these finite-$N$ expressions to our large-$N$
approximation for two paths in a fundamental domain of $\tau$: from the free
theory $\tau\to i\infty$ to the $\mathbb{Z}_{2}$ self-dual point $\tau=i$ with
$\tau_{1}=0$, and from $\tau=i$ to the $\mathbb{Z}_{3}$ self-dual point
$\tau=e^{i\pi/3}$ along the arc $|\tau|=1$. Along the first path, instantons
are very small, so is perhaps not so surprising that our large-$N$
approximation is so accurate, as we have included the exact perturbative
sector. The precise match is more surprising along the second path, where the
range of the plots is much smaller, and the variation is mostly due to
instantons. The greatest discrepancy between large and finite $N$ is at
$\tau=e^{i\pi/3}$, but even that discrepancy is only around $.01\%$ relative
error, and this error decreases as we increase from $N=2$ to $N=3$. We are
thus confident that our approximations are very accurate for all $N$ and
$\tau$.

##  6 Constraints on stress tensor correlator

We now show how the stress tensor correlator can be constrained at large $N$
using the localisation quantities computed in the previous sections. We
consider the four point function of the bottom component $S$ of the stress
tensor multiplet, often also denoted by $\mathcal{O}_{2}$. This operator is a
dimension 2 scalar in the ${\bf 20}^{\prime}$ of $SO(6)_{R}$, and can thus be
represented as a rank-two traceless symmetric tensor $S_{IJ}(\vec{x})$, with
indices $I,J=1,\ldots,6$. For simplicity we will contract these indices with
polarisation vectors $Y^{I}$, with $Y\cdot Y=0$. The four-point function is
then fixed by superconformal symmetry to take the form [51, 52]

$\begin{split}\langle
S(\vec{x}_{1},Y_{1})S(\vec{x}_{2},Y_{2})S(\vec{x}_{3},Y_{3})S(\vec{x}_{4},Y_{4})\rangle=\frac{1}{\vec{x}_{12}^{4}\vec{x}_{34}^{4}}\left[\vec{{\cal
S}}_{\text{free}}+{\cal T}(U,V)\vec{\Theta}\right]\cdot\vec{{\cal
B}}\,,\end{split}$ (6.1)

where $\vec{x}_{ij}\equiv\vec{x}_{i}-\vec{x}_{j}$, and
$U\equiv\frac{\vec{x}_{12}^{2}\vec{x}_{34}^{2}}{\vec{x}_{13}^{2}\vec{x}_{24}^{2}}$
and
$V\equiv\frac{\vec{x}_{14}^{2}\vec{x}_{23}^{2}}{\vec{x}_{13}^{2}\vec{x}_{24}^{2}}$
are the usual conformal invariant cross-ratios. The precise forms of
$\vec{{\cal S}}_{\text{free}}$, $\vec{\Theta}$ and $\vec{\cal B}$ can be found
in [7], and all non-trivial information is given by the $R$-symmetry invariant
correlator $\mathcal{T}(U,V)$, which will be our focus.

We can expand $\mathcal{T}(U,V)$ at large $c=(N^{2}-1)/4$, whose functional
form is fixed by the analytic bootstrap [53, 54, 33] to take the form

$\begin{split}{\cal T}=\frac{8}{c}{\cal T}^{R}+\frac{1}{c^{2}}\big{[}{\cal
T}^{R|R}+B_{0}{\cal T}^{0}\big{]}+\frac{1}{c^{3}}\big{[}{\cal
T}^{R|R|R}+\sum_{i=0,2,3,4}\overline{B_{i}}{\cal
T}^{i}\big{]}+\dots\,,\end{split}$ (6.2)

where the ellipses denote both higher loop terms as well as contact terms
involving higher derivative terms such as $R^{4}$, which at large $c$ and
finite $\tau$ generically come with fractional powers of $c$.161616 The first
few higher derivative corrections to the correlator, including
$R^{4},d^{4}R^{4}$ and $d^{6}R^{4}$, were determined in [14, 33]. Furthermore,
starting from order $O(1/c^{3})$ there is actually no difference between the
loop and higher derivative expansion in the large-$c$ finite-$\tau$ limit,
since both $d^{12}R^{4}$ and ${\cal T}^{R|R|R}$ are the same order in
derivatives. In this case, we simply define the terms written in (6.2) to be
those that are independent of $\tau$, as higher derivative terms generically
depend on $\tau$. We denote the supergravity exchange term by ${\cal T}^{R}$,
which is given by [55]

$\begin{split}\mathcal{T}^{R}=-\frac{1}{8}U^{2}\bar{D}_{2,4,2,2}(U,V)\,,\end{split}$
(6.3)

where $\bar{D}_{a,b,c,d}(U,V)$ are standard functions whose explicit form can
be found e.g. in [56]. The coefficient of ${\cal T}^{R}$ is fixed by the
conformal Ward identity [57]. Next we have the one-loop graviton exchange term
${\cal T}^{R|R}$, which was fixed in position space in [41] from tree-level
data using the AdS unitarity method of [58] up to a contact-term ambiguity

$\begin{split}\mathcal{T}^{0}=U^{2}\bar{D}_{4,4,4,4}\,.\end{split}$ (6.4)

In general, ${\cal T}^{n}$ for integer $n$ denotes contact terms with $2n+8$
derivatives, which are related to the divergences that appear at each loop
level. The two-loop graviton exchange term is denoted by ${\cal T}^{R|R|R}$,
and was conjecturally fixed in [43, 42] up to four contact-term ambiguities,
whose explicit form we will not use. The coefficients of contact-term
ambiguities (i.e. $B_{0}$ and $\overline{B_{i}}$) are not fixed by general
analytic bootstrap constraints, and so require the supersymmetric localisation
constraints discussed here.

We can constrain ${\cal T}(U,V)$ using an exact property of ${\cal
H}_{N}(\tau)$ [8]:

$\begin{split}{\cal H}_{N}(\tau)&=48\zeta(3)c+c^{2}I_{4}[{\cal T}]\,,\\\
I_{4}[{\cal
T}]&\equiv\frac{32}{\pi}\int_{0}^{\infty}\\!dr\\!\int_{0}^{2\pi}d\theta\,r^{3}\sin^{2}\theta\,\frac{1+U+V}{U^{2}}\bar{D}_{1,1,1,1}(U,V){\cal
T}(U,V)\bigg{|}_{\begin{subarray}{c}U=1+r^{2}-2r\cos\theta\\\
V=r^{2}\end{subarray}}\,,\end{split}$ (6.5)

where the $48\zeta(3)c$ term comes from the free theory contribution to (6.1).
A similar integrated constraint between ${\cal T}(U,V)$ and ${\cal
C}_{N}(\tau)$ was derived in [7] and has been studied in the literature so we
will not use its detailed form here. Since we have two integrated constraints,
one from ${\cal C}_{N}$ and one from ${\cal H}_{N}$, that means we can fix at
most two $B$’s at each order in $1/c$ in (6.2). There is just a single $B_{0}$
at $O(1/c^{2})$, which was fixed in [13] using the ${\cal C}_{N}$ integrated
constraint to be

$\begin{split}B_{0}=\frac{15}{4}\,.\end{split}$ (6.6)

For the integrated constraint arising from ${\cal H}_{N}$, we make use of the
following integrals:

$\begin{split}I_{4}[{\cal T}^{R}]&=3-6\zeta(3)\,,\quad I_{4}[{\cal
T}^{0}]=\frac{16}{5}\,,\qquad I_{4}[{\cal
T}^{R|R}]=-\frac{48\zeta(3)}{5}-\frac{83}{10}\,,\\\ \end{split}$ (6.7)

where the first two integrals were computed numerically in [8] using the
explicit expressions in (6.3) and (6.4), while the last integral we compute
now using the explicit expression in [41]. We find that these integrals
applied to the correlator (6.2) using the integrated constraint (6.5) exactly
match ${\cal H}_{N}$ in (1.7) to one-loop order with the one-loop ambiguity
$B_{0}$ fixed as in (6.6) and $C_{0}$ fixed in (4.2). We do not have
sufficient constraints yet to fix the ambiguities at higher loop orders, such
as the four ambiguities at two loops [42, 43], but the localisation data for
${\cal H}_{N}$ is now available for that purpose.

##  7 Discussion

There are four main results of this paper. Firstly, we found recursion
formulae that relate the half-integer terms and part of the integer terms in
the large-$N$ expansion of ${\cal H}_{N}(\tau)$ to ${\cal C}_{N}(\tau)$.
Secondly, we gave striking evidence that the large-$N$ expansion of ${\cal
H}_{N}(\tau)$ can be written as a lattice sum to any order in $1/N$. Thirdly,
we used the large-$N$ expansion of ${\cal H}_{N}(\tau)$ to verify the one-loop
contact term in the stress tensor correlator as originally fixed using ${\cal
C}_{N}(\tau)$ in [13], and gave new constraints at higher loops. Lastly, we
showed how the large-$N$ expansions of ${\cal H}_{N}(\tau)$ and ${\cal
C}_{N}(\tau)$ can be used to accurately estimate these quantities for all
finite $N$ and $\tau$. This estimate will be useful for the numerical
bootstrap at finite $N$ and $\tau$ [45], which was previously limited to low
$N$ due to the difficulty of performing the $(N-1)$-dimensional matrix model
integral.

While the recursion relations we found were sufficient to fix the half-integer
powers of $N$ in ${\cal H}_{N}(\tau)$ to all orders, we have not found a
recursion relation that applies to the full ${\cal H}_{N}(\tau)$, which could
be used to compute ${\cal H}_{N}(\tau)$ at finite $N$ and $\tau$ as was done
for ${\cal C}_{N}(\tau)$ in [16]. Since ${\cal H}_{N}(\tau)$ can be written as
a lattice sum to all orders in $1/N$, it seems likely that a finite $N$
lattice sum expression should exist, analogous to the lattice sum expression
for ${\cal C}_{N}(\tau)$ in (1.5). Importantly that expression is finite even
though it contains the term in the lattice sum with $m=n=0$. Based on the
lattice sum for the large-$N$ expansion of ${\cal H}_{N}(\tau)$ in (4.16), a
natural conjecture for the finite-$N$ expression would be

$\displaystyle{\cal
H}_{N}(\tau)=\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\in{\mathbb{Z}}+\tau{\mathbb{Z}}}}\int_{0}^{\infty}d^{3}t\,{\cal
B}_{N}(t_{1},t_{2},t_{3})\,\exp\Big{(}-\frac{\pi}{\tau_{2}}\sum_{i=1}^{3}t_{i}|p_{i}|^{2}\Big{)}\,,$
(7.1)

where the function ${\cal B}_{N}(t_{1},t_{2},t_{3})$ is a symmetric function
that can be shown to satisfy an inversion relation of the form ${\cal
B}_{N}(t_{1},t_{2},t_{3})=1/Y(t)^{2}\,{\cal
B}_{N}(\hat{t}_{1},\hat{t}_{2},\hat{t}_{3})$ where $\hat{t}_{i}=t_{i}/Y(t)$
and $Y(t)=t_{1}t_{2}+t_{2}t_{3}+t_{3}t_{1}$. The unknown function ${\cal
B}_{N}(t_{1},t_{2},t_{3})$ must also be well-defined when $t_{i}=0$ and the
large-$N$ expansion of (7.1) would have to reproduce the terms with half-
integer as well as integer powers of $1/N$. Perhaps a recursion relation could
be found for ${\cal B}_{N}(t_{1},t_{2},t_{3})$ that would relate it to
$B_{N}(t)$ given in (1.5), just as we found a recursion relation connecting
parts of ${\cal H}_{N}(\tau)$ to ${\cal C}_{N}(\tau)$. This would also allow
us to compute the contribution to ${\cal H}_{N}(\tau)$ that is non-
perturbative in $N$, analogous to that found for ${\cal C}_{N}(\tau)$ in [17].

The identities we found that relate ${\cal H}_{N}(\tau)$ to ${\cal
C}_{N}(\tau)$ suggest that higher mass derivatives of the partition function
$Z_{N}$ might also be related to lower mass derivatives, or at least certain
parts of the large-$N$ expansion as found in this paper. From the localisation
expression for $Z_{N}$, it is clear that the perturbative part of the $2p$-th
mass derivative of $Z_{N}$ can be written in terms of products of $p$ zeta
functions. While the stress tensor correlator was recently shown to contain
single-valued multiple valued zeta functions [59, 60], it appears that these
vanish after taking the integral of the correlator that is related to mass
derivatives of $Z_{N}$171717Similar phenomena were observed in the weak
coupling regime [25, 22].. It would be interesting to find the modular
completions of these products of zeta functions.

The $\tau$-independent terms in the large-$N$ expansions of ${\cal
H}_{N}(\tau)$ and ${\cal C}_{N}(\tau)$ give two constraints at each orders in
the loop expansion of the holographic correlator. At one loop, there was a
single contact-term ambiguity, which was fixed using ${\cal C}_{N}(\tau)$ in
[13], and verified using ${\cal H}_{N}(\tau)$ in this paper. At two loops,
there are four such ambiguities, so we need two additional constraints. It is
possible that derivatives with respect to the squashing parameter $b$ of
$Z_{N}$ on the squashed sphere might give an additional independent
constraint. It would be interesting to find another source of constraints, so
that the 2-loop term can be fixed completely.

The methods developed in this paper to compute the perturbative contributions
to ${\cal H}_{N}(\tau)$ should also apply to other 4d $\mathcal{N}=2$ gauge
theories, since the localisation expressions for all 4d $\mathcal{N}=2$ gauge
theories include Barnes G-functions that lead to Bessel functions at large
$N$. Indeed, our methods were already used to fix the $O(N^{0})$ term in the
$USp(2N)$ gauge theory considered in [61, 31].181818Our methods were shared
with the authors of those papers prior to this publication. It would be
interesting to study the large-$N$ expansion of mass derivatives of $Z_{N}$
for other ${\cal N}=2$ gauge theories.

## Acknowledgments

DD and MBG are grateful for the hospitality of the Pollica Physics Centre and
to conversations with Don Zagier, Kim Klinger–Logan, Ksenia Fedosova and Boris
Pioline, as well as other participants in the programme ‘New Connections
Between Physics and Number Theory’. The work of LFA is supported by the
European Research Council (ERC) under the European Union’s Horizon 2020
research and innovation programme (grant agreement No 787185). LFA is also
supported in part by the STFC grant ST/T000864/1. SMC is supported by a Royal
Society University Research Fellowship, URF\R\221310. CW is supported by a
Royal Society University Research Fellowship, URF\R\221015 and a STFC
Consolidated Grant, ST\T000686\1 “Amplitudes, strings & duality”.

##  Appendix A Matrix model details

In this Appendix, we will discuss some details of the large-$N$
large-$\lambda$ expansion from Section 2. We will first discuss the pole
prescription that gives the leading $O(N^{2})$ result in (2). We will then
discuss how the $O(N^{0})$ term in that expression can be computed using a
similar but more complicated calculation due to the double sum.

###  A.1 $O(N^{2})$ at large $\lambda$

We start with the with the expression for the 4-body contribution in (2.20),
which we repeat here for ease of access:

$\begin{split}I^{(0)}(\lambda)&=\oint\frac{d{s}}{2\pi i}\frac{d{t}}{2\pi
i}\,p({s},{t})c({s},{t})\,,\qquad
c({s},{t})=\sum^{\infty}_{\ell=1}c_{\ell}({s},{t})\,,\\\
p({s},{t})&=\frac{3\lambda^{{s}+{t}+2}\sin(\pi{s})\sin(\pi{t})\Gamma(2{s}+4)^{2}\Gamma(2{t}+4)^{2}\zeta(2{s}+3)\zeta(2{t}+3)}{2^{4{s}+4{t}+3}\pi^{2({s}+{t}+2)}}\,,\\\
c_{\ell}({s},{t})&=\frac{\ell(\ell-{s}-3)(\ell-{t}-3)\Gamma(\ell-{s}-3)\Gamma(\ell-{s}-1)\Gamma(\ell-{t}-3)\Gamma(\ell-{t}-1)}{\pi^{2}\Gamma(\ell+{s}+2)\Gamma(\ell+{s}+3)\Gamma(\ell+{t}+2)\Gamma(\ell+{t}+3)}\,.\end{split}$
(A.1)

We find it convenient to use the following representation for $c({s},{t})$

$\displaystyle c({s},{t})$
$\displaystyle=\sum^{\infty}_{q=0}P^{(3q)}({s},{t})\zeta(4{s}+4{t}+15+2q)$
(A.2)
$\displaystyle=\frac{\zeta(4{s}+4{t}+15)}{\pi^{2}}+P^{(3)}({s},{t})\zeta(4{s}+4{t}+17)+\cdots\,,$

where $P^{(3q)}({s},{t})$ are polynomials of degree $3q$, easily computable by
expanding $c_{\ell}({s},{t})$ for large $\ell$ and performing the sum over
$\ell$. This representation is useful, as it makes explicit poles located at
$\mbox{Re}({s}+{t})<0$ that are easy to miss otherwise. In addition, we will
have poles of the type $\mbox{Re}({s}),\mbox{Re}({t})<0$. We now compute all
contributions to $I^{(0)}(\lambda)$ by degree of complexity.

First there is a contribution arising from the explicit poles of
$p({s},{t})c_{1}({s},{t})$ at ${s}=-1$ and ${t}=-1,-2,-5/2,-7/2,-9/2,\cdots$
(and the same with ${s}\leftrightarrow{t}$). This can be computed analytically
to all orders in $1/\lambda$ and gives

$I^{(0)}_{0}(\lambda)=6-\frac{16\pi^{2}}{\lambda}+3\sum_{n=2}^{\infty}\frac{(-1)^{n}4^{n+2}\pi^{2n-\frac{7}{2}}\lambda^{\frac{1}{2}-n}\Gamma\left(n-\frac{5}{2}\right)\Gamma\left(n-\frac{3}{2}\right)\Gamma\left(n+\frac{1}{2}\right)\zeta^{\prime}(2-2n)}{\Gamma(n-1)\Gamma(2n-2)}\,.$
(A.3)

Naively there is a contribution when ${s}=-m,{t}=-n$, with $m,n=2,3,\cdots$.
It turns out however, that this contribution is zero, as $c(-m,-n)=0$. This
can be seen from the representation for $c({s},{t})$ in terms of zeta
functions. For negative integer ${s},{t}$ the polynomials are such that the
series truncate, and we can directly evaluate $c({s},{t})$ at those points and
see that it vanishes. Then there is a contribution from
${s}=-m-1/2,{t}=-n-1/2$, with $m,n=2,3,4,\cdots$. This contribution can be
written as

$\displaystyle
I^{(0)}_{1}(\lambda)=3\sum^{\infty}_{m,n=2}\frac{(-1)^{m+n}2^{4m+4n-1}\pi^{2(m+n-1)}\zeta^{\prime}(2{-}2m)\lambda^{-m-n+1}\zeta^{\prime}(2{-}2n)}{\Gamma(2m-2)^{2}\Gamma(2n-2)^{2}}c(-m{-}\tfrac{1}{2},-n{-}\tfrac{1}{2})\,,$
(A.4)

where $c(-m{-}1/2,-n{-}1/2)$ can be easily computed for any pair of integers
$m,n$. For ${s}=-m-1/2,{t}=-n-1/2$ we find

$\displaystyle
c({s},{t})=\frac{\Gamma(-{s}-1)\Gamma(-{s})\Gamma(-{t}-1)\Gamma(-{t})}{6\pi^{2}\Gamma({s}+3)\Gamma({s}+4)\Gamma({t}+2)\Gamma({t}+3)}P_{s}({t})\,,$
(A.5)

where the polynomials $P_{s}({t})$ can be computed for ${s}=-5/2,-7/2,\cdots$

$\begin{split}P_{-5/2}({t})&=1,\\\
P_{-7/2}({t})&=\frac{1}{5}\left(-8{t}^{2}-12{t}+5\right),\\\
P_{-9/2}({t})&=\frac{1}{105}\left(-128{t}^{4}-128{t}^{3}+272{t}^{2}+152{t}+105\right),\end{split}$
(A.6)

and so on. Lastly, we have the poles at ${s}\\!=\\!-m$ and ${t}\\!=\\!-n-1/2$
(and the same with ${s},{t}$ exchanged). These are slightly subtle because the
zeta functions have poles at these locations. We find their contribution is
given

$\displaystyle I^{(0)}_{2}(\lambda)$
$\displaystyle=\sum^{\infty}_{m=2}\frac{2^{2m+1}\pi^{2m-\frac{5}{2}}\lambda^{\frac{3}{2}-m}\zeta(3-2m)\Gamma\left(m-\frac{5}{2}\right)\Gamma\left(m-\frac{3}{2}\right)}{\Gamma\left(\frac{9}{2}-m\right)\Gamma(m-2)\Gamma(2m-3)}$
$\displaystyle-3\sum^{\infty}_{m=2}\frac{2^{2m+2}\pi^{2m-\frac{5}{2}}\lambda^{\frac{1}{2}-m}\zeta(3-2m)\Gamma\left(m-\frac{3}{2}\right)\Gamma\left(m-\frac{1}{2}\right)}{\Gamma\left(\frac{7}{2}-m\right)\Gamma(m-1)\Gamma(2m-3)}$
(A.7)
$\displaystyle-3\sum^{\infty}_{m=2}\frac{64(-1)^{m}\pi^{2m-3}\lambda^{\frac{1}{2}-m}\Gamma\left(m-\frac{5}{2}\right)\Gamma\left(m-\frac{3}{2}\right)\zeta^{\prime}(2-2m)}{\Gamma(m-1)^{2}}$
$\displaystyle-3\sum^{\infty}_{m,n=2}\frac{4^{m+n+1}\pi^{2m+2n-\frac{9}{2}}\zeta(3-2m)\lambda^{-m-n+\frac{3}{2}}\Gamma\left(m+n-\frac{5}{2}\right)\Gamma\left(m+n-\frac{3}{2}\right)\zeta^{\prime}(2-2n)}{2\Gamma(2m-3)\Gamma(2n-2)\Gamma\left(-m-n+\frac{9}{2}\right)\Gamma(m+n-2)}\,.$

Now we discuss poles where ${s}+{t}=negative$, which are apparent from the
expression for $c({s},{t})$ as a sum over zeta-functions in (A.2). More
precisely

$c({s},{t})\sim\frac{R_{q}({t})}{{s}-\frac{1}{2}(-(2q+7)-2{t})}\,,\qquad{\rm
with}\quad q=0,1,2,\cdots\,,$ (A.8)

where we will first perform the integral over ${s}$. The residues at those
poles can be calculated, and are given by

$R_{q}({t})=-\frac{4^{-q-1}\Gamma\left(q-\frac{1}{2}\right)\Gamma\left(q+\frac{1}{2}\right)\Gamma\left(q+\frac{3}{2}\right)\Gamma(2(q+{t}+2))}{\pi^{7/2}\Gamma(q+1)\Gamma(2{t}+4)}\,.$
(A.9)

We can now perform the integral over ${s}$ for each $q$, and obtain

$\displaystyle I^{(0)}_{3}(\lambda)$
$\displaystyle=\sum^{\infty}_{q=0}\lambda^{-q-\frac{3}{2}}I_{3,q}^{(0)}\,,\qquad{\rm
where}\qquad I_{3,q}^{(0)}\equiv\oint\frac{d{t}}{2\pi i}Q_{q}({t})\,,$ (A.10)

with

$\displaystyle Q_{q}({t})=$ (A.11) $\displaystyle
3\frac{(-1)^{q+1}4^{q+4}\Gamma\left(q{-}\frac{1}{2}\right)\Gamma\left(q{+}\frac{1}{2}\right)\Gamma\left(q{+}\frac{3}{2}\right)\Gamma(2{t}+4)\zeta(2{t}+3)\Gamma(-2q-2{t}-3)\zeta(-2(q+{t}+2))}{\pi^{-\frac{1}{2}-2q}\Gamma(q+1)}\,.$

There is a subtle question of which poles to include in the integrals for
$I^{(0)}_{3}(\lambda)$. This is akin to the precise choice of the contour. It
turns out we need to choose all poles for negative ${t}$. Let’s consider the
leading order contribution, arising from $q=0$, as this already has the
ingredients we will encounter later on:

$\displaystyle I^{(0)}_{3,0}=3\pi^{2}\oint\frac{d{t}}{2\pi
i}2^{8}\Gamma(-2{t}-3)\Gamma(2{t}+4)\zeta(-2({t}+2))\zeta(2{t}+3)\,,$ (A.12)

summing over all poles with $\mbox{Re}({t})<0$ we can see that the poles in
the region $-4<{t}<0$ cancel each other, and we obtain the following sum

$\displaystyle
I^{(0)}_{3,0}=3\pi^{2}2^{7}\sum_{n=4}^{\infty}\zeta(3-2n)\zeta(2n-4)\,.$
(A.13)

This is a divergent sum, and we should regularise it. Using the integral
representation

$\displaystyle\zeta(n)=\frac{1}{\Gamma(n)}\int_{0}^{\infty}\frac{x^{n-1}}{e^{x}-1}dx\,,$
(A.14)

which is valid for $n>1$, together with the reflection relation, we can obtain
the following expression

$\displaystyle
I^{(0)}_{3,0}=-3\sum^{\infty}_{n=4}\int_{0}^{\infty}dxdy\frac{2^{10-2n}\pi^{4-2n}x^{2n-3}y^{2n-5}\cos(\pi
n)}{\left(e^{x}-1\right)\left(e^{y}-1\right)\Gamma(2n-4)}\,.$ (A.15)

We can now exchange the sum and the integral, and perform the sum over $n$ to
find

$\displaystyle
I^{(0)}_{3,0}=48\int_{0}^{\infty}dxdy\frac{x^{2}\left(2\pi\sin\left(\frac{xy}{2\pi}\right)-xy\right)}{\pi^{2}\left(e^{x}-1\right)\left(e^{y}-1\right)}\,,$
(A.16)

which is a convergent double integral that yields

$\displaystyle I^{(0)}_{3,0}=-{96\zeta(3)}-\frac{8\pi^{4}}{15}+16\pi^{2}\,.$
(A.17)

One may be worried about the regularisation procedure. However, the integral
over ${t}$ can also be computed numerically, choosing a contour slightly to
the left of the imaginary axis, and the result agrees perfectly with our
analytic result. At higher orders we can do exactly the same, but the sums
involved are a bit harder:

$I^{(0)}_{3,q}=3\sum^{\infty}_{n=4+q}\frac{(-1)^{q+1}2^{2q+7}\Gamma\left(q{-}\frac{1}{2}\right)\Gamma\left(q{+}\frac{1}{2}\right)\Gamma\left(q{+}\frac{3}{2}\right)\Gamma(2n{-}2q{-}3)}{\pi^{-1/2}\Gamma(2n-3)\Gamma(q+1)}\zeta(3{-}2n)\zeta(2n{-}2q{-}4)\,.$
(A.18)

Introducing an integral representation for the zeta functions and performing
the sum over $n$ we obtain

$\displaystyle I^{(0)}_{3,q}$
$\displaystyle=-48\pi^{2}\int_{0}^{\infty}dxdy\frac{y^{3}x^{2q+5}\Gamma\left(q-\frac{1}{2}\right)\Gamma\left(q+\frac{1}{2}\right)\Gamma\left(q+\frac{3}{2}\right)H_{q}(x,y)}{(2\pi)^{2q+9}\left(e^{x}-1\right)\left(e^{y}-1\right)\Gamma(q+1)},$
(A.19) $\displaystyle H_{q}(x,y)$
$\displaystyle=x^{2}y^{2}\,_{1}\tilde{F}_{2}\left(2;q+\frac{7}{2},q+4;-\frac{x^{2}y^{2}}{16\pi^{2}}\right)-32\pi^{2}\,_{1}\tilde{F}_{2}\left(1;q+\frac{5}{2},q+3;-\frac{x^{2}y^{2}}{16\pi^{2}}\right)\,.$

For any integer value of $q$ the regularised hypergeometric functions reduce
to something simpler. The integrals over $x$ are always of the form

$\displaystyle\int_{0}^{\infty}dx\frac{x^{a}}{e^{x}-1}$
$\displaystyle=\Gamma(a+1)\zeta(a+1)\,,$ (A.20)
$\displaystyle\int_{0}^{\infty}dx\frac{x\cos\left(\frac{xy}{2\pi}\right)}{e^{x}-1}$
$\displaystyle=\frac{1}{2}\pi^{2}\left(\frac{4}{y^{2}}-\text{csch}^{2}\left(\frac{y}{2}\right)\right),$
(A.21)
$\displaystyle\int_{0}^{\infty}dx\frac{x^{2}\sin\left(\frac{xy}{2\pi}\right)}{e^{x}-1}$
$\displaystyle=\pi^{3}\left(\frac{8}{y^{3}}-\frac{1}{2}\sinh(y)\text{csch}^{4}\left(\frac{y}{2}\right)\right)\,.$
(A.22)

For every fixed $q$, this allows to perform the integral over $x$. The
integral over $y$ can also be performed, using191919In practise, the simplest
way to perform the integrals over $y$ is to multiply by $y^{\epsilon}$, use
the expressions below and then take $\epsilon\to 0$ at the end.

$\displaystyle\int_{0}^{\infty}dy\frac{y^{a}\text{csch}^{2}\left(\frac{y}{2}\right)}{e^{y}-1}$
$\displaystyle=\Gamma(a+1)(2\zeta(a-1)-2\zeta(a))\,,$ (A.23)
$\displaystyle\int_{0}^{\infty}dy\frac{y^{a}\sinh(y)\text{csch}^{4}\left(\frac{y}{2}\right)}{e^{y}-1}$
$\displaystyle=\frac{4\Gamma(a+1)}{3}(2\zeta(a-2)-3\zeta(a-1)+\zeta(a))\,.$
(A.24)

With these results we can compute the integrals for any given $q$ and we can
actually write down the contributions analytically. We find a structure which
resembles very much $I^{(0)}_{2}(\lambda)$. More precisely we get

$\displaystyle
I^{(0)}_{3}(\lambda)=-\sum^{\infty}_{m=2}\frac{2^{2m+1}\pi^{2m-\frac{5}{2}}\lambda^{\frac{3}{2}-m}\zeta(3-2m)\Gamma\left(m-\frac{5}{2}\right)\Gamma\left(m-\frac{3}{2}\right)}{\Gamma\left(\frac{9}{2}-m\right)\Gamma(m-2)\Gamma(2m-3)}$
$\displaystyle-3\sum^{\infty}_{m=2}\kappa_{1}(m)\frac{2^{2m+2}\pi^{2m-\frac{5}{2}}\lambda^{\frac{1}{2}-m}\zeta(3-2m)\Gamma\left(m-\frac{3}{2}\right)\Gamma\left(m-\frac{1}{2}\right)}{\Gamma\left(\frac{7}{2}-m\right)\Gamma(m-1)\Gamma(2m-3)}$
(A.25)
$\displaystyle-3\sum^{\infty}_{m=2}\kappa_{2}(m)\frac{64(-1)^{m}\pi^{2m-3}\lambda^{\frac{1}{2}-m}\Gamma\left(m-\frac{5}{2}\right)\Gamma\left(m-\frac{3}{2}\right)\zeta^{\prime}(2-2m)}{\Gamma(m-1)^{2}}$
$\displaystyle-3\\!\\!\sum^{\infty}_{m,n=2}\\!\\!\kappa_{1}(m)\frac{4^{m+n+1}\pi^{2m+2n-\frac{9}{2}}\zeta(3-2m)\lambda^{-m-n+\frac{3}{2}}\Gamma\left(m+n-\frac{5}{2}\right)\Gamma\left(m+n-\frac{3}{2}\right)\zeta^{\prime}(2-2n)}{2\Gamma(2m-3)\Gamma(2n-2)\Gamma\left(-m-n+\frac{9}{2}\right)\Gamma(m+n-2)}\,,$

with

$\begin{split}\kappa_{1}(m)=2\delta_{m,2}-1\,,\qquad\kappa_{2}(m)=-\frac{(m-3)(2m-1)}{3}\,,\\\
\end{split}$ (A.26)

Note that we can combine the last two contributions into

$\displaystyle I^{(0)}_{2}(\lambda)+I^{(0)}_{3}(\lambda)$
$\displaystyle=\frac{32}{\lambda^{3/2}}-\sum^{\infty}_{m=2}\frac{64(-1)^{m}m(7-2m)\pi^{2m-3}\lambda^{\frac{1}{2}-m}\Gamma\left(m-\frac{5}{2}\right)\Gamma\left(m-\frac{3}{2}\right)\zeta^{\prime}(2-2m)}{\Gamma(m-1)^{2}}$
$\displaystyle+\sum^{\infty}_{m=2}\frac{4^{m+2}\pi^{2m-\frac{1}{2}}\lambda^{-m-\frac{1}{2}}\Gamma\left(m-\frac{1}{2}\right)\Gamma\left(m+\frac{1}{2}\right)\zeta^{\prime}(2-2m)}{\Gamma\left(\frac{5}{2}-m\right)\Gamma(m)\Gamma(2m-2)}\,.$
(A.27)

Finally, the derivative of Riemann zeta function can be simplified using the
identity

$\zeta^{\prime}(-2m)=\frac{(-1)^{m}(2m)!}{2\pi^{2m}}\zeta(2m+1)\,,$ (A.28)

valid for any $m$ strictly positive integer. We can then put all ingredients
together and compute $I^{(0)}(\lambda)$ to any desired order. The lowest few
orders match (2), while higher orders are given in the attached Mathematica
file.

###  A.2 $O(N^{0})$ at large $\lambda$

Recall that there are two contributions to the 4-body term in (2): those that
factorise in $\omega_{1},\omega_{2}$, and those that do not. The factorizable
terms $\mathcal{J}_{0}^{\text{fac}}(\omega_{1},\omega_{2})$ can be computed at
large $\lambda$ just like the 2-body terms by writing them in Mellin space
using (2.13), performing the $\omega_{1},\omega_{2}$ integrals using (2.14),
and closing each contour to the left. We obtain

$\begin{split}I^{(2)}_{\text{fac}}(\lambda)=-\frac{13\sqrt{\lambda}}{6}+\frac{55}{12}+\frac{3}{8\sqrt{\lambda}}-\frac{13\zeta(3)+6}{4\lambda}+\frac{273}{64}\frac{\zeta(3)}{\lambda^{3/2}}+O(\lambda^{-2})\,.\end{split}$
(A.29)

For the non-factorizable terms
$\mathcal{J}_{0}^{\text{fac}}(\omega_{1},\omega_{2})$, the first two lines in
(2) take a similar form as the $O(N^{2})$ terms in the previous section, i.e.
they have $(\omega_{2}^{2}-\omega_{1}^{2})$ in the denominator, so they can be
evaluated using similar methods to get an “easy” contribution

$\begin{split}I^{(2)}_{easy}(\lambda)=\frac{\sqrt{\lambda}}{6}+\frac{17}{12}+\left(-\frac{31}{8}-\frac{\pi^{2}}{6}\right){\frac{1}{\sqrt{\lambda}}}+\frac{\zeta(3)+6}{4\lambda}+\left(\frac{51\zeta(3)}{64}+\frac{5\pi^{2}}{8}\right)\frac{1}{\lambda^{3/2}}+O(\lambda^{-2})\,.\end{split}$
(A.30)

The other terms in $\mathcal{J}_{0}^{\text{fac}}(\omega_{1},\omega_{2})$ have
$(\omega_{2}^{2}-\omega_{1}^{2})^{2}$ in the denominator, so they require a
slightly more complicated approach. As with the other non-factorizable terms,
we use the identity (2.18) to write these “hard” terms as

$I^{2}_{hard}(\lambda)=\oint\frac{d{s}}{2\pi i}\frac{d{t}}{2\pi
i}p({s},{t})c({s},{t}),~{}~{}~{}c({s},{t})=\sum^{\infty}_{\begin{subarray}{c}k=1\\\
\ell=0\end{subarray}}c_{k,\ell}({s},{t})$ (A.31)

with

$p({s},{t})=3\frac{\sin(\pi{s})\sin(\pi{t})\Gamma(2{s}+4)\Gamma(2{s}+6)\Gamma(2{t}+4)^{2}\zeta(2{s}+5)\zeta(2{t}+3)\lambda^{{s}+{t}+3}}{2^{4{s}+4{t}+7}\pi^{2({s}+{t}+4)}}$
(A.32)

and

$c_{k,\ell}({s},{t})=\frac{k(2\ell+1)\Gamma(k-\ell-{s}-2)\Gamma(k+\ell-{s}-1)\Gamma(k-\ell-{t}-2)\Gamma(k+\ell-{t}-1)}{\Gamma(k-\ell+{s}+2)\Gamma(k+\ell+{s}+3)\Gamma(k-\ell+{t}+2)\Gamma(k+\ell+{t}+3)}\,.$
(A.33)

We can apply the method described above by changing variables, and writing
$c({s},{t})=\sum^{\infty}_{p=1}c_{p}({s},{t})$ where each $c_{p}({s},{t})$ is
given by a finite sum

$c_{p}({s},{t})=\sum_{q=1}^{p}\frac{q(2p-2q+1)\Gamma(p-{s}-1)\Gamma(p-{t}-1)\Gamma(-p+2q-{s}-2)\Gamma(-p+2q-{t}-2)}{\Gamma(p+{s}+3)\Gamma(p+{t}+3)\Gamma(-p+2q+{s}+2)\Gamma(-p+2q+{t}+2)}\,.$
(A.34)

In this form it can be shown that $c_{p}({s},{t})$ satisfies a second order
homogeneous recursion relation

$f_{0}(p)c_{p}({s},{t})+f_{1}(p)c_{p+1}({s},{t})+f_{2}(p)c_{p+2}({s},{t})=0\,,$
(A.35)

with

$\displaystyle f_{0}(p)$
$\displaystyle=(p+2)^{2}(p-{s}-1)(p-{s})^{2}(p-{t}-1)(p-{t})^{2}\,,$ (A.36)
$\displaystyle f_{1}(p)$ $\displaystyle=-(p-{s})(p+{s}+3)(p-{t})(p+{t}+3)$
(A.37) $\displaystyle\times$
$\displaystyle\left[\left(2p^{2}+6p+5\right){s}{t}+\left(5p^{2}+15p+11\right)({s}+{t})+2p^{4}+12p^{3}+36p^{2}+54p+29\right]\,,$
$\displaystyle f_{2}(p)$
$\displaystyle=(p+1)^{2}(p+{s}+3)^{2}(p+{s}+4)(p+{t}+3)^{2}(p+{t}+4)\,.$
(A.38)

Supplemented with the leading order terms in the large $p$ expansion, this
recursion relation allows us to compute arbitrarily high orders. The final
answer takes the form

$\displaystyle c_{p}({s},{t})$
$\displaystyle=\frac{1}{(2{s}+2{t}+5)(2{s}+2{t}+7)}\frac{1}{p^{4{s}+4{t}+13}}\left(1-\frac{2{s}+2{t}+\frac{13}{2}}{p}+\cdots\right)$
(A.39)
$\displaystyle+\frac{\pi^{3/2}4^{{s}+{t}+2}\csc(\pi{s})\csc(\pi{t})\Gamma\left({s}+{t}+\frac{7}{2}\right)}{\Gamma(2{s}+4)\Gamma(2{t}+4)\Gamma({s}+{t}+4)}\frac{1}{p^{2(3+{s}+{t})}}\left(1-\frac{{s}+{t}+3}{p}+\cdots\right)\,,$

with two distinct type of terms. Summing over $p$ term by term we obtain a
representation analogous to (A.2). The rest of the computation proceeds as
above and yields

$\begin{split}I^{2}_{hard}(\lambda)=\left(4-\frac{4\pi^{2}}{15}\right)\sqrt{\lambda}+\frac{24\zeta(3)-83}{10}+\frac{21+\pi^{2}}{6\sqrt{\lambda}}+\frac{3}{16}\frac{\left(12\zeta(3)+\pi^{2}\right)}{\lambda^{3/2}}+O(\lambda^{-2})\,.\end{split}$
(A.40)

We then sum (A.29), (A.30), and (A.40) along with the 2-body term (2.15) and
the $-12\zeta(3)$ in (2.3) to get the final answer in (2).

##  Appendix B Generalised Eisenstein series

A Generalised Eisenstein Series satisfies the inhomogeneous Laplace eigenvalue
equation

$\displaystyle(\Delta_{\tau}-s(s-1))\,{\cal
E}(s,s_{1},s_{2};\tau)=E(s_{1};\tau)E(s_{2};\tau)\,.$ (B.1)

Whereas non-holomorphic Eisenstein series arise as the coefficients of the
$R^{4}$ and $d^{4}R^{4}$ interactions in the low-energy expansion of four-
graviton amplitude in type IIB superstring theory, the coefficient of the
$d^{6}R^{4}$ interaction is the generalised non-holomorphic Eisenstein series
${\cal E}(4,{\scriptstyle\frac{3}{2}},{\scriptstyle\frac{3}{2}};\tau)$ [35].
More general examples with half-integral indices in the source term arise in
the low-energy expansion of two-loop eleven-dimensional supergravity
amplitudes compactified on a two-torus [34]. This is the prototype for this
general type of modular functions that arise in ${\cal H}_{N}^{i}(\tau)$. For
generic eigenvalue $s$, we define the generalised Eisenstein series as the
unique modular invariant solution to (B.1) satisfying the boundary condition
that the coefficient of the homogeneous solution $\tau_{2}^{s}$ vanishes.

Within string theory contexts, the eigenvalue spectrum of ${\cal
E}(s,s_{1},s_{2};\tau)$ is sensitive to whether $s_{1}$, $s_{2}$ are integer
or half-integer. It has been argued [48] that when
$s_{1},s_{2}\in{\mathbb{N}}$ the spectrum of $s$ is restricted to
$s\in\mathrm{Spec}(s_{1},s_{2})=\\{|s_{1}-s_{2}|+2,\,|s_{1}-s_{2}|+4,\dots,s_{1}+s_{2}-4,\,s_{1}+s_{2}-2\\}$.
However, for the generalised Eisenstein series with
$s_{1},s_{2}\in{\mathbb{N}}+{1\over 2}$ arising as coefficients of the
$1/N^{r}$ terms in the large-$N$ expansion of the integrated correlator ${\cal
H}_{N}$ [33], the spectrum is found to be restricted to
$s=\\{s_{1}+s_{2}+1,s_{1}+s_{2}+3,\dots,3r+1\\}$. This is discussed further in
the main text.

The generalised Eisenstein series ${\cal E}(s,s_{1},s_{2};\tau)$ can be
decomposed in Fourier modes in the form

${\cal E}(s,s_{1},s_{2};\tau)=\sum_{n,m=0}^{\infty}{\cal
E}^{(n,m)}(s,s_{1},s_{2};\tau_{2})q^{n}\bar{q}^{m}\,,$ (B.2)

where $q=e^{2\pi i\tau}$ and $\bar{q}=e^{-2\pi i\bar{\tau}}$. The non-zero
modes represent the contributions of instantons and anti-instantons, with
instanton number $k=n-m$, where $k>0$ for instantons and $k<0$ for anti-
instantons. Terms with $k=0$ get purely perturbative contributions from the
$m=n=0$ term and instanton/anti-instanton contributions from terms with
$m=n>0$.

The Laurent polynomial of ${\cal E}(s,s_{1},s_{2};\tau_{2})$ (i.e. the ${\cal
E}^{(0,0)}$ term) is given by

$\displaystyle{\cal E}^{(0,0)}(s,s_{1},s_{2};\tau_{2})=$
$\displaystyle\frac{4\pi^{-2s_{1}-2s_{2}}\zeta(2s_{1})\zeta(2s_{2})}{\pi^{s_{1}+s_{2}}(s_{1}+s_{2}-s)(s_{1}+s_{2}+s-1)}(\pi\tau_{2})^{s_{1}+s_{2}}+\frac{4{\pi^{-{1\over
2}-2s_{2}}}\Gamma(s_{1}-{\scriptstyle\frac{1}{2}})\zeta(2s_{1}-1)\zeta(2s_{2})}{(s_{2}-s_{1}+s)(s_{2}-s_{1}-s+1)\Gamma(s_{1})}(\pi\tau_{2})^{1-s_{1}+s_{2}}$
$\displaystyle+$ $\displaystyle\frac{4\pi^{-{1\over
2}-2s_{1}}\Gamma(s_{2}-{\scriptstyle\frac{1}{2}})\zeta(2s_{2}-1)\zeta(2s_{1})}{(s_{1}-s_{2}+s)(s_{1}-s_{2}-s+1)\Gamma(s_{2})}(\pi\tau_{2})^{1-s_{2}+s_{1}}$
(B.3) $\displaystyle+$
$\displaystyle\frac{4\Gamma(s_{1}-{\scriptstyle\frac{1}{2}})\Gamma(s_{2}-{\scriptstyle\frac{1}{2}})\zeta(2s_{1}-1)\zeta(2s_{2}-1)}{\pi(s_{1}+s_{2}-s-1)(s_{1}+s_{2}+s-2)\Gamma(s_{1})\Gamma(s_{2})}(\pi\tau_{2})^{2-s_{1}-s_{2}}+\beta(s,s_{1},s_{2})\,(\pi\tau_{2})^{1-s}\,,$

where the first four terms can be obtained by matching powers of $\tau_{2}$ on
the left-hand and right-hand sides of (B.1) and the last term satisfies the
homogeneous equation and its coefficient is given by [35, 37]202020In these
references, this coefficient is determined by projecting the Laplace equation
(1.10) on $E(s;\tau)$. Alternatively, it may be obtained from a Poincaré
series representation as in [48].

$\beta(s,s_{1},s_{2})={4\pi^{s-1}\over\Gamma(s_{1})\Gamma(s_{2})}\,{\zeta^{*}(s-s_{1}-s_{2}+1)\zeta^{*}(s+s_{1}-s_{2})\zeta^{*}(s-s_{1}+s_{2})\zeta^{*}(s+s_{1}+s_{2}-1)\over(1-2s)\,\zeta^{*}(2s)}\,,$
(B.4)

with $\zeta^{*}(s)=\zeta(s)\Gamma(s/2)/\pi^{s/2}$.

In addition to the Laurent expansion, the $k=0$ sector contains contributions
from instanton/anti-instanton pairs, $(n,n)$ with $n>0$, which have the form

${\cal
E}^{(n,n)}(s,s_{1},s_{2};\tau_{2})=\frac{n^{s_{1}+s_{2}-2}\sigma_{1-2s_{1}}(n)\sigma_{1-2s_{2}}(n)}{\Gamma(s_{1})\Gamma(s_{2})}\Phi_{s,s_{1},s_{2}}(4\pi
n\tau_{2})\,.$ (B.5)

When $s_{1},s_{2}\in{\mathbb{N}}$ the expression $\Phi_{s,s_{1},s_{2}}(4\pi
n\tau_{2})$ is a polynomial in inverse powers of $\tau_{2}$ of degree
$s_{2}+s_{1}-2$. The first two perturbative orders in the $(n,n)$ sector do
not depend on the eigenvalue $s$ and have the form

$\displaystyle\Phi_{s,s_{1},s_{2}}(4\pi n\tau_{2})=\frac{8}{(4\pi
n\tau_{2})^{2}}+\frac{8[s_{1}(s_{1}-1)+s_{2}(s_{2}-1)-4]}{(4\pi
n\tau_{2})^{3}}+O(\tau_{2}^{-4})\,.$ (B.6)

Higher-order coefficients do depend on $s$ and can be computed from the
differential equation (B.1) as performed in [33, 47] or via resurgence
analysis methods [62, 63]. The contribution of the general $(n,m)$ sector with
$n,m\neq 0$ is given by

${\cal
E}^{(n,m)}(s,s_{1},s_{2};\tau_{2})=\frac{n^{s_{1}-1}m^{s_{2}-1}\sigma_{1-2s_{1}}(n)\sigma_{1-2s_{2}}(m)}{\Gamma(s_{1})\Gamma(s_{2})}\Big{(}\frac{1}{4nm\,(\pi\tau_{2})^{2}}+O(\tau_{2}^{-3})\Big{)}+(n\leftrightarrow
m)\,.$ (B.7)

These terms correspond to non-perturbative contributions in the topological
sector with instanton number $n-m$.

For the purely instantonic sector, i.e. the $(n,0)$ term, and purely anti-
instantonic sector, i.e. the $(0,m)$ term, we need to be more careful since in
these sectors we now have the freedom of adding homogeneous solutions to the
differential equation (B.1). An algorithm for constructing a particular
solution, ${\cal E}_{{\rm p}}(s,s_{1},s_{2};\tau)$, to (B.1) was proposed in
[33], and later generalised in [47]. As discussed in the main text, whenever
the eigenvalue $s$ is such that the vector space, $\mathcal{S}_{2s}$, of
holomorphic cusp forms with weight $2s$ is non-trivial, this particular
solution ${\cal E}_{{\rm p}}(s,s_{1},s_{2};\tau)$ does not correspond to a
modular invariant solution to (B.1), i.e. the particular solution does not
correspond to the generalised Eisenstein series.

To construct the generalised Eisenstein series ${\cal E}(s,s_{1},s_{2};\tau)$,
a certain homogeneous solution has to be added to the particular solution, so
that the complete solution is a modular function. Since we have fixed the
boundary condition in the perturbative, i.e. the $(0,0)$ sector, this is only
relevant for the instantonic (and anti-instantoic) pieces, namely we can only
add an homogeneous solution to ${\cal E}_{\rm
p}^{(n,0)}(s,s_{1},s_{2};\tau_{2})$ (and ${\cal E}_{\rm
p}^{(0,m)}(s,s_{1},s_{2};\tau_{2})$).

Based on [48, 38], we can use the novel results [39] to show that even in the
case of present interest the homogeneous solution that must be added is given
by

$\displaystyle\lambda_{\Delta}(s,s_{1},s_{2})\cdot H_{\Delta}(\tau)\,,$ (B.8)

where the coefficient $\lambda_{\Delta}(s,s_{1},s_{2})$ is presented in (4.7).
The function $H_{\Delta}(\tau)$ takes the form,

$H_{\Delta}(\tau)\equiv\sum_{n=1}^{\infty}\frac{a_{\Delta}(n)}{n^{s}}\sqrt{n\tau_{2}}K_{s-\frac{1}{2}}(2\pi
n\tau_{2})\big{(}e^{2\pi in\tau_{1}}+e^{-2\pi in\tau_{1}}\big{)}\,,$ (B.9)

and satisfies the homogeneous equation

$(\Delta_{\tau}-s(s-1))H_{\Delta}(\tau)=0\,.$ (B.10)

Note the important fact that $H_{\Delta}(\tau)$ alone is not modular
invariant. However, the addition of a suitable multiple of $H_{\Delta}(\tau)$
to the particular solution ${\cal E}_{{\rm p}}(s,s_{1},s_{2};\tau)$ is crucial
for obtaining a modular invariant solution to (B.1), i.e. for constructing the
generalised Eisenstein series

${\cal E}(s,s_{1},s_{2};\tau)={\cal E}_{{\rm
p}}(s,s_{1},s_{2};\tau)+\sum_{\Delta\in{\cal
S}_{2s}}\lambda_{\Delta}(s,s_{1},s_{2})H_{\Delta}(\tau)\,.$ (B.11)

##  Appendix C Properties of lattice sums

In general there is no known lattice integral representation for an individual
generalised Eisenstein series. However, we shall now describe the construction
of a lattice representation of the special linear combinations of generalised
Eisenstein series which are relevant for the integrated correlator under
discussion (see [34, 35, 37, 44, 64] for further details). This gives a four-
dimensional lattice representation that generalises the two-dimensional
lattice representation of a non-holomorphic Eisenstein series.

###  C.1 A family of lattice sum integrals

Following [35], we here consider a class of lattice sum integrals constructed
from the special non-holomorphic functions $A_{i,j}(\rho)$, called modular
local polynomials, introduced in section 4.2 and whose precise definition and
properties will be reviewed in appendix C.2. We will consider the lattice sum
integral in the form212121In earlier lieterature, such as [35] the variables
$t_{1}$, $t_{2}$ $t_{3}$ corresponded to ‘inverse Schwinger parameters’
$t_{i}=\hat{t}_{i}/(\hat{t}_{1}\hat{t}_{2}+\hat{t}_{1}\hat{t}_{3}+\hat{t}_{2}\hat{t}_{3})$,
where $\hat{t}_{i}$ ($i=1,2,3$) are the Schwinger parameters for a two-loop
vacuum diagram.

$\displaystyle{\cal
E}^{w}_{i,j}(\tau)=\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\neq
0}}\int_{0}^{\infty}d^{3}t\,[V(t)]^{w-3}\,A_{i,j}(\rho(t))\,\exp\Big{(}-\frac{\pi}{\tau_{2}}\sum_{i=1}^{3}t_{i}|p_{i}|^{2}\Big{)}\,,$
(C.1)

where the parameters $\rho(t)$ and $V(t)$ are defined in (4.18). The $p_{i}$
are discrete loop momenta of a two-loop Feynman diagram on a two-torus of
complex structure $\tau$ and are defined by

$p_{i}=m_{i}+n_{i}\tau\,,$ (C.2)

with $m_{i},n_{i}\in{\mathbb{Z}}$.

The analysis of the inhomogeneous Laplace equation satisfied by ${\cal
E}^{w}_{i,j}(\tau)$ will demonstrate that the parameter $w\in\mathbb{R}^{+}$
corresponds to the total transcendental weight of the sources, i.e. the sum of
the indices of the bilinears in Eisenstein series which appear as source
terms.

In [35] a slightly different lattice sum was considered which gave rise to
ultraviolet divergences arising at the boundary of the integration region from
particular $SL(2,{\mathbb{Z}})$ orbits of the lattice sums variables. With the
present constraints on the lattice sum, the integral (C.1) is completely well-
defined since it is exponentially convergent at all boundaries of the domain
of integration.

###  C.2 Modular local polynomials: definition and properties

We now present the definition and clarify the most important properties of the
non-holomorphic functions $A_{i,j}(\rho)$, which feature prominently in the
formula (C.1) for the modular invariant functions under consideration.

The function $A_{1,0}(\rho)$ arose in the construction of the $d^{6}R^{4}$
coefficient ${\cal E}(4;\frac{3}{2},\frac{3}{2};\tau)$ in [35], and a more
general discussion in [37] presented the general procedure for constructing
particular combinations of ${\cal E}(s,s_{1},s_{2};\tau)$ once the
$A_{i,j}(\rho)$ are known. The following general construction of
$A_{i,j}(\rho)$ was given by Don Zagier (in unpublished notes that are
reproduced in [44] with more details).

The result of this analysis will be that $A_{i,j}(\rho)$ are Laurent
polynomials in $\rho_{2}$ with coefficients which are polynomial in
$\rho_{1}$, and inside the domain $\rho\in\mathcal{H}\backslash\Gamma_{0}(2)$
they satisfy the homogeneous Laplace equation (4.20),

$\left[\Delta_{\rho}-s(s-1)\right]\,A_{i,j}(\rho)=0\,,$ (C.3)

where the eigenvalue $s=3i+j+1$ (for example, $s=4$ for $A_{1,0}$).

The first step in constructing these Laurent series is to define the complex
combinations

$\displaystyle u=\rho^{2}(1-\rho)^{2}\,,\qquad v=\rho^{2}-\rho+1\,.$ (C.4)

The functions $A_{i,j}(\rho)$ are given by

$\displaystyle A_{i,j}(\rho)=D_{-2n}^{(n)}(u^{i}v^{j})\,,\hskip
72.26999ptn=s-1=3i+j\quad\hbox{ with }\quad i,j\geq 0\,.$ (C.5)

The polynomials $u^{i}v^{j}$, holomorphic in $\rho$, can be understood
precisely as the modular local polynomials of modular weight $-2n$ discussed
in [50]. Here $D_{-2n}^{(n)}$ is the iterated derivative operator

$\displaystyle D_{-2n}^{(n)}={(-2i)^{n}n!\over(2n)!}D_{-2}\circ
D_{-4}\circ\cdots\circ D_{-2n+2}\circ D_{-2n}\,,$ (C.6)

where the differential operator $D_{k}=\partial_{\rho}+k/(\rho-\bar{\rho})$
satisfies $\Delta_{k+2}\cdot D_{k}-D_{k}\Delta_{k}=-kD_{k}$, where
$\Delta_{k}=4D_{k-2}\,\rho_{2}^{2}\partial_{\bar{\rho}}$ is the laplacian
acting on weight $k$ modular forms. Explicitly, the operator $D_{-2n}^{(n)}$
can be written as

$\displaystyle
D_{-2n}^{(n)}={(-2i)^{n}n!\over(2n)!}\sum_{m=0}^{n}\begin{pmatrix}n\cr
m\end{pmatrix}\frac{(-n-m)_{m}}{(\rho-\bar{\rho})^{m}}\,\frac{\partial^{n-m}}{\partial\rho^{n-m}}\,,$
(C.7)

where $(x)_{m}$ is the Pochhammer symbol.

The functions $A_{i,j}(\rho)$ constructed in this manner satisfy the Laplace
equation (C.3) in the interior of the fundamental domain
$\mathcal{H}\backslash\Gamma_{0}(2)$. Furthermore they are Laurent expansions
that have the following form,

$\displaystyle
A_{i,j}(\rho)=\sum_{k=0}^{2i+j}A_{i,j}^{(k)}(\rho_{1})\,\rho_{2}^{i+j-2k}\,,$
(C.8)

where $A_{i,j}^{(k)}(\rho_{1})$ is a polynomial of degree $k$ in
$\rho_{1}(1-\rho_{1})$. Since it is invariant under $\rho_{1}\mapsto
1-\rho_{1}$, it may be expressed as a linear combination of Bernoulli
polynomials $B_{2k}(\rho_{1})$ of even index. It is often clarifying to re-
express $\rho_{1},\rho_{2}$ in terms of $t_{1},t_{2},t_{3}$. The function
$A_{i,j}$ is then by construction a homogenous function of the $t_{i}$’s,
invariant under permutations.

### Examples of $A_{i,j}$

From (C.5), it is easy to generate expressions for $A_{i,j}$ with very high
eigenvalues $s=3i+j+1$. Here we display the $A_{i,j}$ that appear in the main
body of the paper – in particular in (4.29)-(4.31). For convenience, we
express $A_{i,j}$ in terms of the $t$-variables using (4.18). Since
$A_{i,j}(\rho(t))$ is a symmetric function of $t_{i}$ we shall make use of a
basis $(\sigma_{1},\sigma_{2},\sigma_{3})$ of symmetric polynomials defined by

$\sigma_{1}\equiv t_{1}+t_{2}+t_{3}\,,\qquad\sigma_{2}\equiv
Y(t)=t_{1}t_{2}+t_{1}t_{3}+t_{2}t_{3}\,,\qquad\sigma_{3}\equiv
t_{1}t_{2}t_{3}\,.$ (C.9)

Below we list examples of $A_{i,j}$ that are relevant for ${\cal
H}_{N}^{i}(\tau)$ up to order $1/N^{3}$,

$\displaystyle\underline{s=4}:$
$\displaystyle\sigma_{2}^{{\scriptstyle\frac{3}{2}}}A_{1,0}=\frac{1}{5}\sigma_{1}\sigma_{2}-\sigma_{3}\,$
(C.10) $\displaystyle\underline{s=5}:$
$\displaystyle\,\sigma_{2}^{2}A_{1,1}=\frac{1}{7}\sigma_{1}^{2}\sigma_{2}+\frac{2}{35}\sigma_{2}^{2}-\sigma_{1}\sigma_{3}\,,$
(C.11) $\displaystyle\underline{s=6}:$
$\displaystyle\sigma_{2}^{{\scriptstyle\frac{5}{2}}}A_{1,2}=\frac{1}{9}\sigma_{1}^{3}\sigma_{2}+\frac{1}{21}\sigma_{1}\sigma_{2}^{2}-\sigma_{1}^{2}\sigma_{3}+\frac{1}{3}\sigma_{2}\sigma_{3}\,,$
(C.13)
$\displaystyle\sigma_{2}^{{\scriptstyle\frac{5}{2}}}A_{0,5}=\sigma_{1}^{5}-\frac{10}{3}\sigma_{1}^{3}\sigma_{2}+\frac{15}{7}\sigma_{1}\sigma_{2}^{2}\,,$
$\displaystyle\underline{s=7}:$
$\displaystyle\,\sigma_{2}^{3}A_{2,0}=\frac{1}{33}\sigma_{1}^{2}\sigma_{2}^{2}+\frac{6}{11}\sigma_{1}\sigma_{2}\sigma_{3}+\frac{16}{231}\sigma_{2}^{3}+\sigma_{3}^{2}\,,$
(C.14) $\displaystyle\underline{s=8}:$
$\displaystyle\sigma_{2}^{\scriptstyle\frac{7}{2}}A_{2,1}=\frac{3}{143}\sigma_{1}^{3}\sigma_{2}^{2}+\frac{28}{429}\sigma_{1}\sigma_{2}^{3}-\frac{6}{13}\sigma_{1}^{2}\sigma_{2}\sigma_{3}+\sigma_{1}\sigma_{3}^{2}-\frac{4}{143}\sigma_{2}^{2}\sigma_{3}\,,$
(C.16)
$\displaystyle\sigma_{2}^{\scriptstyle\frac{7}{2}}A_{1,4}=\frac{1}{13}\sigma_{1}^{5}\sigma_{2}{-}\sigma_{1}^{4}\sigma_{3}{+}\frac{2}{143}\sigma_{1}^{3}\sigma_{2}^{2}{+}\frac{18}{13}\sigma_{1}^{2}\sigma_{2}\sigma_{3}{-}\frac{43}{429}\sigma_{1}\sigma_{2}^{3}{-}\frac{27}{143}\sigma_{2}^{2}\sigma_{3}\,,$
$\displaystyle\underline{s=10}:$
$\displaystyle\sigma_{2}^{\scriptstyle\frac{9}{2}}A_{3,0}=\frac{1}{221}\sigma_{1}^{3}\sigma_{2}^{3}-\frac{3}{17}\sigma_{1}^{2}\sigma_{2}^{2}\sigma_{3}+\frac{72}{2431}\sigma_{1}\sigma_{2}^{4}+\frac{15}{17}\sigma_{1}\sigma_{2}\sigma_{3}^{2}$
(C.18)
$\displaystyle\qquad\quad\phantom{=}-\frac{24}{221}\sigma_{2}^{3}\sigma_{3}-\sigma_{3}^{3}\,,$
$\displaystyle\sigma_{2}^{\scriptstyle\frac{9}{2}}A_{2,3}=\frac{1}{85}\sigma_{1}^{5}\sigma_{2}^{2}-\frac{6}{17}\sigma_{1}^{4}\sigma_{2}\sigma_{3}+\frac{11}{221}\sigma_{1}^{3}\sigma_{2}^{3}+\sigma_{1}^{3}\sigma_{3}^{2}+\frac{2}{17}\sigma_{1}^{2}\sigma_{2}^{2}\sigma_{3}$
$\displaystyle\qquad\quad\phantom{=}-\frac{4}{143}\sigma_{1}\sigma_{2}^{4}-\frac{9}{17}\sigma_{1}\sigma_{2}\sigma_{3}^{2}+\frac{4}{221}\sigma_{2}^{3}\sigma_{3}\,.$

##  Appendix D From lattice sums to generalised Eisenstein series

In order to determine the properties of the modular invariant functions which
can be represented via the lattice sum integrals just discussed, it is very
useful to apply the laplacian $\Delta_{\tau}$ to (C.1) which gives

$\displaystyle\Delta_{\tau}{\cal
E}^{w}_{i,j}(\tau)=\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\neq
0}}\int_{0}^{\infty}d^{3}t\,[V(t)]^{w-3}\,A_{i,j}(\rho(t))\,\Delta_{t}\exp\Big{(}-\frac{\pi}{\tau_{2}}\sum_{i=1}^{3}t_{i}|p_{i}|^{2}\Big{)}\,,$
(D.1)

where
$\Delta_{\rho}=\rho_{2}^{2}(\partial_{\rho_{1}}^{2}+\partial_{\rho_{2}}^{2})$,
and we have used (4.23). We now integrate by parts and note the property

$\Delta_{t}[V(t)^{\alpha}F(t)]=V(t)^{\alpha}\Delta_{t}[F(t)]\,,$ (D.2)

for all $\alpha$. We will also define

$B^{w}_{i,j}(t)\equiv[V(t)]^{w-3}A_{i,j}(\rho(t))\,,$ (D.3)

and note that (C.3) and (D.2) imply

$\Delta_{t}\Big{[}B^{w}_{i,j}(t)\Big{]}=[V(t)]^{w-3}\Big{[}\Delta_{\rho}A_{i,j}(\rho)\Big{]}_{\rho=\rho(t)}=s(s-1)B^{w}_{i,j}(t)\,,$
(D.4)

with $s=3i+j+1$. It follows that

$\displaystyle\Delta_{\tau}{\cal E}^{w}_{i,j}(\tau)$
$\displaystyle=\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\neq
0}}\int_{0}^{\infty}d^{3}t\,B^{w}_{i,j}(t)\,\Delta_{t}\Big{[}\exp\Big{(}-\frac{\pi}{\tau_{2}}\sum_{i=1}^{3}t_{i}|p_{i}|^{2}\Big{)}\Big{]}$
$\displaystyle=\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\neq
0}}\int_{0}^{\infty}d^{3}t\,\Delta_{t}\big{[}B^{w}_{i,j}(t)\big{]}\,\exp\Big{(}-\frac{\pi}{\tau_{2}}\sum_{i=1}^{3}t_{i}|p_{i}|^{2}\Big{)}+{\rm{b.t.}}$
$\displaystyle=s(s-1){\cal E}^{w}_{i,j}(\tau)+{\rm{b.t.}}\,,$ (D.5)

where ‘${\rm{b.t.}}$’ represents boundary terms.

The boundary terms in (D.5) are easily collected by integrating $\Delta_{t}$
by parts using the definition (4.22). We start by noting that since all the
momenta $p_{1},p_{2},p_{3}=-p_{1}-p_{2}$ are non-vanishing we only need to
worry about boundary contributions coming from $t_{i}\to 0$ since the integral
is exponentially suppressed along any direction as $t_{i}\to\infty$. Using the
fact that the integrand is invariant under permutations of
$(t_{1},t_{2},t_{3})$ we arrive at the expression

$\displaystyle{\rm{b.t.}}=$
$\displaystyle\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\neq
0}}3\int_{0}^{\infty}d^{2}t\,\Big{[}t_{1}t_{2}(\partial_{3}-\partial_{1}-\partial_{2})B^{w}_{i,j}(t_{1},t_{2},t_{3})\Big{]}_{t_{3}=0}\exp\Big{(}\\!-\frac{\pi
t_{1}|p_{1}|^{2}}{\tau_{2}}-\frac{\pi t_{2}|p_{2}|^{2}}{\tau_{2}}\Big{)}$
$\displaystyle+\sum_{p\neq
0}\frac{6\pi|p|^{2}}{\tau_{2}}\int_{0}^{\infty}d^{2}t\,\Big{[}t_{1}t_{2}B^{w}_{i,j}(t_{1},t_{2},t_{3})\Big{]}_{t_{3}=0}\exp\Big{(}-\frac{\pi|p|^{2}(t_{1}+t_{2})}{\tau_{2}}\Big{)}\,,$
(D.6)

where we have defined $p=m+n\tau$ with $(m,n)\in\mathbb{Z}^{2}/\\{(0,0)\\}$.
It is easy to check with the definition of $A_{i,j}$ given above that the
first term produces terms bilinear in Eisenstein series, while the second term
will produce terms linear in Eisenstein series.

Although we do not have a closed formula for the boundary term associated with
a specific $B^{w}_{i,j}$, we can make some general observations. It follows
from (C.8) that when $A_{i,j}(\rho)$ is rewritten in terms of $t_{i}$ and
expressed in terms of the symmetric variables
$\sigma_{1},\sigma_{2},\sigma_{3}$ defined in (C.9), it takes the general form

$A_{i,j}(\rho(t))=\sum_{\alpha=0}^{s-1}\sum_{\beta=\frac{s\,{\rm
mod}\,2}{2}}^{\frac{s-1}{2}}\sum_{\gamma=0}^{\lfloor\frac{s-1}{3}\rfloor}c(\alpha,\beta,\gamma)\sigma_{1}^{\alpha}\sigma_{2}^{-\beta}\sigma_{3}^{\gamma}\,,$
(D.7)

where the coefficients $c(\alpha,\beta,\gamma)\in\mathbb{Q}$ implicitly depend
on $i,j$ (and $s=3i+j+1$).

We note that when $s$ is even $\beta\in\mathbb{N}$ while for $s$ odd
$\beta\in\mathbb{N}+\frac{1}{2}$. Furthermore the $A_{i,j}$ are homogeneous
functions of $t_{i}$ since $\rho_{1}$ and $\rho_{2}$ (defined in (4.18)) are,
and hence it follows that

$\alpha-2\beta+3\gamma=0\,.$ (D.8)

From (4.18) we also deduce that the transformation $\rho\to\bar{\rho}$ is
equivalent to $t_{i}\to-t_{i}$, or equivalently

$\displaystyle\sigma_{1}\to-\sigma_{1}\,,\qquad\sigma_{2}\to\sigma_{2}\,,\qquad\sigma_{3}\to-\sigma_{3}\,.$
(D.9)

Hence, using (C.8), we see that

$\displaystyle
A_{i,j}(\bar{\rho})=(-1)^{i+j}A_{i,j}(\rho)=(-1)^{s-1}A_{i,j}(\rho)$
$\displaystyle\Longrightarrow\qquad$ $\displaystyle
A_{i,j}(\rho(-t))=A_{i,j}(\bar{\rho}(t))=(-1)^{s-1}A_{i,j}(\rho(t))\,,$ (D.10)

which from the expansion (D.7) implies

$c(\alpha,\beta,\gamma)=0\,,\qquad{\rm
for}\qquad\alpha+\gamma\not\equiv(s-1)\,~{}{\rm mod}\,~{}2\,.$ (D.11)

We can then focus on the integrand (D.3) and expand it as

$B^{w}_{i,j}(t_{1},t_{2},t_{3})=\sum_{\alpha=0}^{s-1}\sum_{\beta=\frac{s\,{\rm
mod}\,2}{2}}^{\frac{s-1}{2}}\sum_{\gamma=0}^{\lfloor\frac{s-1}{3}\rfloor}c(\alpha,\beta,\gamma)\,\sigma_{1}^{\alpha}\,\sigma_{2}^{\frac{w-3}{2}-\beta}\sigma_{3}^{\gamma}\,,$
(D.12)

from which it is rather easy to compute the boundary contribution (D.6) where
only the terms with $\gamma=0$ and $\gamma=1$ contribute. The boundary terms
are given by

$\displaystyle{\rm{b.t.}}=$
$\displaystyle\,-\\!\\!\\!\\!\sum_{\alpha=(s-1)\,{\rm
mod}\,2}^{s-1}\,\,\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\neq
0}}3\alpha
c(\alpha,\tfrac{\alpha}{2},0)\int_{0}^{\infty}d^{2}t\,(t_{1}t_{2})^{\frac{w-1-\alpha}{2}}(t_{1}+t_{2})^{\alpha-1}e^{-\frac{\pi
t_{1}|p_{1}|^{2}}{\tau_{2}}-\frac{\pi t_{2}|p_{2}|^{2}}{\tau_{2}}}$
$\displaystyle+\sum_{\alpha=s\,{\rm
mod}\,2}^{s-2}\,\,\sum_{\underset{p_{1}+p_{2}+p_{3}=0}{p_{1},p_{2},p_{3}\neq