Multi-Matrix Correlators and Localization

: We study generating functions of 14 BPS states in N = 4 super Yang-Mills at finite N . We do this by extending the Harish-Chandra-Itzykson-Zuber localization formula to models with multiple commuting matrices. This allows us to compute the overlaps of two or more generating functions; such calculations arise in the computation of two-point correlators in the free-field limit. We sketch a proof of the fixed point formula, and discuss some connections with the restricted Schur polynomial operator basis. Our results generalize readily to arbitrary numbers of matrices, opening up the opportunity to study more generic BPS operators.


Introduction
Large operators in large N gauge theories are an important subject of study with relevance to nuclear physics, theories of quantum gravity and strings.Although there has been enormous success in computing the spectrum of anomalous dimensions of light operators in models such as maximally supersymmetric Yang-Mills theory in the planar limit, very little is known about how to tackle generic operators whose dimensions can scale with a power of N .This is an interesting problem for holography [1] and for understanding the structure of conformal field theories more generally [2].One of the difficulties one faces when trying to address these types of problems is that the intuitions from the planar limit are often unjustified for large operators; one must sum over both planar and non-planar diagrams and it is not a priori clear which diagrams dominate in the large N limit.A promising approach is to replace single and multi-trace operators with a different basis that is better behaved at finite N [3], and then perform a systematic expansion around protected states in the large N limit.In the case of maximally supersymmetric Yang-Mills theory, this has been implemented at finite N [4][5][6][7].Even though the expressions found through these techniques at finite N are quite explicit, it is usually difficult to take the large N limit of such quantities.
More recently, there have been works showing that certain generating functions can be used to perform computations in the free-field theory limit [8][9][10][11][12].This technique has been succesfully implemented in the computation of three-point correlators involving large operators made out of a single matrix field [13][14][15], as in the half-BPS sector of N = 4 SYM [16,17], where the dual gravitational description is explicitly realized from the gauge theory.An explicit mapping between BPS states made out of more than one matrix and asymptotically AdS 5 × S 5 geometries is still lacking, though a compelling description in terms of bubbling geometries seems to exist [18][19][20].The study of generating functions for multi-matrix correlators was outlined in [10,21] for certain classes of operators, and more generally in [9].Our goal is to elucidate some of the details regarding the generating functions of 1  4 and 1 8 -BPS operators in N = 4 SYM.We do this by proposing a fixedpoint formula for the overlap of generic coherent state generating functions; this gives us an integral formula that generalizes the Harish-Chandra-Itzykson-Zuber (HCIZ) formula to multiple pairs of matrices.Integrals of this type appear naturally in the study of multimatrix models of commuting random matrices.
This paper is structured as follows.In section 2, we review the generating function techniques, focusing on the case of 1  4 BPS operators in N = 4 SYM.We argue that the form of these operators is protected, so we can restrict to eigenstates of the one-loop dilatation operator.We then evaluate the norm of the generating function for the U (2) theory by explicit integration to motivate our fixed point formula for general N .Finally, we give a prescription for extending the HCIZ formula to the multiple-matrix model using the heat kernel method as outlined in [22]; we will discuss our results for U (2) and U (3) and the insights we may glean from them to extrapolate a general formula for U (N ).In section 3, we connect our results to the construction of restricted Schur polynomials and outline how to generalize to operators associated with Young diagrams with arbitrary number of rows or columns.We will briefly discuss our attempts to arrive at a general formula via the character expansion method.Finally, we conclude with some future directions.

