Fast Conformal Bootstrap and Constraints on 3d Gravity

The crossing equations of a conformal field theory can be systematically truncated to a finite, closed system of polynomial equations. In certain cases, solutions of the truncated equations place strict bounds on the space of all unitary CFTs. We describe the conditions under which this holds, and use the results to develop a fast algorithm for modular bootstrap in 2d CFT. We then apply it to compute spectral gaps to very high precision, find scaling dimensions for over a thousand operators, and extend the numerical bootstrap to the regime of large central charge, relevant to holography. This leads to new bounds on the spectrum of black holes in three-dimensional gravity. We provide numerical evidence that the asymptotic bound on the spectral gap from spinless modular bootstrap, at large central charge $c$, is $\Delta_1 \lesssim c/9.1$.


Introduction
The conformal bootstrap is a rapidly growing array of techniques to solve strongly interacting quantum field theories. It had early success in two spacetime dimensions [1], and more recently has proved successful in higher dimensions [2]. An important impetus has been the development of numerical techniques to analyze the crossing equations, starting with the introduction of the functional bootstrap method in 2008 by Rattazzi,Rychkov,Tonni,and Vichi [3]. This method uses linear (or semidefinite [4]) programming to constrain solutions to the crossing equations, thereby carving out the allowed parameter space of unitary CFTs. It has also been used to solve individual CFTs, including the 3d critical Ising model [5,6]. The method relies on finding extremal functionals, often numerically, which bound the space of CFTs, and in certain cases encode the spectrum of the target theory [7,8].
Despite this exciting progress, the current numerical methods are in their infancy, and there are many potential areas for improvement, both technical and conceptual.
For example, with semidefinite programming, there is no clear way to leverage analytic results, such as the spectrum at high spin, to improve the numerics [9]. This is a tantalizing hint that much more efficient algorithms are waiting to be discovered.
Such algorithms, combined with powerful analytics, will be needed to address certain physics problems, such as solving theories with large global symmetries [10,11], pushing bootstrap into the hydrodynamic regime [12,13], or exploring the boundaries of quantum gravity by combining numerical bootstrap with AdS/CFT in the black hole regime.
In this paper, we develop a new algorithm for numerical bootstrap, and apply it to 3d gravity. The method is based on solving truncated crossing equations. It builds on one of the extremal functional methods of El-Showk and Paulos [14], and is very similar in many respects, but differs in some important details that extend the reach by several orders of magnitude. In principle, these methods are general, but in practice, they are limited to a certain class of bootstrap problems, where one can easily generate approximate solutions to the truncated equations. We will show that this is possible for the modular bootstrap problem relevant to 3d gravity.
The modular bootstrap is a variant of the conformal bootstrap introduced by Hellerman to constrain the spectrum of 2d CFTs [15]. Hellerman used it to prove analytically that any theory of 3d gravity has a primary operator with scaling dimension ∆ c/6.
From a holographic point of view, this bound is perhaps weaker than expected, because all known theories of 3d gravity have BTZ black holes [16] starting at ∆ ∼ c/12. The bound can be systematically improved by numerics, but the difficulty grows rapidly with c [17,18]. Previous numerical calculations were effective only for c 100 and did not reach the semiclassical regime -the bound had not converged in c, and it was not possible to estimate the asymptotics. A stronger bound is worth pursuing, because if it could be pushed down to the black hole threshold, then the numerics could also potentially be used to construct or rule out candidate theories of 3d gravity, or to suggest new analytic methods.
Using our algorithm, we extend the modular bootstrap up to c ∼ 1800. This central charge appears to be in the semiclassical regime, where the bound has nearly converged, and we estimate the asymptotic bound to be ∆ c/9.1. This estimate relies on an extrapolation in c, which appears to be reliable, but strictly speaking it is impossible to rule out a stronger bound as c → ∞. Our bound is not at the black hole threshold, so unfortunately cannot be used to construct candidate theories of 3d gravity or to address the questions above. We have not included the constraints from spin; it is possible that incorporating spin would lead to a stronger bound that saturates the black hole threshold.
The functional bootstrap is an expansion in derivatives, truncated at some order N deriv . Our calculations are performed up to derivative order N deriv ∼ 5000. For a given c, this takes several CPU-hours on a standard laptop. It is impossible to make a direct comparison to existing methods, but a rough extrapolation suggests that semidefinite programming would require at least 10 7 CPU-hours to perform the same calculation.
A basic version of our algorithm, sufficient for modular bootstrap, is straightforward to implement in Mathematica with a few dozen lines of code.
On the down side, our methods do not generalize in any obvious way to bootstrap problems with spin. This has (so far) prevented us from applying the algorithm to the spinning modular bootstrap, or to the 3d critical Ising model.
We begin with a review of the crossing equations and a summary of the truncation in section 2.
We explain under what circumstances the truncated equations can be used to place general bounds. In section 3 we review the functional method, introduce a new method to optimize bootstrap functionals, and show that it is Lagrange dual to the truncation method. Section 4 describes the algorithm, and in section 5 it is applied to modular bootstrap and 3d gravity.
Comparing our results to semidefinite programming requires very accurate bounds.
In a standalone appendix C, we describe a method to find high precision bounds with linear programming, which has some advantages over the standard bisection method.
Although we have stressed the computational aspects, it seems likely that the truncated crossing equations can also be solving analytically in various limits -lightcone [19,20], large N [21], short distance [22], etc. This could be used to derive analytic bounds, to jumpstart the numerics, or to develop variational methods that go to much higher orders. We consider it likely that with more refinement, numerical methods can be extended to functionals with millions of derivatives. Further discussion of the outlook, including the challenges presented by problems with spin, is in section 6.

