1/2 BPS structure constants and random matrices

We study three point functions of half BPS operators in N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 4 super Yang-Mills theory focusing on correlators of two of the operators with dimension of order ∆ ~ N2 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.


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, nontrivial 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 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 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 where C R 1 R 2 R 3 are Richardson-Littlewood coefficients and f R is the norm of the operator 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 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 ∆ ∼ 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 setups 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] ρ(z) = (z 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 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 where the density ρ 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.

Coherent States and Form Factors
Most of our discussion will concentrate on the simplest correlation functions in the N = 4 SYM theory, which are three point functions of half-BPS operators.We will also work mostly with the theory on the cylinder R × 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 • n = 0: 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 In the large N limit with ∆ = k L k √ 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 × S 5 and their bound states.However, for operators with large enough conformal dimension ∆, 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 Y M = 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: where k denotes the number of boxes of the Young diagram associated to the representation R of U (N ) and χ R (π) 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 × 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) A straightforward computation gives an alternate formula for |Λ as a expansion in characters of the unitary group s R , and s R (Λ) is a Schur polynomial.The point of this analysis is that by exploiting the Campbell-Hausdorff formula the free field contractions of the operators O R (Z(x, n)) 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: 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: (2.7) where the initial and final states are created by heavy operators ∆ ∼ 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).

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 : which describes a sphere giant graviton sitting at the origin of global AdS at a position inside S 5 given by e iφ 0 cos θ 0 = λ.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 semiclassical 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: As discussed in [4,11] this is the same state that one obtains from setting Λ to be a rank one projector in (2.6).This state has a natural U (1) gauge symmetry 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: 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 ϕ.Then we can split the coordinates of CP N −1 into ϕ 1 and ϕ 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 ϕ n in terms of an 2N − 2 dimensional spherical slices of radii Finally we can rewrite the radial part of the integral as ˆd (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 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, λ → √ N − 2 λ we finally arrive at the expression which apart from the contour on integration is identical to the integral obtained from the Hubbard-Stratonovich trick used in [9] for determinant operators 1 .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: (2.15) 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 |λ| = cosh ρ 0 is greater than one which simply says that the brane is at a position ρ 0 > 0 in global AdS, and this is true for all half-BPS giant graviton solution.

Form Factors
The next step is compute the following form factor; Our choice of 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 Z by constant matrices.The resulting trace is (2.17) To proceed we use the fact that the integrals over dϕ and d φ are invariant under the action of U (N ), so we can gauge fix φ 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: (2.18) Now instead of evaluating every power the matrix ρ, it is better to work with the resolvent When we evaluate this expression at the saddle point value for ρ, we can immediately recognize that the resolvent R(t) is a generating function for Chebyshev polynomials of the first kind.
where we parametrized the eigenvalue in terms of LLM coordinates λ = e iφ 0 cosh ρ 0 .

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 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 |λ|.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 |λ| = 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 λ 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 λ.Since the wavefunction λ λ is already invariant under shifts in the phase of λ whenever λ = λ * , 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 → √ det ρ t, and then perform the integration over the phase of λ by a contour integral.The result is In the large N limit, we can set det ρ = 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 ρ 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

Two Droplets
The next simplest calculation that we can perform is the case where the matrix Λ 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: We will argue that this state described the wavefunction of a stack of p giant gravitons sitting at position λ in the LLM plane.The first thing to note is that this coherent state has an explicit U (p) gauge symmetry V ∼ V g which we can identify as the gauge symmetry on a stack of D-branes.The expectation value of ∆ 0 on this state is given by p|λ| 2 so that whenever |λ| ∼ √ N and p ∼ 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 Z on this state, we can see that this state breaks the gauge symmetry spontaneously from U (N ) to U (N − p) × U (p), and that the center of mass of the stack of p branes is at the position z = λ 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 where A 11 and A 22 are p×p and (N −p)×(N −p) square matrices, and A 12 is a p×(N −p) matrix (similarly for A 21 ).Then we make an arbitrary choice of frame distinguished by a rectangular matrix this matrix can then be used to build projectors into arbitrary p-dimensional subspaces of 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: Clearly any V in this space is rank deficient, so all of the information about its singular value decomposition is captured a square matrix: In analogy to the rank one calculation, we can write down the integration measure for this space as By a similar calculation we can see that the norm of this state is given by One last fact that we need before continuing is the Jacobian for the coordinate transformation M = A † A for any n × m rectangular matrix A. This change of variables is well known in the context of Wishart distributions; 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 λ, p λ, p = n λ,p ˆA † (3.10)

Large N Limit: Steepest Descent
We will now sketch the evaluation of (3.10) in the large N limit.As before, |λ| 2 will scale with N so that the exponential term is large.The integral can evaluated explicitly using the Andreief identity [13]: 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 λ → √ N − 2p λ the saddle point equations of this integral are of the form 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 λ λ = 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 ≤ N/2.This makes sense, since condition forces the second droplet to be smaller than the first.In the case that p 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 ∼ 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): where 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 1  N −2p as the relevant parameter in the problem; this is the reason why p ∼ 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 |λ, 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 λ as √ p ∼ √ 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 with µ = p N −2p ; this turns out to be simpler after we make the change of variables z = λ λ (1 − x).After this we can normalize the distribution such that ´ρ(z)dz = 2µ and find K = 1.In this convention the eigenvalues are quantized in units of 2g s = 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 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 × 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 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.

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 × p matrix.In this case the moment generating function is The most natural variable to work with is once again z i = λ λ(1 − x i ), which makes the density of eigenvalues be of the form: where z ± are given by the roots of the polynomial inside the square root in (3.14) Expanding 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 th moment we get λ Tr To extract the diagonal part of this form factor we can average over the phase of λ.
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 |λ| 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.

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 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 σ † σ, and the matrix σ is not a normal matrix, so σ and σ † 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: 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 σ satisfying Indeed we will see later that these are the saddle point conditions for the integration over the angular variables for σ.Because of the large exponent, corrections to this are heavily suppressed as long as Λ † = Λ, and so The other saddle points contribution for the norms consists of pairing the eigenvalues λ i with λπ(i)] for all permutations π, which are highly suppressed for non-coincident eigenvalues.
e Tr[Λσ Λσ † ] → π∈Sp e i |s i | 2 λ i λπ(i) × (one-loop) . (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 λ a , for a = 1, . . ., k.As long as the λ 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 λ a is computed by the following matrix model: .
(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.

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 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 multicenter solution.

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 λ 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 Λ † , |Λ , but now we treat the parameters as being independent from one another; we will also force each of the eigenvalues λ i and λi to lie on a unit circle.By multiplying |Λ by the square of the Vandermonde determinant of Λ, and integrating we can recognize that the resulting integration measure is just a Haar measure for a new unitary matrix U = U ΛU † : ˆU(N) The resulting state is clearly proportional to the vacuum state, since there are no U † insertions to feed to the exponential.From this it becomes clear that in order to extract a term proportional to the state |R = S R (Z) |0 one should multiply the integrand by a character S R ( U † ): 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 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 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 in the limit that |R| ∼ |R | ∼ 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: This procedure replaces all of the free-field Wick contractions with unitary integrals which we can evaluate very explicitly.

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: 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 where ν are the eigenvalues of the Hessian of Tr U † u U v ; this determinant factor is well known and it is proportional to ∆(u)∆(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 → v π(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: 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 ; (4.9) To perform the integral over y j we set y j = e 2iβ 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 π = id.After evaluating the integral over y we get This final integral is a simple Fourier integral that can be evaluated by expanding in powers of x i .After normalizing F(t) appropriately we obtain a formula for the generating function of structure constant C RRL : . (4.11) For large representations R i ∼ N with large blocks, the ratio of gamma functions can be replaced by the asymptotic expansion Γ(x) Γ(x−β) x β as x → ∞, and the sum may be replaced by an integral with x = (N − i) /N and α i = R k+1−i /N : here µ 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: where x = r 2 the density ρ R is defined as follows.
Given the Young diagram we rotate it by −3π/4 radians.Then the diagram has (k + 1, k) edges of slopes ∓1 of lengths µ i and κ i for negative and positive slopes respectively.For example the length of the 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 1a.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 ρ R (r) = 1, and ρ 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:

.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 O S k,+k + O S k−k and integrate over φ.This would correspond to the contribution of a spherical harmonic Y L/2 ∼ (z + z + y − ȳ) L .

Off-Diagonal Structure Constant
To compute off-diagonal structure constants we need to change one of the representations R to another representation R whose Young diagram is close to R. In this case most of the integrals will vanish, unless R and R only differ at a single row R l − R 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 l , so we only get one term: where ∆S R R is the ratio of norms . This should be interpreted as the linear response to a fluctuation localized at the edge of a particular Fermi surface within the LLM geometry.

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 λ λ 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): (4.16)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 n k ).Then we use the fact that the integration variables are invariant under permutations to  where N a is the partial sum b≤a n b and N 1 = 0. Putting it all together, we are left with a multi-cut matrix model: 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 /Z p .To make the analogy more precise we can change variables to x i = e −u i after which the partition function becomes: The only difference between this and the matrix model studied in [20] is that the Gaussian term is replaced by (1 − e −u (a) i ) 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 (a) i take the form: (4.20)The resolvents for this type of problem can be taken to be of the form ω (a) (z) = , and the total resolvent is ω(z) = k b=1 ω (b) (z).At large N the eigenvalues condense into k branch cuts and the saddle point equation on the a th cut reads or equivalently This equation defines a Riemann-Hilbert problem for the total resolvent ω(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 Re(z) > 0. Following [20] we can define a set of complex variables 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 ω ∝ 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 turn the terms of the form coth 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: Clearly the eigenvalues λ 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, and the distribution ρ R (r) for r k is the same distribution we found before.

Computing correlators: reduced matrix elements
To compute correlation functions of single trace operators between a set of states |R, p , |R , 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 If we insert a half BPS single trace operator between these states, the scalar fields Z and Z inside the trace are traded by complex valued matrices: We can then perform a color-flavor transformation to express any trace of a power of Z + Z (over N color indices) into a trace over 2p × 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 × p 'radial' matrix and here 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 σ.
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: The critical points of this function are given by the solutions to a pair of matrix equations 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β and U † R s † U † L e iα U L s U R are both diagonal in the same basis.The second equation gives a similar condition for e iα and But since e iα , e iβ 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 π in the first equation for some permutation matrix U π 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 π, τ in S p , which act on e iα and e iβ 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 − (e iα i − e iα j )(e so the one-loop factor becomes simply Each saddle point is weighted by det(τ π), 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 α i and β i , which makes the integrals over the angular variables straightforward For light operators and diagonal structure constants we can simply interchange (Z + Z) L with L L/2 (Z Z) L/2 , which is simpler to work with.The integration over α i − β 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 α i + β 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:

R, p|
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.

Discussion and Future Directions
In this paper we introduced techniques for dealing with large BPS operators in 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 semiclassical 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 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 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 (4.34).For example one can insert a particular Schur polynomial of Z + Z + 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.A more long term goal would be to understand the systematics of less supersymmetric BPS operators in the N = 4 SYM theory.In principle introducing additional scalar matrices in the exponentials gives a way of generating 1 4 and 1 8 BPS [4,3,43].This was carried out for rank two 1  4 BPS operators in [43].Additionally, the generating functions introduced in [4] contain all bosonic 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 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 1  16 BPS states by exponentiating the so-called 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, 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 P 1 instead of 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; 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.

Figure 1 :
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. ) Since the only combination of A 11 and A † 11 that appears in the integral is A † 11 A 11 , we can reduce the computation to an integral over the eigenvalues of A † 11 A 11 .