Higher spin wormholes from modular bootstrap

We investigate the connection between spacetime wormholes and ensemble averaging in the context of higher spin AdS$_3$/CFT$_2$. Using techniques from modular bootstrap combined with some holographic inputs, we evaluate the partition function of a Euclidean wormhole in AdS$_3$ higher spin gravity. The fixed spin sectors of the dual CFT$_2$ exhibit features that starkly go beyond conventional random matrix ensembles: power-law ramps in the spectral form factor and potentials with a double-well/crest underlying the level statistics.


Introduction
The mechanism by which quantum information escapes an evaporating black hole is one of the deepest mysteries of modern theoretical physics. Over the past few years, there has been significant progress on this front that reproduce a unitary Page curve from semi-classical gravity path integrals [1][2][3]. A crucial ingredient in the analysis involves the inclusion of wormhole saddles that interpolate between regions connected by an entanglement cut. In a similar vein, Euclidean wormholes can also connect two separate boundaries. In the context of AdS/CFT, the existence of such wormhole solutions lead to a loss of factorization in the observables of, what is apparently, a direct product of CFTs [4,5].
A promising way out of this conundrum is to formulate versions of holography in which the dual CFT isn't a single theory but an ensemble average of theories. This idea finds its origins in the context of spin-glasses where the effective description emerges from a disorderaverage over Hamiltonians. The partition function is then given by the mean of the partition functions of the ensemble, Z(β) , while the non-vanishing fluctuations or higher moments, Z(β 1 )Z(β 2 ) · · · , encode the wormhole amplitudes of the bulk dual. Such a construction is largely motivated by the fact that topological JT gravity in 2d is precisely dual to an ensemble of random matrices [6]. In one dimension higher, similar ideas have been developed which demonstrate that pure AdS 3 gravity shares common features with random matrix theory (RMT) and perturbative U (1) D × U (1) D Chern-Simons (CS) theory is dual to a theory of D free bosons averaged over Narain moduli [7][8][9][10][11][12][13].
This brings us to a natural question: how fundamental is the notion of ensemble averaged holography? Ultimately, one would like to depart from a semi-classical gravity approximation and understand whether ensemble averaging makes sense in string theoretic constructions of AdS/CFT. String theory in AdS 3 with a single unit of NS-NS flux has been shown to be exactly dual to the symmetric orbifold of T 4 [14]. In this setup, it has been demonstrated that wormhole partition functions do factorize and, therefore, the averaging operation is unnecessary at the free orbifold point [15]. We are then inclined to ask: where does the averaged description break down?
As we lack fundamental principles at this point, it is valuable to explore whether ensemble averaging can be embedded into more general settings and see what lessons these situations can offer. In this work, we explore this possibility in the framework of higher spin AdS 3 /CFT 2 . Theories of massless higher spin fields describe the leading Regge trajectory of string theory in the tensionless limit. These theories are grown-up versions of classical (super)gravity theories and are more tractable than full-fledged string theories. As the symmetries of the gravity theory get enhanced beyond diffeomorphisms, we lose traditional geometrical notions such as horizons and geodesics. The CFT duals, described by coset constructions, have additional higher spin conserved currents that enlarge the chiral algebra to W ∞ [16]. The coset models are however rational CFTs and do not possess the features of sparseness or chaos which are central to describe black hole microstates. We shall therefore focus on irrational CFTs with W N symmetries (with c > N − 1) which are dual to finite tower higher spin fields in AdS 3 , often dubbed pure higher spin gravity [17].
The key object we consider in this paper is the partition function of a Euclidean wormhole in higher spin gravity. This wormhole connects two spacetime boundaries which are tori. As the precise details of the CFT dual are a priori unclear, we employ modular bootstrap techniques (along with some well-informed assumptions from holography) to arrive at the partition function. In particular, we adapt the techniques of [18] to the case where the CFT has W N symmetries instead of Virasoro. The analysis in this work, therefore, provides a concrete realization of the wider applicability of the methods to bootstrap ensembles. The wormhole partition function takes the form of a Poincaré series along with some prefactors encoding contributions from zero-modes and W N descendants. We shall see that the wormhole amplitude appropriately generalizes the pure gravity case and, at first sight, has a form very similar to the Narain moduli average of (N − 1) bosons. The zero-modes and descendant contributions turn out to be the same as the Narain averaged case but the Poincaré series itself is slightly different. This leads to drastic differences in the spectral correlations.
We dissect the wormhole amplitude further by Fourier transforming to sectors of fixed spin. A useful quantity that captures statistics of energy eigenvalues is the spectral form factor: Z(β + it)Z(β − it) . This quantity can be obtained from the correctly projected wormhole amplitude upon analytic continuation. We find that at late times the spectral form factor of a fixed spin sector strikingly displays a power-law ramp ∼ t N −1 , in contrast to the linear one for RMT or the pure gravity/Virasoro case (for N = 2). Although we haven't tracked down the species of RMTs that mimic this behaviour, the faster ramp ties in well with earlier findings that the irrational W N CFTs violate the bound on chaos, under certain approximations [19]. In the non-chaotic N → ∞ limit, where the theory is described by a 't Hooft limit of rational coset CFTs [16], the ramp might be expected to show exponential behaviour similar to integrable fermion models [20,21].
The pair correlation functions of the spectral densities exhibit some novel properties for the higher spin case. We evaluate this quantity directly from the wormhole amplitude using an inverse Laplace transform and, also independently, using the method of resolvents. For even N , we find long-range correlations between the eigenvalues. Whereas for odd N , the spectral correlations turn out to be short-ranged and they localize around delta-function singularities. If at all a random matrix description exists for this, the potentials for the eigenvalues should have some regimes of attraction -we verify this in a toy example. These properties are markedly different from the pure gravity counterpart and the lower-dimensional case of JT gravity. It is undoubtedly imperative to ask whether two-dimensional higher spin gravity, described by a topological BF theory [22][23][24][25], also has these properties in its spectrum. We do not address this question in this paper, hoping to return to it in the near future. This paper is organized as follows. In §2 we evaluate the wormhole partition function in higher spin gravity using modular bootstrap. This section contains an outline of the bootstrap method and the ingredients of W N CFTs we need for the analysis. We study the spectral statistics of wormhole partition function in §3 -this constitutes finding the spectral form factor, the pair correlation function of the spectral density and the potentials describing the level statistics. We conclude in §4 and discuss some avenues for future research. The appendices contain some identities of Bessel functions, additional technical details and consistency checks.

