The Lightcone Bootstrap and the Spectrum of the 3d Ising CFT

We compute numerically the dimensions and OPE coefficients of several operators in the 3d Ising CFT, and then try to reverse-engineer the solution to crossing symmetry analytically. Our key tool is a set of new techniques for computing infinite sums of SL(2,R) conformal blocks. Using these techniques, we solve the lightcone bootstrap to all orders in an asymptotic expansion in large spin, and suggest a strategy for going beyond the large spin limit. We carry out the first steps of this strategy for the 3d Ising CFT, deriving analytic approximations for the dimensions and OPE coefficients of several infinite families of operators in terms of the initial data $\{\Delta_\sigma,\Delta_\epsilon,f_{\sigma\sigma\epsilon},f_{\epsilon\epsilon\epsilon},c_T\}$. The analytic results agree with numerics to high precision for about 100 low-twist operators (correctly accounting for O(1) mixing effects between large-spin families). Plugging these results back into the crossing equations, we obtain approximate analytic constraints on the initial data.


Introduction
Despite the ubiquity of conformal field theories (CFTs) in d > 2 spacetime dimensions, very little is known about their operator dimensions and OPE coefficients away from simplifying limits like large central charge (large-N ) or weak coupling. Unlike in d = 2 dimensions [2], we have no nontrivial exactly-solvable CFTs in d > 2 from which to draw lessons.
In this work, we produce a new numerical picture of the spectrum of the 3d Ising CFT, including about 100 operators, and use it as a guide to explore the theory analytically. In addition to the intrinsic interest of the 3d Ising CFT for its role in second-order phase transitions, our motivation is to develop analytical tools for solving crossing symmetry in general (and eventually apply them to wider classes of theories).
The current most powerful techniques for studying the spectrum of small central charge theories are numerical bootstrap techniques , based on the conformal bootstrap [53,54] and the methods pioneered in [3]. For example, the numerical bootstrap has yielded precise predictions for dimensions of the lowest-dimension scalars σ and in the 3d Ising CFT [12,20,24,31,55]. It is difficult to reproduce these results analytically because the 3d Ising CFT does not admit a (known) controlled expansion in a small coupling constant. 1 But even strongly-coupled theories admit small parameters in kinematic limits. The authors of [59,60] showed that every CFT admits a large-spin expansion, accessible via the lightcone limit of the crossing equations. By studying the lightcone limit, one can prove: Theorem 1.1 (Existence of double-twist operators [59,60]). Suppose a CFT in d > 2 dimensions contains primary operators O 1 , O 2 with twists τ 1 , τ 2 . 2 For each n = 0, 1, 2, . . . , there exists an infinite family of primary operators with increasing spin and twists approaching τ 1 + τ 2 + 2n as → ∞.
Schematically, these operators are (1.1) Of course, composite operators like (1.1) don't make sense in a general strongly-coupled theory. However, theorem 1.1 implies that they do make sense in the large-limit. We denote the family with twist approaching τ 1 + τ 2 + 2n as [O 1 O 2 ] n and refer to such operators as "double-twist" operators (following [60]).
Dimensions and OPE coefficients of double-twist operators have a computable expansion in (generically non-integer) powers of 1/ , where terms in the expansion come from matching operators on the other side of the crossing equation. Recently, there has been significant progress in understanding this expansion [59-62, 1, 63-67]. The large-expansion is asymptotic in general [1], so its usefulness for studying finite-spin operators is not immediately clear. Nevertheless, we might hope that large-spin techniques could enhance numerics or vice versa. Perhaps an analytical solution of the large-spin expansion could help make numerics more efficient, or even replace numerics entirely if crossing symmetry could be solved via the lightcone limit.
With our concrete numerical calculations as a guide, we find the following: • Double-twist operators play an important role in the numerical bootstrap.
• By truncating the asymptotic large spin expansion, and with the help of some new analytical techniques described below, we can describe a large part of the 3d Ising spectrum, including operators with spin as small as = 2 or = 4.
• The large-spin expansion can be used to solve crossing symmetry systematically in a "double lightcone" expansion in z, 1 − z. 3 • The "errors" associated with the fact that the expansion is asymptotic can be precisely characterized (they are "Casimir-regular" terms defined in section 4). Requiring that they cancel gives nontrivial constraints on the spectrum.
Let us describe the structure of this paper in more detail.
In section 2, we perform a (non-rigorous) numerical computation of the 3d Ising spectrum using the extremal functional method [7,14,20]. Importantly, we use a trick from [68] which lets us assign error bars to the resulting operator dimensions and thereby understand which predictions are robust and which ones aren't. The robust predictions turn out to be for lowtwist operators (not just low-dimension operators). Specifically, we find relatively precise predictions for 112 operators in the 3d Ising CFT, of which only 9 do not fall into an obvious double-twist family. The remaining 103 operators give a clear numerical picture of the families [σσ] 0 , [σσ] 1 , [ ] 0 , and [σ ] 0 , up to spin ∼ 40. We give additional details of our computation in appendix A, and list the resulting operators in appendix A. 3. Although many of the results in this work are analytical (and applicable to any CFT), this numerical picture is a crucial guide, helping us ask the right questions and find the right tools to answer them.
We then set out to describe the families [σσ] 0 , [σσ] 1 , [ ] 0 , and [σ ] 0 analytically using the large-spin expansion. To succeed, we must develop two new technologies: • Techniques for summing an infinite family of large-spin operators in the conformal block expansion. (For example, this lets us compute the contribution of a twist family to its own anomalous dimensions.) • Techniques for describing mixing between multi-twist families.
Our key tool is a better understanding of infinite sums of SL(2, R) conformal blocks, which we develop in section 4 (after reviewing the lightcone bootstrap in section 3). By generalizing the conformal block expansion of 1-dimensional Mean Field Theory, we show how to compute exactly, and in great generality, sums of SL(2, R) blocks in an expansion in the crossed channel z → 1 − z. A simple example is 4 The crucial point is that the first term on the right-hand side, 1−z z a , becomes arbitrarily singular at z = 1 after repeated application of the quadratic Casimir of SL(2, R), while the remaining terms do not. We compute general sums of SL(2, R) blocks by exploiting this distinction. Because SO(d, 2) conformal blocks are sums of SL(2, R) blocks, equation (1.2) and similar identities can be used as building blocks for understanding crossing symmetry in general. Using them, we solve the asymptotic lightcone bootstrap to all orders (for both OPE coefficients and anomalous dimensions) in section 5.
In section 6, we explore how well the truncated large-spin expansion describes the families [σσ] 0 and [σ ] 0 . Surprisingly, we find that the first few terms (coming from and T µν in the crossed-channel) fit the numerical data for [σσ] 0 beautifully, even down to spin = 2! 5 To describe [σ ] 0 , we must perform a nontrivial sum over the twist family [σσ] 0 in the σσ → OPE channel. The result is another beautiful fit that works down to spin = 2. In this way, we find analytical approximations for dimensions and OPE coefficients of [σσ] 0 and [σ ] 0 in terms of the data {∆ σ , ∆ , f σσ , f , c T }.
Describing [σσ] 1 and [ ] 0 requires a novel approach because the two families exhibit nontrivial mixing. (For example the OPE coefficient f [σσ] 1 is larger than f [ ] 0 for spins ≤ 26.) In section 7, using our solution of the asymptotic lightcone bootstrap, we show how to define a "twist Hamiltonian" H(h = ∆+ 2 ) whose diagonalization correctly describes this mixing, and matches the numerics well for ≥ 4. In particular, diagonalizing H(h) leads to O(1) anomalous dimensions and variations in OPE coefficients, despite the fact that we have truncated the asymptotic expansion for H(h) to only a few terms. Our tentative conclusion is that by using the appropriate twist Hamiltonian, the large-spin expansion can in practice be extended down to relatively small spins for all double-twist operators in the 3d Ising CFT (and perhaps other theories as well).
In section 8, we ask what the asymptotic large spin expansion is missing. Part of the four-point function is invisible to this expansion, to all orders in 1/ . Demanding that this part be crossing-symmetric gives additional nontrivial constraints on the CFT data. Using our analytical approximations from section 6, we briefly explore some of these constraints. For example, we find conditions that approximately determine c T and f σσ in terms of ∆ σ , ∆ , f , using only the lightcone limit.
We discuss future directions in section 9. 4 The sum over k in (1.2) can be written in terms of 3 F 2 hypergeometric functions. 5 This was conjectured in [62,1]. 2 Numerics and the lightcone limit 2.1 A numerical picture of the 3d Ising spectrum Numerical bootstrap methods have become powerful enough to estimate several operator dimensions and OPE coefficients in the 3d Ising CFT. The strategy is as follows. Consider the four-point functions σσσσ , σσ , and where σ and are the lowest-dimension Z 2 -odd and Z 2 -even scalars in the 3d Ising CFT, respectively. Crossing symmetry and unitarity for these correlators forces the dimensions ∆ σ , ∆ and OPE coefficients f σσ , f to lie inside a tiny island given by [55] ∆ σ = 0.5181489 (10), f σσ = 1.0518537 (41), ∆ = 1.412625 (10), f = 1.532435 (19).
We can then ask: given that (∆ σ , ∆ , f σσ , f ) lie in this island, what other operators are needed for crossing symmetry? Although it is possible in principle to compute rigorous bounds on more operators, it is difficult in practice because we must scan over the dimensions and OPE coefficients of those additional operators.
Instead, we adopt the non-rigorous approach of [68], based on the extremal functional method [7,14,20]. Consider N derivatives of the crossing equation around z = z = 1 2 , which we write as F N = 0, where F N is an N -dimensional vector depending on the CFT data. We assume that OPE coefficients are real and operator dimensions are consistent with unitarity bounds [69]. By the argument of [3], there is an allowed region A N in the space of CFT data such that any point outside A N is inconsistent with F N = 0. 6 For every point p on the boundary of A N , there is a unique "partial spectrum" S N (p): a finite list of operator dimensions and OPE coefficients that solve F N = 0. The number of operators in S N (p) grows linearly with N . 7 If p lies on the boundary of the Ising island and N is large, we might expect that S N (p) is a reasonable approximation to the actual spectrum of the theory. However, it is not obvious how to assign error bars to S N (p). Firstly, the actual theory lies somewhere in the interior of the island, not on the boundary. It is important that the island is small enough that points on the interior are close to points on the boundary. Secondly, S N (p) depends on p, and there is no canonical choice of p.
In [68], we propose the following trick. We sample several different points p on the boundary of the island, and compute S N (p) for each one. As we increase N and vary p, some of the operators in S N (p) jump around, while others remain relatively stable. If an operator remains stable, we can guess that it is truly required by crossing symmetry.
In [68], we used this strategy to estimate the dimensions and OPE coefficients of a few low-dimension operators in the 3d Ising CFT. In figures 1 and 2, we show a more complete computation, giving about a hundred stable operators. To produce figures 1 and 2, we computed 60 different spectra by varying (∆ σ , ∆ , f σσ , f ) and minimizing c T . (We give more details in appendix A.1.) We then superimposed these 60 spectra, and grouped together operators with dimensions closer than 0.03. Each circle represents a group, and the size of the circle is proportional to the number of operators in that group. Thus, large circles correspond to stable operators and small circles correspond to unstable operators. We list the dimensions and OPE coefficients of the stable operators in appendix A.3. Most of the stable operators also appear in figures 7, 9, 12, 13, 14, 17, 18, and 19, where we compare to analytics.  Estimates of Z 2 -even operators in the 3d Ising model. Larger circles represent "stable" operators whose dimensions and OPE coefficients have small errors in our computation. We plot the twist ∆ − versus spin . The grey dashed lines are τ = 2∆ σ + 2n and τ = 2∆ + 2n for nonnegative integer n.

