$1/2$ BPS Structure Constants and Random Matrices

Adolfo Holguin

Department of Physics, University of California Santa Barbara, CA 93106

###### Abstract

We study three point functions of half BPS operators in $\mathcal{N}=4$ super
Yang-Mills theory focusing on correltors of two of the operators with
dimension of order $\Delta\sim N^{2}$ and a light single trace operator. These
describe vacuum expectation values of type IIB supergravity modes in LLM
backgrounds that do not necessarily preserve the same symmetries as the
background solution. We propose a class of complex matrix models that fully
capture the combinatorics of the problem, and describe their solution in the
large $N$ limit. In simple regimes when the dual description is in terms of
widely separated condensates of giant gravitons we find that the models are
solvable in the large $N$ and can be approximated by unitary Jacobi ensembles;
we describe how these distributions are reproduced in the dual bubbling
geometry picture for large droplets. In the case of two eigenvalue droplets
the model is exactly solvable at finite $N$. As a result we compute all half-
BPS structure constants of heavy-heavy-light type, and reproduce the formulas
found via holographic renormalization in the large $N$ limit. We also comment
on structure constants of three heavy operators.

##  1 Introduction

The study of baryonic operators in large $N$ gauge theories is an old subject
[1] that has received renewed attention in the context of holographic field
theories [2, 3, 4]. Such operators are extremely interesting from the point of
view of the large $N$ expansion, since they correspond to heavy non-
perturbative objects that are not very well captured by the conventional t’
Hooft expansion. On physical grounds one expects that such heavy objects
modify the physics at the semi-classical level, and that one should attempt to
approach the problem from a point of view where one treats the dynamics of the
many constituents of the object in terms of a simpler collective coordinates
[5, 6]. This is well understood in string theories; heavy objects can lead to
non-trivial boundary conditions for strings, or in some cases deform the
target space geometry which the string probes. For this reason such an
approach is essential for understanding how gravitational physics arises from
large $N$ models. In examples of the AdS/CFT correspondence [7] these ideas
have sharp realizations in terms of giant gravitons, and in suitable limits,
non-trivial supergravity backgrounds like bubbling geometries and black holes.
Given that maximally supersymmetric Yang-Mills theory is expected to be fully-
fledged theory of quantum gravity, a particularly interesting question to
address is how $\mathcal{N}=4$ SYM resolves the many puzzles of gravitational
theories, in particular the physics of black holes. At the moment such
questions are out of reach and they lie in regions of parameter space where
current non-perturbative techniques such as integrability are expected to
fail. One particular fruitful approach has been to concentrate on observables
which are protected by supersymmetry in order to test and develop tools, and
the half-BPS sector of asymptotically AdS type IIB supergravity and
$\mathcal{N}=4$ SYM is perhaps the simplest non-trivial toy model.

The spectral problem in this sector of the $U(N)$ theory was solved by the
work of Corley, Jevicki, and Ramgoolam [8]; half-BPS operators are Schur
functions and their structure constants are given by multiplicities of
representations of the unitary group

$\left<\mathcal{O}_{R_{3}}(\bar{Z})\mathcal{O}_{R_{2}}(\bar{Z})\mathcal{O}_{R_{1}}(Z)\right>=C_{R_{1}R_{2}R_{3}}f_{R_{1}},$
(1.1)

where $C_{R_{1}R_{2}R_{3}}$ are Richardson-Littlewood coefficients and $f_{R}$
is the norm of the operator $\mathcal{O}_{R}$. Although this solves the
problem in principle and combinatorial algorithms exist which generate these
coefficients for a fixed value of $N$ it is unclear how the asymptotics of
these coefficients are reflected in the corresponding supergravity solutions.
Since these numbers appear naturally in the study of the intersection theory
of Grassmannians, a natural expectation is that there is an alternative
description for such calculations involving only geometric data coming from
the gauge group of the theory. Another issue is that most results in the
existing literature on structure constants of half-BPS operators either focus
on single trace operators, or in operators preserving the same
supersymmetries, or rely heavily on the free fermion description of model.
Holographic computations of one point-functions in half-BPS backgrounds have
also been studied in generality, but explicit calculations are limited to
maximally charged operators which are charged under the same symmetry as the
background and to operators of low dimensions. So an important step towards
understanding very heavy operators in $\mathcal{N}=4$ SYM is to develop tools
that can tackle problems of this kind for generic BPS operators in the large
$N$ limit. Such tools have been developed recently for operators of dimension
$\Delta\sim N$ [2, 3, 9, 10, 11, 12] and in this paper we extend this to
operators of dimensions that scale as $N^{2}$. We show that the computation of
very generic three-point functions of half-BPS operators can be packaged in a
large family complex matrix model of matrices valued on a Grassmannian.
Although we mostly focus on the $U(N)$ theory, our results generalize readily
to orthogonal and symplectic gauge groups. For simple observables, such as
set-ups involving a single stack of AdS giant gravitons, the corresponding
matrix ensemble is a unitary Jacobi ensemble, while for more generic
observables the matrix model cannot be easily reduced to integrals over
eigenvalues. At large $N$, we find that the saddle point equations simplify
the calculation significantly allowing us to either reduce the integrals to
sums over integrals over eigenvalues, where each term in the sum is labelled
by a permutation. The average density of eigenvalues is universal and is given
by the well-known Marchenko-Pastur distribution, which appears as the Poisson
distribution of noncommutative probability theory [13]

$\rho(z)=\frac{\sqrt{(z_{+}-z)(z-z_{-})}}{2\pi z}\;.$ (1.2)

Gaussian matrix integrals have been used extensively in the study of the
combinatorics of half-BPS correlators in the past [14, 15] and in other
contexts [16, 17]; the matrix models we study on the other hand describe large
deviations from the vacuum state. for a similar story for Wilson loops see
[18, 19, 20]. In simple terms they describe wavefunctions of semi-classical
BPS states in $\mathcal{N}=4$ SYM and as such they provide a quantum
mechanical description of half-BPS Coulomb branch configurations. This makes
them ideal candidates for computing quantities that one can match to the dual
geometric description. In fact, we argue that despite the fact that one
expects an exact match for half-BPS observables on both sides of the duality
due to a lack of $g^{2}N$ corrections to the free-field theory answer, the
corresponding supergravity will not in general compute a precise quantity from
the point of view of the conformal field theory but rather a (micro-canonical)
average. This is purely an effect of the large $N$ limit and the full stringy
description should be able to resolve the details of the boundary observable.
These observations have been made in this context before, for instance in
[21], but we clarify how this happens on the field theory side of the
computation. Our main result is a large $N$ formula for all heavy-heavy-light
structure constants of half-BPS operators, for instance

$C_{R\,R\,L/2}=\frac{1}{\sqrt{L}}\int
dzd\bar{z}\,\rho_{R}(z,\bar{z})\,(z+\bar{z})^{L},$ (1.3)

where the density $\rho_{R}$ is determined entirely from Young diagram data in
a well-known way [22]. This is essentially the formula motivated in [21] from
holographic renormalization of low lying operators and Coulomb branch limits.
Our computation provides a check of this one-point function formula for all
single trace primaries and in principle for all LLM geometries without relying
on free-fermion methods. We also compute off-diagonal structure constants
between sufficiently close heavy states suggesting that semi-classical
supergravity calculations should be able to probe the precise microstructure
of bubbling geometries.

This paper is structured as follows. In section 2, we review the BPS coherent
state construction and and discuss the computation of the form factor of a
single trace operator in the background of a giant graviton. In section 3 we
generalize the computation to the case where the number of giant gravitons
scales with $N$. To do this we explain how to reduce the corresponding
integral over a Grassmannian to a more conventional matrix model involving
square matrices and then solve the model at large $N$. The resulting
distribution essentially reproduces the distributions studied in [23], up to a
change of variables. We then study the general problem of multiple stacks of
giant gravitons using steepest decent methods. In section 4 we return to the
problem of computing correlators in the character basis and provide a more
explicit connection from the eigenvalue picture presented by the coherent
state generating functions and the character basis. In section 5 we conclude
by discussing some general lessons and future directions.

##  2 Coherent States and Form Factors

Most of our discussion will concentrate on the simplest correlation functions
in the $\mathcal{N}=4$ SYM theory, which are three point functions of half-BPS
operators. We will also work mostly with the theory on the cylinder
$\mathbb{R}\times S^{3}$, but translating the results to the plane is
straightforward. A convenient parametrization for half-BPS operators is given
in terms of a six dimensional null complex vector $n\cdot n=0$:

$Z(x,n)=n_{I}\phi^{I}(x),$ (2.1)

and any half-BPS operator is obtained by taking gauge invariant combinations
of $Z(x,n)$. One common choice of operators are single and multi-trace
operators

$\mathcal{O}_{\left\\{L_{i}\right\\}}(x)=\prod_{k}\Tr_{N}\left[\left(\vec{n}\cdot\vec{\phi}(x)\right)^{L_{k}}\right].$
(2.2)

In the large $N$ limit with $\Delta=\sum_{k}L_{k}\ll\sqrt{N}$ these provide an
approximately orthogonal basis of operators; this is the usual statement that
in the large $N$ limit planar graphs contribute the most in correlation
functions. This class of operators is naturally associated to the supergravity
modes of $AdS_{5}\times S^{5}$ and their bound states. However, for operators
with large enough conformal dimension $\Delta$, various non-planar effects can
contribute meaningfully or eventually dominate over planar graphs. Even more
strikingly, certain extremal correlators of single traces operators have
enhanced contributions from non-planar diagrams even for small charges [24].
For these reasons it is useful to first perform the computation at finite $N$
with a proper orthogonal basis of states, and then take the large $N$ limit.

By restricting to primary operators, we will often drop the space coordinates
$x_{1,2,3}$ since we will mostly work with constant modes on the $S^{3}$. Due
to non-renormalization properties of half-BPS operators, the two and three
point functions of such operators can be computed in the free field theory
limit $g_{YM}=0$, so that our task reduces to a combinatorial problem of
performing Wick contractions of free fields. This problem was first addressed
in [8] for extremal correlators $n_{1}=n_{2}=n_{3}^{*}$. The main idea is to
construct an orthogonal basis of states for one matrix quantum mechanics with
$U(N)$ gauge symmetry, or equivalently a set of operators that diagonalize two
point functions. The resulting basis is build from characters of the unitary
group and is often referred to as the Schur basis:

$\mathcal{O}_{R}\left(Z(x,n)\right)=\frac{1}{k!}\sum_{\pi\in
S_{k}}\chi^{R}(\pi)\Tr_{\mathbf{N}^{\otimes k}}\left[\pi Z(x,n)^{\otimes
k}\right]=s_{R}\left(Z(x,n)\right),$ (2.3)

where $k$ denotes the number of boxes of the Young diagram associated to the
representation $R$ of $U(N)$ and $\chi^{R}(\pi)$ is the character of the
corresponding representation of $S_{k}$. A different proof that this set of
operators provides a diagonal basis was given in [4]. These operators have a
dual description in terms of giant gravitons, or their bound states, in
asymptotically $AdS_{5}\times S^{5}$ spaces [8, 14, 25].

Instead of performing the explicit contractions for a particular operator, it
was realized that one could instead work with a coherent state for the free
field $Z(x,n)$

$\ket{\Lambda}=\frac{1}{\texttt{Vol}\left[U(N)\right]}\int
dU\;e^{\Tr_{N}\left[U\Lambda U^{\dagger}Z(x,n)\right]}\ket{0}.$ (2.4)

A straightforward computation gives an alternate formula for $\ket{\Lambda}$
as a expansion in characters of the unitary group $s_{R}$,

$\ket{\Lambda}=\sum_{R}\frac{1}{f_{R}}\mathcal{O}_{R}\left(Z(x,n)\right)s_{R}(\Lambda),$
(2.5)

and $s_{R}(\Lambda)$ is a Schur polynomial. The point of this analysis is that
by exploiting the Campbell-Hausdorff formula the free field contractions of
the operators $\mathcal{O}_{R}\left(Z(x,n)\right)$ can all be replaced by an
integral over the unitary group. For the two point functions the resulting
integral is a Harish-Chandra-Itzykson-Zuber integral which has an exact fixed
point formula:

$\bra{\bar{\Lambda}}\ket{\Lambda}=\frac{1}{\texttt{Vol}\left[U(N)\right]}\int
dU\;e^{\Tr[U\Lambda U^{\dagger}\bar{\Lambda}]}=\mathcal{C}_{N}\sum_{\pi\in
S_{N}}\;\det(\pi)\,\frac{e^{\lambda_{i}\bar{\lambda}_{\pi(i)}}}{\Delta(\Lambda)\Delta(\bar{\Lambda})}.$
(2.6)