Partition function of the Euclidean wormhole
In the context of AdS 3 /CFT 2 , Euclidean wormholes have been studied for the pure gravity case in [18,26]. The bulk topology of the 3d Euclidean wormhole is T 2 × I, see Fig. 2.1. The two boundaries of the wormhole are given by two distinct tori, that are connected via the bulk geometry.
The partition function of the wormhole (often referred to as the 'wormhole amplitude') can be obtained from the gravitational path integral using a constrained instanton approach, and this method has been further systematized to higher dimensions [27,56]. In hindsight, it has been realized that the wormhole amplitude can be bootstrapped by imposing modular constraints arising from the boundary tori. However, in this method, it is not just modular invariance alone that fixes the amplitude. Other essential inputs -like smoothness, topological considerations, boundary orientation and charge conservation -have bulk origins and play a key role in determining the partition function.

The modular bootstrap procedure
In this subsection we review the steps involved in the modular bootstrap procedure [26]. It can be seen from Fig. 2.1 that the tori, living at the boundaries of T 2 × I, have no relative Dehn twists and are oppositely orientated with respect to each other (i.e. the outward normals point in opposite directions). This feature implies that modular transformations act oppositely on the tori. We can define a double moduli preamplitude,Z(τ 1 , τ 2 ), that obeys the following invariance constraint:Z (τ 1 , τ 2 ) =Z γτ 1 , γ −1 τ 2 , γ ∈ PSL(2; Z). (2.1) We have suppressed anti-holomorphic dependence to simplify the notation. The action of γ and γ −1 denote the simultaneous modular and inverse modular transformations of τ 1 and τ 2 For future reference, we note that the S-modular transformation is τ → −1/τ , while the T-modular transformation is τ → τ + 1. These two transformations generate the modular group SL(2, Z).
Next, the preamplitude is proportional to the moduli space volume form, V 0 , which arise in the bulk from the zero-mode contributions dictating the relative twist between the two tori. This is a physical effect, henceZ is imbued with this contribution. Furthermore, as the two tori are the boundaries of the same connected bulk, charge conservation constrains the CFTs living on the two tori boundaries to have the same spectrum of primaries. Requiring bulk-smoothness also keeps the conformal dimensions above the BTZ threshold. In the momentum representation, h − c−ccurr 24 = k 2 4 (with h representing the conformal dimension, c the central charge and c curr the number of conserved currents), this implies that 0 ≤ k ≤ ∞. With these inputs, one arrives at an useful ansatz for the preamplitudẽ where, χ k (τ ) is the CFT character. The details of the character depend on the chiral algebra of the CFT which is also the asymptotic symmetry algebra of the bulk theory. The bootstrap constraint (2.1) is sufficient to determine the distribution of primaries, ρ(k,k), upto an overall normalization. Once this is obtained, we plug it back into (2.3) to evaluateZ(τ 1 , τ 2 ). Note that in the above ansatz we have implicitly assumed that the CFT in question is irrational, i.e. we have an infinite number of primaries owing to modular invariance. Furthermore, the character appearing in (2.3) will turn out to be non-degenerate characters of the chiral algebra which have no restrictions coming from null states.
The quantityZ(τ 1 , τ 2 ), however, is not the full wormhole amplitude yet. It misses instances where only one of the two tori gets modular transformed. From the bulk point of view, these are distinct and allowed physical configurations. Therefore, the full partition function should involve a sum over them. Such configurations are generated by γ acting on one of the torus moduli. We finally end up with The sum above is over an infinite number of modular images and it is a priori unclear whether the result is convergent. The convergence depends on the detailed structure of the preamplitudeZ(τ 1 , τ 2 ) itself. We shall return to this point below in Sec 2.4. The modular sum (2.4) is similar to a Poincaré series and, given (2.1), it is invariant under independent modular transformations. This can be seen as follows The partition function (2.4) can also be expressed as a Fourier sum (or q-series) and this naturally projects states onto fixed spin sectors. The BTZ threshold k ≥ 0 then transforms into the spinning BTZ threshold, which for spin s, keeps the energy above the extremality bound E s ≥ 2π |s| − ccurr.
12 [17,54,55]. 1 Before we carry out the procedure outlined above for irrational CFTs with W N symmetries, we describe some essential ingredients that will be useful.