Effectiveness of the large spin expansion
Let us make some comments about these results. Firstly, most of the stable operators fall into families with increasing spin and nearly constant twist τ = ∆ − . We immediately recognize these as double-twist operators -specifically the families [σσ] 0 , [σσ] 1 , [ ] 0 in figure 1, and [σ ] 0 in figure 2. (There are also vague hints of [σ ] 1 .) The fact that these families are stable implies that they play a crucial role in the numerical bootstrap for the 3d Ising CFT. 8 8 Note that even though our numerical calculation uses an expansion of the crossing equation around the Euclidean point z = z = 1 2 , the results are sensitive to the Lorentzian physics of the lightcone limit. The prevailing lore was that, since the conformal block expansion converges exponentially in ∆ in the Euclidean One can compute anomalous dimensions of double-twist operators in a large-expansion using the crossing equation [59-62, 1, 63-67]. The authors of [1] observed that the largeexpansion appears to be asymptotic, but they conjectured that the anomalous dimensions of [σσ] 0 should be well-described by the first few terms in this expansion, coming from the operators and T µν appearing in the σ × σ OPE. The expansion is most naturally organized in terms of the "conformal spin" J defined by regime [70], numerical bootstrap methods should only be sensitive to low-dimension operators. Evidently this is incorrect because certain derivatives probe physics outside the Euclidean regime. Some hints that the numerical bootstrap probes the lightcone limit were given in [71], where an exact extremal functional was constructed that involves the lightcone limit of conformal blocks.
One finds 9 Here, d = 3 is the spacetime dimension and c free is c T for the free boson [72]. We will rederive (2.3) and find its all-orders generalization in section 5. Plugging in (2.1) and the value c T c free ≈ 0.946534(11) (2.5) computed in [20], we find that this prediction fits the numerics beautifully, even at small (figure 7)! This is surprising because the arguments leading to (2.3) only fix anomalous dimensions at asymptotically large . Rigorously speaking, they say nothing about any finite value of .
Nevertheless, inspired by this result, we might try to match the dimensions of [ ] 0 and [σσ] 1 to the leading terms in their large-spin expansions. Unfortunately, the naive analytic predictions disagree wildly with the data. To fit [ ] 0 and [σσ] 1 , we will need a more sophisticated understanding of the large-spin expansion, which we develop over the course of this work.
A clue about what's going on is the fact that the twists of [ ] 0 and [σσ] 1 move away from each other at small . This is reminiscent of the behavior of the eigenvalues of as → 0. If |τ 1 −τ 2 | is small, the eigenvalues repel more at small . Furthermore, the smalleigenvectors become nontrivial admixtures of the large-eigenvectors. In the 3d Ising CFT, it turns out that τ 1 = 2∆ σ + 2 ≈ 3.04 is numerically close to τ 2 = 2∆ ≈ 2.83. This suggests that the repulsion between [ ] 0 and [σσ] 1 is due to large mixing between these families. We will make this notion more precise in section 7 and compute the twists of [ ] 0 and [σσ] 1 in section 7.5. The off-diagonal terms will come from the σ operator in the σ × OPE and behave like −∆σ .

Double-twist operators
Let us review the argument from [59,60] for the existence of double-twist operators. The crossing symmetry equation for a four-point function of scalar operators φφφφ is Here, O runs over primary operators in the φ × φ OPE and ∆, are the dimension and spin of O. The functions g ∆, (z, z) are conformal blocks for the d-dimensional conformal group SO(d, 2).
The lightcone limit is given by z 1 − z 1. 10 Let us replace z → 1 − z so that we have 2) and the lightcone limit becomes z z 1. (We have used g ∆, (1 − z, z) = g ∆, (z, 1 − z).) In this limit, the left-hand side is dominated by the unit operator, z −∆ φ (1 + O(z)). On the right-hand side, no single term dominates the small z limit. However, because we also have small z, we can replace each conformal block by its expansion in small z [73,74], The function k 2h (x) is a conformal block for the 1-dimensional conformal group SL(2, R).

Our equation becomes
The left-hand side of (3.6) has a power-law singularity at small z. However, each individual term on the right-hand side has a logarithmic singularity at small z, 12 A power singularity can only come from the sum over an infinite number of operators on the right-hand side with h → ∞. Also, these operators must have h → ∆ φ as h → ∞ to match fact that z −∆ φ on the left-hand side is independent of z. These are the double-twist operators [φφ] 0 .
One can determine the asymptotic growth of the OPE coefficients f φφ[φφ] 0 by demanding that they reproduce the singularity z −∆ φ . The leading growth is (3.8) The sum in (3.6) is dominated by the regime 2h √ z ∼ 1, 13 where the SL(2, R) block becomes where K 0 (x) is a modified Bessel function. We can then approximate the sum over [φφ] 0 as an integral, which reproduces the required singularity (The factor of 1 2 is because only even spin operators appear in [φφ] 0 .) Matching z −∆ φ only determines the asymptotic density of OPE coefficients f 2 φφ[φφ] 0 at large h. The density (3.8) could be distributed evenly, with one operator per spin, or with one operator every other spin, or in many different ways. We will not see evidence of this freedom when we compare to numerics. The OPE coefficients will always be distributed in the simplest way consistent with the large-spin expansion.
We can determine the anomalous dimensions of double-twist operators by matching additional terms on the left-hand side of (3.6). Let O 0 be the smallest-twist operator in the φ × φ OPE that is not the unit operator (often O 0 = T µν ). Including the contribution of O 0 at small z on the left-hand side of (3.6), we have where we have used (3.7), this time on the left-hand side of the crossing equation. To match the log z term, we can take (3.13) 13 The fact that 2h √ z ∼ 1 is the appropriate regime was shown in [59]. It also follows from the physical arguments of [75]. In the physical picture of [75], this diagram shows the exchange of virtual O 0 -particles between φ-particles separated by a distance log over time log z.
Dividing (3.13) by (3.8) gives the leading large-h expansion of the anomalous dimension δ(h) ∝ h −2h 0 , agreeing with the leading term in (2.3). 14 Again, only the asymptotic density of the combination f 2 φφ[φφ] 0 δ is determined by this computation. An interesting feature of this argument (not realized in [59,60], but pointed out in [61]) is that it most naturally determines a function h(h) instead of τ ( ). We obtain actual operator dimensions by demanding that the spin be an even integer, (3.14) Thinking in terms of h(h) will be even more important when we compute higher-order corrections to (3.13).
It is often useful to draw the contribution of O 0 as a "large-spin diagram" like figure 3 (see, e.g. [76]). Such diagrams are particularly natural in the language of [75], where largespin operators become widely separated particles in a massive two-dimensional effective theory. Figure 3 represents a Yukawa potential between φ-particles induced by exchange of a virtual massive O 0 -particle. The distance between φ's (the width of the figure) is given by χ = log , and the mass of O 0 is the twist m = τ 0 . The Yukawa potential has the form e −mχ = −τ 0 , in agreement with the large-behavior of δ(h). We can also think of figure 3 as having height − log z, so that integration over the vertical position of the O 0 exchange gives a factor − log z, matching the − log z term in the conformal block of O 0 (3.12).
14 δ = γ 2 is half of what is usually called the anomalous dimension.

What about log n z?
Above, we matched the log z terms on the left-hand side of (3.12) to anomalous dimensions on the right-hand side. However, the expansion of z δ contains higher-order terms in log z: What do they map to under crossing? Using (3.8), (3.9), and (3.13), the log n z terms become The z-dependence of (3.16) is what one would expect from an operator of weight nh 0 . Such operators exist: they are the multi-twist operators [O 0 . . . O 0 ] 0 . The log n z behavior is not present in any individual conformal block -instead it must come from a sum over all the operators in the family [O 0 . . . O 0 ] 0 . We will see examples of log 2 z coming from a sum over double-twist operators in sections 6 and 7. We prove that double-twist operators always account for the correct log 2 z terms (i.e. that exponentiation of δ log z works automatically to second order) in appendix C.
We could have immediately guessed this result using large-spin diagrams. Exponentiating the Yukawa potential in figure 3 gives a sum of "ladder diagrams" like figure 4. Reading these diagrams from left-to-right, they look like an exchange of multi-twist operators. If we interpret the figure as having height log z, then integration over the vertical positions of the exchanges gives log n z. In practice, n − 1 "integrations" are achieved by summing over different distributions of derivatives among the operators ∂ · · · ∂O 0 · · · ∂ · · · ∂O 0 (i.e. by summing over all members of the twist family [O 0 · · · O 0 ]), while one integration is encoded in the log z factor in each individual conformal block. This makes it clear why we must sum over all multi-twist operators [O 0 . . . O 0 ] in one channel to recover exponentiation in the other channel.
In this way, crossing symmetry forces multi-twist operators [O 1 . . . O n ] to appear in the conformal block expansion whenever O 1 , . . . , O n do individually. In particular, this implies that multi-twist operators built from the stress tensor and other low-spin operators should appear in the σ × σ OPE in the 3d Ising model. In figure 1, we see some evidence of operators with twist near 2, which would correspond to [T T ] 0 . However, none of them are numerically stable. This is likely because the anomalous dimension δ [σσ] 0 is small (of order 10 −2 ), so higher terms in the expansion of z δ (3.15) are highly suppressed. To get a better picture of these operators, one must study mixed correlators involving σ and T µν together, or perhaps higher-point correlators like σσσσσσ . We return to this point in section 7.1.

The Casimir trick
The derivation of (3.13) makes sense when h 0 < ∆ φ , so that the sum diverges faster than any individual term (log z) at small z. When this happens, the sum must be dominated by large h and can be approximated by an integral. 15 However, [61] argued that the large-spin expansion can be extended to include contributions from operators with h 0 > ∆ φ . For example, there is a calculable correction to δ [σσ] 0 in the 3d Ising CFT coming from , which has h ≈ 0.7 > ∆ σ ≈ 0.52.
To see why, suppose h 0 > ∆ φ . Since each term k 2h (1 − z) is more singular than z h 0 −∆ φ , we cannot conclude that the sum is dominated by large h. However, k 2h (1 − z) obeys a Casimir differential equation with eigenvalue h(h − 1), By repeatedly acting with the Casimir operator D on a power z a , we can make it arbitrarily singular, 16 D n z a = (a − n + 1) 2 n z a−n (1 + O(z)). (3.19) Acting n times on (3.17), we obtain Taking n big, the right-hand side is now dominated by large h when z is small, and we can proceed as before. The resulting correction to f 2 φφ[φφ] 0 δ is again given by (3.7).