Following the ideas of [2, 3, 26], one can reduce the computation of any
correlator in the free theory to a matrix integral by commuting various
generating functions past each other using the Campbell-Hausdorff formula. In
the language of [2], this is equivalent to replacing the fields inside small
operators (such as traces) by their vevs after integrating out the SYM fields
and then performing a saddle point approximation over the auxiliary parameters
(in this case $U$).

More concretely we will be interested in computing form factors such as:

$\bra{\bar{\Lambda},n_{3}}\Tr[Z(t=0,n_{2})^{L}]\ket{\Lambda,n_{1}}\simeq\sum_{R,R^{\prime}}C_{R,R^{\prime}}(\Lambda,\bar{\Lambda})\bra{R^{\prime},n_{3}}\Tr[Z(t=0,n_{2})^{L}]\ket{R,n_{1}},$
(2.7)

where the initial and final states are created by heavy operators $\Delta\sim
N^{2}$ in the large $N$ limit. For relatively simple choices of operators,
such as determinants and traces of fully symmetric tensors [2, 3, 9, 11], the
saddle point analysis can be performed rather explicitly and the correlators
can be matched precisely to their holographic counterparts. For more
complicated operators, such as insertions of many determinant operators, or
operators associated to generic Young diagrams, the saddle point analysis
appears to be less straightforward and the structure of the solutions to the
saddle point equations is not fully understood. The main difficulty lies in
the fact that the resulting matrix models cannot be easily reduced to
integrals over eigenvalues, so that the saddle point equations appear to be
truly matrix equations. We will discuss in the later sections how to overcome
these complications in the regimes relevant to states with nice supergravity
descriptions (i.e. states corresponding to non-trivial geometries with small
curvatures).

###  2.1 Example: AdS giant graviton

Before proceeding to the case of interest, it is convenient to review the
results presented in [2, 11] since many of the parts of the calculations
presented there extend naturally. In the simplest of cases, the heavy
operators can be taken to be of rank one, meaning that they correspond to
Schur polynomials of fully symmetric or fully anti-symmetric representations.
In the case of fully anti-symmetric representations the correct generating
function is the determinant operator [27], for instance :

$\det(\phi_{5}+i\phi_{6}-\lambda)=\det(Z-\lambda),$ (2.8)

which describes a sphere giant graviton sitting at the origin of global AdS at
a position inside $S^{5}$ given by $e^{i\phi_{0}}\cos\theta_{0}=\lambda$. This
operator is a semi-coherent superposition of all subdeterminant operators,
each describing the R-charge eigenfunctions of a sphere giant graviton. The
method for semi-classical computation with this class of operators was
presented in [2] and also [26] so we refer the reader there for details.
Instead we will describe how the analogous computation is done for an operator
describing a semi-classical AdS giant graviton. The reason for this is that in
the end both calculations lead to very similar integrals for the correlators,
but their form is much easier to understand for the computation involving
symmetric tensors.

First we consider the following coherent state:

$\ket{\lambda}=\frac{1}{\texttt{Vol}\left[\mathbb{CP}^{N-1}\right]}\int_{\mathbb{CP}^{N-1}}d\varphi
d\varphi^{\dagger}\;e^{\lambda\varphi^{\dagger}Z\varphi}\ket{0}.$ (2.9)

As discussed in [4, 11] this is the same state that one obtains from setting
$\Lambda$ to be a rank one projector in (2.6). This state has a natural $U(1)$
gauge symmetry

$\varphi\sim e^{i\alpha}\varphi,$ (2.10)

which can be identified with the gauge symmetry on the worldvolume of the
giant graviton, as well as invariance under $U(N)$ gauge transformations of
$Z$. This state is also a coherent superposition of AdS giant graviton
wavefunctions with fixed R-charge. A simple calculation yields:

$\bra{\bar{\lambda}}\ket{\lambda}=\frac{1}{\texttt{Vol}\left[\mathbb{CP}^{N-1}\right]}\int_{\mathbb{CP}^{N-1}}d\varphi
d\varphi^{\dagger}\;e^{\lambda\bar{\lambda}\varphi^{\dagger}P_{1}\varphi}.$
(2.11)

To evaluate this integral we need to do a series of simple coordinate
transformations. Without loss of generality, we can let $P_{1}$ be a rank one
projector into the first component of $\varphi$. Then we can split the
coordinates of $\mathbb{CP}^{N-1}$ into $\varphi_{1}$ and $\varphi_{n}$ with
$n>1$. The reason we emphasize this will become clear when we generalize this
more complicated coherent states. Then, we can parametrize the coordinates
$\varphi_{n}$ in terms of an $2N-2$ dimensional spherical slices of radii
$R=\sqrt{1-|\varphi_{1}|^{2}}=\sqrt{1-r^{2}}$. Finally we can rewrite the
radial part of the integral as

$\int
d\left(r^{2}\right)d\left(R^{2}\right)\;R^{2N-4}\;\delta(R^{2}+r^{2}-1)\;e^{\lambda\bar{\lambda}r^{2}}=\int_{0}^{1}dx\;(1-x)^{N-2}\;e^{\lambda\bar{\lambda}x}.$
(2.12)

This last integral is simply the moment generating function of a particular
unitary Jacobi distribution. To make contact with the calculation involving
determinants, we can rewrite

$\left(\lambda\bar{\lambda}\right)^{N-2}(1-x)^{N-2}=\det\begin{pmatrix}\lambda&&\bar{\lambda}\varphi_{1}\\\
\lambda\varphi_{1}^{*}&&\bar{\lambda}\end{pmatrix}^{N-2}=\det\begin{pmatrix}\rho_{11}&&\rho_{12}\\\
\rho_{21}&&\rho_{22}\end{pmatrix}^{N-2}.$ (2.13)

Although this step is not necessary and the previous integral expression is
simple enough to evaluate explicitly, doing this change of variables makes it
clear that the final answer the large $N$ approximation for AdS giants is the
same as that of sphere giants (up to analytic continuation). After a final
simple re-scaling, $\lambda\rightarrow\sqrt{N-2}\;\lambda$ we finally arrive
at the expression

$\bra{\bar{\lambda}}\ket{\lambda}=C_{N}\int_{\det\rho\;\geq\lambda\bar{\lambda}-2}d\rho\;\det\rho^{N-2}\;e^{\left(N-2\right)\tr_{2}\left[\rho\rho^{\dagger}\right]},$
(2.14)

which apart from the contour on integration is identical to the integral
obtained from the Hubbard-Stratonovich trick used in [9] for determinant
operators 111Our trace convention is so that $\tr_{m}1=1$ as opposed to
$\tr_{m}1=m$ . This will be true generically, once we solve the saddle point
equations for a configuration of AdS giant gravitons, we automatically have
the solution for a configuration of sphere giant gravitons after a simple
analytic continuation. In this case the saddle point equations are simply:

$\displaystyle\rho^{\dagger}$
$\displaystyle=\rho^{-1}\Rightarrow\rho_{12}\rho_{21}=\lambda\bar{\lambda}-1$
(2.15) $\displaystyle\lambda\bar{\lambda}$ $\displaystyle>1.$

The second equation comes from the fact that the exponential needs to be
positive for the saddle point to be a maximum; this implies that the
semiclassical approximation is valid whenever $|\lambda|=\cosh\rho_{0}$ is
greater than one which simply says that the brane is at a position
$\rho_{0}>0$ in global AdS, and this is true for all half-BPS giant graviton
solution.

####  2.1.1 Form Factors

The next step is compute the following form factor;

$\bra{\bar{\lambda}}\Tr_{N}\left[\left(\vec{n}\cdot\vec{\phi}\right)^{L}\right]\ket{\lambda}=\bra{\bar{\lambda}}\Tr_{N}\left[\left(\frac{Z+\bar{Z}+Y-\bar{Y}}{2}\right)^{L}\right]\ket{\lambda}.$
(2.16)

Our choice of $\vec{n}$ is taken from [9] for clarity of presentation and is
arbitrary. To evaluate this quantity, we use the fact that the initial and
final states are coherent states which lets us replace $Z$ and $\bar{Z}$ by
constant matrices. The resulting trace is

$\Tr_{N}\left[\left(\frac{Z+\bar{Z}+Y-\bar{Y}}{2}\right)^{L}\right]=\Tr_{N}\left[\left(\frac{\bar{\lambda}\varphi\varphi^{\dagger}+\lambda\tilde{\varphi}\tilde{\varphi}^{\dagger}}{2}\right)^{L}\right].$
(2.17)

To proceed we use the fact that the integrals over $d\varphi$ and
$d\tilde{\varphi}$ are invariant under the action of $U(N)$, so we can gauge
fix $\tilde{\varphi}$ to be a unit vector $v$ with a one in the first
component. Finally one uses the trick introduced in [2] to exchange the trace
over color indices into a trace over “flavor” indices associated to the in and
out states of the brane:

$\Tr_{N}\left[\left(\frac{Z+\bar{Z}+Y-\bar{Y}}{2}\right)^{L}\right]=2^{-L}\tr_{2}\left[\begin{pmatrix}\lambda&&\bar{\lambda}\varphi_{1}\\\
\lambda\varphi_{1}^{*}&&\bar{\lambda}\end{pmatrix}^{L}\right]\simeq\tr_{2}\rho^{L}.$
(2.18)

Now instead of evaluating every power the matrix $\rho$, it is better to work
with the resolvent

$R(t)=\tr_{2}\left[\left(1-t\rho\right)^{-1}\right]=\frac{1-t\tr_{2}[\rho]}{t^{2}\det\rho-2t\tr_{2}[\rho]+1}.$
(2.19)

When we evaluate this expression at the saddle point value for $\rho$, we can
immediately recognize that the resolvent $R(t)$ is a generating function for
Chebyshev polynomials of the first kind.

$\left<R(t)\right>_{N=\infty}=\sum_{n=0}^{\infty}T_{n}(\cos\phi_{0}\cosh\rho_{0})\;t^{n},$
(2.20)

where we parametrized the eigenvalue in terms of LLM coordinates
$\lambda=e^{i\phi_{0}}\cosh\rho_{0}$.

####  2.1.2 Extracting Structure Constants

Naively one might expect that the saddle point approximation of the resolvent
computes a generating function of some half BPS structure constants. This is
not quite correct for the following reason. First we would need to extract the
contribution to the form factor from a particular set of primary operators. In
this case this is somewhat easy to do, given that the coherent state has a
simple expansion in terms of Schur polynomials for rank one representations

$\ket{\lambda}=n_{\lambda}\sum_{k=0}^{\infty}\frac{\lambda^{N+k-1}}{(N+k-1)!}\mathcal{O}_{(k)}(Z)\ket{0}.$
(2.21)

We should then think of the coefficients in the expansion as the distribution
of lengths for a Young diagram with a single row. The distribution is similar
to a Poisson random variable, so that the average length is of the Young
diagram is of order $|\lambda|$. Another way of seeing this is by using the
Stirling formula for the denominator and extremizing with respect to
$l=N+k-1$; the maximum occurs when $|\lambda|=N+k-1$. This is the reason why
the coherent state calculation in [11] gives the correct answer for structure
constants without the need to project into a particular character. However,
the dependence on the phase of $\lambda$ when we insert an operator will not
be correct due to unwanted contributions coming from off-diagonal terms. To
fix this one should project the intermediate operator into a an $R$-charge
singlet operator. This is done by performing a group average over the phase of
$\lambda$. Since the wavefunction $\bra{\bar{\lambda}}\ket{\lambda}$ is
already invariant under shifts in the phase of $\lambda$ whenever
$\bar{\lambda}=\lambda^{*}$, the only effect of averaging is to project out
off-diagonal terms from the resolvent.

Strictly speaking this averaging should be performed prior to doing the saddle
analysis, since averages do not generally commute. One way of performing this
average is to rescale $t\rightarrow\sqrt{\det\rho}\;t$, and then perform the
integration over the phase of $\lambda$ by a contour integral. The result is

$\overline{R(t)}=\sqrt{\det\rho}\left[\frac{t^{2}-1}{2\sqrt{(1+t^{2})^{2}\det\rho-4t^{2}\lambda\bar{\lambda}}}\right].$
(2.22)

In the large $N$ limit, we can set $\det\rho=1$, and the averaged resolvent
will take the form of a generating function of Legendre polynomials. Since the
saddle point analysis for this case basically involves setting $\rho$ to a
particular value, the averaging procedures commute so the large $N$ limit was
taken first in [9, 11] without any trouble. This is not the case for operators
made out Schur polynomials for large Young diagrams, since the large $N$ limit
leads to a continuous distributions of eigenvalues. Once the eigenvalues
condense, the result of the computation will be highly sensitive to the
analytic properties of the moment generating function, and performing the
averaging and large $N$ limit can lead to contradicting results. The most
natural prescription to remedy this is to perform the projection into a
particular primary operator first by an appropriate averaging, and then take
the large $N$ limit of this quantity.