Multi-matrix Generating Functions
We are interested in studying operators in gauge theories that are made out of more than one matrix-valued scalar field.In particular, we will work with 1  4 -BPS operators in U (N ) N = 4 SYM on the cylinder R × S 3 .At weak coupling, these operators can be built out of symmetrized products of two of the three complex scalar fields of the theory X, Y .Generalizing to more than two matrices is straightforward.This class of operators transforms in the [p, q, p] representations of the SU (4) R symmetry, and the operators are generically of multi-trace form.We will concentrate on scalar primary states at an equal time slice for simplicity.Unlike 1 2 -BPS operators, which can be built explicitly in the free theory, 2.1 Generating 1  4 BPS States Yet another way of generating 1  4 -BPS states can be found by studying operators of the form: If we insist that the coherent state parameters Λ X and Λ Y commute, |Λ X , Λ Y ⟩ is annihilated by the one-loop dilatation operator; it was shown in [12] that this persists to two-loop order.In [24], it was conjectured that the space of BPS states in N = 4 SYM is given by the kernel of the one-loop dilatation operator at all values of the coupling; we will take this as a working assumption and work with the set of states annihilated by the Beisert one-loop dilatation operator: Because the states (2.1) are coherent states of X, Ȳ [9], they form an overcomplete basis of states for any value of N .This has many computational advantages, mostly due to the fact that taking the large N limit is very straightforward, but translating back into a complete orthogonal basis of operators can be complicated.This may be solved by computing the norm of the coherent states.By exploiting the Campbell-Hausdorff formula, we arrive at an integral of the form: Since we can in principle expand (2.1) in terms of an orthonormal basis, we may use this overlap to determine the coefficients relating the multi-trace basis of operators to an orthogonal basis by expanding in a series and matching the coefficients as done in [9].The precise tool relating the multi-trace basis operators and the character expansion in this case is the Weingarten calculus [25]; an example illustrating this technique can be found in [26].The main obstacle we face is evaluating the integral (2.3) for generic coherent state parameters.To our knowledge, these types of integrals have not been studied before, and a closed form expression for them is needed.Our main goal will be to evaluate this class of integrals for any value of N .Although we only explicitly study the case of U (N ) integrals, the methods should apply generally and should generalize to SO(N ) and Sp(N ) groups as well as to quivers.These types of integrals are also a natural object to study in the context of matrix models, since they arise in the study of multi-matrix models of commuting matrices.

The Four-Matrix Model in SU (2)
Before proceeding to the case of general N , we will study the following integral (2.4) for commuting matrices A, B, Ā, B. We will first approximate I 2 by a saddle point approximation; the critical points of the function in the exponential are given by the solutions to the equations For generic enough matrices, this is only satisfied if each of the two terms vanishes individually The only problematic cases occur when a subset of the eigenvalues of B is a permutation of a subset of eigenvalues of −A.From here on, we assume that the eigenvalues are generic enough that this does not happen.This means that, generically, the saddle points are labelled by permutation matrices U π .We are then left with a Gaussian integral around each of the saddle points, which can be evaluated easily; this results in a "one-loop determinant" factor given by: This gives an approximate value for the integral (up to a convention dependent normalization factor): At first sight, it is not clear that this approximation is reliable, since there is no large parameter in the exponential.To gain more intuition, we evaluate I 2 through an explicit computation.
First, we must parameterize our unitary matrix U ; then, we need to compute the Haar measure.We start with the following matrices: We then seek to parametrize our unitary matrix.We know that any arbitrary SU (2) matrix must meet the following conditions: For ease of computation, we choose to parameterize U with Euler angles: We seek to rewrite the Haar measure dU in terms of J(θ, γ, α)dθdγdα, where J(θ, γ, α) is the Jacobian.We may do so by computing the inverse of the unitary matrix and multiplying it by its partial derivatives with respect to the Euler angles.We start by finding the inverse of U : Then we calculate the partial derivatives with respect to γ, α, and θ and multiply by the inverse.We obtain: We calculate the Jacobian matrix using the following basis and ϵ 3 = 0 −e iγ e −iγ 0 : The Jacobian J(θ, γ, α) we seek is the determinant of J : We see that it is only dependent on θ.Our integral becomes: (2.16) Our critical points are θ = 0 and θ = π, so we can remove the absolute value bars.Then we evaluate our integral: (2.17)This is precisely the same result that the saddle point approximation yields.From the intermediate steps, it is clear that there are never any terms that mix the eigenvalues of A and B; if we set either A = 0 or B = 0, we immediately recover the HCIZ formula for U (2).