Higher-order corrections
By including 1/h-corrections in the approximation (3.9), one can compute higher-order corrections to the OPE coefficients f φφ[φφ] 0 (h) and anomalous dimensions δ(h). After applying the Casimir operator enough times, each term in the 1/h-expansion contributes to a singularity at small z, and can thus be calculated by approximating the sum over h as an integral. This gives expansions of the form The authors of [1] showed how to use the Casimir trick to compute the above coefficients.
(Actually, they organize their expansion in terms of the Casimir eigenvalue J 2 = h(h − 1), as in equations (2.2) and (2.3).) In section 5, we will write down an all-orders solution for (3.21).
We have written "∼" to indicate that both sides have the same asymptotic expansion at large h. The arguments above only fix the asymptotic expansion of f φφ[φφ] 0 (h) 2 and δ(h) because it is always possible to throw away a finite number of blocks k 2h (1 − z) and still match the power z h O −∆ φ on the other side of the crossing equation. We can only fix the behavior of f φφ[φφ] 0 (h) 2 and δ(h) for h larger than some h 0 , where h 0 might grow as we include more terms in (3.21). Thus, (3.21) should be interpreted as asymptotic series.
The behavior of f φφ[φφ] 0 (h) 2 and δ(h) at finite h is still important -it contributes to "Casimir-regular" terms defined in the following section.

Sums of SL(2, R) blocks
Our main tool will be a better understanding of infinite sums of SL(2, R) blocks, where h is an increasing series of weights that asymptotes to integer spacing, and p(h) are coefficients that grow no faster than 2 −2h h const. as h → ∞. We start from a simple example, Mean Field Theory (MFT) in 1-dimension, and then generalize it in several ways.

Casimir-singular vs. Casimir-regular terms
Sums of SL(2, R) blocks have two parts that play different roles in the bootstrap. As discussed in in section 3.2, we can make a power z a arbitrarily singular by repeatedly applying the Casimir operator D, D n z a = (a − n + 1) 2 n z a−n (1 + O(z)).

(4.2)
We say that z a for generic a is Casimir-singular. An exception occurs when a is a nonnegative integer, since then (a − n + 1) 2 n vanishes for n ≥ a + 1. In fact, terms of the form z n , z n log z n ∈ Z ≥0 (4.3) do not become arbitrarily singular when we repeatedly apply D. We call such terms Casimirregular.
The lesson of section 3.3 is that Casimir-singular terms can be matched unambiguously to an asymptotic expansion in large h. Furthermore, to compute coefficients in this expansion, we can think of the sum over h as an integral. By contrast, Casimir-regular terms are not determined by a large-h expansion. This is consistent with the fact that a single SL(2, R) block k 2h (1 − z) is Casimir-regular, since it is an eigenvector of the Casimir operator. (We can also see that it is Casimir-regular by noting that its z-expansion (3.7) is a sum of terms of the form (4.3).) For example, suppose Moving the first term on the left-hand side to the right-hand side, we have The Casimir-regular part of the right-hand side has changed, but the large-h expansion of p(h) obviously hasn't.
It will often be useful to work modulo Casimir-regular terms. When we do so, we denote Casimir-regular terms by [. . . ] z .

Matching a power-law singularity
Casimir-singular terms match to a unique asymptotic expansion for coefficients of SL(2, R) blocks at large h. We can find the right expansion by looking at an example. Consider the conformal block expansion of (4.6) Replacing z → 1 − z and writing ∆ = −a, this can be written where .
Many formulae will be much simpler in the variable y (and y, defined similarly) instead of z. Note that y a is Casimir-singular for generic a, while y n and y n log y for nonnegative integer n are Casimir-regular. We will denote Casimir-regular terms by [. . . ] y . The crossing transformation z → 1 − z maps y → 1/y.
Casimir-singular terms can only come from an infinite sum of blocks, and they are sensitive only to the asymptotic density of OPE coefficients. Thus, if we change the weights h entering (4.7), while preserving the same asymptotic density, only the Casimir-regular terms should change. For example, changing −a + → h 0 + , we expect (4.10) The Casimir-singular term y a is independent of h 0 , but the Casimir-regular terms [. . . ] y depend on h 0 . As a sanity check, (4.10) is certainly true when h 0 = −a + n for nonnegative integer n, since we get it by moving the first n terms of (4.7) to the right-hand side.
The coefficients S a (h) will be our building blocks for solving the asymptotic lightcone bootstrap. They encode the all-orders large-h expansion needed to match powers y a . By taking linear combinations, we can match any Casimir-singular term we want. For example, to match an SL(2, R) block k 2h (z) in the crossed channel, we can take a linear combination of S h +m (h) which can be resummed into a 4 F 3 hypergeometric function.
Casimir-regular terms depend on the detailed structure of the weights being summed over. We can determine the Casimir-regular terms in (4.10) as follows. Let us expand k 2h (1 − z) in small y (equivalently small z) inside the sum, Here, we have introduced .
Naively, we might try to switch the order of summation in (4.12), (4.14) However, this cannot be correct. If the result converged, it would be Casimir-regular, a contradiction. Indeed, the summand grows like h 2k−2a−1 , so the terms with k > a diverge. However, let us analytically continue from the region a > k for each term. After some gymnastics, 17 we find We claim that (4.17) gives the correct coefficients for the Casimir-regular terms in (4.10). That is, we have the remarkable exact identity (equation (1.2) from the introduction) One can verify that (4.18) is consistent with the fact that shifting h 0 → h 0 + 1 changes both sides by −S a (h 0 )k 2h 0 (1 − z). We have also extensively checked (4.18) numerically. 18 We slightly generalize (4.18) in equation (4.47). The special case of this formula where h 0 = 0 was proven recently in [78], using hypergeometric function identities from [79].
The key feature of (4.18) is that it expresses an integer-spaced family of conformal blocks in one channel as an expansion in the other channel. Since families of nearly integer-spaced operators are ubiquitous, we can use (4.18) as a building block for understanding crossing symmetry in general.

General coefficients
Consider a sum of SL(2, R) blocks with general coefficients p(h) and integer-spaced weights, 17 We obtained (4.17) in the following shameful way. When b = −1, we can use . (4.16) When b = −k − 1 with k a positive integer, T −k−1 (h) is a polynomial in h and we can write T a (h)T −k−1 (h) in terms of linear combinations of terms of the form Γ( +α) Γ( +β) and use (4.16). We did this for several positive integer k's, guessed an answer for general k, analytically it continued away from integer k, and then checked the result numerically. 18 We expect (4.18) can be derived using Sturm-Liouville theory for SL(2, R) blocks [77].
If p(h) has the same large-h behavior as a sum of S a (h)'s, the structure of (4.19) will be similar to (4.18). To determine the Casimir-singular terms, we match asymptotic expansions, where A is some discrete (possibly infinite) set of values depending on the function p(h), and c a are constants. We then have To compute the Casimir-regular terms, we expand k 2h (1 − z) inside the sum and then naively switch the order of summation, Again, the sums in parentheses are divergent for sufficiently large k. However, we can regulate them by adding and subtracting linear combinations of the known answer (4.18) until the sums become convergent. This gives If we choose K ≥ k, then the sum over h in (4.24) will converge. In fact, the larger we take K, the more quickly the sums converge (since the quantity in parentheses falls off more quickly with h). Note that α k is analytic in k, so we can evaluate its derivative β k .

Non-integer spacing and reparameterization invariance
We often encounter sums over SL(2, R) blocks k 2h (1 − z) where the weights h are not integer-spaced. The Casimir-singular terms depend only on the asymptotic density of OPE coefficients. Thus, for a sequence h that depends sufficiently nicely on , we can compensate for uneven spacing by inserting a factor of ∂h ∂ , giving the same Casimir-singular part as an integer-spaced sum: A way to understand (4.25) is that Casimir-singular terms come from asymptotically large h, where the sum can be treated as an integral. We are then free to redefine the integration variable and include a Jacobian ∂h ∂ . We call this freedom "reparameterization invariance." Let us prove (4.25) for an important class of h . Suppose h is defined implicitly by (4.26) where δ(h) is an analytic function that behaves like a sum of powers h −b as h → ∞. We have Working modulo Casimir-regular terms, we may restrict the sum (4.25) to ≥ L for some large L so that δ(h) is small. Expanding (4.25) in small δ, we find the following identity: (4.28) (One way to motivate why an identity like (4.28) should exist is to pretend the sum over is an integral and consider an infinitesimal change of variables in the integral.) Now summing over , the terms in parentheses are integer-spaced sums of the type in section 4.3. They give Casimir-singular contributions that are independent of h 0 . Thus, only k = 0 contributes in (4.28), modulo Casimir-regular terms. This proves (4.25).
Another way to understand (4.25) is as follows. The non-integer-spaced sum can be written as a contour integral The Casimir-singular part must come from the region of the integral h → ±i∞, since any sum of blocks with bounded h is Casimir-regular. However, in this region the δ-dependent factor in the integrand approaches a δ-independent constant exponentially quickly (assuming δ(h) grows slower than h as h → ±i∞): Thus, the Casimir-singular part is δ-independent and can be obtained by replacing tan(π(h− h 0 − δ(h))) → tan(πh). 19 Sums over general weights h with general coefficients p(h) can be computed using the same strategy as in section (4.3). We obtain Casimir-singular terms from the asymptotic expansion of p(h). We determine Casimir-regular terms by expanding k 2h (1 − z) inside the sum, naively reversing the order of summation, and regulating the resulting sums over h. We give more details in appendix B.

Alternating sums and even integer spacing
We will also encounter sums of SL(2, R) blocks with insertions of (−1) . To understand these, consider the conformal block expansion of Substituting z → 1 − z and ∆ → −a, this can be written Note that (1 + y) a is Casimir-regular. Using the logic of the preceding sections, we conclude that general sums with (−1) insertions are Casimir-regular, where h = h is any sequence of the form discussed in section (4.4).
Let us describe how to compute the Casimir-regular terms in alternating sums. For simplicity, consider the case of integer-spaced weights and general coefficients p(h), The strategy is the same as before: we expand k 2h (1 − z) at small y, switch the order of summation, and regulate the resulting sums by adding and subtracting known answers. We find Again, c a are defined by matching asymptotic expansions analytically continued in a from the region where the sum converges.
We have not found a simple closed-form expression for A − a,b (h 0 ) for general a, b. However, we can evaluate it to arbitrary accuracy as follows. Using similar tricks to before, we can compute the case b = −1: This can be used to regularize the sum for general where K > −a − b − 5/2 is taken large enough that the sum over h converges. The larger we take K, the faster the sum converges. When a or b is a negative integer, the expansion (4.39) truncates and becomes an equality, and we can omit the second line in (4.40).
We will also need to evaluate sums with even-integer spacing. These are an average of alternating and non-alternating sums, Similarly, we define