##  3 Matrix models for general coherent states

###  3.1 Two Droplets

The next simplest calculation that we can perform is the case where the matrix
$\Lambda$ in (2.6) is taken to be a rank $p$ projector. In this case, the
analogous coherent state is an integral over the Grassmannian $Gr(p,N)$. The
expression for the coherent state is simple to write down, but some of the
steps needed to evaluate the resulting matrix integrals require some care; the
main task will be to evaluate the following norm:

$\ket{\lambda,p}=\frac{1}{\texttt{Vol}\left[Gr(p,N)\right]}\int_{Gr(p,N)}dVdV^{\dagger}\;e^{\lambda\Tr_{N}{VV^{\dagger}Z}}\ket{0}.$
(3.1)

We will argue that this state described the wavefunction of a stack of $p$
giant gravitons sitting at position $\lambda$ in the LLM plane. The first
thing to note is that this coherent state has an explicit $U(p)$ gauge
symmetry $V\sim Vg$ which we can identify as the gauge symmetry on a stack of
D-branes. The expectation value of $\Delta_{0}$ on this state is given by
$p|\lambda|^{2}$ so that whenever $|\lambda|\sim\sqrt{N}$ and $p\sim N$, the
average dimension of this state is of order $N^{2}$. By inspection, we can
also deduce that this state is a coherent superposition of Schur operators of
at most $p$ rows. by acting with $\Tr[\bar{Z}]$ on this state, we can see that
this state breaks the gauge symmetry spontaneously from $U(N)$ to
$U(N-p)\times U(p)$, and that the center of mass of the stack of $p$ branes is
at the position $z=\lambda$ on a complex plane.

Before proceeding we need to comment on the choice of coordinates for
Grassmannian, since the details on how to perform these types of integrals are
known but not readily available. We will mostly follow the notation of [28,
29]; for a pedagogical presentation we refer the reader to [13]. First, we can
choose to split any given group element $U$ in terms of block matrices

$U=\begin{pmatrix}A_{11}&&A_{12}\\\ A_{21}&&A_{22}\end{pmatrix},$ (3.2)

where $A_{11}$ and $A_{22}$ are $p\times p$ and $(N-p)\times(N-p)$ square
matrices, and $A_{12}$ is a $p\times(N-p)$ matrix (similarly for $A_{21}$).
Then we make an arbitrary choice of frame distinguished by a rectangular
matrix

$v^{T}=\begin{pmatrix}\;\mathbb{I}_{p}&&\boldsymbol{0}_{(N-p)\times
p};\end{pmatrix}$ (3.3)

this matrix can then be used to build projectors into arbitrary
$p$-dimensional subspaces of $\mathbb{C}^{N}$ by acting on $v$ with unitaries.
This set of projectors precisely gives a parametrization of the affine
Grassmannian $Gr(p,N)$. By an affine Grassmanian we will simply mean the space
spanned by the unitary transformations of $v$:

$Gr(p,N):=\\{V=U\cdot v\;|\;U\in G\\}.$ (3.4)

Clearly any $V$ in this space is rank deficient, so all of the information
about its singular value decomposition is captured a square matrix:

$VV^{\dagger}=U^{\dagger}P_{p}U=\begin{pmatrix}A_{11}^{\dagger}A_{11}&&0\\\
0&&0\par\end{pmatrix}.$ (3.5)

In analogy to the rank one calculation, we can write down the integration
measure for this space as

$dVdV^{\dagger}=dA_{11}dA_{11}^{\dagger}dA_{12}dA_{12}^{\dagger}\times\delta\left(A^{\dagger}_{11}A_{11}+A_{12}^{\dagger}A_{12}-\mathbb{I}_{p}\right)$
(3.6)

By a similar calculation we can see that the norm of this state is given by

$\bra{\bar{\lambda},p\;}\ket{\;\lambda,p}=\frac{1}{\texttt{Vol}\left[Gr(p,N)\right]}\int_{Gr(p,N)}dVdV^{\dagger}\;e^{\lambda\bar{\lambda}\Tr_{p}{V^{\dagger}P_{p}V}}.$
(3.7)

One last fact that we need before continuing is the Jacobian for the
coordinate transformation $M=A^{\dagger}A$ for any $n\times m$ rectangular
matrix $A$. This change of variables is well known in the context of Wishart
distributions;

$dAdA^{\dagger}\propto\det(M)^{n-m}dM,$ (3.8)

where the constant of proportionality is an integral over angular variables
that will not be important for our analysis. So in the end we find

$\bra{\bar{\lambda},p\;}\ket{\;\lambda,p}=n_{\lambda,p}\int_{A_{11}^{\dagger}A_{11}\preceq\;\mathbb{I}_{p}}dA_{11}^{\dagger}dA_{11}\det\left(\mathbb{I}_{p}-A_{11}^{\dagger}A_{11}\right)^{N-2p}\exp\left(\lambda\bar{\lambda}\Tr_{p}\left[A_{11}^{\dagger}A_{11}\right]\right).$
(3.9)

Since the only combination of $A_{11}$ and $A_{11}^{\dagger}$ that appears in
the integral is $A_{11}^{\dagger}A_{11}$, we can reduce the computation to an
integral over the eigenvalues of $A_{11}^{\dagger}A_{11}$.

$\bra{\bar{\lambda},p\;}\ket{\;\lambda,p}=\mathcal{C}_{\lambda,p}\int_{[0,1]^{p}}\prod_{i=1}^{p}dx_{i}\;\Delta\left(x_{i}\right)^{2}\left(1-x_{i}\right)^{N-2p}\;\exp\left(\lambda\bar{\lambda}\sum_{j=1}^{p}x_{j}\right).$
(3.10)

###  3.2 Large $N$ Limit: Steepest Descent

We will now sketch the evaluation of (3.10) in the large $N$ limit. As before,
$|\lambda|^{2}$ will scale with $N$ so that the exponential term is large. The
integral can evaluated explicitly using the Andreief identity [13]:

$\mathcal{C}_{\lambda,p}^{-1}=p!\det_{j,k}\left(\int_{0}^{1}dx\;(1-x)^{N-2p}x^{j+k-2}e^{\lambda\bar{\lambda}x}\right)$
(3.11)

The function inside the determinant is an incomplete Gamma function. Another
way interpreting it is at the moment generating function for the GUE Jacobi
ensemble. Even though we can evaluate simple moments in this distribution
exactly, the final form will always be a large determinant of a Hankel matrix
which we cannot deal with easily. Instead we can perform a saddle point
analysis of (3.10) for large $N$ and large $p$. After rescaling a similar
rescaling $\lambda\rightarrow\sqrt{N-2p}\;\lambda$ the saddle point equations
of this integral are of the form

$\lambda\bar{\lambda}-\frac{1}{1-x_{i}}-\frac{2}{N-2p}\sum_{j\neq
i}\frac{1}{x_{j}-x_{i}}=0.$ (3.12)

The behavior of the saddle point configurations are easy to understand; the
first two terms are the same as in the single eigenvalue problem so that all
the eigenvalues have a tendency to condense around
$\lambda\bar{\lambda}=\frac{1}{1-x_{i}}$ while second term the usual
eigenvalue repulsion term. We should also note that for distribution to be
stable we have to require that $p\leq N/2$. This makes sense, since condition
forces the second droplet to be smaller than the first. In the case that $p\gg
N/2$, this configuration is no longer a small deformation of the original
vacuum configuration and the remaining $N-2p$ eigenvalues become the relevant
degrees of freedom. The case with $p\sim N/2$ has to be treated with
particular care as we will see. To solve these equations at large $N$ with
$p/N$ fixed we simply recast the saddle point equations as a Ricatti equation
for the resolvent matrix $R_{p}(z)$:

$\left[\lambda\bar{\lambda}-\frac{1}{1-z}\right]R_{p}(z)-\frac{K}{1-z}-\frac{p}{N-2p}R_{p}(z)^{2}+\frac{1}{N-2p}R^{\prime}_{p}(z)=0,$
(3.13)

where $K=\frac{1}{p}\sum_{i=1}^{p}\frac{1}{1-x_{i}}=\int\frac{\rho(z)}{1-z}$
is a constant that is determined from $R_{p}(z)$ by imposing self-consistency
conditions of the distribution. From this we can identify $\frac{1}{N-2p}$ as
the relevant $\hbar$ parameter in the problem; this is the reason why $p\sim
N/2$ should be treated with care, since the $1/N$ fluctuations of the
resolvent are no longer under control and one must solve the differential
equation exactly. This can be seen from a simple scaling argument; if $N-2p$
is of order one, the exponential and determinant terms in (3.9) are of order
$e^{p}$, while the Vandermonde term is of order $e^{p^{2}}$. This says that
the dominant effect in this case is the eigenvalue repulsion, so that the
separation of each of the individual eigenvalues is large when compared to the
size of the system. In particular we should expect the corresponding geometry
to have a region with string scale curvature where the supergravity
approximation breaks down. This is expected, since wavefunction
$\ket{\lambda,p}$ no longer has a good semiclassical approximation in the
large $N$ limit. If instead we decide to keep $N-2p$ of order one but this
time scaling $\lambda$ as $\sqrt{p}\sim\sqrt{N}$, then the Vandermonde
contribution is off-set by the exponential term and the determinant term is
still sub-leading, so the distribution of eigenvalues $x_{i}$ is approximately
semi-circular. What this is saying is that now the two droplets are too close
together to be treated as separate (meaning that their size is of the same
order as their separation), and instead the deviation from vacuum is described
by a collection of order one giant gravitons probing the vacuum. In other
words, depending on how we decide to scale $N-2p$ we will obtain qualitatively
different saddle point conditions and the expansion in terms of AdS or sphere
giant gravitons might be more suitable.

The regime we will be interested in is when $N-2p$ is of order $N$, so that
the constant equilibrium configuration is a good approximation to the
eigenvalue distribution. In the large $N$ limit the density of states becomes

$\rho(x)dx=\frac{1}{\pi}\frac{\sqrt{4K\mu(1-x)-\left(1-(1-x)\lambda\bar{\lambda}\right)^{2}}}{1-x}dx,$
(3.14)

with $\mu=\frac{p}{N-2p}$; this turns out to be simpler after we make the
change of variables $z=\lambda\bar{\lambda}\left(1-x\right)$. After this we
can normalize the distribution such that $\int\rho(z)dz=2\mu$ and find $K=1$.
In this convention the eigenvalues are quantized in units of
$2g_{s}=\frac{1}{N-2p}$, where $g_{s}$ is the effective string coupling of
this system.

To compare with the corresponding LLM solution we can express the solution in
coordinates that manifest $\frac{1}{8}$ of the supersymmetries [23]. The point
is to write the 10d metric in terms of a 6d complex basis with coordinates
$x,y,z$. For the vacuum $AdS_{5}\times S^{5}$ solution these coordinates
should be identified with the coordinates of the five-sphere. Translating the
whole metric into these coordinates is a non-trivial task for generic LLM
geometries, but we will only be interested in determining the volume of the
cycle wrapped by the branes. In these coordinates the radius of the three-
sphere wrapped by the giant gravitons for a single droplet solution is given

$\tilde{r}=\frac{\sqrt{(L^{2}-|z|^{2})(|z-a|^{2}-b^{2})}}{|z-a|},$ (3.15)

where the droplet is centered at $z=a$, the size of AdS is $L$ and $b$ is the
radius of the droplet. Notice that this is essentially of the same form as
(3.14) up to some relabellings. The discrepancy between the denominators is
due to the fact that the variable $x$ is actually related to the square of the
radial direction of AdS. The precise mapping between both pictures should
involve some more complicated charge of variables in general, but the analytic
properties of both distributions are the same.

###  3.3 Three Point Functions: Diagonal Case

To compute the correlator of a single trace in the background of these
coherent states we can use the resolvent trick. Using the same kind of color-
flavor transformation the trace over the original color indices can be
replaced by a trace over $p\times p$ matrix. In this case the moment
generating function is

