Bootstrapping the Spectral Function: On the Uniqueness of Liouville and the Universality of BTZ

We introduce spectral functions that capture the distribution of OPE coefficients and density of states in two-dimensional conformal field theories, and show that nontrivial upper and lower bounds on the spectral function can be obtained from semidefinite programming. We find substantial numerical evidence indicating that OPEs involving only scalar Virasoro primaries in a c>1 CFT are necessarily governed by the structure constants of Liouville theory. Combining this with analytic results in modular bootstrap, we conjecture that Liouville theory is the unique unitary c>1 CFT whose primaries have bounded spins. We also use the spectral function method to study modular constraints on CFT spectra, and discuss some implications of our results on CFTs of large c and large gap, in particular, to what extent the BTZ spectral density is universal.


Introduction
Enormous progress in the conformal bootstrap program has been made in recent years based on semidefinite programming . Typically, one aims to bound the scaling dimensions and OPE coefficients of the first few operators in the spectrum based on unitarity and crossing invariance of the 4-point function. Such methods are most powerful in constraining CFTs with simple low lying spectrum, but become less constraining when the spectrum becomes dense.
In this paper we introduce the spectral function method, which allows for constraining not just the gap or the first few OPE coefficients but the distribution of OPE coefficients over a wide range of scaling dimensions. While the method can be applied to CFTs in any spacetime dimension, we focus on two-dimensional unitary CFTs, where the spectral functions are defined by truncating the Virasoro conformal block decomposition of a 4-point function in the scaling dimension of the internal primaries, evaluated at the self-dual cross ratio, or by truncating the Virasoro character decomposition of a partition function in the scaling dimension of the primaries, evaluated at the self-dual modulus. More precisely, consider a scalar 4-point function 1 g(z,z) ≡ φ(0)φ(1)φ(z,z)φ (∞) = Of course, for a compact CFT with a discrete spectrum, f (x) will be composed of step functions. If the CFT is non-compact, then typically f (x) will be a monotonically increasing smooth function that takes value between 0 and 1. We will see that the crossing equation combined with additional assumptions on the spectrum lead to upper and lower bounds on f (x). Numerically, the crossing equation can be utilized by applying to (1.3) linear functionals spanned by the basis ∂ n z ∂ m z | z=z=1/2 , for odd n + m ≤ N . The resulting upper and lower bounds on f (x) (which are rigorous bounds although not optimal) will be denoted f + N (x) and f − N (x). If we make the assumption that the CFT contains only scalar Virasoro primaries, we find that f + N (x) and f − N (x) become closer as N increases, for various values of the central charge c > 1. We conjecture that both converge to the spectral function of Liouville theory, which can be computed by integrating the square of the DOZZ structure constants [22,23] times the Virasoro conformal blocks. Note that this approach can be extended to the 4-point function involving a pair of different primaries, leading to spectral functions that encode the most general structure constants of the CFT.
Convergence of the upper and lower bounds f ± N (x) to the same value f ∞ (x) is related to the completeness of the derivatives of scalar Virasoro blocks in a suitable space of functions. Conversely, this completeness statement implies the uniqueness of the solution to the crossing equations. We propose a numerical test for the completeness and find compelling evidence suggesting that it holds. We moreover obtain numerical approximations to the (conjecturally) unique solution to the crossing equations, which reproduce the DOZZ spectral function with high accuracy. In contrast to the semidefinite methods, this linear approach does not rely on the assumption that the OPE coefficients are real. The linear and semidefinite results above therefore lead us to conjecture that the DOZZ structure constants are the unique solution to the crossing equations for (not necessarily unitary) CFTs with only scalar primaries (of non-negative scaling dimensions) and c > 1.
Interestingly, we find that the bounds on the spectral function f ± N (x) exist for external operator dimensions ∆ φ ≥ c−1 16 ( 3 4 of the Liouville threshold), and converge to a step function when ∆ φ is equal to c−1 16 . When ∆ φ < c−1 16 , the crossing equation cannot be satisfied with only scalar internal primaries, ruling out the possibility of such operators. 2 We will see that all of these are in agreement with the analytic continuation of Liouville 4-point functions.
A caveat in the above uniqueness claim is that we have assumed a non-degenerate scalar spectrum. If degeneracies are allowed, then the operator algebra of Liouville CFT tensored with a topological quantum field theory (TQFT) (or equivalently, a finite dimensional commutative non-unital Frobenius algebra) would also solve the crossing equation. In fact, such a TQFT can always be "diagonalized" by a basis change, and amounts to superselection sectors. We will give partial arguments suggesting that under our assumptions, "Liouville ⊗ TQFT" is the only possibility.
If we further invoke modular invariance, it will turn out that demanding that a unitary CFT contains only primaries of spins in a finite range (s ≤ s max for some finite s max ) implies that the CFT must have a non-compact spectrum with only scalar primaries, and that the spectral density ρ(∆) must be that of Liouville theory, namely This leads us to conjecture that Liouville theory is the unique unitary CFT with c > 1 whose primaries have bounded spins.
The spectral function method also can be applied to modular bootstrap. In this context, we write the torus partition function as (1.5) where χ ∆,s is the Virasoro character associated with a primary of dimension ∆ and spin s and d ∆,s is the degeneracy. The modular spectral function is defined by truncating the Virasoro character decomposition of the partition function (1.6) Once again, upper and lower bounds f ± mod,N (x) can be obtained by acting on the modular crossing equation with linear functionals spanned by the basis (τ ∂ τ ) n (τ ∂τ ) m | τ =−τ =i , for odd n+m ≤ N . In [20] (improving upon [24,25]), an upper bound ∆ mod (c) on the gap in the scaling dimenions was computed numerically as a function of the central charge c. When this bound is saturated, the entire spectrum is fixed by modular invariance, and is determined by the zeroes of the optimal linear functional acting on the Virasoro characters. We will see in examples of small c (between 2 and 8) that under the assumption of maximal dimension gap, f + mod,N (x) and f − mod,N (x) converge with increasing N to step functions, corresponding to the spectral functions of known theories.
For larger values of c, even when the dimension gap is maximized, the convergence of the bounds f ± mod (x) to a sum of step functions is difficult to see numerically, because a good approximation of the optimal linear functional requires larger values of N , and because the step function feature becomes invisible due to an exponentially large spectral density. Nonetheless, for 50 ≤ c ≤ 300, we find empirically that the horizontal average f mod,N (x) of the upper and lower bounds converges rather quickly with N , and the result is in good agreement with the total contribution from thermal AdS 3 and BTZ black hole [26] to the gravity partition function, which results in the modular spectral function Note that this asymptotic spectral function at large c is nontrivial when the dimension x lies in a window of width ∼ √ c around c/6. The agreement with the numerical bounds confirms the validity of the effective field theory of pure gravity in AdS 3 in the canonical ensemble, for temperatures of order 1 in AdS units.
Curiously, BTZ black holes corresponding to operators of scaling dimension ∆ in the range c 12 < ∆ < c 6 never dominate the canonical ensemble, and yet have macroscopic (AdS scale) horizon, provided that ∆− c 12 scales with c. While the naive expectation from effective field theory is that the Bekenstein-Hawking entropy formula should be a valid counting of the microstates of such BTZ black holes, it is unclear to us whether this is a universal property of CFTs with sufficiently large gap. 3 In principle, the modular spectral function bounds at large c should either confirm or disprove such statements. To probe the density of states in the regime ∆ = yc for 1 12 < y < 1 6 and large c would require exponential precision in determining the modular spectral function, which is beyond our current numerical capability. This paper is organized as follows. In section 2 we introduce the spectral function for the scalar 4-point function in a 2D CFT, and explain how to obtain upper and lower bounds f ± N (x) from semidefinite programming. We then specialize to the case where only scalar primaries are present, and demonstrate the convergence of the bounds toward the Liouville spectral function. In section 3, we examine the completeness of scalar Virasoro conformal blocks which would be implied by the aforementioned convergence, and we give numerical evidence that the completeness indeed holds. We then present analytic arguments based on modular invariance that a unitary CFT with c > 1 and Virasoro primaries of bounded spin must be a non-compact CFT with the same spectral density as that of Liouville. This together with the result of section 2 strongly supports the conjecture that Liouville theory is the only CFT with bounded spins. In section 4, we analyze the numerical bounds on the modular spectral function in a number of examples. We conclude with a discussion on the universality of the BTZ spectral density in large-c CFTs with large gaps. 3 Such a universality would in particular require the dimension gap bound ∆ mod (c) to have asymptotic slope 1 12 , namely d∆ mod (c) dc → 1 12 , c → ∞, which is not ruled out by the result of [20] but remains unproven (with no numerical evidence either).
2 Spectral function bounds from semidefinite programming 2.1 A sphere four-point spectral function We begin by considering the conformal block decomposition of the sphere four-point function of a pair of scalar Virasoro primary operators φ 1 , φ 2 of dimensions ∆ 1 , ∆ 2 , (2.1) Here I 12;s is the set of scaling dimensions of spin-s primary operators in the φ 1 φ 2 OPE and C 12;s,∆ = C φ 1 φ 2 O is the OPE coefficient corresponding to the fusion of φ 1 and φ 2 into the primary O with dimension ∆ and spin s. 4 The OPE coefficients are real in a unitary CFT. The conformal block F 12;∆,s takes the form is the holomorphic Virasoro conformal block with external primaries of weight h i and an internal primary of weight h, in a CFT with central charge c. Note that in writing the four-point function this way we have assumed a parity-invariant spectrum. 5 An efficient method for computing Virasoro conformal blocks is Zamolodchikov's recurrence relation [23,27], which we review in Appendix A. It computes the blocks as expansions in the "nome" q(z), defined as Note that as z ranges over the complex plane, q(z) takes value in an eye-shaped region on the unit disc, and the expansion of a conformal block in q converges on the entire unit disc. In the numerical approach, we apply Zamolodchikov's recurrence relation up to a finite depth d q , which generates the correct q-series coefficients up to order q dq . We then truncate the conformal block to this order as an approximation of the exact block.
It follows from the associativity of OPE that the four-point function is crossing symmetric, which amounts to the crossing equation This relation puts highly nontrivial constraints on the spectrum and OPE coefficients of the CFT, some of which were analyzed in [19,21,[28][29][30]. In previous works, one typically either focuses on a limit of the crossing equation in the cross ratio and extracts asymptotic properties of the spectrum, or numerically bounds the scaling dimension and OPE coefficients of the first few operators from the positivity assumption on C 2 12;∆,s . We now introduce a "spectral function" that captures the distribution of OPE coefficients over a range of scaling dimensions of primaries in the φ 1 φ 2 OPE, defined through the conformal block decomposition of the four-point function evaluated at the crossing-symmetric point z =z = 1 2 , truncated on the dimension of internal primary operators: Note that due to the unitarity bound, f (∆ * ) receives no contribution from primary operators with spin s > ∆ * . By definition, obviously, f (∆ * ) is a non-decreasing function that takes value between 0 and 1.
One can place bounds on the spectral function using semidefinite programming as follows. We would like to either maximize or minimize the spectral function subject to the crossing equation expanded around z =z = 1 Note that z = 1 2 corresponds to the nome q = e −π , thus the q-expansion of conformal blocks converges rather quickly at this point. Consider a set of coefficients y 0,0 and y m,n (m + n odd) such that Here Θ(∆ * − ∆) is the step function that takes value 1 for ∆ ≤ ∆ * and 0 otherwise. (∆, s) runs through all possibly allowed values of dimension and spin in the OPE. We could place additional assumptions on the spectrum by restricting the range of (∆, s) in (2.7). For instance, if we are to impose a dimension gap ∆ gap or twist gap t gap , then we have respectively ∆ ≥ max(s, ∆ gap ) or ∆ ≥ s + t gap for the spin-s (non-vacuum) primaries. 6 We shall seek the minimal y 0,0 such that (2.7) holds, which we denote by y min 0,0 . It follows that where we have invoked unitarity by making use of the non-negativity of the squared structure constants, and applied the crossing equation. In other words, y min 0,0 is an upper bound on the value of the spectral function at ∆ * .
Likewise, if we minimize w 0,0 subject to i.e., −w min 0,0 is a lower bound on the value of the spectral function at ∆ * . To obtain these bounds numerically we need to restrict to a finite subset of linear functionals acting on the crossing equation. We will do so by restricting the sums in (2.7) and (2.9) to odd m + n ≤ N ; we refer to N as the "derivative order." The upper and lower bounds on the spectral function derived from the above minimization procedure using linear functionals up to derivative order N will be denoted f + N (∆ * ) and f − N (∆ * ), respectively. While these bounds at every N are rigorous by themselves, the optimal bounds are obtained by extrapolating to the N → ∞ limit.
The numerical implementation of the above procedure is performed using the SDPB package [31], with two practical modifications. Firstly, we will need to truncate the spectrum: while the application of SDPB does not require cutting off the dimension spectrum from above, a sufficiently large but finite truncation on the spin is necessary. In principle, the spin truncation means that we would not be taking into account all inequalities obeyed by the coefficients y 0,0 and y m,n , resulting in stronger-than-correct bounds on the spectral function. 7 Nonetheless, working at a fixed derivative order N , we generally find that the spectral function bounds stabilize to within numerical precision once the maximal spin s max is taken to be sufficiently large (empirically, s max at order N is sufficient). For the application to theories with only scalar primaries in the next few subsections, of course, we do not need to worry about the spin truncation being sufficiently large. In this case, however, we must be especially careful in taking the truncation on the q-series of the conformal blocks to be sufficiently large, as the corrections to the approximate blocks would introduce nonzero spin primary contributions.
Secondly, since SDPB deals with the question of whether there exists a linear combination of polynomials p i (x) that is non-negative for all x ≥ 0, the above minimization problem must be recast in the form of inequalities on polynomial functions of ∆ on a semi-infinite line. For instance, suppose we impose a lower bound ∆ * s on the dimension of spin-s primaries as part of the a priori assumptions on the spectrum, then (2.7) is equivalently written as (2.11) By default, ∆ * s can be set to the unitarity bound. While the first inequality in (2.11) can be implemented in SDPB by a simple shift in the variable ∆, the second inequality which holds for ∆ in an interval is more subtle. It is handled 8 by converting the inequality to one on the semi-infinite line by a change of variable ∆ = (∆∆ * + ∆ * s )/(∆ + 1); now ∆ * s ≤ ∆ < ∆ * amounts to∆ ≥ 0.

Bounding the spectral function in a CFT with only scalar primaries
We now specialize to the case of CFT with only scalar primary operators. We do not specify the normalization of the primaries; as far as the spectral function is concerned, the external primaries are effectively normalized through the 4-point function (thus capturing relative OPE coefficients). This allows us to deal simultaneously with compact and noncompact CFTs. (By a non-compact CFT, we mean one with continuous spectrum and no SL(2, R) × SL(2, R)-invariant vacuum.) As alluded to in the introduction, there is only one known unitary CFT with c > 1 of this type, namely Liouville theory, and we will compare our bounds to the Liouville spectral function which can be obtained by numerically integrating the known OPE coefficients (given by DOZZ formula [22,23], as reviewed in Appendix B) with the Virasoro conformal blocks.
We can write the four-point function involving a pair of primaries φ 1 , φ 2 as g 12 (z,z) = ∞ 0 d∆ C 2 12;0,∆ F 12;0,∆ (z,z), (2.12) and the spectral function as This accommodates both continuous and discrete spectra (in the latter case the integral will receive contributions from delta-functions). To place bounds on f (∆ * ), we simply solve the minimization problem (2.7), (2.9) for s = 0 only. This is implemented with SDPB with a given set of c, ∆ 1 and ∆ 2 , while scanning over a range of ∆ * , at increasing derivative orders N .
First, we consider the case where all external operators in the four-point function have the same scaling dimension (above or below the Liouville threshold, ∆ 0 ≡ 2ξ). Our results for c = 8 are summarized in Figure 1. We observe that as the derivative order N increases, the upper and lower bounds approach one another, narrowing the allowed range of the spectral function. Both bounds appear to be converging upon the spectral function of Liouville theory (whose background charge Q is related to c by c = 1 + 6Q 2 ), which sits in the middle of the allowed window.
There exist solutions to the scalar-only crossing equations when the external operator dimension drops below the Liouville threshold, so long as ∆ φ ≥ 3 4 ∆ 0 . For ∆ φ < 3 4 ∆ 0 , solutions to the crossing equations with only scalar primaries in the OPE are excluded by our numerical analysis. When ∆ φ = 3 4 ∆ 0 , we find that the upper and lower bounds on the spectral function converge quickly to a step function, i.e., f + , already at small derivative order N . This case and an example where ∆ φ lies in between 3 4 ∆ 0 and the Liouville threshold are included in Figure 1.
In fact, for ∆ φ ∈ ( c−1 16 , c−1 12 ), our bounds on the spectral function are entirely consistent with the analytic continuation of the Liouville spectral function to external operator dimensions below the Liouville threshold. Indeed, such analytically continued Liouville correlators arise in the study of certain normalizable BPS correlators in super-Liouville theory [19] as a result of a relation due to Ribault and Teschner between SL(2) WZW model correlators and (d) Figure 1: Upper and lower bounds on the spectral function from linear functionals of increasing derivative order (from green to red), assuming only scalar primaries for c = 8 with ∆ φ /∆ 0 = 3 4 , 7 8 , 1, 24 7 . In all cases, the shaded regions are excluded and the black curve denotes the corresponding spectral function of (analytically continued) Liouville theory.
Liouville correlators [32]. A priori, the crossing invariant Liouville 4-point function involves external primaries of scaling dimension ∆ i = 2α i (Q − α i ), and an integration over internal primaries of scaling dimension ∆ = 2α(Q − α), where both α i and α lie on the half line Q 2 + iR ≥0 . We can analytically continue α i to the real axis, away from Q 2 , provided that no pole in the structure constant C(α 1 , α 2 , α) as a function of α crosses the integration contour Q 2 + iR. This is possible for Q 2 < α 1 + α 2 < Q, but fails for α 1 + α 2 ≤ Q 2 when a pole in α crosses the contour and the 4-point function picks up a residue contribution that violates unitarity. Indeed, α 1 = α 2 = Q 4 corresponds to ∆ φ = 3 4 ∆ 0 , and we find the step function behaviour demonstrated in Figure 1 whenever  12 7 ). The black curve denotes the (analytically continued) DOZZ spectral function. In (c), a small gap of ∆ gap = 0.01 has been imposed to explicitly exclude the vacuum channel which would correspond to a singular conformal block for the mixed correlator.
Next, we study the bounds on the spectral function for the 4-point function involving a pair of primaries φ 1 and φ 2 of different scaling dimensions, of the form (2.1). Note that for a non-compact CFT with only scalar primaries, such spectral functions capture the complete 9 This step function behavior is consistent with the fact that the 4-point conformal block with α 1 +α 2 = Q 2 and internal primary with α = Q 2 is crossing invariant by itself. This conformal block is the same as the holomorphic part of the 4-point function e 2α1φ(z) e 2α2φ(0) e 2α2φ(1) e 2α1φ(∞) in the linear dilaton CFT with background charge Q. Note that in the linear dilaton theory the closure of the OPE demands a non-unitary spectrum.
set of structure constants for three primaries of arbitrary weights. In Figure 2 we plot the upper and lower bounds on the mixed correlator spectral function for c = 8 with external primaries of various dimensions (∆ 1 , ∆ 2 ). Once again, the bounds narrow down the allowed window towards the spectral function of Liouville theory.
Apart from the case of α 1 + α 2 = Q 2 , our numerical upper and lower bounds have not quite converged convincingly to the (analytic continuation of) the Liouville spectral function, due to the computational complexity of computing bounds at high derivative order N . Our results nonetheless suggest such a convergence in the N → ∞ limit, supporting our conjecture that the DOZZ structure constants C(α 1 , α 2 , α 3 ) are the unique solution to the crossing equations for unitary CFTs with c > 1 and only scalar primaries.
Note that the convergence of the bounds on the φφφφ spectral function would determine the φφ OPE up to normalization; if this holds for all ∆ φ , it would then determine, assuming a non-degenerate spectrum, the conformal block decomposition of φ 1 φ 1 φ 2 φ 2 as well. This then determines the most general φ 1 φ 2 OPE, up to normalization. Compatibility with all crossing equations fixes the normalizations of OPE coefficients to be DOZZ up to an overall scale factor which cannot be fixed for a non-compact CFT. 10 Thus, in order to establish our conjecture for the uniqueness of the DOZZ solution for the scalar-only crossing equations in the non-degenerate case, it suffices to consider the OPE of pairs of identical primaries, and then the result for mixed correlators would follow.
One can notice that the bounds appear to change slowly with N in certain regions of the plots. We also observed in the numerical studies of spectral functions in modular bootstrap that the convergence of upper and lower bounds is relatively slow for continuous spectra as compared to discrete spectra (see section 4.2.2) in the cases where we know that the solution to the modular crossing equation is unique. It appears to be quite difficult numerically to push these bounds to higher derivative orders N , due to the need to substantially increase the truncation order d q on the q-expansion of the Virasoro conformal blocks. This is discussed in Appendix D.1. In the next subsection, we consider an alternative method of directly solving the linear system that determines the spectral function assuming that the optimal upper and lower bounds coincide. This method in fact does not rely on the assumption of reality of the OPE coefficients and appears to converge much faster to the DOZZ spectral function.
(3.1) That is to say, on a certain vector space of functions in ∆, the function Θ(∆ * −∆)F 12;0,∆ (1/2, 1/2) can be decomposed on the basis spanned by F 12;0,∆ (1/2, 1/2) and ∂ m z ∂ n z F 12;0,∆ (z,z)| z=z= 1 2 . Since the step functions are themselves complete, our conjecture of the DOZZ structure constants as the unique solution is related to the completeness of this basis on a suitably defined Hilbert 11 space of functions in ∆. 12 While we do not have a proof of this statement, we can analyze the linear problem directly in an attempt to solve for the coefficient y 0,0 (for a truncated system). The stability of the solution and its convergence to the Liouville spectral function will provide strong evidence for the conjecture.
Another way to arrive at (3.1) is the following. In a non-compact CFT with only scalar primaries the crossing symmetry equations, together with a normalization condition g 12 (1/2, 1/2) = 1, (3.1) can be written as We may equivalently express these equations as It is not obvious that the Hilbert space structure is the fundamentally correct one; for example, it might be that the correct notion is denseness in some Banach space. 12 Note that the linear independence of F 12;0,∆ (1/2, 1/2) from ∂ m z ∂ n z F 12;0,∆ (z,z)| z=z= 1 2 as functions of ∆ is guaranteed by the existence of DOZZ structure constants as a solution to the crossing equation.
where the vectors v, p n,m represent the functions for some suitable choices of f v (∆) and f p (∆) (see Appendix D.2 for details), while the inner product is defined by We now hope for completeness of the set p n,m (from hereon we only consider n = m = 0 or n + m odd), assuming that all functions in question have finite norm. We truncate by n + m ≤ N and consider the approximation of v by its orthogonal projection v N = P N v onto P N = span{p n,m } n+m≤N . Note that despite the notation v N , because of the equations (3.3), v N is independent of a particular solution v. It can be computed by evaluating the Gram matrix of p vectors and taking its inverse, where (G N ) n,m n ,m = p n,m , p n ,m , n + m ≤ N, n + m ≤ N.
The spectral function can be computed as the inner product We have an estimate, Note that E N = |(1 − P N )θ ∆ * | 2 /|θ ∆ * | 2 is also independent of a particular solution v and is computable from (3.3).
If (3.1) holds in the norm induced from ·, · , then E N → 0. Conversely, if we show for all ∆ * that lim N →∞ E N = 0, it will imply that any normalizable solution to (3.3) and thus 0.5 to (3.2) is equal to the limit lim N →∞ v N , which is unique if exists. Our strategy would be therefore to evaluate v N and E N numerically and estimate their limits.
We first numerically evaluate f N (∆ * ) = v N , θ ∆ * and find that it converges to the Liouville spectral density in the limit N → ∞. For example, in Figure 3(a) the approximation f N is plotted at successive odd values of N up to N = 25 for c = 8 and ∆ φ = 7/12. We can see that the curves exhibit the expected convergence. Another example where the external operator dimension is far above the Liouville threshold is shown in Figure 3(b), where we studied c = 2, ∆ φ = 55/12, up to N = 33 and d q = 200. While f N (∆ * ) oscillates wildly at smaller N (the case N = 27 is shown for comparison), the oscillation settles down substantially as N is increased.
In Figure 4 we compare f N (∆ * ) with the DOZZ spectral function for c = 8 and c = 30, with ∆ φ at or above the Liouville threshold, as well as an example of a mixed correlator spectral function 13 with two different values of external operator dimensions. In all cases we find good agreement.
To further support the conjecture, we numerically compute the error estimate E N as a function of N . For example, in Figure 5 we show E N as a function of 1/N for ∆ φ = c−1 12 , ∆ * = c 10 , and c = 8. In the figure we also show a linear fit using N ≥ 11. Empirically, we find that the result is consistent with E N ∼ N −1 . We study E N in more detail in appendix D.2.1.
The discussion above depends on the assumption that v has finite norm. This assumption 13 The mixed correlator spectral function was obtained under a technical assumption of a lower dimension bound of ∆0 5 (note that this is below the Liouville threshold ∆ 0 ), see appendix D.2 for details. 0.5 itself depends on the choice of measure. We describe our choice of measure and details of our implementation in appendix D.2. Here we simply note that with our choice, v has finite norm if C 4 12;0,∆ is locally integrable on [0, ∞), and the OPE expansion is convergent in the region |z| < 1. Discrete spectra have infinite norm since C 4 12;0,∆ involves squares of delta-functions, but such spectra are excluded by modular invariance.