Mixed blocks
Correlation functions of operators φ 1 φ 2 φ 3 φ 4 with different scaling dimensions can be expanded in SL(2, R) blocks of the form We include the unconventional factor (1 − z) −r because it simplifies several formulae later on. It also ensures that k r,s 2h (z) is symmetric in r and s, by elementary hypergeometric function identities. Casimir-regular terms for the mixed block (4.43) are of the form y n−r and y n−s for nonnegative integer n.
The mixed block analog of S a (h) is These coefficients satisfy the 1-dimensional MFT equation and its generalization in the spirit of the previous sections 20 Using (4.17), we also find a generalized version of (4.18) giving the explicit Casimirregular terms in an integer-spaced sum of mixed blocks Here, ". . . " represents other operators that are unimportant for this computation. Note that the y variables make the unit operator block very simple. For other operators, expanding in y instead of z is equivalent to shuffling around contributions of descendants.
The h i are weights of primary and descendant operators in the φ × φ OPE. We can match the left-hand side by choosing Here, "∼" means the two sides have the same large-h expansion. We include factors of 2 in (5.4) because the family [φφ] 0 only contains even spin operators. Dividing, we find Once we know δ [φφ] 0 (h), we can obtain the OPE coefficients f φφ[φφ] 0 from (5.3). Expanding in large h gives a series with terms of the form 1/h 2(h i 1 +···+h i k )+n .
In (5.5), we can see explicitly why the large-spin expansion for δ [φφ] 0 (h) is naturally organized in terms of the Casimir eigenvalue J 2 = h(h − 1) as discussed in [61]. The reason is that ratios of S a (h) are also ratios of We have suppressed an important subtlety in (5.1). The OPE φ × φ contains an infinite number of operators with bounded h (for example, the families [φφ] n ) themselves. Thus the sum on the left-hand side, may not converge. For simplicity, suppose all the h i = h are the same. The correct procedure is to perform the sum over i first, before expanding in y, using the methods of section 4.2. This leads to y h a c a y a + y h (A log y + B + O(y)), (5.8) where A and B are regularized versions of the sums over A i and B i . The y a terms are Casimir-singular in y, and will be cancelled by other operators on the right-hand side of (5.1). The remaining y-Casimir-regular (but still y-Casimir-singular) terms y h A log y and y h B contribute to anomalous dimensions and OPE coefficients of [φφ] 0 , respectively. The y-Casimir-singular terms in (5.8) can also include log n y contributions related to higher-order exponentiation of anomalous dimensions, and discussed in section 3.1.1. We will see several examples in section 6.
Thus, the techniques of section 4.2 for summing SL(2, R) blocks have two roles to play. Firstly, they let us match Casimir-singular terms in one channel to h-dependence in the other channel. Secondly, they let us resum operators whose twists have accumulation points.
Naively this leads to an impasse: we must resum [φφ] 0 before finding how it contributes to its own anomalous dimensions δ [φφ] 0 . However, it turns out that [φφ] 0 contributes to its own anomalous dimensions only at order δ 2 [φφ] 0 and higher. (This is related to the fact that Mean Field Theory has no anomalous dimensions.) Thus, both the resummation and the matching to h-dependence will be possible. We will see this explicitly in section 6.1.2.

Why asymptotic?
We have been careful to write "∼" instead of "=" because the relations (5.4) are not necessarily equalities. In fact, taken literally, the expressions on the right-hand side may not even converge to functions of h. Instead, they represent equivalence classes of functions with the same asymptotic expansions at large h. For example, both sides of formally have the same large-h expansion, but they are different. In fact, the sum on the right diverges. We must interpret (5.9) in terms of large-h equivalence classes.
The asymptotic nature of the large-h expansion for double-twist operators makes mathematical and physical sense. Mathematically, a given Casimir-singular term only determines an asymptotic density of coefficients on the other side of the crossing equation. Any change in the density at finite h contributes to Casimir-regular terms. Thus, we cannot fix the actual function of h without simultaneously considering all Casimir-regular terms.
Physically, it is ambiguous which twist family (if any) we should assign a given operator to. For instance, should we assign T µν to the family [σσ] 0 , or should the family should start at spin-4 or higher? Twist families only make sense as infinite collections of operators with unbounded spin. We shouldn't necessarily expect to write analytic expressions that interpolate between their OPE coefficients and dimensions at finite . On the other hand, we might expect a convergent large-h expansion for an object that packages together all operators in the theory, and does not try to distinguish them into twist families.
When our theory has extra structure, twist families may become well-defined even at finite spin. For example, in a large-N expansion, we have a well-defined classification of operators into single-trace, double-trace, etc.. Consequently, large-h equivalence classes in large-N theories should have distinguished representatives. See, for example, in [80]. Similar remarks hold in weakly-coupled theories.

General double-twist families
Let us be more explicit and derive all-orders expansions for OPE coefficients and anomalous dimensions of double twist families [φ i φ j ] n for all n ≥ 0. For generality, we study mixed four-point functions φ 1 φ 2 φ 3 φ 4 of scalars with possibly different external dimensions.
We use a slightly unconventional definition for SO(d, 2) blocks, (z, z) are the mixed scalar blocks of [81] with coefficient c = 1. 21 Using identities from [81], one can show that our G r,s h,h (z, z) is symmetric under r ↔ s. The extra factors ((1 − z)(1 − z)) −r = v −r simplify the crossing equations in the y, y variables and make the symmetry between r and s manifest. For brevity, we omit r, s when they are zero.

Sums over n and
The coefficients S r,s a (h) give a simple result when summed over a single family of SL(2, R) blocks. However, in d-dimensions, double-twist operators come in doubly-infinite families, 21 Our blocks differ from those of [24] . 22 The ordering f 12O f 43O differs from the f 12O f 34O ordering in [24] because our blocks differ by (−1) times positive factors. We have reabsorbed this (−1) by using f 34O = (−1) O f 43O . A useful way to remember the correct sign is to note that φ 1 (0)φ 2 (z)|O|φ 2 (1)φ 1 (∞) is the norm of a state in radial quantization, where |O| is a projector onto the conformal multiplet of O. Thus, it should be positive, which implies that it should have coefficient f 2 12O in the conformal block expansion.
labeled both by and n such that h ≈ h 0 + n. The d-dimensional analog of S r,s a (h) will be coefficients C (n)r,s a (h 0 , h) that, when summed over both and n, produce a simple result, We can obtain the C (n)r,s a (h 0 , h) by expanding SO(d, 2) blocks in terms of SL(2, R) blocks and using what we know about the coefficients S r,s a (h). A simple example is in 2-dimensions, where SO(2, 2) blocks are just products of SL(2, R) blocks, 23 (for simplicity we take r = s = 0). Then we have In general, SO(d, 2) blocks have an expansion of the form 24 The coefficients A r,s n,j (h, h) can be determined, for example, by solving the SO(d, 2) Casimir equation order-by-order in y. Alternatively, we can obtain them from the decomposition of d-dimensional blocks into 2-dimensional blocks [82]. The first few coefficients are then the y h 0 terms on both sides of (5.13) will agree, by equation (4.46). We can then choose the n > 0 coefficients to cancel higher-order terms in y. This gives a recursion relation that determines all the higher C (n) 's.
As a cross-check, recall that d-dimensional MFT has conformal block expansion with coefficients given by [83] C MFT n, To be consistent with (5.13), we must have We have checked this explicitly for n = 0, 1, 2. Although C MFT n, (∆ 1 , ∆ 2 ) has a simple formula, we have not found a closed-form expression for C (n)r,s a (h 0 , h) in general dimensions.

Small y expansion of the left-hand side
On the left-hand side of the crossing equation, we should expand the blocks G h,h (z, 1 − z) in small y. As a starting point, the SL(2, R) blocks have an expansion In the special case r = s, this becomes . (5.27)

Matching the two sides
Using (5.25), the left-hand side of the crossing equation (5.12) is Let us assume that the terms y k+h 1 +h 2 match the families [φ 1 φ 2 ] n with n ≤ k on the righthand side, while y k+h 3 +h 4 match [φ 3 φ 4 ] n with n ≤ k. (We return to this assumption in section 7.) As before, define λ ij[kl]n by Using (5.13) and working order-by-order in y, we find The sum O∈2×3,m≥0 runs over operators O in the φ 2 × φ 3 OPE and their descendants organized by weights under SL(2, R) L . The prime indicates that we must regularize the sum, as discussed above and demonstrated in sections 6 and 7.
We are free to add contributions that do not change the Casimir-singular part of the sum over blocks. As we learned in section 4.5, sums with a (−1) insertion are Casimir-regular. Thus, we can safely add the two contributions, and this single formula produces the correct Casimir-singular terms in both cases. 26 The two terms in (5.33) are illustrated in figure 5.
In the special case h 1 + h 2 = h 3 + h 4 , (5.28) develops log y-dependence (because P r,s m,k has a pole at r = s), and we instead find a formula for products of OPE coefficients and anomalous dimensions, λ 12 [12]n (h)λ 43 [12]n (h)δ [12]n (h) + λ 12 [34]n (h)λ 43 [34]n (h)δ [34] (5.37) More explicitly, they are given by Specializing further, we will need the case where the pairs of operators φ 1,2 and φ 3,4 are actually the same. Since now only a single family [12] n reproduces y k+h 1 +h 2 and y k+h 1 +h 2 log y in (5.28), we must drop the [34] n terms in (5.35) before setting 12 = 43. This gives The identity operator is the leading contribution to (5.41). Its coefficients are those of Mean Field Theory, analytically continued to Finally, when all the operators are equal, (5.40) and (5.41) become We will often replace 1 + (−1) → 2 and simply remember that only even-spin operators appear in the OPE φ 1 × φ 1 .

Checks
Knowing CFT data up to weight h max unambigiously determines the large-h corrections up to order h −2hmax , or equivalently J −τmax . To get this information, we could alternatively use the technology of [1]. It is straightforward to check that the first few J −τ O corrections to anomalous dimensions agree: where δ( , h) has a large-h expansion that includes powers of h and factors of (−1) , The proof in section 4.4 then works, provided we replace where in the derivative ∂δ ∂h we treat (−1) as constant.

Application to the 3d Ising CFT
Let us now apply these results to the 3d Ising CFT. We would like to see how well the truncated large-h expansion describes the spectrum at finite h. The more operators we can describe precisely, the better the prospects for hybrid analytical/numerical approaches like those discussed in section 9.1. We will find that a few terms in the expansion match numerics surprisingly well, even down to relatively small spins.
We will organize our expansions in terms of S a (h)'s. This simplifies several computations (in particular it makes it simpler to compute Casimir-regular terms). However, one could just as well use powers of the SL(2, R) Casimir J 2 , as in [61,62,1,[65][66][67]. A sum of S a (h)'s is a partial resummation of a series in J 2 .
We will work our way upwards in twist, first understanding [σσ] 0 in section 6.1, then [σ ] 0 in section 6.2, and finally [σσ] 1 and [ ] 0 in section 7.5. Because 2h σ is so small, the family [σσ] 0 is particularly important. Its contribution to other large-h expansions is competitive with those T µν and . Thus, we will use our formulae for OPE coefficients and dimensions of [σσ] 0 in several subsequent computations. We expect this approach should also work well for the O(N ) models. It is an interesting question whether it works in a general CFT.