$\displaystyle\mathcal{F}(t)$
$\displaystyle=\tr_{p}\left[\left((2-t(\lambda+\bar{\lambda})\right)\left((1-t\lambda)(1-t\bar{\lambda})\mathbb{I}_{p}-t^{2}\lambda
A_{11}^{\dagger}A_{11}\right)^{-1}\right]$ (3.16)
$\displaystyle=2\sum_{i=1}^{p}\frac{1-\frac{t}{2}(\lambda+\bar{\lambda})}{t^{2}\lambda\bar{\lambda}(1-x_{i})-t(\lambda+\bar{\lambda})+1}.$

The most natural variable to work with is once again
$z_{i}=\lambda\bar{\lambda}(1-x_{i})$, which makes the density of eigenvalues
be of the form:

$\rho(z)=\frac{\sqrt{(z_{+}-z)(z-z_{-})}}{\pi z},$ (3.17)

where $z_{\pm}$ are given by the roots of the polynomial inside the square
root in (3.14)

$z_{\pm}=1+\frac{2\mu}{\lambda\bar{\lambda}}\pm\frac{2\sqrt{\mu(\mu+\lambda\bar{\lambda})}}{\lambda\bar{\lambda}}.$
(3.18)

Expanding $\mathcal{F}(t)$ as a function of $z$ gives an expression for the
moment generating function in terms of the moments of the Marchenko-Pastur
distribution (3.17). After extracting the $L^{\text{th}}$ moment we get

$\displaystyle\bra{\bar{\lambda}}\Tr[\frac{(Z+\bar{Z}+Y-\bar{Y})^{L}}{2^{L}}]\ket{\lambda}=$
(3.19)
$\displaystyle\sum_{k=0}^{\infty}\binom{\frac{L+k}{2}}{\frac{L-k}{2}}2^{-(L-k)}\left(\frac{\lambda+\bar{\lambda}}{2}\right)^{k}\,m_{\frac{L-k}{2}}-\sum_{k=0}^{\infty}2^{-(L-k)}\binom{\frac{L+k}{2}-1}{\frac{L-k}{2}}\left(\frac{\lambda+\bar{\lambda}}{2}\right)^{k+1}\,m_{\frac{L-k-1}{2}},$
$\displaystyle
m_{l}=\int_{z_{1}}^{z_{+}}dz\,\rho(z)\,z^{l}=\sum_{k=1}^{l}\frac{1}{k}\binom{l}{k}\binom{l}{k-1}\left(\frac{z_{+}-z_{-}}{2}\right)^{k}=\;_{2}F_{1}\left(1-l,-l;\,2;\
,\frac{z_{+}-z_{-}}{2}\right)$

To extract the diagonal part of this form factor we can average over the phase
of $\lambda$.

$\mathcal{F}_{\lambda\,\bar{\lambda}\,L}=\frac{1^{L}+(-1)^{L}}{2\sqrt{L}}\sum_{k=0}^{L}\binom{\frac{L+k}{2}}{\frac{L-k}{2}}\binom{k}{\frac{k}{2}}\left(\frac{1^{k}+(-1)^{k}}{2}\right)\,|\lambda|^{k}\,\;m_{\frac{L-k}{2}}\,\left(\frac{L}{L+k}\right).$
(3.20)

This correlator should be interpreted as encoding part of the angular
distribution of the bubbling geometry associated to the condensate of
eigenvalues and should compute the one-point function of a scalar operator on
an LLM geometry with two circular droplets. It would be interesting to compute
these correlators holographically, for instance with the methods developed in
[21]. By performing an additional contour integral over $|\lambda|$ we can
obtain three point functions for operators with fixed scaling dimensions as
opposed to coherent states as is done in [9]. One can similarly perform
computations that extract off-diagonal form-factors between heavy states with
different dimensions.

###  3.4 General Matrix Model: Eigenvalue Picture

Now we will proceed to the general case and sketch how to extract specific
operators associated to a particular Young diagram of at most $p$ rows with
order $N^{2}$ boxes, and we will use this to compute diagonal form factors. To
do this we need to consider the generating function

$\bra{\bar{\Lambda}_{p}}\ket{\Lambda_{p}}=n_{p}\int
d\sigma^{\dagger}d\sigma\det\left(\mathbb{I}-\sigma^{\dagger}\sigma\right)\;e^{N\tr_{p}[\Lambda\sigma\bar{\Lambda}\sigma^{\dagger}]}.$
(1.1)

This formula is the analog of (3.9) for generic coherent state parameters.
Computing this integral is not an easy task for generic eigenvalues, since the
argument of the exponential is no longer just a function of
$\sigma^{\dagger}\sigma$, and the matrix $\sigma$ is not a normal matrix, so
$\sigma$ and $\sigma^{\dagger}$ cannot be simultaneously diagonalized. In this
section we will give a heuristic argument for a saddle point approximation to
this integral. We will give a more concrete proof in the next section where we
will need a more careful analysis of these type of integrals. Intuitively one
would like to say that the integral is dominated by points where the all the
matrices in the trace are diagonal. One way to see why this could be true is
that the exponent can be written in the form:

$\tr_{p}\left[\Lambda\bar{\Lambda}\sigma\sigma^{\dagger}-\sigma\bar{\Lambda}[\Lambda,\sigma^{\dagger}]-\Lambda\bar{\Lambda}[\sigma,\sigma^{\dagger}]\right].$
(1.2)

Since the first term is manifestly positive, the exponent is the largest when
this term is maximized which happens whenever all the matrix components are
concentrated on the diagonals. Any deviation from this contributes to the
second and third terms, which are not necessarily positive. So it is natural
to expect that a good approximation to the integral is obtained by integrating
over the set of $\sigma$ satisfying

$[\lambda,\sigma^{\dagger}]=[\sigma,\sigma^{\dagger}]=0.$ (1.3)

Indeed we will see later that these are the saddle point conditions for the
integration over the angular variables for $\sigma$. Because of the large
exponent, corrections to this are heavily suppressed as long as
$\Lambda^{\dagger}=\bar{\Lambda}$, and so

$\bra{\bar{\Lambda}_{p}}\ket{\Lambda_{p}}\simeq\int_{[0,1]^{p}}dx_{i}\Delta_{p}(x_{i})^{2}\,(1-x_{i})^{N-2p}\,e^{N\sum_{i}|\lambda_{i}|^{2}x_{i}}+O(e^{-N^{2}}).$
(1.4)

The other saddle points contribution for the norms consists of pairing the
eigenvalues $\lambda_{i}$ with $\bar{\lambda}_{\pi(i)]}$ for all permutations
$\pi$, which are highly suppressed for non-coincident eigenvalues.

$e^{\Tr[\Lambda\sigma\bar{\Lambda}\sigma^{\dagger}]}\rightarrow\sum_{\pi\in
S_{p}}e^{\sum_{i}|s_{i}|^{2}\lambda_{i}\bar{\lambda}_{\pi(i)}}\times\left(\texttt{one-
loop}\right).$ (1.5)

We work out the appropriate one-loop determinant for this saddle point
approximation in section 4. This structure precisely explains the saddle point
structure found in [2] for determinant operators; the solutions to the matrix
form of the saddle point equations always involve summing over permutations of
initial and final giant graviton states. For widely separated eigenvalues the
Vandermonde determinant does not contribute meaningfully to the saddle point
approximation and the state is well approximated by a collection of widely
separated giant gravitons. A more interesting regime is whenever we have
$n_{a}$ coincident branes at a point $\lambda_{a}$, for $a=1,\dots,k$. As long
as the $\lambda_{a}$ are sufficiently separated interactions between different
droplets can be neglected and the eigenvalues $x_{i}$ are distributed along
$k$ cuts whose distribution is approximately given by the Marchenko-Pastur
distribution. More precisely, the norm of the coherent state with $k$ lumps of
eigenvalues centered around $\lambda_{a}$ is computed by the following matrix
model:

$\mathcal{Z}(\lambda_{1},\dots,\lambda_{k})=\int_{[0,1]^{p}}\prod_{a=1}^{k}\,\prod_{i_{a}}^{n_{a}}dx_{i_{a}}^{(a)}\,\Delta_{n_{a}}\left(x^{(a)}\right)^{2}\,\left(1-x_{i_{a}}^{(a)}\right)^{N-2p}e^{N\lambda_{a}\bar{\lambda}_{a}x_{a}}\times\prod_{c>b}\,\prod_{i_{b},j_{c}}\left(x^{(c)}_{j_{c}}-x^{(b)}_{i_{b}}\right).$
(1.6)

The derivation of this class of models was presented in [20], and we outline
the details in the next section. As long as the cuts are not exponentially
close to one another, the last eigenvalue repulsion term is far enough from
zero that it does not affect the saddle point. The analysis for the three
point function is then relatively straightforward; the moment generating
function can be block diagonalized and each block is dealt just like the
single-cut case. In the regimes when two cuts approach each other this
approximation is no longer valid and one has to solve the corresponding
monodromy problem exactly as the fluctuations around the stationary eigenvalue
distribution will not be suppressed. We expect that these corrections
reproduce the supergravity picture of [23], with the eigenvalue distribution
being related to the volume of the three-sphere on which the giant gravitons
are wrapped (see for example equations (5.110) and (5.116) in [23]). We expect
that the spectral curve for the matrix model is precisely encoded by the dual
LLM geometry written in the coordinates advocated by [23], since our
distribution of eigenvalue for the single cut case is essentially identical to
their equation (5.117). This aligns with the results coming from numerical
tests performed in [30] and we advocate for a similar viewpoint; $1/N$ effects
will generically give some amount of granularity to the edges of eigenvalue
droplets, specially if multiple droplets are close to one another when
compared to the characteristic size of each eigenvalue. In order to resolve
these details one would need to solve the full interacting saddle point
equations. We will not do this, since in this regime there will not be a
reliable geometric picture for the state. In other words, the saddle point
approximation we described above breaks down for states that describe half-BPS
geometries with string scale curvature, and instead we should think of the
state as being a deformation of a smooth geometry with some branes inserted.

###  3.5 Coulomb Branch Limit

One last interesting limit that we can consider is the limit in which the
droplets are widely separated from each other and from the origin. In this
limit, the Marchenko-Pastur distribution reduces to a delta function

$\rho(x)\rightarrow\sum_{i}\delta\left(\lambda_{i}\bar{\lambda}_{i}-\frac{1}{1-x_{i}}\right).$
(1.7)

This is exactly the Coulomb branch limit discussed in [21]. In this sense the
operators we study here can be understood as quantum mechanical analogs of
Coulomb branch vacua of the theory. This makes the relation between the
geometry of the moduli space of vacua and asymptotically anti-de-sitter spaces
clear; in the large $N$ limit, the moduli space should get quantum correction
which deform its geometry into a bubbling geometry and only in the dilute gas
approximation can we approximate such a geometry by a multi-center solution.

##  4 Matrix Models for the Character Basis

Now we will concern ourselves with computing three point functions where the
initial and final state are specified by specific Young diagrams, as opposed
to a collection of eigenvalues $\lambda_{i}$. The idea will be to make a
somewhat unconventional choice of integration contour for the coherent states
parameters. First we start with a pair of states
$\ket{\tilde{\Lambda}^{\dagger}},\,\ket{\Lambda}$, but now we treat the
parameters as being independent from one another; we will also force each of
the eigenvalues $\lambda_{i}$ and $\tilde{\lambda}_{i}$ to lie on a unit
circle. By multiplying $\ket{\Lambda}$ by the square of the Vandermonde
determinant of $\Lambda$, and integrating we can recognize that the resulting
integration measure is just a Haar measure for a new unitary matrix
$\mathcal{U}=U\Lambda U^{\dagger}$:

$\int_{U(N)}dU\oint
d\lambda_{i}\,\Delta(\lambda_{j})^{2}\rightarrow\int_{U(N)}d\,\mathcal{U}.$
(4.1)

The resulting state is clearly proportional to the vacuum state, since there
are no $\mathcal{U}^{\dagger}$ insertions to feed to the exponential. From
this it becomes clear that in order to extract a term proportional to the
state $\ket{R}=S_{R}(Z)\ket{0}$ one should multiply the integrand by a
character $S_{R}(\,\mathcal{U}^{\dagger})$:

$\ket{R}=\frac{f_{R}}{\texttt{Vol}\left[U(N)\right]}\int
d\,\mathcal{U}\;e^{\Tr[\,\mathcal{U}Z]}\;S_{R}(\,\mathcal{U}^{\dagger})\ket{0}.$
(4.2)

A similar trick was used in [31] to study expectation values of Wilson loop
for arbitrary representations in large $N$ Chern-Simons theory, with a
slightly different generating function. For our choice of generating function
the exponential factor can be expanded in terms of unitary characters using
Schur-Weyl duality and the resulting integrals are easily evaluated using
elementary orthogonality relations. On the other hand, for sufficiently large
representations we will be able to perform the integral using steepest descent
after we perform all contractions of the $\mathcal{N}=4$ SYM fields. This
exponential generating function is also useful for computing Wick contractions
with other operators since we can exploit the Campbell-Hausdorf formula. With
all this in mind the coefficient $f_{R}$ is easily determined to be

$f_{R}=n_{R}!\,\frac{\Tr_{R}[\,1\,]}{\chi_{R}(\text{id})}=n_{R}!\,\frac{\text{Dim}_{R}(N)}{d_{R}}\,;$
(4.3)

which is the norm of the corresponding state; this is done by expanding the
exponential and matching the terms as in [4].

We will want to compute quantities such as