Some ingredients of W N CFTs
For W N CFTs, the symmetry algebra is generated by modes of the stress tensor and the conserved higher spin currents, W (s) m for 2 ≤ s ≤ N (see e.g. [28,29] for reviews). The irrational regime corresponds to the value of the central charge being larger than the number of conserved currents, c > N −1. In what follows, we shall require the characters of non-vacuum primaries on the torus. These are given by Here, η(τ ) is the Dedekind eta-function and we have used a Liouville-like parametrization for the conformal dimension. As usual, these characters contain the contribution of left/right moving descendants of the primary state, (W (s 1 ) −2 ) k 2 · · · |h . The case with N = 2 reduces to the usual Virasoro CFTs. The lightest primary states in irrational W N CFTs scale with the central charge; this fact was found using unitarity constraints and modular bootstrap in [31].
For carrying out the bootstrap procedure for the wormhole partition function, we shall also need the fusion kernel S kk for the following S-modular transformation Note that the above relation is somewhat non-standard due to the presence of additional (−iτ ) factors; the origin of these factors is the moduli space volume, V 0 , that we will encounter momentarily. We can simplify the above relation by using explicit expressions for the characters (2.6) and the S-modular transformation of η(τ ) Our task is to extract S kk . We multiply both sides by e For the LHS of (2.8), we expand the exponential, e − πik 2 2τ , for small k and integrate over x term-by-term. The final result can be resummed into a modified Bessel function: Comparing the last two equations, we obtain the fusion kernel to be For N = 2 this becomes S kk = πk J 0 (πkk ), in agreement with the Virasoro case considered in [18]. We remark at this point that one can parametrize the conformal dimensions of a W N CFT in a manner similar to Toda theories, where the momenta k of vertex operators live in the root lattice of sl(N ) [28,29]. An analogous fusion kernel can be derived which would be a function of the (N − 1) momenta. However, we shall not require that representation for our present purposes and (2.11) will be sufficient.

Evaluating the wormhole amplitude
In this subsection we evaluate the Euclidean wormhole partition function using the modular bootstrap approach. As explained in §2.1 we start by determining the preamplitude and then sum over modular images.

Fixing the preamplitude
In order to evaluate the preamplitude we need to determine the density of primaries, ρ(k,k), by imposing the bootstrap equation (2.1). The preamplitude is understood to have its origins via a three dimensional higher spin gravity path integral, with two torus boundaries characterized by complex structures, τ 1 , τ 2 . In the case of pure gravity (or Virasoro CFTs at the boundaries) the path integral is proportional to the moduli space volume form, V 0 , arising from the constrained saddle nature of the wormhole solutions [56]. This volume form can be written down in terms of the zero modes that control various features of the wormhole geometries, including the length and boundary twists, and turns out to be Im(τ 1 )Im(τ 2 ) [26]. For W N theories at the wormhole boundaries, the current algebra OPEs can be realized via the Miura construction which involves (N − 1) free bosons [28,29]. This is a generalization of the single linear dilaton realizing the Virasoro algebra. It is then reasonable to expect that the net contribution of the zero-modes is A first principle derivation of this will follow from analyzing the moduli space field ranges of pure higher spin theory on T 2 × I. In the pure gravity case, this volume specifically arises from the symplectic measure induced by the Chern-Simons term which involves only a single time derivative. In the SL(N, R) case, this will involve (N − 1) fields corresponding to the Cartan directions, which should give the factor in (2.12).
Let us now consider the preamplitude. We can check that the ansatz (2.3) by construction is invariant under simultaneous T and T −1 modular transformations. Therefore, the nontrivial bootstrap constraints arise only for τ → −1/τ , and its inverse. Explicitly, it enforces the following condition where we also used, Im (−1/τ ) = Im (τ )/|τ | 2 . Next, upon using (2.7) in the LHS we arrive at (2.14) Using the explicit form of the fusion kernel (2.11) we obtain We multiply both sides by k J N −2 (πk q ) andk J N −2 (πk q ). Then using the orthogonality relationship of Bessel functions (see eq. (A.1)) we integrate over k ,k and that results in This equality can only be satisfied if the above ratio is a constant. Therefore We have chosen a convenient normalization constant C, such that simplifications occur later. 2 Plugging this back into (2.3) and explicitly evaluate the preamplitude to bẽ Therefore, upto the overall constant, the higher-spin preamplitude is simply the pure gravity result raised to the (N − 1)'th power. We now turn to the sum over modular images of this quantity.