Truncating the Primal Bootstrap
The functional bootstrap is an optimization problem, and as such, it has a primal and dual formulation. The primal bootstrap refers to the crossing equations themselves.
The dual approach uses linear functionals. We will rely on both pictures, starting on the primal side. In this section and section 3, we consider a broad class of spinless bootstrap problems. We will specialize later to modular bootstrap.

Setup
Consider a crossing equation of the form ∆, where is spin, and W ∆, is the difference of a conformal block and a crossed block, Unitarity requires a ∆ > 0, and places lower bounds on ∆. Here we use the notation of the correlator bootstrap, with z andz conformal cross ratios, but the equations apply also to modular bootstrap, by replacing (z, 1 − z) → (τ, −1/τ ) with τ the torus modulus.
Following [3], let us act on (2.1) by linear functionals (The methods may also apply to more general bases; see discussion section.) Choose a set {F i } i=1...P of P such functionals, and let where for each , the f i are linearly independent. Applying F i to the crossing equation (assuming the F i are sufficiently well behaved to bring the derivatives inside the sum).
These are the primal bootstrap equations. There are P equations for an infinite number of parameters: the spectrum of (∆, ) and the coefficients a ∆ . After absorbing a positive prefactor into the a ∆ , these equations are well approximated by polynomials [2].
We will restrict our attention to bootstrap problems without spin dependence. In this case we set z =z (or τ = −τ in the modular bootstrap), and drop the spin label . This class of problems includes spinless modular bootstrap, 1d correlators, the sl (2) subsector of correlator bootstrap in higher dimensions.