$C_{RR^{\prime}L}=2^{-L}\times\frac{\bra{R^{\prime}}\Tr_{N}\left[\left(Z+\bar{Z}+Y-\bar{Y}\right)^{L}\right]\ket{R}}{\sqrt{L\bra{R^{\prime}}\ket{R^{\prime}}\bra{R}\ket{R}}},$
(4.4)

in the limit that $|R|\sim|R^{\prime}|\sim N^{2}$. To compute the quantity in
the numerator we can substitute the equation (4.2) and perform the Wick
contractions using the Campbell-Hausdorff formula:

$\displaystyle\bra{R^{\prime}}\Tr_{N}\left[\left(Z+\bar{Z}\right)^{L}\right]\ket{R}=\frac{f_{R^{\prime}}f_{R}}{\texttt{Vol}\left[U(N)\right]^{2}}\int
d\,\mathcal{U}d\,\mathcal{V}\;S_{R}(\mathcal{U}^{\dagger})S_{R^{\prime}}(\mathcal{V}^{\dagger})\;e^{N\Tr[\mathcal{U}\,\mathcal{V}]}\Tr\left[\left(\mathcal{U}+\mathcal{V}\right)^{L}\right].$
(4.5)

This procedure replaces all of the free-field Wick contractions with unitary
integrals which we can evaluate very explicitly.

###  4.1 Diagonal Structure Constant

As before it will be easier to work with the moment generating function for
the matrix $U+V$ instead of dealing with each individual trace. We now proceed
by diagonalizing both $U$ and $V$, after which we are left with an integral of
HCIZ-type:

$\mathcal{F}(t)\simeq\int
d\tilde{U}d\mu(u)d\mu(v)S_{R}(u^{*})S_{R}(v^{*})e^{N\Tr[\tilde{U}^{\dagger}u\,\tilde{U}v]}\Tr[\left(1-t(u+\tilde{U}v\tilde{U}^{\dagger})\right)^{-1}].$
(4.6)

This integral is quite challenging to evaluate exactly, mainly due to the
appearance of the unitaries $U$ inside of the trace of the resolvent. At large
$N$, the integral over $U$ can be evaluated the method of steepest descent;
the saddle point equations for the matrix $U$ are solved by permutation
matrices which allows us to replace the integral over $U$ by a sum over
permutations times a one loop determinant

$\int d\tilde{U}\rightarrow\frac{1}{N!}\sum_{\pi\in
S_{N}}\times\prod_{j}\frac{1}{\nu_{i}},$ (4.7)

where $\nu$ are the eigenvalues of the Hessian of $\Tr[U^{\dagger}u\,Uv]$;
this determinant factor is well known and it is proportional to
$\Delta(u)\Delta(v)$. Now, due to the permutation invariance of the measure of
integration for $u,v$, this sum over permutations can be performed by changing
variables $v_{i}\rightarrow v_{\pi(i)}$ in each of the terms in the sum, so
that in the end we are left with an integral over eigenvalues all lying inside
a unit circle:

$\mathcal{F}(t)\simeq\oint\prod_{i}\frac{du_{i}dv_{i}}{u_{i}v_{i}}\,e^{Nu_{i}v_{i}}\det\left(\bar{u}_{j}^{N+R_{k}-k}\right)\det\left(\bar{v}_{j}^{N+R_{k}-k}\right)\sum_{l=1}^{N}\frac{1}{1-t(u_{l}+v_{l})}$
(4.8)

Now the main obstacle is that we have a pair of determinants in the integrand.
To solve this issue we can expand the determinants as sums over permutations,
and exploit the symmetry of the measure under index relabelling to reduce the
number of sums. It is also more convenient to work with the variables
$x_{i}=u_{i}v_{i}$ and $y_{i}=u_{i}/v_{i}$;

$\displaystyle\mathcal{F}(t)$ $\displaystyle\simeq\frac{1}{N!}\sum_{\pi\in
S_{N}}(-1)^{\pi}\oint\prod_{i}\frac{dx_{i}}{x_{i}}\,e^{Nx_{i}}\,\sqrt{\bar{x}_{i}}^{N+R_{i}-i}\sqrt{\bar{x}_{\pi(i)}}^{N+R_{i}-i}$
(4.9)
$\displaystyle\times\oint\prod_{j}\frac{dy_{j}}{2y_{j}}\sqrt{y_{j}}^{R_{j}-j}\sqrt{\bar{y}_{j}}^{R_{\pi(j)}-\pi(j)}\,\sum_{i=1}^{N}\,\sum_{L=0}^{\infty}t^{L}x_{i}^{L/2}\left(2+y_{i}+\bar{y}_{i}\right)^{L/2}.$

To perform the integral over $y_{j}$ we set $y_{j}=e^{2i\beta_{j}}$, since the
coordinate $y_{i}$ winds around the unit circle twice. In order to get a non-
zero value, all of the integrals over $y_{j}$ should be non-zero and since the
moments only multiply by one particular value of $y_{j}$ for each term in the
sum we can conclude that the integral is only non-zero for $\pi=id$. After
evaluating the integral over $y$ we get

$\mathcal{F}(t)\simeq\frac{1}{N!}\oint\prod_{i}\frac{dx_{i}}{x_{i}}\,e^{Nx_{i}}\,\bar{x}_{i}^{N+R_{i}-i}\;\sum_{i}\frac{1}{\sqrt{1-4t^{2}x_{i}}}.$
(4.10)

This final integral is a simple Fourier integral that can be evaluated by
expanding in powers of $x_{i}$. After normalizing $\mathcal{F}(t)$
appropriately we obtain a formula for the generating function of structure
constant $C_{RRL}$:

$\mathcal{F}(t)=\sum_{L=0}^{\infty}t^{L}\,\frac{1^{L}+(-1)^{L}}{2\sqrt{L}}\times\left\\{\binom{L}{L/2}N^{-L/2}\sum_{i=1}^{N}\frac{\Gamma(N+R_{i}-i+1)}{\Gamma(N+R_{i}-i-L/2+1)}\right\\}.$
(4.11)

For large representations $R_{i}\sim N$ with large blocks, the ratio of gamma
functions can be replaced by the asymptotic expansion
$\frac{\Gamma(x)}{\Gamma(x-\beta)}\simeq x^{\beta}$ as $x\rightarrow\infty$,
and the sum may be replaced by an integral with $x=\left(N-i\right)/N$ and
$\alpha_{i}=R_{k+1-i}/N$:

$\mathcal{F}(t)=\sum_{L=0}^{\infty}t^{L}\,\frac{1^{L}+(-1)^{L}}{2\sqrt{L}}\times\left\\{\binom{L}{L/2}\sum_{i=1}^{k}\int_{\mu_{i}+\alpha_{i}}^{\mu_{i+1}+\alpha_{i}}dx\;x^{L/2}\right\\},$
(4.12)

here $\mu_{i}$ are filling fractions that measure the number of rows of size
greater than or equal to $R_{k+1-i}$ in units of $N$ and we assume that there
are $k$ non-zero blocks in the Young diagram for the representation $R$. Then
the structure constants are:

$\displaystyle C_{R\,R\,L/2}$
$\displaystyle=\frac{1^{L}+(-1)^{L}}{2\sqrt{L}}\times\left\\{\binom{L}{L/2}\int_{0}^{\infty}dx\;\rho_{R}(x)\times
x^{L/2}\right\\}$ (4.13) $\displaystyle=\frac{1}{\sqrt{L}}\int
dr\,r\,d\phi\,\rho_{R}(r)\,(z+\bar{z})^{L},\;\;\;z=re^{i\phi/2},$

where $x=r^{2}$ the density $\rho_{R}$ is defined as follows.

(a)

(b)

Figure 1: (a) A Maya diagram associated to a large Young diagram. The
rightmost positive slope edge is mapped to the Fermi sea, while negative slope
edges represent large gaps of unfilled states. (b) A sketch of the LLM
geometry associated to the Maya diagram.

Given the Young diagram we rotate it by $-3\pi/4$ radians. Then the diagram
has $(k+1,k)$ edges of slopes $\mp 1$ of lengths $\mu_{i}$ and $\kappa_{i}$
for negative and positive slopes respectively. For example the length of the
$l^{th}$ block is $R_{i}=N\sum_{i}^{k+1-i}\kappa_{i}=N\alpha_{i-k-1}$. We
color the edges with a negative slope black, and edges with positive slopes
white and unfold edges of the diagram into an infinitely long colored strip as
in 1(a). This strip is to be identified with a radial slice of the LLM plane
[22]. The variable $r$ is taken to start at the right-most colored edge; this
is $r=0$ in the LLM plane, and the asymptotic region is taken to be towards
the left of the strip. For every black region we have $\rho_{R}(r)=1$, and
$\rho_{R}(r)=0$ for every white region. Substituting this into (4.13) will
reproduce the integral expression in (4.12).

This quantity can be matched precisely to the formula found in [21] for the
one point point function of a chiral operator computed using holographic
renormalization:

$\langle\mathcal{O}_{S^{k\pm k}}\rangle_{LLM}=\frac{1}{\sqrt{k}}\int
dr\,rd\phi\;\rho(r)\,r^{k}e^{\pm ik\phi}.$ (4.14)

For our background this quantity vanishes since there is a conserved
$U(1)_{R}$ charge in the background. To match this to the expression above we
take the uncharged combination $\mathcal{O}_{S^{k,+k}}+\mathcal{O}_{S^{k-k}}$
and integrate over $\phi$. This would correspond to the contribution of a
spherical harmonic $Y^{L/2}\sim(z+\bar{z}+y-\bar{y})^{L}$.

###  4.2 Off-Diagonal Structure Constant

To compute off-diagonal structure constants we need to change one of the
representations $R$ to another representation $R^{\prime}$ whose Young diagram
is close to $R$. In this case most of the integrals will vanish, unless $R$
and $R^{\prime}$ only differ at a single row $R_{l}-R^{\prime}_{l}=k_{l}$.
Since $R$ has large blocks, this can only happen when $R_{l-1}>R_{l}$, meaning
after an edge. The only non-zero integral over the $y$ variables comes from
the $y_{l}$ associated to the row $R^{\prime}_{l}$, so we only get one term:

$\displaystyle C_{R\,R+k_{l}\,L}$ $\displaystyle=e^{-\Delta
S_{R\,R^{\prime}}}N^{-L/2}\frac{1^{L-k}+(-1)^{L-k}}{2\sqrt{L}}\times\frac{L!}{(\frac{L-k}{2})!(\frac{L+k}{2})!}\,\frac{\Gamma(N+R_{l}+k_{l}/2-l+1)}{\Gamma(N+R_{l}+k_{l}/2-L/2-l+1)}$
(4.15) $\displaystyle\simeq e^{-\Delta
S_{R\,R^{\prime}}}\;\frac{1^{L-k}+(-1)^{L-k}}{2\sqrt{L}}\times\frac{L!}{(\frac{L-k}{2})!(\frac{L+k}{2})!}\
\alpha_{l}^{L/2},$

where $\Delta S_{R\,R^{\prime}}$ is the ratio of norms
$\frac{\bra{R}\ket{R}}{\sqrt{\bra{R^{\prime}}\ket{R^{\prime}}\bra{R}\ket{R}}}$.
This should be interpreted as the linear response to a fluctuation localized
at the edge of a particular Fermi surface within the LLM geometry.

###  4.3 Comparing to the eigenvalue picture: fixing the number of rows

The method outline in this section allows us to compute expectation values of
light operators in a particularly radially symmetric bubbling geometry. A
natural question to address is how this connects to the eigenvalue coherent
state picture. In other words, given a particular configuration of droplets,
how can we determine which radially symmetry modes make up the state. Clearly
a single droplet made out of $p$ giant gravitons can only be made out of Young
diagrams with $p$ rows. This is because the overlap
$\bra{\tilde{\lambda}}\ket{\lambda}$ has a character expansion, and setting
$N-p$ of the eigenvalues to zero makes characters associated to
representations with more than $p$ rows vanish. To project to a particular
diagram made out of $p$ rows we repeat the same trick where we integrate the
remaining eigenvalues over a unit circle with an appropriate measure. After
regrouping the integration variables we will end up with a pair integrals over
$U(p)$:

$\displaystyle\bra{R^{\prime},p}\ket{R,p}$ $\displaystyle\propto\int
d\mu(s^{\dagger})d\mu(s)\,\det\left(\mathbb{I}-s^{\dagger}s\right)^{N-2p}\int_{U(p)}dU\int_{U(p)}dV\,S_{R}(U^{\dagger})S_{R^{\prime}}(V^{\dagger})e^{N\Tr[UsVs^{\dagger}]}$
(4.16) $\displaystyle\propto\frac{\delta_{RR^{\prime}}}{Dim_{N}(R)^{2}}\int
d\mu(s^{\dagger})d\mu(s)\,\det\left(\mathbb{I}-s^{\dagger}s\right)^{N-2p}S_{R}(s^{\dagger}s)$
$\displaystyle\propto\frac{\delta_{RR^{\prime}}}{Dim_{N}(R)^{2}}\int_{[0,1]^{p}}\prod_{i=1}^{p}dx_{i}\,\Delta_{p}(x)\,\left(1-x_{i}\right)^{N-2p}x_{i}^{p+R_{i}-i}.$