Differences from [1]
Let us comment briefly on the (inconsequential) differences between the above calculation and the series (2.3) computed in [1]. Firstly, we have not included descendants of , T , namely terms of the form W with m ≥ 1, whereas [1] included descendants at first order in z. This is because it doesn't make sense to include level-1 descendants of , T without also including the double-twist operators [ T ] 0 , [T T ] 0 , which contribute at the same order in the large-h expansion. Also, because we organize everything as a series in y instead of z, the contributions of descendants will differ somewhat (though the sum over all of them will be the same). In addition, we have partially resummed the J series into sums of S a (h)'s.
All these alternatives represent different choices of subleading terms in a series that we are truncating anyway. Fortunately, they turn out to be inconsequential at the truncation order and precision at which we are working. A plot of (2.3) looks essentially identical to figure 7. However, S a (h)'s will begin to differ from powers of J when a = h O + m − 2h σ is larger (i.e. for higher-twist primaries and descendants in the crossed-channel). This is because S a (h) has poles at h = a + 1, a, a − 1, . . . , whereas J −2a does not. (In Sturm-Liouville theory for SL(2, R) blocks [77], these poles come from the region near y ∼ 1, outside the validity of the small-y expansion. Thus, they are artifacts of our expansion in small-y in the crossed-channel.) These differences reflect the fact that we are comparing different truncations of an asymptotic expansion outside the regime of validity of those truncations.

Contributions of [σσ] 0 to itself
We should also include higher-spin members of the family [σσ] 0 in (6.1), (6.2). Their contributions for = 4, 6, . . . are small because (6.8) 27 An alternative approach to computing corrections to anomalous dimensions from an infinite family of operators is given in [66,67].  .2)). The leading nonzero term has m = 2, corresponding to two vertical lines, or a "box diagram." The terms with k = 0 contribute to λ σσ[σσ] 0 and δ [σσ] 0 as follows where "above" represents terms already present in (6.1) and (6.2), and β even k = ∂ ∂k α even k . The sums start at m = 2 because S a (h) has a second-order zero at a = 0. (Equivalently, the terms proportional to log m y are Casimir-regular in the other channel when m = 0, 1.) We illustrate the contributions (6.9) in figure 8. Equation (6.9) might look complicated because λ σσ[σσ] 0 and δ [σσ] 0 are defined in terms of themselves. However, δ [σσ] 0 is small, so (6.9) is easily solved by iteration starting with the approximations (6.1), (6.2). Above, we show the result from plugging in (6.1), (6.2) and setting 0 = 4. The corrections in (6.9) are so small that we mostly omit them in what follows. By contrast, similar corrections for [σ ] 0 begin at m = 1, and for [ ] 0 they begin at m = 0. In these cases, one must sum the whole family [σσ] 0 to get accurate results.
We compare analytics and numerics for f σσ[σσ] 0 in figure 9. There is an interesting wrinkle in interpreting the numerics. Although the numerical spectra include operators O with twists τ [σσ] 0 , they also sometimes include spurious higher-spin currents J at the unitarity bound with small but nonzero OPE coefficients. Because τ [σσ] 0 is close to the unitarity bound, these spurious operators can "fake" the contribution of O in the conformal block expansion. 28 The J are artifacts of the extremal functional method. They should disappear at sufficiently high derivatives, but working at higher derivatives is not currently feasible. Instead, we remove them by hand and add their OPE coefficients to the correct operators O . In other words, we use (f 2 σσO +f 2 σσJ ) 1/2 as our numerical prediction for f σσ[σσ] 0 . Indeed, the numerical errors in in this modified quantity are smaller than the errors in f σσJ , and the results agree beautifully with the analytical prediction. We show numerical data both before and after the modification in figure 9. The leading contribution to the OPE coefficients λ [σσ] 0 comes from σ-exchange in the 28 Higher spin currents are disallowed in interacting CFTs [84][85][86][87]. We get a factor of (−1) in (6.12) and (6.13) because σ and switch places.
This agrees with numerics within 1% for all spins ≥ 2. In the next section, we compute additional corrections from the family [σ ] 0 and improve the agreement.

[σ ] 0
The leading correction to OPE coefficients and anomalous dimensions of [σ ] 0 comes from exchange of σ in the σ → σ channel (figure 10), To go further, we must include the contribution of the family [σσ] 0 in σσ → . Doing so will provide a nontrivial test of the tools we have developed.
Because we will discuss both channels simultaneously, let us write the crossing equation in a way that emphasizes the important terms: Our first goal is to compute the sum over [σσ] 0 on the right-hand side, = y-Casimir-singular + α(y) log y + β(y) + O(y).
The terms α(y) log y and β(y) have the correct form to contribute to anomalous dimensions and OPE coefficients of [σ ] 0 on the left-hand side of (6.14). However, the Casimir-singular terms do not, and must be cancelled in some other way. We work through an explicit example in section 6.2.1.
Before performing the sum over [σσ] 0 , let us understand what part we will need. Consider O = [σσ] 0, on the right-hand side of (6.14), and suppose is large so that δ [σσ] 0, = h O − 2h σ is small. The y-dependence of the O-block maps to the following h-dependence of λ σ [σ ] 0 on the left-hand side: The k = 0 term vanishes because S r,s r+a has a simple zero at a = 0. The first nontrivial correction has k = 1 ( figure 11). Thus, the leading correction to λ σ [σ ] 0 and δ [σ ] 0 in the sum over [σσ] 0 comes from expanding to first-order in the anomalous dimension δ [σσ] 0 : Here, ". . . " represents non-log y terms that do not contribute to λ σ [σ ] 0 and δ [σ ] 0 . We will treat T µν separately, so the family [σσ] 0 starts at h 0 = 2h σ + 0 with 0 = 4.
The quantities λ σσ[σσ] 0 , λ [σσ] 0 , and δ [σσ] 0 can be obtained from (6.1), (6.2), and (6.11). For simplicity, we approximate their product by the first two leading terms at large h, coming from the corrections to δ [σσ] 0 due to and T , This approximation has the correct asymptotics and also matches numerics within 1% for all ≥ 4. This is sufficient accuracy for our purposes, since we are already computing a small correction to [σ ] 0 .
For the O = T term, we have where The function lim a→0 Γ(−a) 2 S a (h) is finite, but when we insert it in a sum over blocks, both the Casimir-singular and Casimir-regular terms are naively infinite. However, 1/a 2 and 1/a poles cancel between them, leaving a finite result: where (6.25) Here, ψ (m) (z) ≡ d m+1 dz m+1 log Γ(z) is the polygamma function, ψ(z) = ψ (0) (z), and γ = −ψ(1) is the Euler-Mascheroni constant. (Even though we do not have a simple formula for β even 0 [S a ](h 0 ) in general, the limit B 0 (h 0 ) is computable in closed form and given by (6.25).)

Cancellation of Casimir-singular terms
Equation (6.23) again has the form anticipated in (6.15), where the log 2 y term in (6.23) is y-Casimir-singular. Combining (6.22) and (6.23), this term is in the σσ → channel.
We claimed earlier that the Casimir-singular terms in (6.15) should be canceled by other contributions, and it is instructive to see how this works explicitly. The expression (6.26) has the correct form to match the exchange of [σ ] in the σ → σ channel, where log 2 y comes from expanding y δ [σ ] 0 to second order in δ [σ ] 0 . We could have guessed this by reinterpreting figure 11 as the second order term in the exponentiation of figure 10 (in the bottom-to-top channel).
The important terms in δ 2 [σ ] 0 come from squaring the contribution of σ-exchange. From (6.12) and (6.13), we have Using (4.47), the relevant sum over blocks is (We have written the y −s part of the Casimir-regular terms because we will need them shortly.) Again, 1/a poles cancel between the Casimir-regular and Casimir-singular part, leaving a finite result. It follows that which exactly matches (6.26).
Thus, the other channel indeed cancels the Casimir-singular term in (6.23). This phenomenon, which has been explored previously in [1,65], is a special case of a more general result. The y-Casimir-singular part of the exchange of double-twist operators in one channel matches the y-Casimir-singular part of the exchange of double-twist operators in the other channel. Another way to say this is that box diagrams like figure 10 give the same Casimirsingular parts when interpreted from bottom-to-top or from left-to-right. 29 We prove this claim in appendix C. 30 The Casimir-regular term proportional to y hσ−h in (6.30) determines the leading correction to f [σσ] 0 coming from [σ ] 0 exchange. Including also level-one descendants of σ, which contribute at similar order in the 1/h-expansion to [σ ] 0 , we have where 0 = 2 is the lowest spin appearing in the [σ ] 0 family. As we show in figure 12, (6.31) agrees with numerics for all spins with accuracy ∼ 10 −3 .

Putting everything together
Combining the Casimir-regular terms from (6.20) and (6.23), we have From the above, we can read off the contributions to λ σ [σ ] 0 and δ [σ ] 0 from exchange of the family [σσ] 0 . Including also the corrections from exchange of and T µν , we have (6.34)

Comparison to numerics
We plot the twists τ [σ ] 0 = ∆ σ + ∆ + 2δ [σ ] 0 in figure 13 and OPE coefficients f σ [σ ] 0 in figure 14, comparing the formulae (6.33) and (6.34) to numerical results. In both cases, analytics matches numerics to high precision (∼ 10 −4 ) at large h, and moderate precision (< 10 −2 ) for all h. The agreement is particularly impressive because the corrections are large compared to Mean Field Theory, in contrast to the case of [σσ] 0 . Correctly summing the family [σσ] 0 is crucial for achieving this. One might guess that the asymptotic large-h expansion simply breaks down earlier for these operators -that it just doesn't work for 40. This turns out to be false. In this section, we give a procedure that extends the validity of the large-h expansion down to smaller values of h.
The key idea is to relax the assumption from section 5.3 that the double-twist operators [ij] n on one side of the crossing equation map only to terms of the form y h i +h j +k on the other side. Instead, we will compute a fully y-dependent asymptotic expansion in h and identify operators by diagonalizing an effective "twist Hamiltonian." Let Suppose that, using crossing symmetry, we can find the combination One way to extract the twist Hamiltonian is as follows. Given the elements M ijkl (y, h), we form the matrix where for brevity, we've defined ∂ ≡ ∂ ∂ log y . , which contribute at similar order in y, have small OPE coefficients in the σ × σ and × OPEs. This assumption is supported by numerics (which likely means that it follows from unitarity). However, we do not know how to derive it using the information in this work. Instead, we should enlarge our system of crossing equations to include additional external operators. For example, by studying the matrix . . · · · · · · · · · · · · · · · · · · . . . ∂ m+p M ijij · · · · · · ∂ m+q M ijkl · · · · · · . . . . . . . . . · · · · · · · · · · · · . . . . . . . . . . . . · · · · · · · · · . . . To summarize, we have where M 0 = M (y 0 , h) and M 0 = ∂M (y, h)| y=y 0 . The OPE coefficients Λ(h) can be obtained as follows. Let It must be the case that where Q 1 , Q 2 are orthogonal matrices. To determine the Q 1,2 , consider the combination The right-hand side has the form of a singular value decomposition (SVD), so Q 1 , Q 2 can be obtained by from an SVD of U −1 1 U 2 . Finally, we solve for Λ (and hence Λ) using either equation in (7.11). 32 Note that this procedure gives us λ ij[kl]n . To determine the actual OPE coefficients f ij[kl]n , we must multiply by Jacobian factors (5.29), which are different for each eigenvalue of the twist Hamiltonian h [kl]n .