Proof of Localization for U (2)
One important issue to understand is why the integral I 2 has an exact saddle point approximation, while the naive saddle point approximation for I N fails to be exact for N > 2. The idea is to try to follow the proof [22] for the HCIZ integral and see exactly how the analysis differs for multi-matrix models.One key observation is that the integral I N (A, B, Ā, B) is an eigenfunction of a holomorphic Laplacian: Before continuing, it is worthwile to explain what we mean by holomorphic in this context, and why this is important.One way in which the integral I N appears is as the Jacobian factor for a Gaussian matrix integral over a pair of commuting normal matrices As it stands, this expression is formal unless we specify a contour of integration for the eigenvalues of A, Ā and B, B. A choice of contour corresponds to a choice of polarization in the space of eigenvalues; this makes the eigenvalues of A and Ā canonically conjugate.This is quite natural from the interpretation of the integral as the norm of a coherent state in matrix quantum mechanics, where the collective coordinates a i and b i are holomorphic phase space coordinates.The barred coordinates are then conjugate momentum variables.Thus, the correct Laplacian operator has to be constructed from the metric of the space of commuting normal matrices.This is exactly the quantization discussed in [16].We rescale the matrices by a constant factor t; the resulting equation implies that the integral is related to a holomorphic heat kernel on the space of commuting matrices: As we take t to zero, the integral will be very well approximated by the saddle point approximation.We see that the integral approaches a delta function; we can use the kernel itself to propagate this initial condition to a finite t.This would imply that the integral comes from a sum over the real saddle points of the integral.If the kernel is a plane wave, then the integral localizes, which is to say that the steepest descent contour gives exactly a Gaussian integral centered around each saddle point.This occurs if the heat equation for the kernel corresponds to a Schödinger equation for an integrable system, since we can in principle change the variables into a set of action-angle variables, where the wavefunction is a plane wave.Whenever the kernel cannot be written this way, true localization fails, and instead, the integral is given by a sum over thimbles, with the kernel giving a parametrization of the integration contour.Now we review some coordinate transformations for the Laplacian in the space of normal matrices.Given that the squared distance of two n × n normal matrices A and A ′ is d(A, A ′ ) = Tr |A − A ′ | 2 and invariant under unitary transformations A → U AU † , the metric is: (2.21) We know that the Laplace-Beltrami operator is: We now perform a coordinate transformation A = U ΩaU † to rewrite the matrix in terms of its n eigenvalues a i and n(n − 1) angular variables θ α .Setting dH = iU † dU , we may rewrite our invariant distance as: where we have defined: where The square root of the metric tensor's determinant is then: where ∆ is the Vandermonde determinant of the eigenvalues a i .Because the new metric tensor is block diagonal with an eigenvalue sector and a unitary sector, its inverse is block diagonal as well with an eigenvalue sector and a unitary sector, and thus the Laplacian may be separated into two operators, one for each sector.In its entirety, the Laplacian is: (2.26) We now consider the space of two N × N commuting normal matrices A and B. After diagonalization the metric for this space becomes: Then the square root of the metric tensor's determinant becomes: and we rewrite We may rewrite the holomorphic Laplacian as [16]: Notice that all of the eigenvalue dependance is on the first two terms and because the integral averages over the angular variables of A and B it is annihilated by the last term, so we will omit it from now on.So now our problem is reduced to finding eigenfunctions for this operator.As is common in matrix quantum mechanics, one can often reabsorb the measure factor µ into the definition of the eigenfuntion, so we will express I N (A, B, Ā, B) in terms of an auxiliary function Ψ N (A, B, Ā, B): After this rescaling, the Laplacian operator becomes a sum of two terms, one being the flat space Laplacian and the other an effective potential: So far our discussion applies to general rank of matrices.Focusing on N = 2, we can easily check that the potential term vanishes.This is because µ is linear in a i and b i .In this case, the problem reduces to finding eigenfunctions for the Laplace operator in flat space: (2.32) The solutions to this equation are plane waves: This ansatz does not respect the symmetry properties of the integral under simultaneous permutations of a i and b i , so the correct solution is a symmetrized sum of plane waves: After dividing by the measure factor, we reproduce the expected answer.At this point, it becomes clear that the heat kernel proof works for the SU (2) integral, since the kernel is Gaussian and the saddle point approximation as t → 0 can be propagated forward to obtain the integral for finite t.Intuitively, we should be able to localize the integral much like the single matrix case, because the wavefunction Ψ 2 is an eigenfuction of an integrable (free) Hamiltonian.In the case where B = B = 0, the measure factor µ reduces to a Vandermonde determinant, which is also annihilated by the flat space Laplacian; thus we see that the usual HCIZ integral is associated with a wavefunction of a free fermion or boson.But this is not the case for B ̸ = 0 and N > 2, since the potential term does not vanish.Note that this does not necessarily mean that the integral is not localizable.An example that comes to mind are the integrals of the Harish-Chandra type for the symplectic groups, which are associated with the wavefunctions of integrable Calogero models.While these integrals are known to localize by the heat kernel methods [27], in this case, it was noted that the naive localization argument nevertheless still fails [28], and that one must include additional instanton solutions to the WKB approximation.Returning to our guess for Ψ N , it becomes clear that a closed form for I N (A, B, Ā, B) must include µ as its denominator.This is because µ is the natural integration measure for the eigenvalues a i , b i .The fact that the eigenvalue is Tr Ā2 + B2 also suggests that the denominator should have an exponential factor: (2.35) By a symmetry argument it is also plausible that the numerator is also given by a determinant.We seek the missing factor in the numerator; we know that the Bethe ansatz is unlikely to provide a solution to our differential equation, because our effective potential appears to contain three-body interactions.Finding such a formula would amount to finding eigenfunctions of ∇ A,B along the lines of [16,29], but in our case we are only interested in the ground state wavefunction in the effective potential.While a complete analysis is beyond the scope of this paper, we may still lay out a prescription for finding an analytical solution to our modified integral.We previously argued that such a solution must have µ as its denominator; following [16], the problem can be simplified by rewriting the equation for I N and removing the effective potential in the Laplacian at the cost of adding first order derivative terms.Thus we may rewrite (2.35) as: where we have set f π = c π i e a i āπ(i) +b i bπ(i) .We see then that we have: Expanding, we arrive at: (2.38) We compute the terms containing the second derivatives and find: Thus we can simplify our differential equation: We see then that we have a multivariable PDE; while there are numerical methods to approximate an analytical solution, they are all computationally laborious and ill-suited to solving PDEs containing more than two independent variables.Nevertheless, given the nature of the HCIZ integral, we note that there is another method to compute an analytical expression for the missing factor that, while tedious, is still tractable.That is to evaluate the integral explicitly and use this to determine the missing factors.We solve the integral for U (3) in Appendix A. We will reproduce the result directly below: The expressions for w i are listed in (A.29); we see immediately that the a i and āj eigenvalues do not mix with the b i and bj eigenvalues.We make note of the SU (2) symmetry between a i and b i , and āj and bj ; we note that if we remove the b i and bj , we recover the expression for the two-matrix Harish-Chandra integral evaluated over the U (3) group.
Given that we may remove b i and bj , replace a i and āj with u i = a i b i and ūj = āj bj , and still recover the same expressions for w i , we could naively expect to recover an expression similar to the original Harish-Chandra integral [22]: if we use u i and ūj in lieu of a i , āj , b i , and bj ; set Ω as the normalization constant; and specify that when multiplying the Vandermondes in the denominator, one must take the dot product of (u i − u j ) and (ū i − ūj ).But upon this substitution, we find that we only recover: We are missing the factors of χ(a i , b i , āj , bj ) in the numerator.We can solve for the missing factors by examining the case of U (3); given that the effective potential is the same general form for N > 2, we may extrapolate the missing factor for a general U (N ) formula from the U (3) results.
In the case of U (3), we can expand the original HCIZ integral: until we arrive at a series that takes the form of (A.39) with modified w i to reflect the omission of the b i and bj .We examine the effects of replacing a i and a j with u i and u j at each step; we then compare the results to an expansion of (2.42).The additional terms that replacing a i and a j with u i and u j yields must sum up to the missing factors.We briefly sketch out the start of such an expansion.We note that for SU (3), we have: det exp a i āj + b i bj = −e a 3 ā1 +a 2 ā2 +a 1 ā3 +b 3 b1 +b 2 b2 +b 1 b3 + e a 2 ā1 +a 3 ā2 +a 1 ā3 +b 2 b1 +b 3 b2 +b 1 b3 + e a 3 ā1 +a 1 ā2 +a 2 ā3 +b 3 b1 +b 1 b2 +b 2 b3 − e a 1 ā1 +a 3 ā2 +a 2 ā3 +b 1 b1 +b 3 b2 +b 2 b3 − e a 2 ā1 +a 1 ā2 +a 3 ā3 +b 2 b1 +b 1 b2 +b 3 b3 + e a 1 ā1 +a 2 ā2 +a 3 ā3 +b 1 b1 +b 2 b2 +b 3 b3 = −e s 1 + e s 2 + e s 3 − e s 4 − e s 5 + e s 6 (2.45)We have set: If we expand each term as a Taylor series, we may rewrite our determinant as: We note that: We see immediately that if we remove the b i and b j , the factors listed above cancel out factors in the Vandermonde determinants of the original HCIZ integral; however, once we add in b i and b j , some of the factors in (2.48) no longer cancel factors in µ.This suggests that there are missing saddle points and that the missing factors should add the terms needed to restore the overall factor of µ in the numerator.We note that our results should generalize to an arbitrary number of matrices; we would simply modify µ to account for the additional matrices and add the relevant derivative terms to (2.40), as well as modify u i and ūj to account for the new eigenvalues.
We leave off here and save this computation for future works.