Constraints from modular invariance
Strong constraints on the primary spectrum, especially in the scalar-only case, follow from modular invariance alone. In fact, there is a simple argument that shows any 2D CFT with c > 1 and primary operators with bounded spin must have a spectrum identical to that of Liouville theory: that is, the spectrum is non-compact, has scalar primaries only, and has a spectral density that is uniformly distributed in Liouville momentum P = 2(∆ − c−1 12 ). Suppose the primaries have spins no greater than s max . We can write the reduced torus partition function in the following way: is the degeneracy of primary operators in the spectrum with conformal weights (h,h) and f s ( . For now we assume that the CFT is compact, and the vacuum character is degenerate and so s max ≥ 1. The non-compact CFTs may be viewed as limiting cases, where the spectral density diverges and we divide the partition function by an infinite normalization factor which removes the vacuum contribution. Now consider the following change of variables chosen so that the modular S transformation exchanges x and y. We can then write the modular crossing equation in terms of these variables as Of course, the functions above have branch cuts at x = y −1 , but since the sum over spins is finite by assumption, the analytic continuation around the branch is straightforward. Furthermore, f s (y) is an analytic function for Re(y) > 0. To proceed, we fix x = re −iα with r > 0 and 0 < α < π 2 and y = with → 0 + , so that the modular crossing equation becomes where in the first line we dropped the phase factors e 2πis √ r e iα − 2 (which are close to 1 due to the boundedness of s) in front of f s ( ); this is a valid approximation since f s ( ) is positive for all s. In the second line we again invoked modular invariance (this particular equality is realized as the modular crossing equation with τ 1 = 0, τ 2 = ). In the case that the CFT is compact, the right-hand side is dominated by the contribution of the vacuum, in particular By comparing to the → 0 limit of the left-hand side, which is dominated by the term with maximal spin s e 2πs √ r (sin α 2 +i cos α we arrive at a contradiction and deduce that unitary 2D CFTs with primary operators of bounded spin must have non-compact spectra: namely, there is no SL(2, R) × SL(2, R)invariant vacuum and the dimension of the lowest-lying primary operator obeys ∆ min > 0.
In fact, this same logic allows us to conclude that the dimension of the lowest-lying operator must obey ∆ min ≥ c−1 12 . In the → 0 limit, we have where ρ(∆) is the density of primary operators in the spectrum with dimension ∆ (of any spin). The two sides of the equation are clearly incompatible if the minimum scaling dimension for which ρ(∆) is nonzero is smaller than 2ξ. Furthermore, by non-negativity of the spectral density, the right-hand side can grow no faster than − 1 2 as → 0 + . On the other hand, the absolute value of (3.18) grows like e 2πsmax √ r sin α 2 in this limit. Modular invariance thus demands that s max = 0: that is, a unitary 2D CFT with primary operators of bounded spin must in fact have only scalar primary operators in addition to having a non-compact spectrum. Moreover in this case the modular crossing equation becomes 20) which demands that f 0 (x) is a constant. Thus the required spectral density is nothing other than that of Liouville theory, namely ρ(∆) = ρ Liouville (∆) ∝ (∆ − 2ξ) − 1 2 Θ(∆ − 2ξ), 14 completing the argument. In particular, the dimension of the lowest-lying operator must be exactly ∆ min = 2ξ.
By this result, our conjecture that the DOZZ structure constants are the unique solution to the crossing equations for a unitary 2D CFT with central charge c > 1 and only scalar primaries, as supported by substantial numerical evidence in sections 2.2 and 3.1 leads us to conjecture that Liouville theory is the unique unitary c > 1 CFT with Virasoro primaries of bounded spin. 15

Degenerate spectrum and TQFT
In our analysis of the crossing equation so far, we have implicitly assumed that the scalar primaries are labeled by a continuous parameter, namely the scaling dimension ∆ φ , without further degeneracy. If this assumption is relaxed, one can construct more examples of (noncompact) c > 1 CFTs with only scalar primaries, by taking the tensor product of Liouville CFT with a topological quantum field theory (TQFT); the latter has a finite dimensional Hilbert space on the circle and its structure constants are governed by those of a commutative Frobenius algebra [34]. 16 We conjecture that this is the only possibility.
Let us assume that the scalar primaries are labeled by their scaling dimension ∆ and an extra index i, and denote the structure constants by where we have explicitly factored out the DOZZ structure constants. Our numerical results in the previous sections on the spectral function of mixed correlators of the form φ 1 φ 2 φ 2 φ 1 14 To normalize the reduced partition function of Liouville theory to 1, the constant of proportionality is √ 2.
15 Note that we have not made use of the torus 1-point function. A priori, the modular invariance of the torus 1-point function puts nontrivial constraints on the structure constants with a pair of primaries identified. For the purpose of establishing our conjecture regarding the uniqueness of Liouville, once the OPE coefficients are pinned down to those of DOZZ by the crossing equation, the torus 1-point functions are already modular invariant [33]. 16 To be precise, we do not need to require the TQFT to have a vacuum state (or the algebra to be unital).
indicate that for a CFT with degenerate scalar-only primary spectrum, is independent of ∆. In fact, we can strengthen this result slightly. Let us consider a mixed correlator φ i φ j φ k φ where φ i , φ have scaling dimension ∆ 1 , φ j , φ k have scaling dimension ∆ 2 , and the crossing equation By taking the part of (3.23) that is odd under z → 1 − z,z → 1 −z, our earlier claim of the uniqueness of scalar-only solution to the crossing equation implies that is independent of ∆. On the other hand, for the even part of (3.23) under z → 1−z,z → 1−z, the numerical analysis described in appendix D.2.1 is consistent with the conjecture that {∂ n z ∂ m z F 12;0,∆ | z=z= 1 2 , n, m ∈ Z ≥0 , n + m even} form a complete basis on the space functions of ∆ on the positive real axis defined by the same norm as in section 3.1, which implies is independent of ∆.
It is likely that by analyzing a system of crossing equations for multiple scalar correlators involving φ i , φ j , φ k , φ of generally different scaling dimensions, one could establish that the spectral function for φ i φ j φ k φ (with only scalar Virasoro primaries in the OPEs) is proportional to that of Liouville CFT, which would be equivalent to the statement that is independent of ∆, extending (3.24). We leave the numerical bootstrap of the spectral function with four generic external weights to future work. We now argue that if (3.25) holds, then our conjecture follows.
To each pair-of-pants decomposition of a genus g Riemann surface, represented by a trivalent graph, we may associate a sum of product of A ijk 's, with indices contracted and scaling dimensions identified along each edge of the graph, which we denote by Z g . (3.25) implies the crossing relation between graphs with fixed weights on the edges, and by applying crossing one can always turn the trivalent graph into one that does not contain tadpole subgraphs. 17 (3.25) further implies that Z g is independent of the scaling dimension on every edge that connects a pair of distinct vertices, and thus the genus g partition function of the CFT is equal to Z g times the Liouville partition function. It then follows from modular invariance that Z g is independent of the pair-of-pants decomposition, and depends on the genus g only.
To proceed, pick a finite set of scaling dimensions ∆ a , a = 1, . . . , M and let N ∆ be the number of degenerate primaries of dimension ∆, which we will assume to be finite. Set N = max N ∆a and extend the ranges of the discrete labels to run up to N for all ∆ a by setting the previously undefined structure constants to zero. Then the totality of A ijk (∆ a , ∆ b , ∆ c ) gives an element A in C = S 3 (⊕ a R N ). The space C is equipped with an action of a O(N ), corresponding to changes of basis for the discrete labels. Z g regarded as polynomials generate the algebra of a O(N ) invariants on C. It follows that A is equivalent to any other A with the same values of Z g by a a O(N ) reparametrization. In particular, since A is such that values of Z g on it are independent of the internal labels, A is equivalent to A 0 in which all A ijk (∆ a , ∆ b , ∆ c ) are replaced by a ijk = A ijk (∆ 1 , ∆ 1 , ∆ 1 ). It then follows that we can choose N ∆a = N .
By taking various finite sets of scaling dimensions sharing the dimension ∆ 1 , we find that N ∆ = N is independent of ∆ and thus the density of states is given by N copies of Liouville density. Furthermore, for any such finite set we have up to a reparametrization of finite labels. Note that such reparametrizations depend on the choice of our finite set of ∆ a 's, being defined only up to automorphisms of a ijk . As we show below, these automorphisms are scarce. Compatibility between different ∆ a then completely fixes them after we fix the reparametrization for ∆ 1 . Thus we can pass to the full continuous set of scaling dimensions, and conclude that the CFT in question is a tensor product of Liouville with a TQFT defined by the structure constants a ijk (or the partition functions Z g ).
In fact, we can always find a basis in which the structure constants a ijk are diagonalized. To see this, note that the crossing equation for a ijk implies that the matrices M i with entries (M i ) jk = a ijk are mutually commuting N × N symmetric matrices, and thus can be simultaneously diagonalized by some O(N ) matrix R, namely Λ ijk = mn R jm R kn a imn are diagonal in jk. Then Λ ijk = m R im Λ mjk is still diagonal in jk and completely symmetric, and thus Λ ijk = δ ijk λ k . Multiplying by a diagonal matrix with ±1 entries if necessary, we can set λ m > 0. (If some λ m = 0 they do not contribute to the correlators and we can obviously add or remove such λ's at will.) It is straightforward to check that the automorphisms of Λ ijk are just the permutations preserving the λ's. The partition functions are Z g = n λ g−1 n . This diagonalization implies that the algebra defined by a ijk is given by ⊕ n G λn . Here G λ = Re with (e, e) = 1 and e 2 = λe. Forming the tensor product Liouville ⊗ G λ corresponds to rescaling all OPE coefficients by λ. The overall scale of OPE coefficients cannot be fixed in the absence of the vacuum, and thus we can regard all these theories as isomorphic to Liouville. Therefore, the TQFT structure amounts to superselection sectors.

The modular spectral function 4.1 The minimization problem
We now consider the decomposition of the reduced torus partition function of a compact, unitary CFT (assumed to be parity-invariant) 18 with no conserved currents into non-degenerate Virasoro characterŝ where I s is the discrete spectrum of dimensions of primary operators, d ∆,s = d( ∆+s 2 , ∆−s 2 ) = d( ∆−s 2 , ∆+s 2 ), and the reduced characters are given bŷ

(4.2)
Analogously to the four-point spectral function introduced in section 2.1, we define a "modular spectral function" by truncating the Virasoro character decomposition of the reduced par-tition function up to a cutoff dimension ∆ * , evaluated at the self-dual modulus τ = −τ = i, As with the four-point spectral function, it is straightforward to place bounds on f mod (∆ * ) due to modular invariance using semidefinite programming. DefiningẐ ∆,s (τ,τ ) =χ ∆+s 2 (τ )χ ∆−s 2 (τ )+ χ ∆−s 2 (τ )χ ∆+s 2 (τ ) andẐ 0,0 (τ,τ ) =χ 0 (τ )χ 0 (τ ), the modular crossing equation demands that where we have redefined τ = ie z ,τ = −ie −z . We then seek to minimize y 0,0 subject to the inequalities (4.5) for arbitrary coefficients y m,n . In the first line we have singled out the inequality involving the vacuum primary. In the second line, we made the extra assumption of a gap ∆ * s in the spin-s sector of the spectrum, as will be useful in later applications. As before, the minimal such y 0,0 gives an upper bound on the modular spectral function, since  Similarly, the minimal w 0,0 subject to the constraints provides a nontrivial lower bound on the modular spectral function Working up to a finite derivative order m + n ≤ N , we denote the corresponding upper and lower bounds obtained in this way f + mod,N (∆ * ) and f − mod,N (∆ * ) respectively.

Some consistency checks 4.2.1 Extremal spectra with maximal gap
In [20], an upper bound ∆ mod (c) on the gap in the scaling dimension of primary operators due to modular invariance of the torus partition function was computed numerically as a function of the central charge. Given a dimension gap ∆ gap (≤ ∆ mod (c)), an upper bound on the degeneracy of primaries at dimension ∆ gap can be obtained provided ∆ gap > c−1 12 . When this upper bound on the degeneracy at the gap is saturated, the entire modular invariant spectrum is determined by the locations of the zeros of the optimal linear functional (optimized with respect to the degeneracy bound) acting on the Virasoro characters. Such (candidate) CFT spectra were dubbed 'extremal.' Furthermore, it is expected that for each given c > 1, there is a unique modular invariant spectrum (imposing positivity but not the integral condition on the degeneracy of primaries) whose dimension gap saturates the upper bound ∆ mod (c) [35].
In [20], a number of examples of CFTs with spectra that saturated the bound on the dimension gap were identified at small values of the central charge. Here, we study the bounds on the modular spectral function at these values of the central charge assuming the maximal dimension gap. We will find that the resulting bounds indeed pin down the extremal modular spectral functions. To compute the bounds on the modular spectral functions in these cases, we impose (4.5,4.7) with ∆ * s = max(s, ∆ mod (c)). For c = 2, the dimension gap bound of ∆ mod (2) = 2 3 is realized by the spectrum of the SU (3) WZW model at level one. This theory admits a description in terms of free bosons with T 2 target space at the Z 3 -invariant point in its complex structure and Kähler moduli spaces, with partition function where k 2 L,R = G mn α (n m + B mk w k ± G mk w k )(n n + B nl w l ± G nl w l ) (4.10) . The bounds on the modular spectral function collapse precisely to this extremal modular spectral function when the maximal gap is imposed, as shown in Figure 6.
For c = 4, the dimension gap bound of ∆ mod (4) = 1 is realized by the spectrum of the SO(8) WZW model at level 1, which also admits a description in terms of 8 free fermions with diagonal GSO projection. This theory occupied the kink on the curve ∆ mod (c). The partition function of this theory is given by Once again, in Figure 6 we see that the bounds on the modular spectral function collapse to that of the extremal spectrum.
For c = 8, there is a nontrivial bound on the dimension gap in the spectrum of scalar primaries, ∆ s=0 mod (8) = 2. This bound on the scalar gap is saturated by the spectrum of the E 8 WZW model at level one. This theory, which occupied the first kink on the bounding curve ∆ s=0 mod (c), admits an equivalent description in terms of 8 compact bosons at the holomorphically factorized point in the moduli space; the holomorphic factor is described by the Narain compactification on Γ 8 , the root lattice of E 8 . The partition function is where j(τ ) is the elliptic j-invariant. Figure 6 shows that the bounds on the modular spectral function (derived using ∆ * s = ∆ s=0 mod δ s,0 + s) collapse to that of the extremal spectrum.

Only scalar primaries
As an additional example to illustrate the convergence of the bounds on the modular spectral function, we revisit the case of a CFT with only scalar primary operators. In section 3.2, we showed that a unitary c > 1 CFT with primaries of bounded spins must have a non-compact spectrum with only scalar primaries and a density of states equal to that of Liouville theory. The modular spectral function of Liouville theory is given by Assuming a scalar-only spectrum, and a dimension gap 2ξ (that is, we impose (4.5,4.7) for s = 0 only with ∆ * 0 = 2ξ), the numerical results of upper and lower bounds on the modular spectral function are shown in Figure 7. Note that while the bounds do appear convergent toward the Liouville modular spectral function (as they must), the rate of convergence is The bounds for c = 2 with an assumed dimension gap one-half of (left) and equal to (right) the maximal gap allowed by modular invariance. Bottom: The bounds on the modular spectral function for c = 4 with the maximal dimension gap (left) and for c = 8 with the maximal gap in the spectrum of scalar primaries (right). In all cases, the dotted lines denote the modular spectral function for the corresponding extremal spectrum.
rather slow compared to the previous examples of discrete extremal spectra at small c. On the other hand, such a slow convergence with N is qualitatively similar to our bounds on the 4-point spectral function in the scalar-only case, as analyzed in section 2.2, where we also expect a continuous spectrum, and also to the non-compact example in the next subsection.
We thus expect that the slow convergence is associated with continuity of the spectrum. This is natural from the point of view of extremal functionals -the extremal spectrum should converge to the continuous one as N → ∞, but at any finite N the extremal spectrum has finitely many operators, which therefore should condense. On the other hand, in the discrete case we typically need only a small number of operators to accurately determine the partition function or the correlation function in the neighborhood of the crossing/modular symmetric point.

No scalar primaries
In [20], following the observation that the bound on the gap of the dimension of scalar primaries diverged as c → 25 − , it was shown that for c ≥ 25 there exist (non-compact) modular-invariant spectra with no scalar primary operators. This is due to the fact that the modular invariant function  To compute the bounds on the modular spectral function in this case, we impose (4.5,4.7) for s > 0 with ∆ * s = max(s, 2ξ−1). As shown in Figure 8, the bounds on the modular spectral function do indeed appear to be converging to (4.15) as the derivative order of the linear functional is increased, suggesting that the no-scalar spectrum is unique for c = 25. Note that for c = 25 the dimension gap c−13 12 coincides with the unitarity bound (since we assume no scalars). We expect that for c > 25 the uniqueness holds only under the assumption of the dimension gap c−13 12 .

CFTs at large c with large gap
In [20], the upper bound on the dimension gap ∆ mod (c) due to modular invariance of the torus partition function was computed numerically for central charge up to c ∼ O(10 2 ). As c is increased, the convergence of the upper bound with increasing derivative order N slows, and accurate determinations of the optimal bound on the gap require a careful extrapolation to the limit N → ∞. Nonetheless, a conjecture on the monotonicity of the slope of the optimal bounding curve d∆ mod (c) dc leads one to conclude that the asymptotic slope is less than 1 9 . Potentially, the asymptotic slope could be as small as 1 12 , a possibility that is natural from the holographic perspective (see the discussion in the next section) but with no direct evidence from the analysis of the modular crossing equation.
Thus at large c it has been difficult to determine ∆ mod (c) accurately, and furthermore the exponential growth of operator degeneracies (combined with the need to go to very large N to get a good approximation of the optimal linear functional at large c) makes it practically impossible to resolve the discreteness of the spectrum by bounding the modular spectral function even when the bound ∆ mod (c) is saturated. Nonetheless, we can study the bounds on the modular spectral function assuming a gap ∆ gap close to ∆ mod (c), at values of c where the value of ∆ mod (c) can be reliably computed by numerical extrapolation of ∆ (N ) mod (c) to N = ∞. Figure 9 shows plots of the bounds on the modular spectral function for c = 50, 100, 300 with assumed dimension gap ∆ gap close to the bound ∆ mod (c). (c) Figure 9: Upper and lower bounds on the modular spectral function in the case that the dimension gap is close to the maximal value allowed by modular invariance for c = 50, 100, 300. The dotted black curve shows the modular spectral function of perturbative pure gravity due to the thermal AdS 3 and Euclidean BTZ saddles in the gravitational path integral.
The plots reveal several interesting features that we believe are universal at large c assuming a sufficiently large gap ∆ gap (> c− 1 12 ). Firstly, for ∆ * c 6 , the upper and lower bounds on the modular spectral function converge to f mod (∆ * ) = 1 2 : that is, modular invariance demands that the vacuum character accounts for exactly half of the partition function at the self-dual temperature. Furthermore, the bounds appear to be convergent upon a smooth function that interpolates between 1 2 and 1 in a window of size ∼ √ c about ∆ * = c 6 .
It is generally expected that 2D CFTs at large c with large gap should be holographically dual to a semiclassical theory of pure gravity in AdS 3 . To the best of our knowledge, this statement has not been precisely formulated: how large does the gap need to be? If we merely demand that the dimension gap 19 grow linearly in c, corresponding to a Planckian mass gap in the bulk theory, but with a coefficient less than 1 12 , then the entropy need not follow Bekenstein-Hawking in the entire range ∆ ≥ c 6 . This is the range of masses for which BTZ black holes dominate the canonical ensemble at its Hawking temperature. On the other hand, one might expect that CFTs with gap close to ∆ mod (c) (if they exist) are holographic duals to suitable non-perturbative completions of pure gravity in AdS 3 , in the sense that observables such as the spectral density are correctly captured by the perturbative expansion around known saddle points of the gravitational path integral in the bulk up to exp(−c) corrections.
This suggests that we compare the bounds on the modular spectral function to that of pure gravity, which, up to a priori unknown non-perturbative corrections, is computed by the contributions from thermal AdS 3 and the Euclidean BTZ black hole saddle points, which are known to be perturbatively 1-loop exact. We derive this modular spectral function in Appendix C, see in particular (C.3,C.5). The bounds shown in Figure 9 indeed appear to be converging upon the pure gravity result (1.8) for dimensions above the assumed gap.
Note that for ∆ * c 6 , that the bounds on the modular spectral function with large gap converge to 1 2 can be explained by the fact that in the semiclassical limit, the gravitational path integral evaluated at the self-dual temperature is dominated by the contributions of two saddles mentioned above, which are exchanged by the modular S transformation, and thus the vacuum contribution accounts for 1 2 . Interestingly, at large c, the upper and lower bounds on the vacuum contribution already converge to 1 2 when the dimension gap is slightly above c−1 12 , not necessarily close to ∆ mod (c). This is illustrated in Figure 10, where we plot the bounds on the contribution of the vacuum to the modular spectral function as a function of the dimension gap for c = 8, 50, 100. Note that the vacuum contribution to the spectral function determines the partition function Z(τ,τ ) itself at the self-dual temperature (τ = −τ = i). From the bulk perspective, that the vacuum accounts for 1 2 of the modular spectral function amounts to the statement that the thermal AdS 3 and BTZ saddle points are the two dominant saddle points in the gravitational path integral, while all other saddle points are exponentially suppressed. 20 In [36], it was shown that in theories where the light spectrum is appropriately sparse (a condition of which the maximal gap is an extreme case), the microcanonical entropy is given universally by the Cardy formula to leading order for all operators with dimension ∆ ≥ c 6 . Note that the sparseness criterion is satisfied by our assumption on the gap, but our statement regarding the modular spectral function and thereby the spectral density extends to the regime of ∆ slightly below c 6 (see further discussion in the next section). We conjecture that in the large c limit, assuming a gap sufficiently close to ∆ mod (c), 21 the bounds f ± mod (c) converge onto the modular spectral function of pure gravity described above, up to order exp(−c) corrections. Indeed, we note that for c ∼ O(10 2 ), the horizontal average of the bounds f mod,N (∆ * ) is already well approximated by the pure gravity modular spectral function at moderate values of N , as shown in Figure 11.

On the universality of the BTZ spectral density
The BTZ black hole in AdS 3 has a striking feature that is unlike black holes in other spacetime dimensions (in asymptotically either AdS or flat spacetime): a Planckian mass BTZ black hole has a macroscopic horizon radius (rather than, say, Planckian radius), provided that the mass is an order 1 fraction above the BTZ threshold in Planck units. The microstates of such a black hole would be dual to an operator in the CFT of dimension ∆ > (1+ ) c 12 , where is an order 1 fraction that does not scale with c. When there is a sufficiently large mass gap in the spectrum, standard effective field theory reasoning in the bulk would suggest that the the entropy of the BTZ black hole, which captures the degeneracy of microstates, should be computed from the Bekenstein-Hawking formula based on the Einstein-Hilbert action, as any local higher-derivative corrections to the Einstein-Hilbert action of pure gravity in three dimensions can be absorbed by field redefinition. This would predict a degeneracy or spectral density to leading order in the large c limit, for ∆/c > 1 12 . For ∆ > c 6 , i.e., above twice the BTZ threshold, this universal behavior of the spectral density was demonstrated in [36] to be a consequence of the sparseness of the spectrum and modular invariance. From the gravity perspective, this is also the regime in which the Euclidean BTZ black hole solution is the dominant saddle point of the Euclidean pure gravity path integral, i.e., the BTZ black hole dominates the canonical ensemble at its Hawking temperature, and therefore the spectral density must be (5.1) in order to produce the correct free energy above the self-dual temperature.
The regime c 12 < ∆ < c 6 is much more interesting. Here the BTZ black hole does not dominate the canonical ensemble. Its contribution to the gravitational free energy is nonperturbatively suppressed compared to the thermal AdS 3 contribution. A priori, since we do not know the most general non-perturbative contributions to the pure gravity path integral, we cannot draw any reliable conclusion on the spectral density in this regime. This also puts doubt on the validity of the Bekenstein-Hawking formula, despite the macroscopic size of the horizon. If the Bekenstein-Hawking formula is violated in this regime, it then indicates some sort of breakdown of the effective field theory reasoning based on locality.
What conclusion can we draw from our numerical bounds on the modular spectral function for c ∼ O(10 2 )? We saw that in a window of size √ c around c 6 , the modular spectral function is constrained to well approximate the AdS 3 + BTZ answer, in agreement with the expectation from the known perturbative contributions to the Euclidean gravity path integral. Unfortunately, our numerical results do not have sufficient resolution to allow for distilling a contribution of order exp(−c), thus preventing us from concluding whether the Bekenstein-Hawking entropy of BTZ correctly accounts for the spectral density in the regime ∆ = yc, for 1 12 < y < 1 6 . In fact, if the latter is true, then the asymptotic slope of the modular bound on the dimension gap, lim c→∞ d∆ mod (c)/dc, must be equal to 1 12 , but this has not been shown. Thus, the fate of the small-yet-large BTZ black holes below twice the BTZ threshold remains a mystery.
where the nome q(z) is defined as If we define then H(λ 2 i , h|q(z)) satisfies Zamolodchikov's recurrence relation where h m,n are the weights of degenerate Virasoro representations, and R m,n ({λ i }) are The product of (r, s) is over

B Liouville CFT and DOZZ structure constants
The Liouville CFT is parameterized by the central charge c = 1 + 6Q 2 , where Q = b + b −1 , and a cosmological constant −µ < 0. It is governed by the action To study Liouville theory on the sphere, one typically works with a flat reference metric g mn supplemented with the boundary condition The field φ(z,z) is not a primary operator under holomorphic coordinate transformations z → w(z). In this case one must take care to regulate the action and introduce boundary terms to ensure that the action is finite and invariant under conformal transformations.
The Hilbert space consists of a continuous spectrum of scalar primary operators V α with α ∈ Q 2 + iR ≥0 and conformal dimension ∆ = 2α(Q − α). Operators with α outside this range, such as the identity operator, do not correspond to normalizable states and thus do not belong to the Hilbert space. Making use of a somewhat nonstandard convention (the reason for which will become clear soon), we normalize the primaries so that in the asymptotic regime where the Liouville potential vanishes (the φ → −∞ limit) they take the form 22 V α ∼ S(α) − 1 2 e 2αφ + S(α) where S(α) is the reflection amplitude The torus partition function is the same as that of a single non-compact free scalar. The sphere two-point function of primary operators is Note that with our choice of conventions the two-point function is canonically normalized. 22 In the literature on the Liouville CFT, usually considered are operators with the asymptotic form V α ∼ e 2αφ and which do not have standard two-point functions.
The sphere three-point function is given by the DOZZ structure constants [22,23] .
(B.6) The special functions are given by the following Note in particular that the upsilon function satisfies is a real function of P . To extend Υ b (x) beyond the range of its definition, one notes the following identities which can be proven by considering an integral representation of log Γ(x). The function Υ b (x) has simple zeroes at x = 0, x = Q as well as x = mb + n b when m and n are both nonpositive integers, and when m and n are both positive integers. It is instructive to rewrite the Liouville three-point function coefficient as a manifestly real function of the Liouville momenta P i = −i(α i − Q 2 ), since P takes non-negative values for operators in the physical Hilbert space: where we have used that the reflection amplitude can also be written as where C 0 = 1, C 1 = −2, C 2 = 1. Note that the other known saddle points of the gravitational path integral, related by SL(2, Z) transformations, are always non-perturbatively suppressed for purely imaginary τ .
Thus the perturbative pure gravity modular spectral function, obtained from the thermal AdS 3 and Euclidean BTZ saddle points in the gravitational path integral, can be written as We are interested in the behaviour of this function for ∆ * in a window of size ∼ √ c about c 6 in the semiclassical limit. From the asymptotic form of the Bessel function, it is easy to see that for y ∼ O(1), we have , we end up with the modular spectral function where we have kept only the leading terms in the semiclassical approximation. This is the same as the spectral function one would obtain from applying the "naive" Cardy formula (5.1).

D Details of the numerical computations D.1 Details of the solution of the semidefinite problem
Here we provide some details of the numerical computations of the bounds on the spectral functions, implemented using the SDPB package [31]. In practice, there are several truncations that must be made. First, we must restrict to a finite basis of linear functionals acting on the crossing equation, of total derivative order N . We must also approximate the Virasoro conformal blocks; we can only compute the blocks to a finite order d q in the elliptic nome q(z) from Zamolodchikov's recurrence relation (reviewed in Appendix A). Finally, recall that the upper and lower bounds on the spectral function are derived as the minimal coefficients such that a certain set of positivity conditions (for instance (2.7) or (4.5)) can be satisfied by a linear combination of derivatives of the conformal blocks or characters evaluated at the crossing symmetric point. In practice, we can only impose the positivity conditions on the blocks or characters of a finite set of spins in the spectrum; we denote the maximal spin considered by s max .
The truncation on spin means that for a fixed derivative order N , we will not have taken into account all inequalities that the coefficients in (2.7) or (4.5) must satisfy to constitute a bound on the spectral function, leading to bounds that are in principle too strong. Meanwhile, the truncation to finite d q introduces a controlled error into the computation of the (derivatives of the) conformal blocks evaluated at the crossing-symmetric point. Thus to derive bounds at a fixed N , we must ensure that both s max and d q are sufficiently large so that a bound exists and is stable against further increasing these parameters to within our numerical precision. It is worth emphasizing that while the truncations to finite s max and d q are controlled approximations, when these parameters are sufficiently large the bounds derived using a fixed derivative order N are rigorous. Of course, the optimal bounds are obtained in the N → ∞ limit.
Let us begin by discussing the bounds on the sphere four-point spectral function in the case that there are only scalar primaries in the spectrum. Of course in this case we need not worry about the spin truncation. We should note that in practice, the inequalities we feed into semidefinite programming are not quite of the form (2.7,2.9), for the simple reason that the Virasoro blocks are not polynomials in the dimension of the internal primary. To illustrate the procedure, we write where P 12 (∆; q,q) is a binomial in q,q with ∆-dependent coefficients. Derivatives of the blocks can then be cast in terms of ∂ n z ∂ m z (qq) ∆ 2 −ξ P 12 (∆; q,q) where P n,m (∆) is a polynomial in ∆ and where ∆ i are the locations of the poles kept at the given order of approximation in the computation of the Virasoro block. Importantly, the prefactor (16e −π ) ∆ Q −1 (∆) is nonnegative for unitary values of the internal dimension. The positivity conditions (2.7), (2.9) then amount to the following (y 0,0 − Θ(∆ * − ∆))P 0,0 (∆) + In the case that there are only scalar primaries in the spectrum, we must take particular care to ensure that d q is sufficiently large, for the reason that the four-point function when decomposed into Virasoro blocks truncated at a finite order in q would appear to have contributions from primaries of nonzero spin. Upper and lower bounds on the spectral function in this case can only be found numerically at a fixed derivative order N when d q is sufficiently large. Empirically, for central charges and derivative orders in the ranges considered in section 2.2, we find that d q = 4N is sufficient to compute stable bounds on the spectral function. To illustrate the convergence of the bounds as the derivative order is increased, Figure 12 shows the upper bound on the spectral function f N + (∆ * = 7 12 ) as a function of N −1 for c = 8 with the external operator dimensions at the Liouville threshold. It is clear that we have not been able to access N sufficiently large so that extrapolation to the N → ∞ limit can be reliably performed. Since the q-truncation order is the bottleneck for the speed of the numerical computations, this limits the range of derivative orders we are able to consider in computing bounds on the spectral function. For this reason, it is convenient to consider bounds obtained by further truncating the basis of linear functionals to ∂ n z ∂ m z | z=z= 1 2 with m + n ≤ N and either m ≤ 1 or n ≤ 1. This basis leads to weaker bounds at a fixed N , but renders bounds at larger N accessible. For instance, as shown in section D.3, we are able to compute bounds on the spectral function up to N = 25 with d q = N + 9 using linear functionals in this reduced basis. However, we caution that it is not clear that the N → ∞ limit of the bounds obtained using this reduced basis of linear functions converges to that of the full-basis bounds.
We now turn to the bounds on the modular spectral function. The implementation of the positivity conditions (4.5) and (4.7) proceeds similarly as in deriving bounds on the four-point spectral function; here, for each spin one simply factors out q fv(∆) to be locally integrable and decaying sufficiently quickly. If C 4 12;0,∆ gives a convergent OPE expansion for |z| < 1, it should decay at infinity at least as 16 −2∆ . The decay condition is therefore automatically satisfied if fp(∆) fv(∆) grows slower than 16 2∆ . We need also to ensure that p n,m have finite norm. The norm is where q = e −π , and we get that the ratio fv(∆) fp(∆) should grow slower than (16q) −2∆ . We then choose this ratio to be where λ ∈ (0, 2π). We will set λ = π.
Let us describe the details of the calculation. We want to compute the Gram matrix p n,m , p n ,m . For this, we need to be able to compute integrals with the Virasoro conformal blocks. Recalling (D.1), the inner product p n,m , p n ,m is given by It is a standard fact that such integrals can be evaluated in terms of incomplete gamma functions. However, we want to do this efficiently, since P and Q are high-degree polynomials. 23 The computation of products P n,m (∆)P n ,m (∆) can be optimized by means of fast Fourier transform. After the product is computed, it suffices to compute the integrals Reduction to incomplete gamma function is immediate if we shift ∆ by ∆ i . However, in this case we need to expand (∆ + ∆ i ) n which produces n terms, and n can be large. Instead we write Here k is bounded by 4, so we get a compact sum. The parameter n does enter into the incomplete gamma function, but it satisfies a recursive relation which allows one to compute it as n increases, effectively making the complexity of computation of this integral O(1) for every value of n.
Having found the Gram matrix, it is immediate to find the coefficients of the expansion of v N in the basis p n,m ; they are given by the first row of the inverse of the Gram matrix. Computation of P N θ ∆ * proceeds similarly, except that now we need to know all the inner products p n,m , θ ∆ * . These can be computed just as above, shifting everything by ∆ * . In practice, we find that the basis of p n,m is ill-conditioned and thus we need to know the Gram matrix to a high precision. This typically demands a large q-truncation order d q . For example, in Figure 3(a) and Figure 4 in section 3.1, values of d q between 60 and 100 were used. In Figure 3(b), however, we needed to go to d q = 200. In general, the required value of d q grows with N , similarly to what we observed in semidefinite problems. On the other hand, the linear method computes faster than the semidefinite one, which allows us to study much higher values of N . 23 For example, in order to do calculations for c = 8, ∆ φ = c−1 12 , N = 25, it is necessary to compute the q-expansion of the conformal blocks to the order q 100 (see below). At this order deg P 25,0 = 679 and Q has 327 zeroes.
There is a small subtlety in the computation of v N for the mixed correlator, due to the fact that Q(∆) has a double zero at ∆ = 0, which in the case ∆ 1 = ∆ 2 is not cancelled by zeroes of P n,m . In this case (e.g. in the Figure 4) we have tried two approaches. The first approach is introducing lower bound on the intermediate scaling dimension ∆ gap below the Liouville threshold. The second approach is to modify equation (D.8) as (D.14) While the obtained results differ slightly at small N , already at N = 13 both provide equally good approximations for the Liouville spectral functions in Figure 4.

D.2.1 Numerical checks of completeness
Here we consider the question of completeness of the systems with respect to the measure described above. In both cases we attempt to approximate the step functions θ ∆ * P N θ ∆ * , where P N is the projection onto the subspace spanned by elements of either system with n + m ≤ N , and compute the residual errors (D. 17) We do this for a range of ∆ * in the case of the mixed correlator with c = 8, ∆ 1 = ∆ 0 , ∆ 2 = 12 7 ∆ 0 . The results are shown in Figure 13 for B even and in Figure 14 for B odd , consistent with the completeness of both bases. In the plots of E N as a function of N −1 , we have rescaled E N (∆ * ) by an N -independent factor for each sample value of ∆ * (denoted by E N ) so that for all ∆ * the slope of the linear fit with N −1 (which appears to be valid asymptotically for large N ) is approximately 1.

D.3 Bounds from a reduced basis of linear functionals
Here we consider the bounds on the scalar-only spectral function using the following reduced basis of linear functionals simplicity of the reduced basis it is now possible to access bounds at higher N within the same computing time. This provides a useful arena to study the convergence of the bounds at large derivative orders, with the caveat that the N → ∞ limit of the bounds obtained from the reduced basis are likely weaker than the optimal bounds from the most general linear functionals. If the latter is the case, one may eventually need to relax the restriction on min(m, n) in (D.18).