Choice of external states
We can understand the twist-Hamiltonian prescription as follows. The four-point function is the amplitude for creating a state with φ i (x 1 )φ j (x 2 ) and annihilating it with φ k (x 3 )φ l (x 4 ). States created by pairs of local operators are not eigenstates of the twist-Hamiltonian H. Our task is to compute the change of basis between pair states φ i (x 1 )φ j (x 2 )|0 and H-eigenstates (the OPE coefficients f ij[ab]n ), and to find the eigenvalues h [ab]n . For this, we need matrix elements of y H between enough states to span the Hilbert space.
Although generically any eigenstate O will appear in the span of φ i (x 1 )φ j (x 2 )|0 (when global charges allow it), it should be easier to study O precisely if we use states that have large overlap with O. Specifically, we expect to get a better picture of the [φ i φ j ] n operators if we study matrix elements that include φ i (x 1 )φ j (x 2 )|0 . Similarly, one might learn about multi-twist operators [O 1 · · · O n ] by performing very high-precision studies of four-point functions. However, it may be more efficient to study matrix elements of O 1 (x 1 ) · · · O n (x n )|0 , i.e. to study higher-point correlators.

Analogy with the renormalization group
The difference between the twist-Hamiltonian approach and the approach of section 5.3 is analogous to the difference between RG-improved perturbation theory and fixed-order calculations. In fixed-order perturbation theory at L loops, one finds powers of logarithms log 2 x, . . . , log L x (where x is some kinematic variable) whose coefficients are related by exponentiation to coefficients at lower loop order. In RG-improved perturbation theory, we exploit this fact by choosing a scale x 0 and deriving a differential equation for the xdependence near x/x 0 = 1. The log 1 x/x 0 terms at L-loops give L-th order corrections to anomalous dimensions, beta functions, etc..
In the context of large-spin operators, the role of L-loops is played by L-twist operators in the crossed-channel. To see exponentiation of anomalous dimensions, we must in principle sum all multi-twist operators. Instead, in analogy with RG-improved perturbation theory, we assume exponentiation works and find anomalous dimensions by working at some scale y 0 . L-twist operators also give corrections to anomalous dimensions, given by the Casimirregular terms after summing their conformal blocks. These are analogous to log 1 x/x 0 terms in L-loop perturbation theory. To compute them, we must understand the detailed structure of the L-twist operators.

Crossing symmetry for the twist Hamilonian
To compute M (y, h), we need the following lemma.  where the operator C is defined as follows. Let and extend C linearly to arbitrary sums of powers and logs of y, y. Here, C Proof. By linearity, it suffices to consider f (y, y) = c(y)y a for some function c(y). Let us assume , h c(y), (7.16) and show that the sum (7.13) has Casimir-singular part c(y)y a . Since Casimir-singular terms uniquely determine an asymptotic h-expansion for coefficients of blocks, the claim follows.
As before, let ∂ = ∂ ∂ log y . The SO(d, 2) blocks have expansion Applying the differential operator in parentheses to (7.16), we get where i runs over twist families in the 1 × 2 and 3 × 4 OPEs. As in section 5.3.3, we must add (−1) (3 ↔ 4) for consistency with the symmetry properties of λ 43O i .
The contribution of an individual block to (7.21) is, Using (5.37), this has a smooth limit as As a special case, the unit operator contributes  contributions to the anomalous dimensions of [ ] 0 and [σσ] 1 come from exchange of the stress tensor T µν , and mixing is negligible. However, the operators [σσ] 0 have twist only slightly larger than T µν , so all of their contributions become important at slightly smaller h. 34 As illustrated in figure 15, exchange of large-spin [σσ] 0 operators (namely operators where the vertical distance between σ lines in figure 15 is large) looks like a product of off-diagonal terms that transition between [ ] and [σσ], coming from σ-exchange in the σ σ four-point function. This is part of the exponentiation of a twist Hamiltonian with structure The off-diagonal terms are unimportant at very large h. (We should compare the square of the off-diagonal terms to the diagonal terms.) However, they become important at slightly smaller h. In fact, because 2h ≈ 2h σ + 1, they cause the eigenvalues to repel significantly.

Computing the twist Hamiltonian
To find the twist Hamiltonian for [ ] 0 and [σσ] 1 , we must compute M σσσσ , M , and M σσ . For example, We will take the first few terms in an asymptotic expansion in large-h, so we should truncate powers of y (so that only low-twist operators contribute) before applying C. In the G σσσσ correlator, we will include terms up to order y h . Let us describe the low-twist part of the correlators G σσσσ , G , and G σσ in more detail.

G σσσσ
We have where ". . . " represents terms of higher order than y h . 35 Let us split the sum over [σσ] 0 into a finite part which we treat exactly and an infinite part which we expand in small anomalous dimensions δ [σσ] 0 , We can make δ [σσ] 0 arbitrarily small by choosing 0 large enough. Taking 0 = 6 will be sufficient for our purposes. Thus, the finite sum in (7.31) will contain the stress tensor and the spin-4 operator [σσ] 0,4 . For these contributions, we use the expansion of k 2h (1 − z) up to first order in y, Meanwhile, expanding the infinite sum in δ [σσ] 0 , we obtain The quantities λ σσ[σσ] 0 and δ [σσ] 0 are given in (6.1) and (6.2). We can compute the sums over using the methods of appendix B, where c  Only the m ≥ 2 cases will survive the C operation (because log m y is Casimir-regular for m = 0, 1). However the m = 2 term is already quite small, so it will be sufficient to truncate the series here. The first few c The S 2h −2hσ term is important because it gives a contribution to M σσσσ (y, h) proportional to y 2h , which contributes to mixing with [ ] 0 . In general, we find terms of the form Here, we can see the exponentiation discussed in section 3.1.1 at work. Summing over the family [σσ] 0 gives terms of the form y h 1 +···+h i +k , where h i are half-twists of other operators in the theory. These give contributions to the twist Hamiltonian proportional to h 1 +· · ·+h i +k, which must be matched by multi-twist operators [O 1 . . . O i ] k . This is illustrated in figure 16.
Plugging in the values (6.4) and multiplying by y 2hσ y −2hσ , the infinite sum is where ". . . " represent terms higher order in y or log y. We stress that while we have written the above coefficients numerically for brevity, they all have analytic formulae. For example, the coefficient of log 2 y y 2h is given by 1 2 c (2) 2h −2hσ in (7.36). We have written "≈" instead of "=" because the above formula is based on the approximations (6.1) and (6.2) for λ σσ[σσ] 0 and δ [σσ] 0 . Because those formulae match the numerical data to high accuracy, the same is true of (7.37). However, a more sophisticated approximation for λ σσ[σσ] 0 , δ [σσ] 0 would include contributions from operators O i other than , T , giving rise to additional terms like y h i in (7.37). 36

Comparison to numerics
We compare analytics to numerics in figures 17, 18, and 19. In figure 17, we show two sets of curves: the solid lines correspond to y 0 = 0.1, and the dotted lines correspond to y 0 = 0.02. As expected, the smaller value of y 0 introduces errors that behave approximately like log 4 yh −4hσ . The value y 0 = 0.1 gives beautiful agreement with numerics for all spins 2, so we take y 0 = 0.1 in the remaining plots.
The results show several interesting features. Firstly, we have correctly modeled the large mixing between the two families. For example, the fact that f [σσ] 1 is larger than f [ ] 0 for 26 is reproduced nicely.
We also find that M (y 0 , h) ceases to be positive-definite at h ≈ 3.4. This suggests that we cannot continue one of the twist families below this value. Indeed, in the numerical data, the family [ ] 0 ends at spin 4, which is the lowest spin such that h > 3.4. It is surprising that one can predict such a detailed fact about the low-spin spectrum using the first few terms in an asymptotic expansion at high spin. It may be a happy coincidence. Zeros in the determinant of M (y 0 , h) are responsible for the rapid oscillations and poles at the leftmost edges of figures 17, 18, and 19. 37 Additionally, we truncate M (y, h) at order y 2 , which also removes the effects of higher twist families. 38 It should be possible to surmount these difficulties with a more sophisticated analysis. If we include higher-twist families [σσ] n≥1 and [ ] n≥0 in the twist Hamiltonian, there is less downside to working at larger y 0 . On the other hand, we expect these higher families to have larger mixing with other families like [T T ], [T ], etc.. So it may be necessary to study a larger system of correlators at the same time. Alternatively, we might try to restore exponentiation of (7.42) by approximating the contribution of multi-twist operators [O · · · O] in some way. By matching Casimir-singular terms on one side of the crossing equation to asymptotic large-h expansions on the other, we can systematically solve the crossing equations orderby-order in y, y. In particular, we can reproduce a conformal block on one side with a particular large-h expansion on the other side. Our techniques for summing over twist families remove the difficulties associated with accumulation points in twist space. 39 If this order-by-order solution to crossing is systematic, where are the nontrivial constraints on the spectrum?
Note that the asymptotic large-h expansion misses terms that are Casimir-regular in both channels. That is, terms that are both Casimir-regular in y and Casimir-regular in y. 39 See [66,67] for another approach to this problem. If we write the crossing equation as then these are terms of the form y m y n log p y log q y with p, q ≤ 1. We call such terms "biregular." We have already seen examples of biregular terms in computations: for example, the y 2hσ log y and y 2hσ log y terms in the sum over [σσ] 0 in (7.37) are bi-regular, as we can see by multiplying by y −2hσ as on the right-hand side of (8.1). These are certainly nonzero, but they map to zero in the large-h expansion in either channel because S a (h) has a double zero at a = 0.
It is somewhat subtle to define the biregular part of a correlator separately from the Casimir-singular part. For example, y δ is Casimir-singular, but its expansion in small δ contains nonzero Casimir-regular terms (1 and δ log y). Indeed, no individual term in the sum (8.1) is biregular. However, biregular terms can appear when we evaluate the sum by expanding in the anomalous dimensions of double-twist operators.
To make sense of this, we propose the following prescription. Let us define an "asymptotic solution" to crossing symmetry as a set of CFT data where the dimensions and OPE coefficients of multi-twist operators have the correct asymptotic large-h behavior to reproduce all Casimir-singular terms on the other side of the crossing equation. Given Claim. If S is an asymptotic solution to crossing, then the "biregular limit" is finite. Furthermore if S is a true (not just asymptotic) solution to crossing, then L S = 0.
One can define similar biregular limits to extract biregular terms of the form y m y n log p y log q y with p, q ≤ 1. Demanding that biregular terms are crossing-symmetric gives nontrivial constraints on the spectrum.