Connection with Restricted Schur Polynomials and Collective Coordinates
A natural question to ask is: what sort of basis of operators do the coherent states (2.1) actually generate?This is quite non-trivial, since there are in principle many different ways of orthogonalizing the two point function of 1  4 -BPS operators at finite N .As a concrete example, we can take the simplest multi-matrix coherent state we obtain from choosing the coherent state parameters to be rank one projectors for arbitrary N .These states will describe semi-classical configurations of single quarter BPS giant gravitons.Similar generating functions were introduced in [10,30]; we will clarify the relationship between them and the coherent states studied here.This should help in generalizing to higher rank cases corresponding to bound states of AdS giants.The idea is to consider the following state: To evaluate this, we need a formula for the moments of φ † i φ j with respect to the flat Haar measure on CP N −1 .The measure can be rewritten as follows: In other words, we can trade the integral over projective space for a regular Gaussian integral at the cost of introducing an additional contour integral over an auxiliary parameter s.Then the moments have a simple expression in terms of the projection operators P (k) [10]: Borrowing the results of [10], we may rewrite the coherent state as a sum of the so-called restricted Schur polynomial operators χ (k 1 +k 2 ),(k 1 ) (k 2 ) (X, Y ): Now we would like to understand the analogue of this formula in the general case.First, we need to recall the definition of the restricted Schur polynomials: Here, R is a Young diagram associated with an irreducible representation of S n+m ; the labels (r, s) correspond to an irreducible representation of S n × S m contained in R. The object P R,(r,s) αβ can be understood as follows: starting with S m × S n ⊂ S m+n , we can find representations r × s sitting within R. Generically, the representation r × s can appear more than once inside of R, so one needs to keep track of how one embeds r × s into R.If the multiplicity of (r, s) is n (r,s) and its dimension d (r,s) , then a generic element of S n+m will be block diagonalized into n (r,s) d (r,s) × n (r,s) d (r,s) blocks.The matrix indices α, β keep track of this information, where α and β range from 1 to n (r,s) .The P R,(r,s) αβ are then intertwining operators between each of these blocks.More formally, we can label each of the embeddings of r × s by an index γ and consider the space R γ ⊂ R. The restricted Schur polynomial is then given by: where Γ R (σ) is the matrix representing σ [4].The most complicated part of the restricted Schur polynomials is the evaluation of Tr Rγ [Γ R (σ)], which involves building R γ explicitly.By expanding the exponential and evaluating the unitary integrals, we obtain: where W g(σ, N ) is the Weingarten function.Explicit combinatorial formulas for Weingarten functions are well known from the work of Collins (see [25] for an elementary introduction); before delving into specific details, we should contrast this with the situation where one of the Λ X,Y is zero.In this case, the resulting sum can be recast as a diagonal sum of the products of unitary characters; right now, we have a complicated sum of traces.For a moment, let us consider the situation for a single matrix.The resulting sum is: The last line is obtained from the character expansion of the integral, which was computed in [9].Then for two matrices, we have: Clearly this has a similar structure to the definition of the restricted Schur polynomials (3.6), but the restricted characters have been replaced with ordinary symmetric group characters instead.We can always formally rewrite each of the terms in the series as a sum over restricted characters by decomposing the trace over R into a sum of traces over each of the R γ : This allows us to rewrite the integral as a sum of restricted Schur polynomial operators.However, this sum does not capture every restricted character; the problem comes from elements σ ∈ S n+m that are not elements of S n × S m .Generically, the representation R of S n will be expressed as a sum of the irreducible representations R α of S n × S m with multiplicities n α .This means that σ is not fully block diagonal on the space n α R α , and indeed, restricted traces on each of these blocks lead to different orthogonal states in the free theory.This is expected, since the coherent state only generates operators that have vanishing one-loop anomalous dimension, and the dilatation operator acts non-trivially on generic restricted Schur polynomial operators.This means that the operators obtained from diagonal traces over n α R α have vanishing one-loop anomalous dimension, while offdiagonal operators should be associated with excited states with open strings.We expect that operators that are approximate eigenstates of the dilatation operator in the large N limit take the form of open string modifications of these coherent states, and are likely more closely related to the Gauss graph operators as in [31].
An interesting direction to take would be to construct a generating function for all restricted Schur polynomials.Clearly, something as simple as (2.1) cannot work.This can be traced back to the fact that the sum over S n+m has many redundancies owing to the fact that we can conjugate by an element of S n × S m while leaving the traces fixed.This is the statement that we can permute the n X's and m Y 's among themselves while simultaneously permuting the Λ X,Y 's.As explained in [32], there is an equivalence relation between elements of S n+m in such a way that which happens exactly when σ can be conjugated into τ by an element of S n × S m .In other words, the construction of restricted Schur polynomials is equivalent to constructing generalized class functions on restricted conjugacy classes, which means that the coherent state generating function (2.1) cannot differentiate between different restricted Schur polynomials by itself for the simple reason that the Weingarten function is a class function.If we want to replace the characters in (2.1) with restricted characters, we must either change the domain of integration or integrate against an appropriate measure factor that is sensitive to this information.This is equivalent to finding an analytic formula for restricted characters, which may be recast as a Schrödinger problem over the space of commuting matrices [18,29].The point is that the norm of the quarter-BPS coherent state is related to the heat kernel over the space of commuting matrices, or equivalently to the Green's function of the Schrödinger equation.In practice, the coherent states still form an overcomplete basis of operators that can be used for computations, even if we do not currently know how to project into a particular primary state; the leading contribution at large N will come from the saddle point approximation.
In Appendix B, we checked a few low order terms by explicit calculation and found that the quarter-BPS coherent state generating function is given by a sum of product of restricted Schur polynomials much like the half-BPS coherent state.As we explained, this should hold for all the terms in the series, but checking higher order terms is difficult.This can be taken as further evidence that the coherent states span all possible BPS states and makes manifest many of the ideas in [18], since free field theory correlators can be encoded in integrals of polynomial functions of the collective coordinates.In the saddle point approximation, the integration measure is given by (3.12)After an appropriate choice of contour for which λ and λ are canonically conjugate variables, this reproduces the strong coupling ansatz in [18].It should also be clear that this measure describes the ground state wavefunction.Although our analysis is strictly on the weak coupling regime, we note that the dilatation operator in the SU (2) sector can only act by permutations; this sector is closed, so it is very plausible that the quarter-BPS states are not renormalized.Even at weak coupling, the vacuum structure becomes quite nontrivial, modifying the Coulomb branch analysis for small collective coordinates, since the eigenvalues behave as strongly coupled bosons at low energies.This modifies the topology of the moduli space of vacua near the center of mass of the eigenvalues, even for half-BPS configurations.Since this sector preserves as much supersymmetry as N = 2 SYM, it is quite plausible that the g Y M corrections are under control for sufficiently small modifications of BPS operators.At strong coupling, the expectation is that such states describe rotating strings propagating on bubbling geometries; the energy of these strings will follow the dispersion relation of a centrally extended BPS state with the central charge M being related to the length of the string in the bubbling geometry times the string tension.This is purely a kinematic effect; all of the dynamics should be encoded in the central charge M .Near the core of the geometry, the naive Coulomb branch analysis certainly breaks down, but S-duality considerations suggest that the string tension is not renormalized [33], so the corrections to curvature on the moduli space should come from finite N effects.This should be captured by non-planar corrections involving the exchange of gravitons between the background and a probe string.Such geometries can be engineered by integrating against the wavefunctions that break the SO(6) R symmetry of the vacuum.It is not hard to come up with such wavefunctions (for instance a non-symmetric Gaussian perturbation), and there is a schematic mapping between the eigenvalue distribution at large N and bubbling geometry [19].The relation (3.12) should be corrected with additional 1/N effects, since the naive saddle point approximation fails to give the exact overlap, but these effects should only be relevant when we try to probe eigenvalues are placed in non-generic configurations.These should be thought of as microstate configurations for coarse-grained eigenvalue droplet configurations associated to superstar geometries.For 1 8 BPS states, the analysis is more subtle; we have to take into account the effects from fermions since we would be working in a SU (2|3) subsector.One should be able to ignore the effect of the fermions for large enough semiclassical operators.This makes this class of coherent states ideal for studying near-BPS limits around large operators without having to deal with the mixing of multi-trace structures.