Truncation
Let us take P even, and truncate to P/2 + 1 states, including the identity term with a 0 ≡ 1, ∆ 0 ≡ 0. This results in the truncated primal equations: This is now a closed system of P polynomial equations for P unknowns.
This seems a rather brutal truncation of the crossing equations, so at this point, it is not clear that it is useful. Nonetheless we will show that, sometimes, solutions of (2.6) can be used to place bounds on the space of all unitary CFTs.
In the numerical bootstrap, the goal is to maximize some physical quantity, such as the first scaling dimension ∆ 1 , over all unitary solutions to the P equations (2.5). In practice, the optimal solution has a finite number of states in the sum. If it happens that the optimal solution has exactly P/2 + 1 states, then it can be found by solving the closed equations (2.6). We will show that this is the case for 2d modular bootstrap at zero angular potential.
This leads to our main analytic result: In modular bootstrap, at any order P and central charge c, there is a solution to the truncated bootstrap equations (2.6), with scalar gap ∆ It follows that general bounds on all unitary CFTs can be found by solving (2.6). To derive this, we will show that a solution to (2.6) can be used to build an extremal functional, identical to that obtained in the usual numerical bootstrap. Note that this claim does not apply to all solutions of (2.6), but to one special solution, which must be carefully sought among many.
This is based on the same idea used for numerical bootstrap of 1d correlators in [14].
Other approaches to the correlator bootstrap using truncated crossing equations have also been discussed in [23][24][25].

Comments on Monotonicity
The largest dimension participating in (2.6), ∆ P/2 , grows with P . Therefore the truncation order P acts as a UV cutoff. Qualitatively, there is a parallel between changing P → P − 2 and a renormalization group transformation: As P is decreased, the effect of high-dimension operators must be absorbed into the low-dimension spectrum and OPE coefficients, in order to preserve the crossing equations. A corollary of the result stated above is monotonicity: Under reducing the cutoff P → P − 2, the gap ∆ This point of view is useful to keep in mind when interpreting solutions to the truncated equations, because we should not expect operators with µ ∼ P/2 to have their physical scaling dimensions; they must also encode the effects of the UV states that are removed by the truncation. We will see this effect in the numerics.

Dual Bootstrap
In this section, we derive the claims of section 2, by relating the primal equations to the functional bootstrap. We will first recast the functional bootstrap as a nonlinear optimization problem, then dualize to the primal equations.

Review
The functional method [3] relies on the following observation. Consider the action of a linear functional on W ∆ , The α i parameterize the functional. If we find a functional that is positive on the vacuum, and on all dimensions above some gap, then every term in the sum (2.1) is positive for dimensions in this range. Therefore, to satisfy crossing, every unitary CFT must have an operator at or below ∆ 1 .
The strongest bounds come from minimizing ∆ 1 . For our purposes, a useful way to state the optimization procedure is as follows. Fix P , the gap ∆ 1 , and the normalization condition α P = 1, and: 1 Here and below, the indices run over the ranges (3.3) is a linear program with an infinite number of constraints, or can be recast as a semidefinite program [4]. If the constraints are feasible, and the objective f (0) is positive at the maximum, then the initial choice of ∆ 1 is ruled out. This procedure is repeated, for different values of ∆ 1 , until the theory is marginally excluded, i.e., f (0) = 0 (to some desired precision). This is the edge of the exclusion region, and the optimal functional at this point is called the extremal functional. The dual bootstrap refers to the problem of finding the extremal functional.