For generic (non-rectangular) diagrams, the integral does not have a simple
solution. For large diagrams we can use the saddle point approximation to find
the density of eigenvalues $x_{i}$. The procedure to evaluate this class of
integrals was outlined in [20] and also [32, 33, 19]. For a Young diagram with
$k$ large blocks with $n_{a}$ rows, we split the variables $x_{i}$ into $k$
groups of size $n_{a}$,
$x_{i}=(x_{1}^{(1)},\dots,x_{n_{1}}^{(1)},x_{(1)}^{2},\dots,x_{n_{2}}^{2},\dots,x_{1}^{k},\dots,x_{n_{k}}^{(k)})$.
Then we use the fact that the integration variables are invariant under
permutations to rewrite

$\displaystyle\int\prod_{a=1}^{k}\prod_{i_{a}=1}^{n_{a}}dx_{i_{a}}^{(a)}\left(x_{i_{a}}^{(a)}\right)^{R_{a}-N_{a}}\left(x_{i_{a}}^{(a)}\right)^{n_{a}-i_{a}}$
(4.17)
$\displaystyle\int\prod_{a=1}^{k}\prod_{i_{a}=1}^{n_{a}}\,dx_{i_{a}}^{(a)}\,\times\frac{1}{n_{a}!}\sum_{\pi\in
S_{n_{a}}}\left(x_{i_{a}}^{(a)}\right)^{R_{a}-N_{a}}\left(x_{\pi(i_{a})}^{(a)}\right)^{n_{a}-i_{a}}$
$\displaystyle\int\prod_{a=1}^{k}\prod_{i_{a}=1}^{n_{a}}dx_{i_{a}}^{(a)}\left(x_{i_{a}}^{(a)}\right)^{R_{a}-N_{a}}\Delta_{n_{a}}(x^{(a)}),$

where $N_{a}$ is the partial sum $\sum_{b\leq a}n_{b}$ and $N_{1}=0$. Putting
it all together, we are left with a multi-cut matrix model:

$\mathcal{Z}(R)=\int_{[0,1]^{p}}\prod_{a=1}^{k}\,\prod_{i_{a}}^{n_{a}}dx_{i_{a}}^{(a)}\,\Delta_{n_{a}}\left(x^{(a)}\right)^{2}\,\left(1-x_{i_{a}}^{(a)}\right)^{N-2p}\left(x_{i_{a}}^{(a)}\right)^{R_{a}-N_{a}}\times\prod_{c>b}\,\prod_{i_{b},j_{c}}\left(x^{(c)}_{j_{c}}-x^{(b)}_{i_{b}}\right)$
(4.18)

A matrix model quite similar to this one studied in [20] for computing vevs of
giant Wilson loops in Chern-Simons theory on Lens spaces
$S^{3}/\mathbb{Z}_{p}$. To make the analogy more precise we can change
variables to $x_{i}=e^{-u_{i}}$ after which the partition function becomes:

$\displaystyle\mathcal{Z}(R)$
$\displaystyle=\int\prod_{a=1}^{k}\prod_{i}du_{i}^{(a)}\prod_{j<i}\left(2\sinh\frac{u_{j}^{(a)}-u_{i}^{(a)}}{2}\right)^{2}\left(2\sinh\frac{u_{i}^{(a)}}{2}\,\right)^{N-2p}$
(4.19) $\displaystyle\times
e^{-(L_{a}+N/2-p)\,u_{i}^{(a)}}\prod_{b>a}\prod_{i,j}\left(2\sinh\frac{u_{j}^{(b)}-u_{i}^{(a)}}{2}\right).$
$\displaystyle L_{a}$
$\displaystyle=R_{a}-\frac{1}{2}\sum_{b=1}^{a-1}n_{b}+\frac{1}{2}\sum_{b>a}n_{b}$

The only difference between this and the matrix model studied in [20] is that
the Gaussian term is replaced by $(1-e^{-u_{i}^{(a)}})^{N-2p}$, and $u$
differs by a sign. In this representation the saddle point analysis is quite
straightforward. The equations of motion for the eigenvalues $u_{i}^{(a)}$
take the form:

$\frac{1}{N}\;\sum_{j\neq
i}\coth\frac{u_{j}^{(a)}-u_{i}^{(a)}}{2}+\frac{1}{2N}\sum_{b\neq
a,\,l}\coth\frac{u^{(b)}_{l}-u_{i}^{(a)}}{2}=\frac{(1-2p/N)}{2}\left(\coth\frac{u^{a}_{i}}{2}-1\right)-L_{a}/N.$
(4.20)

The resolvents for this type of problem can be taken to be of the form
$\omega^{(a)}(z)=\frac{1}{N}\sum_{i}\coth\frac{z-u_{i}^{(a)}}{2}$, and the
total resolvent is $\omega(z)=\sum_{b=1}^{k}\omega^{(b)}(z)$. At large $N$ the
eigenvalues condense into $k$ branch cuts and the saddle point equation on the
$a^{th}$ cut reads

$\omega^{(a)}(z+i\,0)+\omega^{(a)}(z-i\,0)+\sum_{b\neq
a}\omega^{(b)}(z)=\frac{(1-2p/N)}{2}\left(\coth\frac{u^{(a)}_{i}}{2}-1\right)-L_{a}/N,$
(4.21)

or equivalently

$\omega(z\pm
i\,0)=-\omega^{(a)}_{\mp}(z)+\frac{(1-2p/N)}{2}\left(\coth\frac{z}{2}-1\right)-L_{a}/N.$
(4.22)

This equation defines a Riemann-Hilbert problem for the total resolvent
$\omega(z)$. Notice that for large enough values of $z$, the potential rapidly
approches a constant value of $-L_{a}/N$, and for small values of $z$ there is
exponential barrier pushing the eigenvalues away from $z=0$. This means that
the eigenvalues will sit far from $z=0$, and will be uniformly distributed
along the cut. To solve (4.22) we need to find a function of the resolvents
that is regular for $\real(z)>0$. Following [20] we can define a set of
complex variables

$\displaystyle X_{0}$
$\displaystyle=e^{-(1-2p/N)\left(\coth\frac{z}{2}-1\right)/2}e^{\omega}=W\,e^{\omega}$
(4.23) $\displaystyle X_{a}$ $\displaystyle=e^{-L_{a}/N}e^{-\omega_{a}};$

then the equations (4.22) are equivalent to a monodromy condition for the
$X_{I}$ as we around each of the cuts. The solution to this problem is given
by a polynomial of degree $(k+1)$, $f(Y,W)$, with roots at the $X_{I}$ and the
spectral curve is the zero locus of this polynomial. The precise eigenvalue
distribution is then found by solving $f(Y,W)=0$ for $Y$ and taking the
discontinuity of the resolvent $\omega\propto\log Y$ at each cut.

In the thermodynamic limit, the branch cuts become widely separated and the
saddle point equations simply significantly. Substituting in the ansatz
$u^{(b)}_{j}-u^{(a)}_{i}\gg 1$ in (4.20) turn the terms of the form
$\coth\frac{u^{(b)}_{j}-u^{(a)}_{i}}{2}$ into a constant. The interaction term
between eigenvalues in the same cut becomes a step function for large
$u^{(a)}$ beginning at the first eigenvalue on the cut which given an equation
for each of the eigenvalues on a given cut:

$\displaystyle\frac{(1-2p/N)}{2}\left(\coth\frac{u^{a}_{k}}{2}-1\right)=\frac{1}{N}\left(L_{a}-(n^{(a)}-k)+\frac{1}{2}\sum_{b>a}n_{b}\right)\Rightarrow$
(4.24) $\displaystyle
2g_{s}\left(\frac{N-2p}{1-x_{k}}\right)=g_{s}\left[R_{a}-\sum_{b=1}^{a}n_{b}+\sum_{c>a}n_{c}+\left(2k+\sum_{b>a}n_{b}+R_{a}\right)+2\right].$

Clearly the eigenvalues $\lambda_{i}=1/(1-x_{i})$ are uniformly distributed
along each cut, and they are quantized in units of $2g_{s}=1/(N-2p)$. So again
we see that this is indeed the correct parameter controlling the fluctuations
around the saddle. These equations are also exactly analogous to the saddle
point equations for the single giant graviton case,

$\frac{1}{1-x_{k}}=|\lambda_{k}|^{2}=r_{k}^{2}$ (4.25)

and the distribution $\rho_{R}(r)$ for $r_{k}$ is the same distribution we
found before.

####  4.3.1 Computing correlators: reduced matrix elements

To compute correlation functions of single trace operators between a set of
states $\ket{R,p},\ket{R^{\prime},p}$ we need to perform a series of elaborate
coordinate substitutions to simplify the matrix integral calculations. First,
the states are obtained by a certain projection of the reduced coherent
generating function

$\displaystyle\ket{R,p}$ $\displaystyle\propto\int_{U(p)}dU_{1}\int
dV_{1}dV_{1}^{\dagger}\;e^{\tr_{p}[U_{1}V_{1}^{\dagger}ZV_{1}]}S_{R}(U_{1}^{\dagger})\ket{0}$
(4.26) $\displaystyle\bra{R^{\prime},p}$
$\displaystyle\propto\int_{U(p)}dU_{2}\int
dV_{2}dV_{2}^{\dagger}\;\bra{0}e^{\tr_{p}[U_{2}V_{2}^{\dagger}\bar{Z}V_{2}^{\dagger}]}S_{R^{\prime}}(U_{2}^{\dagger}).$

If we insert a half BPS single trace operator between these states, the scalar
fields $Z$ and $\bar{Z}$ inside the trace are traded by complex valued
matrices:

$\displaystyle Z\rightarrow V_{2}U_{2}V_{2}^{\dagger}$ (4.27)
$\displaystyle\bar{Z}\rightarrow V_{1}U_{1}V_{1}^{\dagger}.$

We can then perform a color-flavor transformation to express any trace of a
power of $Z+\bar{Z}$ (over $N$ color indices) into a trace over $2p\times 2p$
matrices. Then after performing all of the contractions between the two
exponentials we end up with a pair of unitary integrals over $U_{1,2}$ and an
integral over the $p\times p$ ‘radial’ matrix $\sigma=V_{1}^{\dagger}V_{2}$.

$\displaystyle\bra{R^{\prime},p}\Tr\left[\left(Z+\bar{Z}+Y-\bar{Y}\right)^{L}\right]\ket{R,p}=$
(4.28) $\displaystyle\int
dU_{1}dU_{2}\,d\sigma^{\dagger}d\sigma\Bigg{\\{}\det\left(1-\sigma^{\dagger}\sigma\right)^{N-2p}e^{\tr_{p}\left[U_{1}\sigma
U_{2}\sigma^{\dagger}\right]}\,S_{R}(U_{1}^{\dagger})S_{R}(U_{2}^{\dagger})$
$\displaystyle\times\tr_{2p}\left[\begin{pmatrix}U_{1}&\sqrt{U_{1}}\,\sigma\sqrt{U_{2}}\\\
\sqrt{U_{2}}\,\sigma^{\dagger}\sqrt{U_{1}}&U_{2}\end{pmatrix}^{L}\right]\Bigg{\\}},$

and here $\sqrt{U_{1,2}}$ refers to the unitary matrix whose eigenvalues are
square roots of the eigenvalues of $U_{1,2}$. To solve this integral at large
$N$ we will diagonalize $U_{1,2}$ and perform a singular value decomposition
of $\sigma$.

$\displaystyle\sigma$ $\displaystyle=U_{L}\,s\,U_{R}$ (4.29) $\displaystyle
U_{1}$
$\displaystyle=\mathcal{U}^{\dagger}e^{i\boldsymbol{\alpha}}\,\mathcal{U}$
$\displaystyle U_{2}$
$\displaystyle=\mathcal{V}^{\dagger}e^{i\boldsymbol{\beta}}\,\mathcal{V}.$

Although this at first glance looks daunting one should realize that the
translational invariance of the Haar measure allows us to reabsorb most of the
unitary integrals into the integration over $U_{L}$ and $U_{R}$. This reduces
the calculation to a pair of integrals involving simple complex exponential
associated to the eigenvalues of $U_{1,2}$ which encode the angular dependence
of the correlator, one integral over the radial eigenvalues, and a pair of
difficult unitary integrals over $U_{L,R}$. First we address the integration
over $U_{L}$ and $U_{R}$. At large $N$ we can perform a saddle point
approximation for this integral; since the only term that is relevant for this
is the exponential we are left with the task of finding the critical points of
the following function:

$S(U_{L},U_{R})=\tr_{p}\left[U_{L}^{\dagger}e^{i\boldsymbol{\alpha}}\,U_{L}\,s\,U_{R}\,e^{i\boldsymbol{\beta}}U_{R}^{\dagger}\,s^{\dagger}\right].$
(4.30)

The critical points of this function are given by the solutions to a pair of
matrix equations

$[s^{\dagger}U_{L}^{\dagger}e^{i\boldsymbol{\alpha}}\,U_{L}\,s\,,U_{R}\,e^{i\boldsymbol{\beta}}U_{R}^{\dagger}]=[s\,U_{R}\,e^{i\boldsymbol{\beta}}U_{R}^{\dagger}\,s^{\dagger},U_{L}^{\dagger}e^{i\boldsymbol{\alpha}}\,U_{L}]=0.$
(4.31)

These equations essentially imply that these two pairs of matrices are
simultaneously diagonalizable. For example, the first equation is unitarily
equivalent to the condition that $e^{i\boldsymbol{\beta}}$ and
$U_{R}^{\dagger}s^{\dagger}U_{L}^{\dagger}e^{i\boldsymbol{\alpha}}\,U_{L}\,s\,U_{R}$
are both diagonal in the same basis. The second equation gives a similar
condition for $e^{i\boldsymbol{\alpha}}$ and
$U_{L}\,s\,U_{R}\,e^{i\boldsymbol{\beta}}U_{R}^{\dagger}\,s^{\dagger}\,U_{L}^{\dagger}$.
But since $e^{i\boldsymbol{\alpha}}$, $e^{i\boldsymbol{\beta}}$ and $s$ are
all diagonal the only way that this can be achieved is if $U_{L,R}$ are
permutation matrices. This is clear because making the ansatz $U_{R}=U_{\pi}$
in the first equation for some permutation matrix $U_{\pi}$ immediately forces
$U_{L}$ to be a permutation matrix and vice versa. This means that we have one
critical point for every pair of permutations $\pi,\tau$ in $S_{p}$, which act
on $e^{i\boldsymbol{\alpha}}$ and $e^{i\boldsymbol{\beta}}$ by permuting their
eigenvalues independently from each other. To get the correct answer we should
also compute the one-loop determinant around each of the saddles. This boils
down to computing the Hessian of (4.30), which is the same for each each
critical point up to a sign associated to the determinant of the permutation.
To quadratic order $S(e^{M},e^{N})$ is

$\displaystyle S(m,n)$
$\displaystyle\simeq\tr_{p}[e^{i\boldsymbol{\alpha}}\,s\,e^{i\boldsymbol{\beta}}\,s^{\dagger}]-\frac{1}{2}\sum_{i,j}\Bigg{\\{}(e^{i\alpha_{i}}-e^{i\alpha_{j}})(s_{i}s_{i}^{*}\,e^{i\beta_{i}}-s_{j}s_{j}^{*}e^{i\beta_{j}})|m_{ij}|^{2}$
(4.32)
$\displaystyle+(e^{i\beta_{i}}-e^{i\beta_{j}})(s_{i}s_{i}^{*}\,e^{i\alpha_{i}}-s_{j}s_{j}^{*}e^{i\alpha_{j}})|n_{ij}|^{2}$
$\displaystyle-(e^{i\alpha_{i}}-e^{i\alpha_{j}})(e^{i\beta_{i}}-e^{i\beta_{j}})s_{i}s_{j}^{*}\left(m_{ij}n_{ij}^{*}+n_{ij}m_{ij}^{*}\right)\Bigg{\\}}+\dots,$

so the one-loop factor becomes simply

$\det\left(\partial_{ij}\partial_{kl}S(m,n)\right)=\Delta_{p}\left(e^{i\boldsymbol{\alpha}}\right)\Delta_{p}\left(e^{i\boldsymbol{\beta}}\right)\Delta_{p}\left(s^{\dagger}s\right)\Delta_{p}\left(e^{i(\boldsymbol{\alpha}+\boldsymbol{\beta})}s^{\dagger}s\right).$
(4.33)

Each saddle point is weighted by $\det(\tau\pi)$, so after a coordinate
transformation every term in the sum will give the same value. As before, the
one-loop factor when combined with the denominators of the determinantal
expressions for the Schur polynomials will cancel the Vandermonde determinants
in the integration measure of $\alpha_{i}$ and $\beta_{i}$, which makes the
integrals over the angular variables straightforward

$\displaystyle\bra{R^{\prime},p}\Tr\left[\left(Z+\bar{Z}+Y-\bar{Y}\right)^{L}\right]\ket{R,p}=$
(4.34)
$\displaystyle\int\frac{d\alpha_{i}}{2\pi}\frac{d\beta_{i}}{2\pi}\int_{[-1,1]^{p}}ds_{i}\,\Bigg{\\{}\frac{\Delta_{p}(s^{\dagger}s)}{\Delta_{p}(e^{i(\boldsymbol{\alpha}+\boldsymbol{\beta})}s^{\dagger}s)}\left(1-|s_{i}|^{2}\right)^{N-2p}\exp\left(e^{i(\alpha_{i}+\beta_{i})}|s_{i}|^{2}\right)$
$\displaystyle\det_{l,k}\left(e^{-i(p+R_{l}-l)\alpha_{k}}\right)\det_{l,k}\left(e^{-i(p+R^{\prime}_{l}-l)\beta_{k}}\right)\times\binom{L}{L/2}\tr_{p}\left[(s^{\dagger}s)^{L/2}\,e^{i(\frac{\boldsymbol{\alpha}+\boldsymbol{\beta}}{2})L}\right]\Bigg{\\}}.$

For light operators and diagonal structure constants we can simply interchange
$(Z+\bar{Z})^{L}$ with $\binom{L}{L/2}(Z\bar{Z})^{L/2}$, which is simpler to
work with. The integration over $\alpha_{i}-\beta_{i}$ will again force the
pair of sums over permutations coming from the two determinants to collapse to
a single sum, and the integral over $\alpha_{i}+\beta_{i}$ cancels the factor
of $(s_{k}s_{k}^{*})^{L/2}$. Finally the integral over $s_{i}$ factors out
completely and the correlator comes only from the integrals over the angular
variables:

$\displaystyle\bra{R,p}\Tr\left[\left(Z+\bar{Z}+Y-\bar{Y}\right)^{L}\right]\ket{R,p}$
(4.35)
$\displaystyle=\mathcal{Z}(R)\sum_{i}\frac{1^{L}+(-1)^{L}}{2}\times\left\\{\binom{L}{L/2}N^{-L/2}\sum_{i=1}^{p}\frac{\Gamma(p+R_{i}-i+1)}{\Gamma(p+R_{i}-i-L/2+1)}\right\\}.$

This is the same answer as (4.13) except that we are missing one term coming
from the droplet made out of $N-p$ spectator branes. The reason that we miss
this term in this calculation is that the reduced generating functions project
out all Young diagrams with more than $p$ rows, which essentially freezes
$N-p$ of the branes that make up the background. This agrees with our
interpretation of the off-diagonal structure constant as the linear response
of a particular Fermi sea level. Since we projected out all the contributions
from Young diagrams with more than $p$ rows, ripples of the first Fermi sea
surface are projected out from this computation.

##  5 Discussion and Future Directions

In this paper we introduced techniques for dealing with large BPS operators in
$\mathcal{N}=4$ SYM theory and revisited the computation of one-point
functions in the background of states corresponding to the bubbling geometries
of [34]. Our method is based on the semi-classical techniques introduced in
[4, 2, 3] and further developed in [11, 35, 36, 10, 37] and provides an
independent derivation of the results of [15]. One advantage of our methods is
that they do not rely on diagonalizing any the field operators of the theory,
or performing any kind of consistent truncation of the model, which makes weak
coupling computation of non-BPS observables possible. We also fleshed out the
relation between coherent states and characters by providing an explicit
integral transform between both pictures. This gives somewhat complementary
methods that can be used to compute correlators of operators describing
somewhat generic LLM backgrounds. Although our results are expected and
perhaps unsurprising, we hope that the techniques developed here serve as a
starting point for performing a systematic large $N$ expansion for large
operators.

One question worth asking is whether we learned anything about the statistics
of (BPS) OPE coefficients in $\mathcal{N}=4$ SYM. We are hesitant to claim
that our results display any kind of chaotic behavior predicted by the
eigenstate thermalization hypothesis [38], since half-BPS operators in
$\mathcal{N}=4$ SYM are completely captured by a free theory. We instead
attribute the appearance of random matrix statistics to the averaging
necessary to describe large semi-classical states, and to effects due to the
large charge limit. In fact our computations suggest that true structure
constants (as opposed to fixed charged three point functions) only receive
contributions from a single term, and that the would be distribution of
eigenvalues in these cases are essentially constant distributions. On the
other hand, correlators involving large operators that break the $R$-symmetry
spontaneously but with fixed charge are generically given by complicated
averages. In the large $N$ limit with the charges of the operators scaling
with $N^{2}$ this leads to the appearance of random matrix behavior in the OPE
coefficients. Interestingly similar distributions were found to emerge from
the large charge limit of extremal correlators in rank one SCFTs [39].

We conclude by commenting on some future directions of work.

### Holographic computation of off-diagonal three-point functions

Our methods allow us to make predictions for the value of the off-diagonal
structure constant between two LLM geometries. This quantity is only non-zero
when the two geometries differ by a small fluctuations. It would be ideal to
try develop semi-classical techniques for computing such things using the
gravitational path integral. Since states are highly degenerate (there is one
for each Young diagram!), we expect that there is a set of commuting charges
that differentiates between different states in gravity. These charges would
correspond to some kind of higher spin asymptotic symmetries of LLM solutions.
Then, computing off-diagonal three-point functions would correspond to
dressing the semi-classical saddle by a wavefunction charged with the
appropriate asymptotic charge as suggested by [9]. This would be a very simple
toy model of the soft hair proposal [40] in the sense that one is able to
probe very precise details of the interior of the geometry from simple
boundary manipulations. Such a technique would also have to go far beyond our
current methods of holographic renormalization [41]. We take the fact that the
formulas for one-point functions in LLM backgrounds are quite simple as
indicative of the existence of a different method for computing holographic
one point functions in them. Perhaps a careful WKB analysis [42] might be able
to reproduce the result for operators of arbitrary charge without the need for
a non-linear Kaluza-Klein reduction.

### Three Heavy Operators

One obvious extension of our calculations would be to consider the structure
constants of three really heavy operators half-BPS operators. There are in
principle no conceptual obstructions for performing such computations, since
we only have to include an additional exponential generating function.
Technically speaking the saddle point analysis is more involved and we expect
the form of the structure constants to be much richer. More precisely many of
the simplifications that occur for the one-point functions of single trace
operators are basically due to the fact that for small operators the term in
the tensor decomposition only show up with multiplicity one. This is why we
are able to reduce the number of sums in (LABEL:reduced_correlator). For
example one can insert a particular Schur polynomial of $Z+\bar{Z}+Y-\bar{Y}$
between two coherent states. In this case we can commute the coherent states
past the Schur polynomial resulting in a matrix model generalizing (3.9).
Obtaining an approximate formula even for extremal three point functions of
very large operators would be a rather non-trivial since would encode
information about the statistics of Littlewood-Richardson coefficients. On the
other hand, such correlators predict the existence of supergravity solutions
which interpolate between different LLM geometries [2]. Since these types of
three point functions are protected, we expect that one should be able to
match both results precisely in the large $N$ limit, so obtaining an
approximate form for such correlators might give some intuition about how to
construct such geometries. Perhaps a simpler problem is to consider
correlators of two LLM geometries and a giant graviton, or between three giant
gravitons. These quantities give predictions for giant graviton nucleation
amplitudes.

### Extremal Correlators

Another immediate generalization would be to study higher point extremal
correlators involving various combinations of single trace, giant graviton,
and LLM geometries. These correlators are also protected, so we can hope to be
able to perform a holographic check for these quantities as well. In practice
this will likely involve a very careful group averaging procedure to be able
to overcome any possible ambiguities that often arise in extremal correlators.
For instance one can conceivably study correlators in which many branes
nucleate into a bubbling geometry, or where large droplets all fuse into a
single droplet, or a large droplet splits into smaller ones. Such processes
are somewhat reminiscent of baby universe creation and annihilation processes.
Whether such a geometric picture can be realized from the bulk point of view
is unclear.

### Worldsheet and spin chain interpretation of one-point functions in LLM
background