Discussion
In this paper, we studied multi-matrix coherent states for bosonic matrices that generate 1 and 1 8 BPS states in N = 4 SYM.We showed that the norm of these coherent states admits a fixed point formula generalizing the Harish-Chandra-Itzykson-Zuber formula for gauge group U (2), and provided evidence of an expansion in terms of restricted Schur polynomials for U (N ).This gives in principle a way of generating expressions for BPS states for any value N in N = 4 SYM.One technical obstacle we face is that our construction does not give an alternative construction of the so-called restricted Schur polynomial operators [4].This is related to the expectation that there is a hidden symmetry under which different operators are charged.One idea is that determining the Casimir charges should be enough to differentiate between different operators, but this problem is quite non-trivial even in the 1  2 BPS sector [34].It is also unclear how to implement this idea efficiently at large N since the number of Casimirs needed to distinguish between different operators grows with the complexity of the operators.Despite this obstacle, our results are important for computing correlators of 1  4 and 1 8 BPS operators dual to bound states of giant gravitons [35] and generic bubbling geometries [19].Understanding the precise map between the overcomplete 'eigenvalue basis' of coherent states and specific orthogonal bases of operators remains an important problem.We conclude with a few more immediate directions for future work.

BPS States and Black Hole Microstate Operators
One of the more interesting generalizations would be to the case of 1  16 BPS operators.By now, there is ample evidence that there exists a class of 1  16 BPS operators describing the microstates of supersymmetric black holes in AdS 5 × S 5 [36][37][38][39].Recently, there have been some studies of these types of states for small values of N [40,41]; see [42] for a more general discussion.our results imply that one should be cautious in extrapolating results about operators for the U (2) theory, since multi-matrix states appear to be qualitatively different for other values of N .We expect that most of the interesting qualities of such operators are missing from the U (2) and SU (2) theory.It would be nice to develop more systematic techniques to build these types of operators.In principle, there are no obstructions to generalizing our techniques to this setup, with the working assumption that finding states with vanishing one-loop anomalous dimension is enough [24].The idea would essentially be to build a superfield coherent state [43]: where Ψ(z, θ) is the C 2|3 superfield discussed in [38,43], and Φ is an auxiliary superfield of coherent state parameters.The combined effect of the exponentiation and integration over the unitary matrices is to generate all possible gauge invariant tensor contractions.One should expect that the operators generated by this generating function are generalizations of the SU (2|3) restricted Schur polynomials constructed in [44].Generically, the terms in the expansion of (4.1) will not be of multi-graviton form, so they are natural candidates for microstates of supersymmetric black holes.In practice, the main disadvantage of an expression like (4.1) is that it might not be practically useful, in the sense that the expansion necesarily involves an infinite number of matrix fields associated with covariant derivatives acting on the fields.One way of avoiding this difficulty is to use generating functions such as the ones studied in [12].Alternatively, one can view the auxiliary superfield Φ as a full-fledged dynamical collective coordinate.One would then hope that integrating out the SYM fields leads to an effective matrix quantum mechanics describing (near)-BPS black hole microstates, with the lightcone coordinate z acting as a time variable.