Functionals parameterized by zeroes
The zeroes of the extremal functional correspond to solutions of the primal equations (2.5). This is a general fact of linear programs, which has been used to find the spectrum of various CFTs, including the critical Ising model in two [8] and three [6] dimensions. In certain spinless bootstrap problems, including modular bootstrap, the extremal functional has a very simple pattern of zeroes. Empirically, for even P , there are single zeroes at ∆ = 0 and ∆ = ∆ 1 , and double zeroes at P/2 − 1 additional states. An example is plotted in figure 1. This set of P/2 + 1 zeroes corresponds to a solution of the truncated primal equations (2.6).
We will make an ansatz for the functional by assuming that this pattern of zeroes continues at higher P . We will confirm a posteriori that this assumption holds for modular bootstrap, at least for all c and P that we consider.
Given this assumption about the pattern of zeroes, we now replace the linear optimization problem (3.3) by a nonlinear optimization. Choose a spectrum ∆ µ for µ = 1 . . . P/2. From this spectrum, construct the unique functional obeying This is P − 1 linear equations for the P − 1 variables α a , so the solution is unique. In other words, we can view the functional as parameterized by its roots. Next, the linear The optimization step can also be eliminated altogether by solving directly for the extremal functional. This has f (0) = 0 and is an extremum with respect to ∆ µ . So we simply append to (3.5) the equations where in the latter equation, µ = 2 . . . P/2, and the gradient is taken along solutions to (3.5). That is, α a = α a (∆ µ ) as determined by (3.5).
This is now a closed system of equations for the extremal functional. There is no guarantee that the positivity conditions will be obeyed in this formulation, but this is easily checked after the functional has been constructed, and holds in the applications we will consider. A similar system of equations, but for the linear version of the problem and involving also the primal variables, appears in [14].
To summarize, we have formulated the problem of finding the extremal functional as 3P/2 − 1 equations (3.5) and (3.7) for the P/2 dimensions ∆ µ and P − 1 functional parameters α a . In the derivative basis, the equations are polynomial (after absorbing a positive ∆-dependent factor into the coefficients a µ ).
This reduces the problem of finding the extremal functional to the problem of finding zeroes of a system of polynomial equations. The degree of the polynomials grows with P , so this is non-trivial. Moreover, we have replaced a convex optimization by a nonconvex root-finding problem, where there may be runaways and spurious solutions. In practice, these obstacles can be overcome if we can generate a good guess to use in the first step of Newton's method. We will describe how to do this for modular bootstrap in section 4.

Duality
As the nomenclature suggests, the truncated versions of the primal and dual bootstrap are, in fact, Lagrange duals. With our nonlinear formulation of the dual problem, this duality holds under the assumption that the double roots ∆ µ of the extremal functional have multiplicity exactly two -not higher.
The duality immediately proves the main results stated in section 2. In the dual bootstrap, monotonicity of the gap under changing P is obvious -if we allow a more general functional, the bound on ∆ 1 can only get stronger. When translated into primal language, this monotonicity property becomes more surprising.
The duality that we describe is already very well known in the linear formulation of the problem. We will rederive it in terms of the nonlinear version above.
We will start with the dual bootstrap -i.e., the problem of finding extremal functionals (3.6) -then dualize to the primal bootstrap (2.6). The constrained optimization (3.6) is equivalent to maximizing the Lagrangian over the independent variables α a , ∆ µ and the Lagrange multipliers a µ , λ ν (with µ = 1 . . . P/2, ν = 2 . . . P/2.) Using the normalization condition α P = 1 and rearranging terms, To dualize, we reinterpret this Lagrangian, viewing ∆ µ , a µ , λ ν as independent variables, and α a as a Lagrange multiplier. Assuming f i (∆ µ ) = 0 so that no roots have multiplicity higher than two, the equation of motion for ∆ µ imposes λ ν = 0. Therefore, the primal optimization problem is These are the same sums as appear in the primal crossing equation (2.6), so the constraints consist of the first P −1 primal equations. Finally, let us also impose extremality, f (0) = 0, so the assumed gap is marginally excluded. This is equivalent to the P th primal equation.
The conclusion is that the problem of finding the extremal functional is precisely dual to the problem of solving the P truncated crossing equations, (2.6), for the spectrum and coefficients. In appendix A, we derive the same statement by direct analysis of the extremal functional equations (3.5)-(3.7), and in the process, give a formula for the truncated OPE coefficients a µ in terms of the extremal functional. It is also guaranteed that the bounds we derive are optimal, in the sense that linear programming at the same value of P will always produce exactly the same bounds. It is impossible to derive a stronger bound (at a given P ), because doing so would contradict the primal solution that we have found. In other words, the method we have described produces solutions which are primal-dual feasible, and therefore optimal.

Modular Bootstrap Algorithm
So far, our discussion has applied to a general class of spinless bootstrap problems, with a certain pattern of zeroes in the extremal functional. We now focus on the modular bootstrap exclusively. In this section we first review the (spinless) modular bootstrap [15,17], then describe our algorithm in detail.