Constraints from low-twist operators
This suggests an interesting way to derive approximate constraints on the data of the 3d Ising CFT. From our work in sections 6 and 7.5, we have approximate expressions for OPE coefficients and dimensions of the twist families [σσ] 0 , [σσ] 1 , [ ] 0 , and [σ ] 0 in terms of a finite set of initial data, namely {∆ σ , ∆ , f σσ , f , c T }. By plugging these expressions back into the correlator and demanding that biregular limits vanish, we obtain constraints on the initial data. 40 Because we have not found exact asymptotic solutions to crossing, we must approximate the limits L S in some way. We also do not have analytic approximations for the lowest spin members of the families [ ] 0 and [σσ] 1 , so we will restrict ourselves to limits involving [σσ] 0 and [σ ] 0 .
In our expressions (6.1) and (6.2) for the dimensions and OPE coefficients of the [σσ] 0 family, we treated the and T operators exactly. The biregular terms are approximately given by the log y, log y terms from expanding in small anomalous dimensions of the remaining operators [σσ] 0, ≥4 , Similarly, by demanding that the leading biregular terms cancel in the sums over [σσ] 0 and [σ ] 0 in σσ , we find the conditions and The regularized sums α and β are defined in appendix B.1.) On the right, we show the values of these quantities using the approximations in section 6 and the numericallydetermined {∆ σ , ∆ , f σσ , f , c T }. In all cases, the contributions to the limits L S , L S , L S cancel to reasonable precision.
The L S , L S , L S are interesting because their dominant contributions come from the lowest-twist operators in the theory, namely [σσ] 0 , [σ ] 0 , and indirectly σ, , T . This is based on our empirical observation that the contributions of these operators to the large-spin expansion give approximations that work well for all the operators in the twist families [σσ] 0 , [σ ] 0 . Thus, we can explore them without fully understanding the larger-twist spectrum.
By sampling values near the actual Ising point, we find that L S is much more sensitive to c T and f σσ than the other quantities ∆ σ , ∆ , f . The tangent plane to L S (c T , f σσ ) at the Ising point is given by Demanding that L S vanish gives a relationship between c T and f σσ .
In figure 20, we plot all three limits L S , L S , L S as a function of c T and f σσ , with the other quantities ∆ σ , ∆ , f held fixed at the values (2.1). The three quantities vanish nearly simultaneously at the correct values of c T and f σσ . Thus, requiring that L S , L S , L S vanish gives a way to fix c T and f σσ analytically in terms of ∆ σ , ∆ , f , to accuracy ∼ 10 −2 -10 −3 , using only the lightcone limit! 9 Discussion

Lessons for the numerical bootstrap
Traditional numerical bootstrap methods clearly probe the lightcone limit. This might explain why one must typically study a large number of derivatives around the Euclidean point z = z = 1 2 before the bounds saturate: many derivatives are needed to reach the lightcone limit, and the bounds may not saturate until the lightcone limit has been explored.
However, the Euclidean regime is also important. Because of the convergence properties of the conformal block expansion, the Euclidean regime effectively receives contributions from a small number of operators [70,88,89], and one can make surprising progress by demanding that these contributions (almost) cancel among themselves [15,90,91].
This suggests the following hybrid analytical/numerical approach 1. First solve the lightcone limit analytically using the techniques in this work. The result will be an asymptotic expansion in h, as a function of a small amount of initial data.
This step is likely easiest for theories with a relatively sparse spectrum in twist space.
Since the spectrum becomes less sparse at high-twist, we expect mixing effects in the twist-Hamiltonian to become more important in this regime. We may not find an accurate picture of the high-twist spectrum without studying a large system of crossing equations (enough to build all the necessary multi-twist operators). However, figures 1 and 2 suggest that we may not need a perfectly accurate high-twist spectrum to make progress -we just need some high-twist spectrum with approximately the right density in h-space.
2. Choose some lower cutoff h ≥ h 0 , and compute the dimensions and OPE coefficients of multi-twist operators above this cutoff using the asymptotic expansion in step 1.
Larger h 0 will mean more accurate expressions. However, smaller h 0 will leave fewer operators to solve for in step 3.
3. Plug the large-spin operators from step 2 back into the crossing equation and solve for the remaining operators in the Euclidean regime using traditional numerical bootstrap methods, the techniques of [15,90,91], 41 or some other method. We suspect that many fewer derivatives may be needed. It would also be interesting to see if this hybrid method reduces the need for high-precision arithmetic.
Unfortunately, this approach sacrifices the rigor of traditional numerical bootstrap methods because the large-h expansion is asymptotic. One must take h 0 large enough that the results saturate. (Though working with larger h 0 likely requires more derivatives.) It is encouraging that h 0 ≈ 4 is good enough for most of the results in this work. Another disadvantage is that some theories might require a large amount of initial input to compute the large-spin spectrum. For example, in this work we used {∆ σ , ∆ , f σσ , f , c T } to parameterize the large-spin spectrum of the 3d Ising CFT. We must scan over each of these parameters to explore the space of theories. In a larger system of crossing equations, we would have even more parameters.
On the other hand, the possibility of working with fewer derivatives, at lower precision, with larger systems of correlators, and perhaps without imposing unitarity (using the methods of [15,90,91]) 42 makes this hybrid approach worth exploring.

Moving towards analytics
Although we have made progress in reverse-engineering a solution to crossing symmetry analytically, numerics were crucial throughout. Let us catalog the ways in which we used numerics and discuss whether/how they can be replaced with analytics.
• Because the large-h expansion is asymptotic, numerics were crucial in determining how many terms to keep in the expansion to get reasonable results. We could also see explicitly which operators were well-described by a truncated asymptotic expansion and which ones were not. For example, [σσ] 0, =4 fits well to the analytic predictions (6.1, 6.2), while [σσ] 1, =0 does not fit the prediction in figure 17. We used this information implicitly in several ways. For example, in section 8.2, we used that the analytic predictions (6.1, 6.2) fit well all the operators in the family [σσ] 0 .
To understand these issues without numerics, it will be important to prove rigorous bounds on error terms in the large-h expansion. It would also be interesting to understand convergence properties of the twist expansion in a way analogous to the dimension expansion [70,88,89,94].
Ideally, perhaps there is a way to identify the correct representative of a given largeh equivalence class. Consistency with causality and the chaos bound [95] may be relevant, since it requires delicate cancellations between high-spin operators in a certain kinematic limit, see e.g. [96,78]. 43 It also implies bounds on dimensions of 41 The helpfulness of including higher spin Z 2 -odd operators in the "severe truncation" method of [15] has been observed previously [92]. 42 Not imposing unitarity could also help in studies of boundary and defect crossing equations, which in some cases haven't been formulated in a way that takes advantage of positivity (even in unitary theories) [91,13,93]. 43 We thank Douglas Stanford for this suggestion.
• We used numerics to discover that the contribution of other multi-twist families like [T T ] and [ T ] to the four-point functions σσσσ , σσ , is small. Consequently, we could ignore these families when diagonalizing the twist Hamiltonian for [σσ] 1 and [ ] 0 in section 7.5.
We might guess that [ T ] 0 and [T T ] 0 should be unimportant because the anomalous dimension of T is small, so only the leading term in the exponentiation of y h T −2hσ matters. However, a better treatment of this issue would involve studying correlators with T as an external operator in addition to σ and . In fact, to get a full picture of the small-twist spectrum of the 3d Ising model, we should study correlators including all the operators in [σσ] 0 as external operators. This will likely require new techniques, since the mixing matrices will be infinite-dimensional.
• We also used numerics to help choose the value y 0 at which to evaluate the twist Hamiltonian in section 7.5. The results should become less sensitive to y 0 when we study all the twist families that could contribute to M (y, h). This includes additional double-twist families like [ T ] 0 and [T T ] 0 discussed above, as well as higher-twist towers like [σσ] 2 and [ ] 1 . To completely recover exponentiation, we must also understand n-twist families with n ≥ 3. Although this may be possible with four-point functions, in practice it might require studying higher-point functions, as discussed in section 7.2.
• Although we parameterized most of the low-twist spectrum in terms of a small amount of initial data {∆ σ , ∆ , f σσ , f , c T }, it would be difficult to fix this data in practice without already knowing the answer (2.1). The biregular limits in section 8 give constraints. It will be important to understand whether they can be solved systematically.
The Euclidean regime is also important and currently the best techniques for exploring it are numerical. Perhaps the hybrid approach suggested above can help. It may also be interesting to study how recent Mellin-space approaches to the bootstrap [56][57][58] interact with the results of this work.

More future directions
A central question is: why do the truncated large-h expansions for [σσ] 0 , [σ ] 0 , etc., work so well even at small h? Perhaps our all-orders asymptotic solutions are close to an exact answer. Our work in section 7 suggests the following ansatz for the conformal block expansion: where H(h) is the twist-Hamiltonian and Λ(h) is a matrix of OPE coefficients. We have shown how to compute the large-h asymptotics of using crossing symmetry. (Asymptotics as h → ∞ along the real axis are related to asymptotics as h → ±i∞ for the class of functions we consider.) However, perhaps one could compute the full function on the left-hand side of (9.2) using a crossing kernel for SO(d, 2) conformal blocks [101]. This would remove the difficulty of working with asymptotic expansions.
One could then try to solve crossing symmetry via an iterative procedure: 1. Start with a few known operators like σ, , and T .

Compute H(h) and Λ(h) and diagonalize H(h).
3. Plug the results back in to compute corrections to H(h) and Λ(h) from multi-twist operators.
4. Repeat until the spectrum converges.
It will be interesting to explore this program in the future.
While we have focused on multi-twist operators (in particular double-twist operators), it is also interesting to consider other types of operator families like logarithmic Regge trajectories in conformal gauge theories. Such trajectories can still be described using the techniques in this work, by writing log = ∂ ∂ | =0 . 44 We expect that the techniques of section 7 give the right language for studying interesting phenomena like mixing between large-spin single-and double-trace operators in non-planar N = 4 SYM.
It would also be interesting to apply our techniques to large-N theories. Summing up the effects of graviton exchange in the bulk is important for understanding the emergence of geometry in AdS/CFT. While Virasoro symmetry makes this relatively simple in 2d CFTs, it is a difficult task in d > 2. Our all-orders results for large-spin operators may help make headway on this problem.
Finally, our new data for the 3d Ising CFT may have interesting applications to condensed matter and statistical physics. In [102], we used the low-dimension operators in table 2 to plot the Euclidean four-point function and check some inequalities from the lattice Ising model. (In this work, we can see explicitly some of the non-Gaussianity discussed in [102], from the large mixing between [ ] 0 and [σσ] 1 .) It would be interesting if some of the new operator dimensions and OPE coefficients in this work could be checked with Monte-Carlo techniques, the -expansion, or experiment.
A Numerical calculation of the 3d Ising spectrum A.1 More details on the numerics As explained in section 2.1, our strategy is to compute a partial spectrum S N (p) for several different points p on the boundary of the allowed region A N , and then choose the operators that are stable under varying p. To get to the boundary of A N , we can minimize or maximize any quantity. It is not actually necessary that (∆ σ , ∆ , f σσ , f ) themselves lie on the boundary of (2.1), as long as they don't lie outside (2.1) and some other quantity is minimal or maximal.
Extremizing an OPE coefficient is technically easier than extremizing an operator dimension because it can be done in a single optimization step. 45 In [68], we chose to maximize f σσ . In this work, we minimize c T (equivalently maximize f σσT ) as in [20]. The answers are essentially identical. We describe how to extract a partial spectrum by extremizing an OPE coefficient in section A.2.
To get a sense of the errors in the extremal functional method, we must choose a variety of points on the boundary of A N . The space of CFT data is infinite-dimensional, so random sampling is imposible. We must simply try different things, and hope some results will be invariant.
Let tan θ σ = f /f σσ . We minimize c T at the following 20 points in (∆ σ , ∆ , θ σ )-space: We do not specify the norm f 2 + f 2 σσ -it is an output of the spectrum computation, together with a list of other operators.
We assume that σ and are the only relevant scalars in the theory. In addition, we impose gaps in the Z 2 -even scalar sector (above ) and spin-2 sector (above T µν ) of the following form When we impose a gap in the spin-2 sector, we also impose the stress-tensor Ward identity The resulting spectra are mostly independent of the gaps, with one exception: in the extremal functional method, spurious operators often appear at the gaps. Some examples are the higher spin operators at the unitarity bound discussed in section 6.1. Similarly, the spectra computed using the above assumptions often (not always) have scalars of dimensions 3 or 3.5 or spin-2 operators of dimension 4 or 5. (In addition, there are occasionally Z 2odd scalars with dimension 3 due to the gap in that sector.) By varying the gaps, these operators become "unstable" in the sense that their dimensions depend on the boundary point p. Hence, in practice we can ignore them compared to the stable operators. Their OPE coefficients are usually quite small, so they don't affect the crossing equations much. We have removed these spurious operators by hand in figures 1 and 2.
We minimize c T by setting up a semidefinite program and solving it with the solver SDPB [31]. We work with Λ = 43, corresponding to 1265 derivatives of the crossing equations, and our SDPB parameters are given in   Here, ∆, , R run over the dimension, spin, and symmetry representations of exchanged operators. Each F i ∆, ,R (z, z) is a matrix whose entries are combinations of conformal blocks, and the λ ∆, ,R are vectors of OPE coefficients. For example, for a four-point function of identical scalars φφφφ , k = 1 and