Three Point Correlators, Bubbling Geometries, and Twisted Holography
Although eventually we would like to study black holes, it is important to build intuition from simpler examples.One class of such examples is the BPS bubbling geometries [19] generalizing LLM geometries [17].Although the droplet description of such states in supergravity is compelling, a precise mapping between the weak coupling BPS states is not fully developed1 .The coherent states (3.7) have a more natural connection to such geometries [15].A worthwhile exercise would be to study correlators of single trace chiral primaries in the background of heavy coherent states corresponding to both giant gravitons or bubbling geometries; see [45] for some finite N results.The holographic renomalization techniques of [46] are also applicable in these cases, but it would be interesting to develop more efficient computational techniques in supergravity along the lines of [47].A good toy model for this would be to study these types of questions in Twisted Holography [48].
In that set-up, the eigenvalue droplets should be related to bubbling on a 2 dimensional complex base with holomorphic coordinates X, Y , with the vacuum configuration being a droplet with the topology of S 3 ; this is the deformed conifold description of SL(2, C).More generic eigenvalue configurations should lead to other non-compact Calabi-Yau threefolds such as in [49,50].The expectation is that the geometry on both sides of the duality is encoded by a spectral curve [51] on which stacks of branes are supported, which should appear as the spectral curve of the corresponding matrix model.Understanding this could help in clarifying the dictionary between collective coordinates and bubbling geometries in the AdS 5 case.