Setup
Consider the partition function of a 2d CFT at zero angular potential, with β > 0, and c the central charge. In the second line we have organized the sum into characters of the algebra, Virasoro × Virasoro. 2 Crossing symmetry in this context is modular invariance: Z(β) = Z(1/β). It is convenient to define the reduced partition where ρ µ is the degeneracy of primaries with dimension ∆ µ . This function is also and, for the vacuum, (The extra factors in the vacuum block are to account for the null state at level one, , with χ V ir h (τ ) the standard chiral Virasoro character.
We choose the basis of derivative functionals for k = 1 . . . P . This is a convenient basis because the resulting f k are simply odd-index Laguerre polynomials, L 2k−1 . Indeed, acting with F k on the crossing equation (4.7) gives where and, for ∆ > 0, The vacuum state is different due to the extra factors in (4.6), with x 0 = − c−1 12 . At this point, we have written the primal bootstrap equations in the same format as section 2. The truncated version, a closed system of P equations for P unknowns, is P/2 µ=0 a µ f k (∆ µ ) = 0 , (4.14)

Algorithm
We now turn to numerical algorithms based on the reformulation of the crossing problem described in sections 2-3. We focus on the simplest bootstrap question, which is to 3 Derivation: The generating function for Laguerre polynomials is Therefore acting with 1 p! (∂ t ) p | t=0 on the right-hand side gives L p (2w)e −w . With the identification β = 1+t 1−t , this is the same as F k [W ∆ (β)]. The even-index Laguerre's do not appear because the corresponding differential operator vanishes on the modular-odd function W ∆ (β).
place an upper bound on the scalar gap ∆ 1 , given the central charge c. The algorithms can be modified to accommodate similar problems.
We can consider two main approaches: • Dual : Solve the functional equations (3.5) and (3.7), for the spectrum and the functional.
• Primal : Solve the P truncated primal crossing equations, (2.6), for the spectrum and degeneracies.
Note that standard optimization methods for functional bootstrap require an additional scan over ∆ 1 . In our case, ∆ 1 is a free parameter, so this step is not necessary.
For modular bootstrap, both the dual and primal algorithms are much faster than semidefinite programming. In the rest of this paper we will consider only the primal algorithm, which was somewhat faster in our preliminary testing, and conceptually simpler.
The purpose of the algorithm is to solve the equations (4.14) for some large value of P . The strategy is to start at small P , solve the equations by Newton's method, then gradually increase P , using the previous results to generate a good initial guess for Newton's method. The same strategy was used for 1d correlator bootstrap in [14] (specifically, our algorithm is similar to the 'upgrading + error correction' method).
Besides the differences between 1d correlator bootstrap and 2d modular bootstrap, our algorithm has two significant improvements: better numerical stability, and a more elaborate, more accurate method to guess the initial points. This will allows us to increase P in large increments.
For numerical stability, we first introduce normalization factors to rescale the P × P/2 matrix f k (∆ µ ). Rewrite (4.14) as The matrixf k (∆ µ ) is scaled such that the vacuum column is all 1's, and every other column has an absolute maximum entry of ±1.
We now apply Newton's method to the equations (4.15) for the P unknowns (ã µ , ∆ µ ). The Laguerre polynomials, and their derivatives, are calculated by 3-term recurrence relations. This is faster and more numerically stable than direct evaluation.
Our implementation of Newton's method is completely standard so we will not review it.
An important element in the algorithm is a guess generator, to provide the initial point in Newton's method. The method is largely insensitive to the initial guesses for a µ , so for these coefficients we always guessã µ = 1. Convergence is very sensitive to the starting ∆ µ , so these must be chosen carefully. We start by solving the equations for P = 4, 6. For these initial points, it is easy to find a good guess by hand (or using a linear/semidefinite optimizer such as SDPB [26]). Then we proceed upwards in P recursively, using previous results to generate a guess for ∆ µ .
Our guess generator is based on the observation that the spectrum, at different values of P , is self-similar. That is, the curve ∆ µ as a (discrete) function of µ has almost exactly the same shape for all P , once P is large enough. This is illustrated by the plots in fig. 2, which are the exact spectra at P = 110 and P = 200 for the c = 12 modular bootstrap. The shapes in the two plots are almost identical; only the scaling of the axes and the state spacing differs. This motivates a guess where we rescale the lower-P spectrum, and re-sample at the appropriate points. This generates a good guess, and we use precisely this method for P 200. However, we can do better, and take larger P -steps, by correcting for the fact that the spectrum is not quite self-similar.
We do this by parameterizing the shape curve, fitting the parameters as functions of P , and extrapolating. The full algorithm is given in appendix B.
The resulting guesses are very accurate. For example, in the c = 100 modular bootstrap, using data from P = 120, 122, 124, . . . , 200, the extrapolation to P = 280 produces a spectrum ∆ guess µ that differs from the exact spectrum by under one part in 10 3 -10 5 . Results at higher P are comparable. This means that we never need to run more than a handful of Newton steps to converge to an accurate spectrum.
This completes the algorithm. The only computationally expensive step is inverting the Jacobian in Newton's method, which scales roughly as P 3 . The computations are performed in Mathematica with high precision arithmetic. (right). The horizontal axis is the state number, µ = 1, 2, . . . , P/2. Note that the curves are approximately the same shape. This observation is used to generate very accurate guesses for the initial point in Newton's method.

3d gravity and summary of existing bounds
Three-dimensional gravity in the semiclassical limit is holographically dual to a 2d CFT with a large central charge, c 1. Quantum fields in the bulk, aside from the graviton, are dual to low-dimension operators in the CFT. It follows that the CFT should also have a sparse spectrum of low-dimension operators, in that the number of states below any fixed ∆ * should be finite as c → ∞.
'Pure' 3d gravity, which consists only of the graviton plus black holes, would have no new primaries between the vacuum and black hole threshold. That is, ∆ 1 ∼ c 12 . It is unknown whether pure gravity exists as a quantum theory. A putative partition function was suggested by Witten [27], but explicit calculations by Maloney and Witten [28] gave a different answer that was incompatible with a holographic interpretation.
It is possible that this calculation can be modified to produce a consistent quantum theory [29], but at this point, the realm of possibilities is poorly understood.
Even the simplest question is unanswered: What is the largest gap ∆ 1 compatible with diffeomorphism invariance? Large diffeomorphisms in 3d gravity correspond to modular transformations in the dual CFT, so this is a question for the modular bootstrap.
The functional approach to modular bootstrap was initiated in [15], where, working analytically at P = 2, it was proved that ∆ 1 ≤ c 6 + 0.474. This is known as the Hellerman bound. At large c, the Hellerman bound is a factor of 2 above the black hole threshold. It has also been shown analytically that, for ∆ c 6 , 2d CFTs with a large gap have the same entropy as BTZ black holes [30]. However, the range c 12 ≤ ∆ ≤ c 6 remains enigmatic. This range certainly does not have a universal spectrum, even in holographic theories, but it is an open question whether there must be any states at all within this window.
Numerical evidence strongly suggests that the asymptotic bootstrap bound can be improved to somewhere between c/12 and c/6 [17,18,31]. (For related work, see [32][33][34][35][36][37][38][39][40][41][42][43][44][45][46] needed to find the bound in the semiclassical, large-c limit. The difficulty is that for large c, the bound on ∆ 1 converges slowly with P , so large c requires large P . This is a general obstacle to bootstrapping quantum gravity in the regime ∆ ∼ c, which presents a significant challenge to probing black hole physics with the bootstrap. (On the other hand, many interesting questions in AdS/CFT already arise for states with ∆ c, and in these cases numerical bootstrap has been very successful.) The largest truncation order considered in the literature on modular bootstrap is P = 92 [18], using the semidefinite program solver SDPB [26]. We apply our algorithm up to P ∼ 2250 (derivative order 4500). At small P , we have checked that our bounds agree exactly with those produced by SDPB.

Bound as a function of c
The upper bound on ∆ 1 /c is plotted in fig. 3. The dots are numerical data, and the solid line is the extrapolation to P = ∞ with an exponential fit, ∆ 1 (P ) = ∆ ∞ 1 + ae −bP . To estimate the asymptotic bound, we fit the data in the range 300 ≤ c ≤ 700 to a function of c. In this range, the bound has converged very close to its P → ∞ limit. The best strict (non-extrapolated) bound on ∆ 1 /c is c = 1800, P = 2310, for which In addition to the gap, ∆ 1 , the method also provides a spectrum with over a thousand operators. As an example, the spectrum for c = 500, P = 2310 is plotted in fig. 4. As discussed in section 2.3, only the operators µ P should be trusted, but in practice, roughly the lower half of the spectrum with µ P/4 is converged.  This coincidence will be explained analytically in [47]. Note that we did not impose invariance under τ → τ + 1, so it is surprising to find the j-function appearing.
Extrapolating the runtime of SDPB to large P may be unreliable, but naively fitting to a power law and setting P = 2000 gives the estimate ∼ 10 9 s.
We can also compare to the extremal functional method of El-Showk and Paulos [14]. They solve a different problem (1d correlators) so it is not a direct comparison, but show every 15th scaling dimension to reduce clutter. The smaller inset figure shows the first 30 scaling dimensions with deviation away from linearity. In [47], the linear regime of the spectrum with ∆ µ = c−4 8 + µ (dashed line) will be explained analytically.
they report runtime of 45 minutes at P = 150. For modular bootstrap, our algorithm runs in seconds at the same P . The algorithm is very similar to that described in [14], so we do not know the exact origin of the speedup, but it presumably comes from some combination of improvements to numerical stability, the guess generator described in precisions from 16 to 500 digits, chosen to ensure that the matrix inversion step in Newton's method is well behaved.

Discussion
In summary, we have described a reformulation of the modular bootstrap equations, amenable to both analytic and numerical analysis. On the analytic side, the main logic was to show that (i) extremal functionals can be efficiently parameterized by their roots, and (ii) this is dual to solving the truncated crossing equations.
For modular bootstrap, our algorithm is orders of magnitude faster than semidefinite programming or any other previous method, extending the maximum truncation order from P ∼ 100 to P ∼ 2500. Other potential advantages include: • It is numerically stable, so intermediate steps do not require very high precision.
• The method might apply to non-derivative or non-polynomial functionals. (One choice has been considered [48].) We restricted to derivative functionals only because in this basis, we have confirmed explicitly (in examples) the pattern of zeroes invoked in section 3. The method applies to any basis of functionals with a similar pattern of zeroes.
• The algorithm can be sped up significantly by hot starting, i.e., inputting a guess for the spectrum of the target CFT. This means that analytic results can potentially be leveraged for faster numerics, similar to the way we used results at small P for faster convergence at large P .
• Although the strict bounds we derived apply only to unitary CFTs, the truncated primal equations may also have solutions corresponding to non-unitary CFTs, with potential applications to weakly first order phase transitions [2,49]. This would be similar to the approach taken in [23,24].
• The truncated bootstrap equations (2.6) can be studied analytically. This may lend additional insights.
There is also a major drawback to our approach, at least in our current implementation. We considered only bootstrap problems without spin dependence. More generally, for example in the modular bootstrap with angular potential, or the correlator bootstrap of the 3d Ising model, there are more constraints on the extremal functional, labeled by . This typically leads to a different pattern of zeros, not the simple pattern we observed here. In primal language, the problem is that with spin, the solution of the primal optimization problem at truncation order P typically does not have exactly P/2 + 1 states. 5 This makes it impossible to truncate to a closed, primal-only problem, as we have done here. Instead one can write primal-dual equations for the extremal functional, involving (α a , a µ , ∆ µ ) simultaneously, as in [14]. In principle, we could use Newton's method to solve this enlarged problem, but in practice, we have not found any efficient way to generate initial guesses. The pattern of states does not change predictably as the truncation order is increased.

A Direct derivation of optimization duality
In this section, we re-derive the duality of section 3.3 by direct analysis of the extrema of the functional optimization (3.6). This is useful in practice for translating between dual and primal solutions. Denote and form the block matrix where i = 1 . . . P , µ = 1 . . . P/2, ν = 2 . . . P/2, and β = 1 . . . P − 1. The roots (3.5), using (3.1), impose where we have used the normalization condition α P = 1. Therefore, given ∆ µ , the functional with the correct zeroes is where V −1 denotes the inverse of the square submatrix V aβ (recall a is restricted to the range 1 . . . P − 1). The functional evaluated on the identity term is According to (3.6), we seek to maximize this over ∆ µ . Ultimately we are interested in a maximum with f (0) = 0, so we can assume the maximum is actually obtained, and therefore occurs at a critical point. The gradient, with ∇ µ ≡ d d∆µ , is Therefore we have two options: f iµ = 0, so some ∆ µ is actually a triple rather than double root, or V −1 βa f a (0) = 0 (β = P/2 + 1 . . . P − 1) (A.10) We assume the latter, so that no roots have multiplicity higher than two. This is precisely the first P − 1 equations of the primal bootstrap, (2.6), with a µ the OPE coefficient. Now let us impose extremality, i.e., that we are at the boundary of the excluded parameter space. That is, we let ∆ 1 be an independent variable and impose 0 = f (0) = α a f a (0)+f P (0) = −V P β V −1 βa f a (0)+f P (0) = −V P β c β +f P (0) = V P µ a µ +f P (0) . (A.14) This is the P th equation of (2.6).
Therefore, the solution of (3.6), at extremality, is equivalent to a solution of the P primal bootstrap equations (2.6) (assuming no roots of multiplicity greater than two).
We have also shown that, given the extremal functional, the OPE coefficients in the

B Generating guesses for Newton's method
The logic of our guess generator is described in section 4.2. We fit to the curves in fig. 2, extrapolate in P to generate a new curve, then sample at the appropriate discrete points. In this appendix we describe the complete algorithm. The choice of parameters and fitting functions that we use is ad hoc, but produces good results.
Suppose we have already found the spectrum for various truncations P = P A , A = 1, . . . , K − 1. Denote these spectra by ∆ P A µ . Define ∆ P A (x) to be the continuous function that numerically interpolates this spectrum at integer x, i.e., ∆ P A (µ) = ∆ P A µ for µ = 1, 2, . . . , P A /2 . (B.1) We use a degree-22 Hermite interpolation (in Mathematica: Interpolation[..., InterpolationOrder->22]). Now suppose we want to guess the spectrum at truncation order P K . Let This is defined at discrete P A . Fit this this, as a function of P A , to the continuous function F µ (y) = a 1 + a 2 y + a 3 log(y) + a 4 log(y) 2 + a 5 log(y) 3 . (B. 3) The parameters a i are chosen such that F µ (P A ) ≈ F P A µ , with least-squares fitting. Finally, our guess for the spectrum at P = P K is ∆ guess µ = F µ (P K ) . (B.4)

C High precision bounds with linear programming
In this appendix we describe a simple method to find high-precision bounds on the gap, ∆ 1 , with linear programming. This is used when we compare our results to SDPB, but it is more general. For example, it is faster than the usual method to find bounds on the 3d Ising model (the 'kink' in [5]), if the bounds are needed to high accuracy.
The standard method is to run a program solver such as SDPB in feasibility mode, to rule in or out each gap ∆ 1 , and to zoom in on the optimal ∆ 1 using bisection.
Instead, we run the linear program (3.3). Call the objective at the optimum M (∆ 1 ).
To find the optimal bound, i.e., the gap ∆ 1 which is marginally excluded, we use the secant method to solve M (∆ 1 ) = 0.
In the first few steps, this method is slower, because the program solver takes longer to solve an optimization problem than a feasibility problem. But the secant method converges exponentially faster than bisection, so for high accuracy bounds on ∆ 1 , this method is faster.