A.2 Extracting spectra and OPE coefficients from SDPs
For simplicity, we first consider minimizing the 0 function with respect to the constraints (A.3). Thus, we are simply asking when it is possible to find real λ ∆, ,R such that (A.3) is true. (We comment about the case where we minimize something nontrivial later.) In [24,31], it was shown how to reformulate this question as a Polynomial Matrix Program (a special type of semidefinite program) of the following form: 46 Find y ∈ R N such that where α = (1, y) ∈ R N +1 . The notation "M 0" means "M is positive semidefinite." The M n j (x) are matrices with polynomial entries In our case, the dual objective function b vanishes because we are minimizing the 0 function.
The relation between the matrices M n j (x) and the functions F i ∆, ,R (z, z) entering the crossing equation is as follows. Firstly, j corresponds to tuples ( , R). Secondly, n corresponds to tuples (i, a, b), where i labels crossing equations and a, b are positive integers labeling derivatives. We have where ∆ j,min is the minimum dimension for j = ( , R) (e.g., the unitarity bound for an operator with spin and representation R). χ j (∆) is a positive function of ∆, written explicitly in [31]. The accuracy of the approximation (A.6) can be made arbitrarily good by increasing the polynomial degree of M n j (x). Consequently, the first N + 1 derivatives of the crossing equations (A.3) are (approximately) equivalent to where λ ∆, ,R = χ j (∆ j,min + τ ) v j,τ , ∆ = ∆ j,min + τ.
(A.8) 46 We follow the notation of [31] for the rest of this appendix.
For each j, the τ -sum in (A.7) ranges over a discrete set of nonnegative real numbers. Equation (A.7) can be rewritten as where V j,τ is a sum of outer products of v's and is thus positive semidefinite. The vectors v j,τ can be recovered from V j,τ via Cholesky decomposition. 47 Thus, if we can find τ 's and V j,τ 's such that (A.9) holds, then (A.8) gives a set of dimensions and OPE coefficients that solve N + 1 derivatives of the crossing equations and are consistent with unitarity.
It is simple to find the appropriate τ 's. The solver SDPB returns a vector y ∈ R N which can be assembled into α = (1, y) ∈ R N +1 satisfying the constraints (A.7). Taking the inner product of (A.9) with α, we find 0 = j,τ Tr(V j,τ (α · M j (τ ))).
(A. 10) This implies that each term in (A.10) vanishes individually, since each term is nonnegative. However, this is only possible if α · M j (τ ) is a degenerate m j × m j matrix. Thus, τ must be a nonnegative zero of det(α · M j (x)).
The function det(α·M j (x)) is constrained to be positive for x ∈ [0, ∞). Thus, its positive zeros for must be double zeros. In numerical computations, it never actually attains the value zero, but instead dips very close to to the x-axis. Thus, it's more convenient to compute the τ 's as local minima of det(α · M j (x)), which can be computed as zeros of its derivative (together with a possible zero at x = 0, which must be checked separately).
The matrices V j,τ can be obtained by solving the linear algebra problem (A.9). However, they are also already encoded in the primal solution computed by SDPB. Let d j = max n [deg(M n j (x))]. The primal solution is a vector u ∈ R P where P counts the number of tuples (j, r, s, k) with 1 ≤ r ≤ s ≤ m j , k = 0, . . . , d j . 48 We can assemble u into symmetric matrices where E rs is a symmetrized unit matrix with components (A.12) 47 The matrix V j,τ typically has low rank, which means that numerically it may have very small negative eigenvalues. Thus, instead of using a Cholesky decomposition, we simply compute its eigenvectors and throw out those with very small eigenvalues. For computations in this work, it suffices to keep only the first eigenvector. We expect low rank matrices because a higher-rank matrix would mean that the determinant of the functional has a higher-order zero at a fixed ∆, which is non-generic. An exception occurs if an operator is isolated in ∆-space, which is why one can obtain stronger constraints by imposing that the matrix associated to the operator is rank-1 as in [55]. We thank Slava Rychkov and Alessandro Vichi for discussions on this point. 48 We use the notation u instead of [31]'s x to avoid confusion with the sample points x We estimate errors as standard deviations in our sample set of 60 partial spectra. It is important to remember that these error estimates are non-rigorous (in contrast to the bounds on ∆ σ , ∆ , f σσ , and f , in (2.1), which are rigorous). In fact, the tables show that this method of assigning errors is imperfect. In table 6, for example, the precision of the OPE coefficients f σ O varies significantly at large . For instance, the error for the OPE coefficient of [σ ] 0, =27 is 0.3%, while the error for [σ ] 0, =28 is 2%. It is surprising that these should be so different. Perhaps a wider scan of the boundary of the allowed region A N would equalize the errors somewhat. Regardless, the reader should take the error estimates with a grain of salt.
Because we have chosen different conventions for conformal blocks, our OPE coefficients are normalized differently from those in [68]. Specifically, the leading terms in the conformal block expansion are given by , [68] . To compute the Casimir-singular terms, we must match asymptotic expansions To compute Casimir-regular terms, we expand k 2h (1 − z) in small y and naively switch the order of summations,   We must now regularize the sum over h. Using ∂h ∂ = 1 + ∂δ ∂h 0 , one can show Now form the asymptotic expansions with coefficients c (m) a and sets A m . (When m = 0, these reduce to c a and A above.) Note The derivative ∂ h 0 decreases degree in by 1. Thus, the combination Note that f k ( , h 0 ) as we've defined it is analytic in k, so we can form the derivative β k [p, δ](h 0 ). The above result generalizes easily to the case of alternating or even sums, where we must simply replace A → A − or A → A even and modify the sum over appropriately.

B.1 Special cancellations between singular and regular parts
We sometimes encounter sums where both the Casimir-singular and Casimir-regular part naively diverge, but the divergences cancel to leave a finite quantity. This occurs in sums over un-mixed blocks with coefficients lim →0 Γ(− ) 2 S (h) and in sums over mixed blocks with coefficients Γ(− )S r,s −r (h). In such sums, the naive Casimir-singular parts are proportional (23)       We define regularized quantities α k by replacing S 0 (h) → S (h) and S r,s −r (h) → S r,s −r (h), adding the quantities in parentheses in (B.9) or (B.10) to α k , and then taking the limit → 0. Examples of this procedure are given in (6.22) and (6.29). In general, we must apply it whenever the asymptotic large-h expansion of p(h) contains terms of the form lim →0 Γ(− ) 2 S (h) or Γ(− )S r,s −r (h).

C Box diagrams
We claim that the Casimir-singular part of the box diagram in figure 21 is the same whether we read the diagram from left-to-right or bottom-to-top. Consequently, the Casimir-singular part of the sum of box diagrams over all possible internal legs is crossing-symmetric. 51 We can regard any CFT as a 2d theory with SL(2, R) × SL(2, R) symmetry. If our claim holds in 2d, it holds in general dimensions. A benefit of working in 2d is that tensor structures are extremely simple, so we can prove the claim for external operators of any spin (not just scalars). Some conventions in 2d theories are different from those in the main text. We have the two and three-point functions φ(z 1 , z 1 )φ(z 2 , z 2 ) = 1 (C.1) In unitary theories with these conventions, operators satisfy the reality property φ(z, z) † = (−1) φ(z, z).

(C.2)
That is, even-spin operators are real and odd-spin operators are imaginary. The three-point coefficients have the same reality properties as the product of operators. They satisfy the symmetry property f abc = f bac (−1) a+ b + c and similarly for other permutations.
Consider a four-point function φ 1 · · · φ 4 where the operators φ i have weights (h i , h i ) (not necessarily equal). We have the conformal block expansion φ 1 (z 1 , z 2 ) · · · φ 4 (z 4 , z 4 ) = On the left-hand side, this matches to the contribution of the double-twist operators [12] n . We have λ 12 [12]n (h)λ 34 [12] where we have used f abc = f cab . Thus, Isolating the contribution from a, c, the sums factorize into holomorphic and antiholomorphic parts. The antiholomorphic sum is That is, we have the truly remarkable identity We should sum over unordered pairs b, d, with weight 1/2 for the case b = d because only even-spin operators are present. This is equivalent to dropping the second term and summing over ordered pairs b, d. Overall, the Casimir-singular part in y, y is a,b,c,d The summand is invariant under (1, a, d, y, m) ↔ (3, b, c, y, n) and multiplication by the phase (−1) 1 + 2 + 3 + 4 coming from the prefactor in the crossing equation (C.4). Thus, the above result is crossing-symmetric.
There are some special cases of this calculation that we must take care with. Each special case will require us to modify the calculation slightly. We should check that the results are still crossing-symmetric after these modifications.

(C.24)
The combination above is f 12 [12]n (h)f 34 [12]n (h) γ [12] Recall that we are interested in computing only the Casimir-singular terms with respect to the crossed-channel. The Casimir-regular terms are proportional to y k−h 3 +h 2 , y k−h 1 +h 4 (C. 26) Only the parts proportional to y k−h 3 +h 2 log y or y k−h 1 +h 4 log y are Casimir-singular. To determine them, we can use ∂ ∂h k r,s 2h (z) = log y k r,s 2h (z) +  The y-Casimir-regular terms in the 12 → 34 channel start at n = 1 (i.e. not at leading order in y), because of the fact that k 2h has leading order y h , with a coefficient that is independent of h.
In the case where we have [db] = [12] = [34], the Casimir-singular terms are given by expanding to second order in the anomalous dimension γ [12]n .
n,h f 12 [12]n (h) 2 k 2(h 1 +h 2 +n+γ [12]  where 3, 4 = 1, 2 or 3, 4 = 2, 1, depending on which correlator we're studying. The Casimirregular terms start at order n = 1, by the same logic as before. Another type of special case is when h a + h c = h 1 + h 4 or h a + h c = h 2 + h 3 or both. Naively, this leads to poles in the holomorphic part of the Casimir-singular terms. However, the Casimir-singular and Casimir-regular parts combine to give a finite result in these cases. and then take the limit. This exactly corresponds with the prescription for h double poles (C.32), and again these cases are related by crossing.