Bound States of Giants and Branes at Angles
One reasonable goal would be to understand coherent states associated with pairs of eigenvalues, which are built from simpler integrals over the Grassmanian Gr(2, N ).The main difficulty of such a character expansion already appears in this simpler case: where the λ α are not 2 × 2 diagonal matrices.The norm of this state has a rather explicit integral form over 2 × 2 matrices over a compact domain; thus an explicit evaluation should be feasible.We expect that this state has a non-trivial expansion in terms of restricted Schur polynomials [10].Since the domain of integration is simpler than that of the general case, the integral might lead itself to a saddle point analysis.One might be able to explicitly compute Lefshetz thimbles for this case and determine whether localization fails or not, or whether complex saddle points are needed.It would be interesting to construct multi-matrix analogues of the generating functions found in [21] for determinant operators.This would help in constructing the precise operators dual to intersecting giants [52], and particularly in understanding their integrable boundary states [53].
A The Four-Matrix Model in U (3) We now consider the following integral: We know that we can parameterize our U (3) matrix as: U = e iλ 3 α e iλ 2 β e iλ 3 σ e iλ 5 θ e iλ 3 a e iλ 2 b e iλ 3 c e iλ 8 ϕ , (A.2) where λ i denotes the ith generators of U (3).We list the relevant SU (3) generators below [54]: , we multiply U by an additional phase e iψ .This yields the parameterization of U (3) that we will use to compute the argument of the exponential in Eq. (A.1).Simplifying and expanding, we arrive at: × e 2i(σ+a) cos θ sin 2β sin 2b (A.4)We plug this into the integral.We may use the method outlined in [54] to compute the Haar measure; we arrive at: dU = sin 2β sin 2θ sin 2b sin 2 θdαdβdσdθdadbdcdϕdψ (A.5) We cite the angle limits from [54]: We seek to integrate over σ and a first.We observe that: Returning to our integral, we note that rather than integrate over σ, we may perform a change of variables and integrate over σ + a and change our integration limits to [0, 2π].We would integrate over a twice, but the integral over a is trivial; as long as we note that the integral over the angular variable a is normalized to 1, we may proceed with integrating.Then the relevant integral is: ) where I 0 is the modified Bessel function of the first kind of order zero.There are no factors of a left in the integrand; as long as π 0 da is normalized by a factor of π, we may consider the previous integral to have simultaneously integrated over both σ and a.
We now seek to integrate over θ.First, we collect the relevant terms and group the coefficients for simplicity's sake.We have: Our integral over θ thus becomes: v 1 e v 2 +v 3 cos 2 θ+v 4 sin 2 θ I 0 (v 5 cos θ) sin 2θ sin 2 θdθ, (A.10) where I θ is the integral over θ in I from eq. (A.1) and v i denote the grouped coefficients.
We now observe that we may rewrite sin 2θ as 2 sin θ cos θ.Then our integral becomes: Setting x = cos θ, our integral becomes: We find the Taylor expansion of e (v 3 −v 4 )x 2 : Now, we expand the Bessel function in series.We find that: Putting everything together, we arrive at: We now seek to integrate over β.Before we start, we first examine the combinations v 2 + v 4 and v 3 − v 4 : We note that v 1 = sin 2β sin 2b, which means we can repeat the process of rewriting sin 2β as 2 sin β cos β, but absorbing cos β behind the derivative instead.Then we rewrite the expression above as: Once again, we expand e w 2 z using the Taylor series: Then our integral becomes: This is a hideous series, and we would be forgiven for thinking that we should define a new index that matches 2(2k + p − m − g).But q does the job, if more subtly, and so we 3) exp Tr ĀU AU † + Tr BU BU † , B = b1 b2 b3 .