Modular sum
We are now in a position to use the preamplitude (2.18) to obtain the full partition function. As outlined earlier, we need to perform a sum over one-sided modular images of the preamplitudẽ Z. As the Z 0 (τ )'s in (2.18) are modular invariant, the final result for the higher spin Euclidean wormhole partition function is with Z 0 (τ ) defined in (2.18). This equation is one of the key results of this paper. Note that the wormhole amplitude above does not depend on the central charge. Let us now focus on the sum over modular images (2.20) As alluded to earlier, this Poincaré sum will be performed for fixed spin sectors, i.e. we will decompose the wormhole partition function as follows where, s 1 and s 2 denote the specific spin configuration. To this end, we rewrite τ 1 = z 1 + iz 2 , and, τ 2 = w 1 + iw 2 . The integer valued spins, s 1 and s 2 , arise as Fourier conjugate variables to z 1 and w 1 respectively. Since R is invariant under independent modular transformations (and especially the T transformation) on either of the moduli, the Fourier series exists. More explicitly this is Therefore the fixed spin contribution to the wormhole partition function is given by The principal object here,R s 1 ,s 2 , can be expressed using the inverse Fourier relatioñ In order to proceed, we shall exchange the order of the Fourier integrals and the modular sum. However, before we perform the Fourier integrals, it is useful to separate out the PSL(2, Z) summation into a part which involves only the T transformations and a part that involves at least a single S modular transformatioñ T s 1 ,s 2 involves a single sum over integers, while S s 1 ,s 2 , which is more difficult to evaluate, involves summing over co-primes. Fortunately the universal low energy contribution to Z(τ 1 , τ 2 ) comes only from the part involving T transformations alone and we present it here For the first integral, we may join all the summations by changing variables inside the integral of each summand, w 1 → w 1 + n. We thereby get rid of n dependence inside the integrand at the expense of an extended contour of integration, (−∞, ∞). Therefore we evaluate (2.27) Taking advantage of the infinite w 1 range, the integrals can be decoupled by changing variables w = z 1 + w 1 , which results in The z 1 integral results in a Kronecker delta that enforces s 1 = s 2 , whereas the w integral furnishes an integral representation of the modified Bessel K function. We finally get We note that when N = 2 the above reduces precisely to the gravity answer. Details about the S transformation part, which we denote with S s 1 ,s 2 can be found in Appendix B. We see, from (B.9), for the case of N = 3, while the exponential suppression in S s 1 ,s 2 at low temperatures is the same as in T s 1 ,s 2 , the former is polynomially suppressed while the latter is polynomially enhanced. Therefore, as stated before, we can see that low temperature behaviour is dominated by T s 1 ,s 2 . Putting everything together, the fixed spin contribution to the wormhole partition function takes the form

Comparison to other ensembles
Now that we have obtained the partition function of the Euclidean wormhole in higher spin gravity (2.19), it is useful to contrast the result with other ensembles considered in the recent past.
Let us first consider wormholes in pure gravity, which correspond to irrational CFTs with Virasoro symmetry at the boundaries. The result for the wormhole amplitude was derived in [18,26] and is given by The higher spin result (2.19) along with (2.20), agrees with the above upon setting N = 2.
The overall constant 1/2π 2 is fixed using the JT gravity limit; unfortunately, we do not have an analogous 2d higher-spin gravity computation at our disposal to fix the overall constant C in (2.19). Note that the Poincaré series in (2.31) does not converge and needs to be suitably regulated [26]. The higher spin amplitude, on the other hand, has convergence built in. In fact, the 'zeta-function' used in [26] to regularize the Poincaré sum of (2.31) is exactly the same as the one appearing in the higher spin wormhole partition function (2.20). 3 A more interesting comparison is with the wormhole partition function in 'perturbative' U (1) D × U (1) D Chern-Simons theory. This theory is the bulk dual of D free bosons averaged over the Narain moduli space [7,8]. The wormhole amplitude can be obtained from the connected piece of the averaged genus-2 partition function, with a diagonal period matrix, Ω = Diag(τ 1 , τ 2 ). The result is [8] Z Narain (τ 1 , τ 2 ) = Z 0 (τ 1 ,τ 1 ) D Z 0 (τ 2 ,τ 2 ) D γ∈PSL(2,Z) This amplitude is also the result of the modular bootstrap problem for U (1) D × U (1) D chiral algebra on the boundaries, along with the moduli space volume set as V 0 = 1 [18]. The Z 0 prefactors in (2.32), which count the zero modes and descendant states, agree with (2.19) upon setting D = N − 1. This fact is isn't a mere coincidence since N − 1 free bosons form a realization of the W N algebra.This realization is the Miura transformation and it works at arbitrary central charge [29]. 4 Curiously enough, the Poincaré sum in (2.32) is also of the same form appearing in the higher spin amplitude (2.20). However, for D = N − 1, each term of the above sum is the square-root of the one appearing in the higher spin amplitude (2.20).
For the Poincaré sum of the Narain ensemble (2.32), we can use D/2 = N − 1 in (2.29). to get the fixed spin sector contribution For z 2 = β 1 and w 2 = β 2 the result agrees with [26, first line of (3. 20)]. The factor (z 2 w 2 ) D/2 cancels out exactly with factors from Z 0 (τ 1 ,τ 1 ) D Z 0 (τ 2 ,τ 2 ) D ; this leads to an absence of a ramp in the spectral form factor.
Despite the minor differences in the wormhole amplitude, the Narain and higher spin ensembles have very different energy eigenvalue statistics. To analyze this in detail, it is beneficial to extract the spectral density correlations from Z hs (τ 1 , τ 2 ). This is the topic of the next section.