One particularly puzzling issue is to interpret the result of the computation
of a three point function of a non-protected single trace primary and two very
large half-BPS operators. For large operators of dimension of order $N$ the
interpretation as a worldsheet $g$-function was advocated in [2], and various
checks were performed. This makes sense since operators of dimension of order
$N$ can change the boundary conditions of the string worldsheet in the bulk.
For operators of dimension of order $N^{2}$ this interpretation is inadequate,
since we expect that the correct description is instead in terms of a string
moving in an LLM background, rather than simply an open string ending on a
stack of branes. From the worldsheet point of view one would expect that the
worldsheet CFT flows as the background is deformed. Understanding what this
would mean from the point of view of the spin chain picture would be quite
interesting.

###  $\frac{1}{4}$, $\frac{1}{8}$, and $\frac{1}{16}$ BPS operators

A more long term goal would be to understand the systematics of less
supersymmetric BPS operators in the $\mathcal{N}=4$ SYM theory. In principle
introducing additional scalar matrices in the exponentials gives a way of
generating $\frac{1}{4}$ and $\frac{1}{8}$ BPS [4, 3, 43]. This was carried
out for rank two $\frac{1}{4}$ BPS operators in [43]. Additionally, the
generating functions introduced in [4] contain all bosonic $\frac{1}{8}$
operators. The main difficulty lies in constructing explicit expressions for
restricted Schur polynomials with which one can project to particular BPS
operators of the weakly coupled theory. In other words, finding the
integration rules for the analog of the Schur polynomials in (4.2) with
multiple matrices would be a big step in this program. This problem is even
more stark for $\frac{1}{16}$ BPS operators where one seems needs to need to
introduce an infinite number of matrices corresponding to covariant
derivatives of the scalar fields. In analogy to the construction in [4], one
can generate $\frac{1}{16}$ BPS states by exponentiating the so-called
$\frac{1}{16}$ BPS letter introduced in [44]. Understanding precisely what
kind of excitations are relevant for studying supersymmetric black holes would
be a first step towards a boundary derivation of the results of [45].

### Twisted Holography

A natural setting where our techniques can be readily applied is in the
context of the Twisted Holography program. For instance sphere giant gravitons
and non-conformal vacua have already been studied [26, 46] and their geometric
picture is quite clear. Extending this to include the analog of AdS giant
gravitons and more generic Schur polynomial operators seems like a
straightforward task. It would be nice to develop the more conventional view
of giant gravitons wrapping compact cycles on the deformed conifold
$SL(2,\mathbb{C})$ in the B-model by developing a global version of Twisted
Holography, in analogy to global AdS holography. This should be related to
studying the chiral algebra on $\mathbb{P}^{1}$ instead of $\mathbb{C}$. The
point of this exercise is that when we insert an operator associated to a
Young diagram $R$, there is a well-studied recipe for producing a spectral
curve, and hence a bubbling Calabi-Yau manifold [20]. This story is well-
understood for the $A$-model on the deformed conifold, so it is not
implausible that a similar story exists for the $B$-model. This might help
bridge the conceptual gap about the relation between the full physical
holography and the twisted version.

### Semiclassics and Large $N$

Inevitably one will need to include $g^{2}N$ corrections in a systematic way
when dealing with non-protected operators. A good place to start would be
consider simpler models where the techniques developed in [47, 48, 49] can be
combined with the large $N$ expansion. Performing a near BPS expansion seems
like a natural starting point since one expects the coupling constant to be
enhanced by additional kinematic effects [50] making a reliable extrapolation
to strong coupling a possibility for certain observables. A good target would
to understand the correlation functions of more complicated ‘baryonic’
operators in the Wilson-Fisher fixed point of the $O(N)$ model [51] in the
large $N$ limit;

$\mathcal{B}\sim S_{i_{1},\dots,i_{N}}\,\phi^{i_{1}}\dots\phi^{i_{N}}$ (5.1)

where $S$ is a symmetric tensor. The simplest example of these operators where
studied in [51], although it would be interesting to extend these kinds of
results to operators associated to other representations of $O(N)$. These kind
of operators can also be exponentiated with the help of real gaussian
integrals over Grassmannians, so the main difficulty would be to study the gap
equations in the presence of these baryons.

## Acknowledgments

We would like to thank David Berenstein, Sanjaye Ramgoolam, and Antal Jevicki
for discussions. We are pleased to thank Kwinten Fransen, Shota Komatsu,
Robert de Mello Koch, and Wayne Weng for comments on an earlier draft. A.H. is
supported in part by funds from the University of California.

## References

  * [1] Edward Witten.  Baryons in the 1/n Expansion.  Nucl. Phys. B, 160:57–115, 1979.
  * [2] Yunfeng Jiang, Shota Komatsu, and Edoardo Vescovi.  Structure constants in $\mathcal{N}$ = 4 SYM at finite coupling as worldsheet g-function.  JHEP, 07(07):037, 2020.
  * [3] Gaoli Chen, Robert de Mello Koch, Minkyoo Kim, and Hendrik J. R. Van Zyl.  Absorption of closed strings by giant gravitons.  JHEP, 10:133, 2019.
  * [4] David Berenstein and Shannon Wang.  BPS coherent states and localization.  JHEP, 08:164, 2022.
  * [5] Jean-Loup Gervais and B. Sakita.  Extended Particles in Quantum Field Theories.  Phys. Rev. D, 11:2943, 1975.
  * [6] Jean-Loup Gervais, A. Jevicki, and B. Sakita.  Collective Coordinate Method for Quantization of Extended Systems.  Phys. Rept., 23:281–293, 1976.
  * [7] Juan Maldacena.  The large-n limit of superconformal field theories and supergravity.  International journal of theoretical physics, 38(4):1113–1133, 1999\.
  * [8] Steve Corley, Antal Jevicki, and Sanjaye Ramgoolam.  Exact correlators of giant gravitons from dual N=4 SYM theory.  Adv. Theor. Math. Phys., 5:809–839, 2002.
  * [9] Peihe Yang, Yunfeng Jiang, Shota Komatsu, and Jun-Bao Wu.  D-branes and orbit average.  SciPost Phys., 12(2):055, 2022.
  * [10] Hai Lin.  Coherent state operators, giant gravitons, and gauge-gravity correspondence.  Annals Phys., 451:169248, 2023.
  * [11] Adolfo Holguin and Wayne W. Weng.  Orbit Averaging Coherent States: Holographic Three-Point Functions of AdS Giant Gravitons.  11 2022.
  * [12] Warren Carlson, Robert de Mello Koch, and Minkyoo Kim.  Generating functions for giant graviton bound states.  Journal of High Energy Physics, 2023(1), jan 2023.
  * [13] Giacomo Livan, Marcel Novaes, and Pierpaolo Vivo.  Introduction to Random Matrices.  Springer International Publishing, 2018.
  * [14] David Berenstein.  A toy model for the ads/cft correspondence.  Journal of High Energy Physics, 2004(07):018, 2004.
  * [15] Yastoshi Takayama and Asato Tsuchiya.  Complex matrix model and fermion phase space for bubbling AdS geometries.  JHEP, 10:004, 2005.
  * [16] Nadav Drukker and David J. Gross.  An Exact prediction of N=4 SUSYM theory for string theory.  J. Math. Phys., 42:2896–2914, 2001.
  * [17] Vasily Pestun.  Localization of gauge theory on a four-sphere and supersymmetric Wilson loops.  Commun. Math. Phys., 313:71–129, 2012.
  * [18] Nadav Drukker, Simone Giombi, Riccardo Ricci, and Diego Trancanelli.  On the D3-brane description of some 1/4 BPS Wilson loops.  JHEP, 04:008, 2007.
  * [19] Jaume Gomis, Shunji Matsuura, Takuya Okuda, and Diego Trancanelli.  Wilson loop correlators at strong coupling: From matrices to bubbling geometries.  JHEP, 08:068, 2008.
  * [20] Nick Halmagyi and Takuya Okuda.  Bubbling Calabi-Yau geometry from matrix models.  JHEP, 03:028, 2008.
  * [21] Kostas Skenderis and Marika Taylor.  Anatomy of bubbling solutions.  JHEP, 09:019, 2007.
  * [22] Robert de Mello Koch.  Geometries from Young Diagrams.  JHEP, 11:061, 2008.
  * [23] Bin Chen, Sera Cremonini, Aristomenis Donos, Feng-Li Lin, Hai Lin, James T. Liu, Diana Vaman, and Wen-Yu Wen.  Bubbling AdS and droplet descriptions of BPS geometries in IIB supergravity.  JHEP, 10:003, 2007.
  * [24] Eric D’Hoker, Daniel Z. Freedman, Samir D. Mathur, Alec Matusis, and Leonardo Rastelli.  Extremal correlators in the AdS / CFT correspondence.  pages 332–360, 8 1999.
  * [25] Akikazu Hashimoto, Shinji Hirano, and N. Itzhaki.  Large branes in AdS and their field theory dual.  JHEP, 08:051, 2000.
  * [26] Kasia Budzik and Davide Gaiotto.  Giant gravitons in twisted holography.  6 2021.
  * [27] David Berenstein.  Giant gravitons: a collective coordinate approach.  Phys. Rev. D, 87(12):126009, 2013.
  * [28] M. Adler and P. van Moerbeke.  Integrals over grassmannians and random permutations, 2001.
  * [29] Eduardo Duenez.  Random matrix ensembles associated to compact symmetric spaces, 2003.
  * [30] David Berenstein and Randel Cotta.  A Monte-Carlo study of the AdS/CFT correspondence: An Exploration of quantum gravity effects.  JHEP, 04:071, 2007.
  * [31] Takuya Okuda.  A Prediction for bubbling geometries.  JHEP, 01:003, 2008.
  * [32] Takuya Okuda and Diego Trancanelli.  Spectral curves, emergent geometry, and bubbling solutions for Wilson loops.  JHEP, 09:050, 2008.
  * [33] Jaume Gomis and Takuya Okuda.  Wilson loops, geometric transitions and bubbling Calabi-Yau’s.  JHEP, 02:083, 2007.
  * [34] Hai Lin, Oleg Lunin, and Juan Martin Maldacena.  Bubbling AdS space and 1/2 BPS geometries.  JHEP, 10:025, 2004.
  * [35] Davide Gaiotto and Ji Hoon Lee.  The Giant Graviton Expansion.  9 2021.
  * [36] Adolfo Holguin and Shannon Wang.  Giant gravitons, Harish-Chandra integrals, and BPS states in symplectic and orthogonal $\mathcal{N}$ = 4 SYM.  JHEP, 10:078, 2022.
  * [37] Hai Lin.  Coherent state excitations and string-added coherent states in gauge-gravity correspondence.  Nucl. Phys. B, 986:116066, 2023.
  * [38] Mark Srednicki.  Chaos and Quantum Thermalization.  3 1994.
  * [39] Alba Grassi, Zohar Komargodski, and Luigi Tizzano.  Extremal correlators and random matrix theory.  JHEP, 04:214, 2021.
  * [40] Stephen W. Hawking, Malcolm J. Perry, and Andrew Strominger.  Soft Hair on Black Holes.  Phys. Rev. Lett., 116(23):231301, 2016.
  * [41] Massimo Bianchi, Daniel Z. Freedman, and Kostas Skenderis.  Holographic renormalization.  Nucl. Phys. B, 631:159–194, 2002.
  * [42] Matan Grinberg and Juan Maldacena.  Proper time to the black hole singularity from thermal one-point functions.  JHEP, 03:131, 2021.
  * [43] Warren Carlson, Robert de Mello Koch, and Minkyoo Kim.  Generating functions for giant graviton bound states.  JHEP, 01:104, 2023.
  * [44] Chi-Ming Chang and Xi Yin.  1/16 BPS states in $\mathcal{N}=$ 4 super-Yang-Mills theory.  Phys. Rev. D, 88(10):106005, 2013.
  * [45] Jan Boruch, Matthew T. Heydeman, Luca V. Iliesiu, and Gustavo J. Turiaci.  BPS and near-BPS black holes in $AdS_{5}$ and their spectrum in $\mathcal{N}=4$ SYM.  3 2022.
  * [46] Kasia Budzik and Davide Gaiotto.  Twisted holography without conformal symmetry.  11 2022.
  * [47] Simeon Hellerman, Domenico Orlando, Susanne Reffert, and Masataka Watanabe.  On the CFT Operator Spectrum at Large Global Charge.  JHEP, 12:071, 2015.
  * [48] Gil Badel, Gabriel Cuomo, Alexander Monin, and Riccardo Rattazzi.  The Epsilon Expansion Meets Semiclassics.  JHEP, 11:110, 2019.
  * [49] Gabriel Cuomo.  OPE meets semiclassics.  Phys. Rev. D, 103(8):085005, 2021.
  * [50] David Eliecer Berenstein, Juan Martin Maldacena, and Horatiu Stefan Nastase.  Strings in flat space and pp waves from N=4 superYang-Mills.  JHEP, 04:013, 2002.
  * [51] Simone Giombi and Jonah Hyman.  On the large charge sector in the critical $o(n)$ model at large $n$, 2020\.