Spectral statistics from the wormhole amplitude
In a quantum chaotic system, it is universally expected that there is repulsion amongst energy eigenvalues [32,33]. It has recently emerged that the Euclidean wormhole encodes analogous features for black hole microstates. The objective of this section is to quantify these features for the higher spin wormhole and understand the details of the dual ensemble description. We shall do this by studying the spectral form factor and the spectral density 2-point function.

The spectral form factor
Random matrix theory captures very universal features of quantum chaotic systems. In this context, the spectral form factor (SFF) serves as a useful tool towards diagnosing quantum chaos. The SFF is defined as follows The factors, Z(β 1 ) and Z(β 2 ), denote the partition functions with inverse temperatures β 1 and β 2 which are analytically continued to β + it and β − it. The SFF is a simpler proxy for Lorentzian two-point correlation functions, i.e. the SFF depends only on the details of spectrum and has the dependence on matrix elements of operators stripped off. The product Z(β + it)Z(β − it) probes the discreteness of the spectrum, and this aspect is realized by considering the late time behaviour of the quantity [39]. Furthermore, in chaotic systems the SFF exhibits an universal profile [34] -this consists of an initial dip, followed by a linear ramp and then a plateau at very late times. In the above definition the average is taken over ensemble irrational CFTs with W N symmetries. 5 On the other hand, in random matrix theory, an integral over matrices plays the role of averaging. The linear ramp in particular, is related to the spectral rigidity of random matrix eigenvalues. The averaging operation also leads to a loss of factorization. In our case the non-factorization is geometrically realized via the three dimensional wormhole geometries, and the connected piece of the 2-point function (3.1) is built into the bootstrapped wormhole amplitude (2.19).
In order to extract the SFF however, one needs to focus on a superselection sector of the wormhole partition function. 6 In our context this corresponds to first focusing on the contribution of primaries of Z(τ 1 , τ 2 ) within a fixed spin-sector. We can obtain this directly by stripping off the descendant counting functions from the spin-decomposed amplitude Z s 1 ,s 2 (τ 1 , τ 2 ) in equation (2.30) The superscript P indicates the contribution from primary states. The SFF can now be obtained by analytic continuation, This object encodes correlations of energy eigenvalues across sectors of fixed spin. We focus on the low-temperature regime. This corresponds to the region in parameter space where energies are close to the threshold energy of the spin-s BTZ black hole Let us write down the Euclidean result first. We take the boundary tori to be rectangular and set the modular parameters as τ 1 = iβ 1 and τ 2 = iβ 2 . At low-temperatures (β 1,2 → ∞) the dominant contribution arises from the sum over T -modular transformations. Taking into account the temperature dependence from the Z 0 prefactors (sans the descendant contributions) and using large argument approximation of the Bessel function (A.2), we have the following 5 The details of the averaging of this is unknown at the moment. However, we imagine the wormhole amplitude can be reproduced by a suitable averaging over the moduli of the conformal manifold, along the lines of what has been done for free theories [7,8]. 6 We thank Kristan Jensen for emphasizing this point.
As a consistency check, this result agrees with pure gravity case for N = 2 [26, eq (3.27)]. It is now straightforward to perform the analytic continuation to obtain the SFF (3.3). We obtain Therefore, at late times we have a power-law ramp t N −1 within a given fixed spin sector. This generalizes the pure gravity case (N = 2) for which the ramp is linear.
It is worthwhile to compare the above results with the Narain ensemble (see [35] for an exhaustive study). The wormhole partition function at low temperatures is, from (2.33) where, E s = 2π(|s| − D 12 ). Upon setting β 1,2 = β ± it, the above expression does not contain any time dependence and, therefore, the ramp is absent in the SFF. In a sense, this conclusion is well expected as the CFTs being averaged over are free theories and they do not exhibit chaotic properties. 7 The above features are for the connected piece of the SFF, which is given by Euclidean wormhole amplitude. There is also a disconnected component in the SFF, Z(β+it) Z(β−it) , which gives a decay and dictates early time behaviour. In the context of RMT, this decay is universal and depends only on the symmetry class of the model, for instance it is t −3/2 for Gaussian ensembles, and, t −1/2 for Wishart-Laguerre ensembles [36]. Note that the SFF starts its life from Z(β) 2 . By the time the disconnected piece decays the rise coming from the connected piece begins to take over -this is also predicted by RMT. This leads to a clear transition between the dip and the ramp. Writing down the connected SFF in spectral decomposition we note that at late times, the randomness of the energies favor only small E 1 − E 2 due to phase cancellations. This makes the ramp sensitive to nearest level distributions, which is encoded in the connected density-density correlator, ρ(E 1 )ρ(E 2 ) . For one-matrix models, this can be calculated in terms of the eigenvalue correlation kernel, K(E 1 , E 2 ), which can determine any arbitrary joint probability distribution [36]. For large random matrices, when E 1 − E 2 is small, K(E 1 , E 2 ) is given by an universal function known as the sine-kernel [37]. Using this sine-kernel, one can show that the connected SFF will start growing (as a ramp) at late times; an array of one-matrix model examples has been reviewed in [38]. The dominant behaviour of the ramp is always linear, whilst there are non-linear corrections that become important at later times.
The linear ramp of the SFF arises from the Fourier transform of the divergent contribution to the density-density correlator, ρ( [39]. Since for higher spins we find a power-law ramp, t N −1 , we expect a stronger divergence in the spectral density correlator. In the next subsection, we explore this expectation in detail.

Pair correlation function of spectral densities
We now want to extract the two-point function of spectral densities in a specific spin-sector. As in (3.8) this two-point function is related to the wormhole amplitude in the following manner: The density-density correlator can either be obtained from the discontinuities of the double resolvent or by directly double inverse Laplace transforming the wormhole partition function. In this section we use the latter method, and discuss the former in Appendix C. We focus on the low temperature regime, where we observed a power-law ramp in the SFF. The partition function is given in (3.5). The expression for ρ s 1 (E 1 )ρ s 2 (E 2 ) is (3.10) where, we have defined the following to lighten the notation Let's consider the inverse Laplace transform (ILT) wrt β 1 first. This is (3.12) The integral can be performed by expanding the E 1 independent factor as a power series in β 1 /β 2 and then integrating term by term. The result is Here, (a) n = Γ(a + n)/Γ(a), is the Pochammer function. We now need to inverse Laplace transform wrt β 2 . The details turn out to be quite different depending on whether N is even or odd. So let's consider these cases separately.

Even N
The density-density correlator can be obtained from (3.13) as We write U(β 2 ) as in the first line of (3.13) and perform the ILT term-by-term, i.e. the integral transform acts on powers of β 2 . The basic identity we can use here is the following which makes sense only if ν isn't a positive integer. Upon resumming the integrated terms, we finally have the following result for the pair correlator In Appendix C, we verify the above formula using the method of resolvents for N = 2, 4. For N = 2 we get the result (C = (2π 2 ) −1 for N = 2 in our conventions) This agrees with the pure gravity case. We observe that the two-point function (3.16) displays long-range eigenvalue attraction for odd N/2 and repulsion for even N/2.

Odd N
For odd N (or half-integer N/2), we can no longer use the identity (3.15) to perform the second inverse Laplace transform. Instead, we can rewrite the hypergeometric function appearing in (3.13) in terms of Laguerre polynomials [40]. For N = 2P + 1 we have the following identity The inverse Laplace transform (3.14) can then be carried out by acting on each term of the quantity in square brackets above. For a given value of P , this is a finite number of terms. The ILT acts on integer powers of β 2 in the following manner +i∞ −i∞ with E 21 = E 2 − E 1 . The pair correlation function can then be written as This shows that the result is a linear combination of derivatives of the Dirac delta function. Unlike the case for even N where we have long-range correlations, we see that the pair correlation function localizes at the contact-term singularities. As examples, we have the following for W 3 and W 5

Is there a matrix model description?
We would now like to gain a better understanding of the spectral correlations. For simplicity, we shall focus on the even N case for which we obtained the spectral density 2-point function to be (3.16). In the context of Gaussian unitary ensembles (GUE), a similar behaviour of the pair correlation function is seen, ρ(λ 1 )ρ(λ 2 ) ∼ 1/(λ 1 − λ 2 ) 2 . This arises from the 1d Coulomb repulsion between the eigenvalues, V (λ i , λ j ) = − log |λ i − λ j |. Since we observe a slightly generalized spectral density correlator for W N CFTs, it is worthwhile to ask: what matrix ensembles or potentials V (λ i , λ j ) do they correspond to? In what follows, we shall consider a simplified version of this problem and glean the lessons.
To keep the discussion self-contained, let us recall the well studied GUE case first. The matrix integral can be written in terms of the eigenvalues after diagonalizing the matrices and performing a unitary change of basis. It reads Here, L is the size of the matrix and ρ(λ) is the unit normalized eigenvalue density. The logarithimic potential is essentially the (exponentiated) Vandermode determinant that arises from the Jacobian while changing integration variables from matrix elements to eigenvalues.
Our next step is to derive the pair correlation function using (3.24) as the starting pointcf. [39,41]. The quadratic fluctuations about the saddle at large L is where, δρ(λ) = ρ(λ) − ρ saddle (λ). Fourier transforming these density fluctuations as δρ(λ) = du 2π e iuλ δρ(u) and performing the integrals, we get The λ 1,2 integrals in (3.26) are performed by changing variables, r = λ 1 − λ 2 , and then using the standard identity for the Fourier transform of log |r|. 8 From this we can find the propagator and revert back to λ-space (3.27) The form of this 2-point function is analogous to the N = 2 or Virasoro case (3.17), for small energy separations.
We would now like to reverse engineer the potential for eigenvalue interactions from a spectral density correlator, e.g. given (3.27), we want to reconstruct (3.25). A bare-bones version of (3.16) that retains the dependence on eigenvalue differences and the overall sign is (3.28) where, N is an even integer. For N = 2 this reduces to the GUE case (3.27). The inverse Fourier transform of (λ 1 − λ 2 ) −N , with respect to conjugate variable u, is given by (−1) N/2 |u| N −1 .
We can then write down the action of the density fluctuations in the Fourier space In the eigenvalue space, if the eigenvalue potential is given by, V (λ 1 − λ 2 ) = V (r), then we also have the following analogue of (3.25) The two equations above describe the same quantity and should be equivalent. In order to extract V (r), we introduce the Fourier transforms as in the previous paragraph, δρ(λ j ) = du j 2π e iuλ j δρ(u j ). The integral over λ 2 yields a δ(u 1 + u 2 ) factor, cf. (3.29), and we are left with the following equation for V (r) du e iur V (r) = 1 |u| N −1 . (3.31) Finally, we find that the potential takes the form This is the generalization of the V (r) = − log |r| behaviour of the GUE case. Here a N and b N are numerical coefficients and their exact N -dependence isn't particularly illuminating. However, the details are such that the potentials take characteristic shapes -see Fig. 3.1. Quite strikingly, when N/2 is even we find a double-well shape. This is mostly attractive, with two degenerate basins located away from the origin at values that increase with N . This implies that most eigenvalue differences are confined in the wells and their vicinity. On the other hand, we see that when N/2 is odd, the potential resembles an inverted double-well. Therefore, it is mostly repulsive with a weak attraction/confinement at very small eigenvalue differences. 9 For N > 2 we get "double-well" or "double-crest" type potentials, depending on whether N/2 is even or odd.
These features are in sharp contrast to the GUE matrix models. Although we considered a fairly simple form of the pair correlation function (3.28), the potential (3.32) that gives rise to such behaviour is quite different. At this moment, it is unclear to us whether suitable matrix models (or deformations thereof) can give rise to these interactions between the eigenvalues.

Discussion
In this paper we analyzed Euclidean wormholes in 3d higher spin gravity. We used modular bootstrap methods, generalizing the case of wormholes in pure gravity [18]. The wormhole amplitude gives the connected part of Z(τ 1 )Z(τ 2 ) of a suitably ensemble averaged, irrational, W N CFT. This amplitude captures the eigenvalue statistics of black hole microstates. We observed that the spectral form factors have a power-law ramp (∼ t N −1 ), and, the inferred eigenvalue dynamics exhibits many interesting features, including strong and weak attractive behaviours, as well as localization, which depend on the value of N . These features are novel, and perplexing at the same time -they have not been seen previously for pure 3d gravity [18,26] or in the lower dimensional case of JT gravity [6]. A drawback of our analysis is that the twist zero-mode volume, V 0 of equation (2.12), was argued based on symmetries of the dual CFT. It would be desirable to have a derivation of the same from the higher spin gravity perspective.
At very late times, the spectral form factor saturates to a plateau of height Z(2β); this follows from the very definition of the SFF and is universal. However, in order to reproduce this from the wormhole amplitude we require to take into account contributions beyond the low-temperature regime. In particular, one needs to explicitly evaluate the Kloosterman sums that include the S-transformed images of the preamplitude (see (B.9) for the N = 3 case). Reproducing the value of the plateau (which should equal the torus partition function Z(2τ )) will be a good consistency check of the full bootstrapped amplitude. It is to be noted that this remains an open problem even for the pure gravity case [26].
Our results point towards very interesting features for the matrix ensembles that can produce the pair correlation of spectral densities. Typically random matrices exhibit eigenvalue repulsion. However, when N is a multiple of 4, we found that the eigenvalues show attractive behaviour. Even though such an effective attraction is unexpected in quantum chaos discussions, there are some rare examples like [20]. Furthermore, some integrable systems show an exponential ramp in the SFF; an example is the q = 2 SYK model [21,43]. Therefore, the power-law ramp for the irrational W N CFTs interpolates between the chaotic case (with a linear ramp) and the fully integrable one. It would then be valuable to study these higher spin ensembles further, even with the goal of understanding generic eigenvalue dynamics. With this in mind and inspired by recent developments in JT gravity [6,[44][45][46][47], it will be fascinating to translate the analysis here in terms of a matrix model, if at all it exists.
One can imagine more general wormhole backgrounds in higher spin gravity which contain higher spin charges, in the same spirit of the higher spin black holes constructed in the past [48]. Studying the corresponding wormhole amplitude (in the grand canonical ensemble of non-zero higher spin chemical potentials) will reveal the statistics of the higher-spin charges. Unfortunately, there are technical obstacles in carrying this out; even for the case of a single torus boundary these partition functions are not known only perturbatively [49][50][51][52][53]. Furthermore, the modular properties of the partition functions are not clearly known which is a hindrance to the bootstrap method employed here.
The results of this work and that of [54] provide information about the one-and two-point functions Z(τ ) and Z(τ 1 )Z(τ 2 ) for the low-temperature regime of pure higher spin gravity in 3-dimensions. This is the near-horizon regime of near-extremal black holes in which an AdS 2 -throat appears [55]. For the case of higher spins, the 2d gravity description is provided by a topological BF theory [22][23][24][25]. It would be reassuring to derive Z(τ ) , Z(τ 1 )Z(τ 2 ) and higher point correlators which translate to BF theory on the disk, double-trumpet and geometries with multiple boundaries respectively. A related question is: how does topological recursion generalize for 2d BF theory? Given the topological nature of BF theory it is very likely that a recursive machinery will exist that would fruitfully allow the evaluation of partition functions of n-boundary wormholes in a genus expansion. Alternatively, one can hope to obtain the n-boundary amplitude by generalizing techniques of Liouville gravity, developed in [57], to the case of Toda gravity. These amplitudes would enable the determination of higher moments of the spectral densities ρ(E 1 )ρ(E 2 )ρ(E 3 ) · · · . Relatedly, it can be investigated to what extent higher moments of the spectral density are fixed/constrained by the first few moments -this is an incarnation of the truncated moment problem.

A Bessel function identities
We list a couple of properties of (modified) Bessel functions, J ν (x) and K ν (x), that have been useful in our analysis. 2. At large arguments the modified Bessel function has the following behaviour

B Further details on the Poincaré sum
We start with the Poincare series in the Fourier space indexed by spins Using the properties of P SL(2, Z), it can be shown that γ can be decomposed into where, T generates τ → τ + 1, and all other variables are integers. The first part which involves only T transformations, give rise to T s 1 ,s 2 , which is written down explicitly in (2.26) and evaluated in (2.29). The second part of γ gives rise to S s 1 ,s 2 . The integer c is positive, while d ∈ (Z/cZ) * . This means, d is an integer such that, 1 ≤ d ≤ c − 1 is co-prime with respect to c. dw 1 e 2πiz 1 s 1 +2πiw 1 s 2 (Im(τ 1 )Im(T n γ c,d T m · τ 2 )) |τ 1 + T n γ c,d T m · τ 2 | 2 3) The effect of the multiple actions by the T transformations, present in S s 1 ,s 2 , can be dealt with in the same way, as was done with T s 1 ,s 2 . Namely, via variable redefinitions, we absorb into z 1 and w 1 , these translations at the expense of extending the contour of integrations to the entire real line. As a result we obtain These integrals can be performed in the complex z 1 and w 1 planes, by closing the contours appropriately. Either of the integrals present us with two N − 1 order poles, one in the UHP while another one is its reflection into the LHP. Therefore, we need to evaluate a (N − 2)-th order derivative to find the residue. This complicates the expressions, thus we focus on the case with N = 3. For notational simplicity we now denote d by d, For the w 1 integral we find two poles of order two. They are symmetrically placed with respect to the real w 1 axis. For the case, s 2 > 0 the UHP pole contributes, whereas for s 2 < 0, the residue contribution arises just from the pole located at w * = − b + i aw 2 + (d + i cw 2 )z a + c z .
The location of the pole can be shown to be in the LHP since both w 2 and z 2 are non-negative. Without loss of generality, we choose the case, wherein we close the contour via the LHP, this contribution gives πd 4 e 2πs 2 (w 2 +(cw 2 −id)(b+dz)) 1+bc+cd z 2 (d 2 z 2 + w 2 ((1 + bc + cdz 1 ) 2 + c 2 d 2 z 2 2 )) It is important to note that in the low temperature regime, apart from the exponential suppression the polynomial suppression goes as 1/(β 2 β 1 ) 2 . The computation for other values of N can be done in a similar manner.

C Density correlators from the resolvent
The resolvent is a useful quantity in the context of matrix models, whose discontinuities have information about densities and their correlators. The single resolvent has information about the density of states, while the double resolvent encapsulates the pair correlation function of spectral densities. The double resolvent is defined as . (C.1) We can obtain the density-density correlator from its double discontinuities in the complex E 1 , E 2 planes Here, P(x) denotes the principal value. It then follows that where we have used definitions of (3.11). Once again the odd and even cases of N require separate treatment due to very similar reasons as in the inverse Laplace transform method.
Here we present a few even cases, and show that the answers agree with the ones obtained directly using inverse Laplace transforms. This provides a useful consistency check of the results.