Bootstrapping the O(N) Vector Models

We study the conformal bootstrap for 3D CFTs with O(N) global symmetry. We obtain rigorous upper bounds on the scaling dimensions of the first O(N) singlet and symmetric tensor operators appearing in the $\phi_i \times \phi_j$ OPE, where $\phi_i$ is a fundamental of O(N). Comparing these bounds to previous determinations of critical exponents in the O(N) vector models, we find strong numerical evidence that the O(N) vector models saturate the bootstrap constraints at all values of N. We also compute general lower bounds on the central charge, giving numerical predictions for the values realized in the O(N) vector models. We compare our predictions to previous computations in the 1/N expansion, finding precise agreement at large values of N.


Introduction
Conformal field theories (CFTs) offer delightful examples of quantum field theories that are strongly coupled, yet contain enough symmetry and structure that they may turn out to be tractable if the right techniques are found. Until recently, the idea of exploiting this structure in order to find complete non-perturbative solutions to theories was only carried out successfully in 2D, most notably in the seminal work of [1]. However, over the last several years great progress has been made at developing the conformal bootstrap [2] approach to CFTs in D > 2, where a large number of nontrivial bounds have been found which follow very generally from the constraints of crossing symmetry and unitarity [3][4][5][6][7][8][9][10][11][12][13][14][15]. The results obtained so far have been particularly striking in 3D, where it was found in [12] that the CFT described by the critical 3D Ising model occupies a special place in the space allowed by crossing symmetry and unitarity. Moreover it appears possible that a robust numerical solution to this theory can be obtained using bootstrap techniques [16].
In this paper we will extend the work of [12] to study 3D CFTs with an O(N) global symmetry using the conformal bootstrap. We will focus on theories containing a scalar field φ i (x) in the vector representation of O(N). The most notable theories falling into this class are the critical O(N) vector models [17,18], which describe second-order phase transitions in a variety of real-world systems at small values of N [19], and are also solvable in a 1/N expansion at large values of N (see [20] for a review). Moreover, the O(N)-singlet sector of this theory is thought to be holographically described by a higher-spin gauge theory in AdS 4 [21]. Previously, bootstrap ideas have been applied to the O(N) vector models in the 1/N expansion, for example in work by Lang and Rühl [22][23][24][25][26][27], Petkou [28,29], and more recently Maldacena and Zhiboedov [30,31]. Our approach allows bootstrap constraints to be studied (albeit numerically) at any value of N.
Our primary goal, following previous numerical studies of the bootstrap [3][4][5][6][7][8][9][10][11][12][13][14][15], will be to place general upper bounds on the scaling dimensions of the first nontrivial scalar operators (both O(N) singlets S and O(N) symmetric tensors T ij ) in the φ i × φ j operator product expansion (OPE). We will also place general lower bounds on the central charge c, defined as the coefficient appearing in the two-point function of the stress-energy tensor. We will then compare these bounds to the best previous results based on Monte Carlo simulations, analytical estimates, and the 1/N expansion. In all cases where results are known (including small values of N), we find that the O(N) vector models saturate our bounds, and moreover sit near special locations in the space allowed by crossing symmetry. Inputting the previously measured values of ∆ φ , this allows us to give sharp predictions for ∆ S , ∆ T , and c for different values of N.
In order to efficiently implement the bootstrap in the presence of a global symmetry, we will use techniques based on semidefinite programming, developed for 4D CFTs in [10]. Here we will show how to adapt this technique for CFTs in an arbitrary number of spacetime dimensions. This requires approximating conformal blocks as rational functions of the exchanged operator dimension ∆, and we will show that such a rational approximation follows directly from a recursion relation expressing a general conformal block as a sum over poles ∼ 1/(∆ − ∆ * ) occurring at special (non-unitary) values of the dimension ∆ * where the conformal multiplet contains a null state. This conformal block representation generalizes an idea of Zamolodchikov, first applied to Virasoro blocks in 2D [32,33], to arbitrary space-time dimensions. This paper is organized as follows. In section 2 we review the formulation of the conformal bootstrap for CFTs containing an O(N) global symmetry, as well as convex optimization techniques for placing bounds on operator dimensions and OPE coefficients. In section 3 we show how to find rational representations for conformal blocks in arbitrary space-time dimensions, presenting a new recursion relation for conformal blocks as a sum over poles in ∆. In section 4 we present our bounds and a comparison with the O(N) vector models. We conclude in section 5.
2 Conformal Bootstrap with O(N ) Global Symmetry

Statement of Crossing Symmetry
Let us briefly review the formulation of the conformal bootstrap for 3D CFTs with an O(N) global symmetry. Further details can be found in [8][9][10]. We focus on theories containing a scalar primary operator φ i of dimension ∆ φ , transforming as a fundamental under O(N). The operator product of φ i with itself takes the schematic form where S + denotes O(N) singlets of even spin, T + denotes O(N) symmetric tensors of even spin, and A − denotes O(N) anti-symmetric tensors of odd spin.
By pairing up φ's and performing the OPE, a four-point function can be decomposed into conformal blocks as where each sum runs over primary operators O of dimension ∆ and spin ℓ appearing in 3) The four-point function itself should be independent of how we perform the OPE. Swapping (1, i) ↔ (3, k), we find two different conformal block expansions of a single fourpoint function which must agree with each other. Writing out this condition and isolating the coefficient of each tensor structure that appears, we obtain three equations which can be grouped into a vector "sum rule" The coefficients λ 2 O appearing in (2.4) are unknown a-priori, with the exception of the unit operator which has λ O = 1 if φ i is canonically normalized. However, we do know that the OPE coefficients λ O must be real in unitary theories, which means that λ 2 O is positive. Further, in D-dimensional unitary theories the operator dimensions must satisfy the lower bounds [34][35][36][37][38] ∆ ≥ where saturation occurs for free scalars (ℓ = 0) or conserved currents (ℓ > 0).

Bounds from Convex Optimization
From here, we follow the general strategy of [3] for putting bounds on CFT data. Let us begin by isolating the unit operator in (2.4), The procedure is as follows: 1. Make an assumption about the CFT spectrum, for instance that all singlet scalars have dimension above some ∆ * .
2. Try to find a linear functional α such that 3. If such a functional exists, the assumption (1) is ruled out, since applying α to Eq. (2.8) gives a contradiction. If not, we cannot conclude anything about our assumption.
Step 2 requires us to solve a convex optimization problem: we must search over the vector space F of linear functionals, subject to linear constraints of the form α(V ) ≥ 0.
Each linear constraint restricts us to a half-space of F , and together these half-spaces carve out a convex subset C ⊂ F . We would like to determine whether this subset is non-empty.
In the case at hand, we will take our functional α to be linear combinations of derivatives with respect to the cross-ratios z, z around the crossing symmetric point z = z = 1/2, The parameter k controls the dimension of the space of linear functionals that we search over. Note that truncating this search space leads to valid, though possibly suboptimal, bounds. As we increase k, our bounds get better and better, converging to an optimal bound as k → ∞.

Formulation as a Semidefinite Program
A key difficulty in our convex optimization problem is that we have an infinite number of constraints on α -one for each O which could appear in the OPE. We must impose α(F R,∆,ℓ ) ≥ 0 for all representations R, dimensions ∆, and spins ℓ obeying our assumptions about the spectrum. An efficient way to deal with this infinity is to reformulate our problem as a semidefinite program, which can include constraints of the form: Systems of these inequalities can be solved efficiently using interior point methods.
To write our constraints α(V ) ≥ 0 in this form, it suffices to find an approximation where χ ℓ (∆) are positive functions, and P R,ℓ,i (∆). The dimensions ∆ satisfy bounds ∆ ≥ ∆ min,ℓ , so writing ∆ = ∆ min,ℓ + x yields a set of inequalities in the form (2.11).
In [10], special analytic expressions for conformal blocks in even dimensions [39,40] were used to derive approximations of the form (2.12), which proved sufficient for applying semidefinite programming to even dimensional CFTs. These approximations worked surprisingly well, but it was unclear how to generalize the techniques to CFTs in odd (or fractional) space-time dimension.
In the next section, we will show that the existence of approximations (2.12) in any space-time dimension follows naturally from representation theory of the conformal group. This is sufficient for formulating our optimization problem as a semidefinite program, which can then be solved using one of the many freely available semidefinite program solvers. We give details of our implementation using the solver SDPA-GMP in Appendix B.

Why Rational Approximations Exist
To compute CFT bounds using semidefinite programming, we need precise, systematic approximations for conformal blocks g ∆,ℓ in terms of positive functions times polynomials in ∆, or equivalently positive functions times rational functions of ∆ with positive denominator. The existence of such approximations follows from conformal representation theory. Recall that the conformal block g ∆,ℓ is a sum over states in radial quantization where O is a conformal primary of dimension ∆ and spin ℓ, and α runs over O and all of its descendants. It will be convenient for our discussion to use the radial coordinates of [41], where the points x i are arranged as in Fig. 1, and the coordinate ρ is given in terms of cross-ratios by The states P µ 1 · · · P µn |O have eigenvalue ∆+n under dilatation and can be decomposed into traceless symmetric tensor representations of the rotation group. It follows that the sum over states (3.1) can be written where r = |ρ|, η = cos θ = (ρ + ρ)/2|ρ|, ν = D/2 − 1, and the C ν j (η) are Gegenbauer polynomials. The coefficients B n,j can be computed straightforwardly by solving the conformal Casimir equation term by term around r = 0. Here, D is a second-order differential operator in the crossratios u, v, representing the action of the quadratic Casimir of the conformal group on a four-point function, Further details can be found in [41].
Solving the Casimir equation to second order in r, we find [41] .
Let us make some comments about these coefficients. First, B n,j is a rational function of ∆. This is a simple consequence of the Casimir equation (3.4), but it follows more directly from the expression for g ∆,ℓ as a sum over states. Each term in the numerator and denominator of (3.1) can be computed from the action of the conformal algebra on |O . This action is polynomial in ∆, so B n,j is rational in ∆.
Secondly, the denominators of B n,j are positive, as long as ∆ obeys the unitarity bound (2.7). This is because poles in B n,j occur when a state |α becomes null. However, the absence of null or negative-norm states in a conformal multiplet is precisely what defines the unitarity bound.
These facts are exactly what we need. We can obtain good approximations to conformal blocks by truncating the expansion (3.3) at some high order r ∆+N . Each term in the resulting expression will be a rational function of ∆, times an overall factor of r ∆ . The error is naïvely of order r N , which is quite small at the crossing-symmetric point r = 3−2 √ 2 ≈ 0.17. Indeed, this naïve estimate is correct, since it can be shown that the coefficients B n,j are uniformly bounded as a function of ∆ for all n, j.
The cost of keeping more terms in the series expansion (3.3) is that the degree of the resulting rational approximation grows. As we'll see in the next subsection, this growth is slow and under control. Further, additional tricks can improve the rational approximation significantly without increasing the degree. We discuss these in Appendix A.

Poles in ∆ and Recursion Relations for g ∆,ℓ
We can be more precise about the structure of our rational approximations by exploiting an idea of Zamolodchikov, originally applied to Virasoro blocks in 2 dimensions [32,33].
Recall that poles in g ∆,ℓ as a function of ∆ occur precisely at special values ∆ = ∆ * where some state |α = P n |O in the conformal multiplet of |O becomes null. When this happens, all the descendants of |α become null as well, and together they form a nontrivial sub-representation of the (now reducible) multiplet of |O . Since the pole in ∆ gets contributions from all states in this sub-representation, its residue is proportional to a conformal block for a primary with the same dimension and spin as |α : where ∆ α = ∆ * + n, and c α is a coefficient which is independent of conformal cross-ratios.
Poles in ∆ determine g ∆,ℓ up to a function which is analytic on the entire complex plane. Thus, we can write is an entire function of ∆. The block g ∆,ℓ has an essential singularity of the form r ∆ as ∆ → ∞. Stripping this off, we have Note that the entire function h is now independent of ∆, since it has no singularities as ∆ → ∞.
With the recursion relation (3.10), we can be more explicit about the form of our positivetimes-polynomial approximation (2.12). Truncating (3.10) to a finite number of poles ∆ i , it's clear that , (3.11) and r = 3 − 2 √ 2 is evaluated at the crossing-symmetric point.
Let us now determine the data entering (3.10). The function h (∞) ℓ can be computed easily by solving the conformal Casimir equation to leading order in ∆, In principle, the pole positions ∆ i and coefficients c i are determined by the conformal algebra. It would be interesting to compute them directly. In practice, by inspecting the solution to the conformal Casimir equation, we find that there are three series of poles, described in Table 1, whose coefficients are as follows: k = 1, 2, . . . , ⌊ℓ/2⌋ Table 1: The positions of poles of g ∆,ℓ in ∆ and their associated data. There are three types of poles, corresponding to the three rows in the table. The first two types exist for all positive integer k, while the third type exists for positive integer k ≤ ⌊ℓ/2⌋. The coefficients c 1 (k), c 2 (k), c 3 (k) are given in Eqs. (3.13).
The recursion relation Eq. (3.10), with poles listed in Table 1, reveals another fact that proves useful for our implementation: when we truncate the series expansion for h ∆,ℓ at order r N , the degree of the resulting rational approximation grows like N. Not only is each coefficient B n,j a rational function, but the number of new factors which enter the denominator as we increase N → N + 2 is either 2 or 3 (depending on whether N ≤ ℓ). In practice, this means that we can compute the expansion (3.3) to extremely high order without incurring too much of a performance hit from dealing with large degree polynomials.
Once we know the residues c i (k), our recursion relation (3.10) provides an extremely efficient way to compute conformal blocks, either analytically in a series expansion, or numerically (for any number of derivatives around any point (r, η)). It would be very interesting to generalize these ideas to other four-point functions -for example, external scalars φ i with different dimensions ∆ i , or external operators with spin. Although much progress has been made computing conformal blocks for these types of operators using series manipulations [39,40,42,12,41], derivative relations [43,44], and conformal integrals [45], we believe the residues c i are the most convenient and directly useful data for numerical bootstrap applications.

Results and Comparison to O(N ) Vector Models
In this section we will present our results from the bootstrap. We focus our attention on computing bounds on the lowest singlet dimension ∆ S , the lowest symmetric tensor dimension ∆ T , and the central charge c, as a function of the external scalar dimension ∆ φ .
Let us begin by summarizing what is known about these quantities in the O(N) vector (5)   models. In these theories, results are typically phrased in terms of critical exponents, rather than scaling dimensions. Concretely, ∆ φ is related to the critical exponent η via ∆ φ = 1/2 + η/2. The dimension of the O(N) singlet operator S is related to the critical exponent ν via ∆ S = 3 − 1/ν. Finally, the dimension of the O(N) symmetric tensor operator T is related to the crossover exponents φ c and η c , through the relations ∆ T = 3 − φ c /ν = 1 + η c . In Table 2 we show the most accurate determinations of these dimensions that we have found in the literature for N = 1, 2, 3, 4, 5, 6 (with N = 1 being the 3D Ising model).
These quantities have also been computed in the large N limit. ∆ φ is known to order 1/N 3 while ∆ S has only been computed to order 1/N 2 (see [20] and references therein). The crossover exponent connected to ∆ T was also computed to order 1/N 2 in [53]. The results are: Finally, let us mention that the leading correction to the central charge c was computed in [28,29] to be where c free = D/(D − 1) is the central charge of a free scalar field.

Bounds on O(N ) Singlets
Now we determine a general bound on the O(N) singlet operator dimension ∆ S following the procedure described in section 2.2. We assume there is a gap in the CFT spectrum so that all singlet scalar operators have dimension greater than ∆ S , all symmetric tensor scalars have dimension greater than 1, and the dimensions of all the other operators are constrained only by the unitarity conditions. Note that due to the assumption on symmetric tensor scalars this is not the most general bound. However, we found that this mild assumption improves numerical stability while not significantly affecting the bound on ∆ S -moreover the assumption is certainly satisfied for O(N) vector models, as can be seen from previous determinations of the operator dimensions (see Table 2).
The boundaries for the allowed values of ∆ S as a function of ∆ φ are shown in Fig. 2. These bounds are determined by a bisection search in ∆ S to within 10 −3 . The parameter k of section 2.3, controlling the number of derivatives in the functional α, is set to k = 10 everywhere. For a given N, only the values of ∆ S below the corresponding solid line are allowed.
In Fig. 2 we see that the bounds on ∆ S grow monotonically from ∆ S = 1 at ∆ φ = 0.5, the point corresponding to the non-interacting theory. At a certain value of ∆ φ , each boundary line exhibits a change in the slope. This type of behavior was already discussed for the Ising model bound in Ref. [12], where it was found that the change in slope occurs at the values of ∆ φ and ∆ S corresponding to Ising model point. The changes in slopes are somewhat less sharp for the O(N) theories; however, they occur in the vicinity of the O(N) vector model points, shown as points with error bars in Fig. 2. These points always seem to lie very close to the boundary; in fact it is possible to rigorously reduce the previous error bars using the bound we obtained.
For large values of N, we can also compare the bound with the results for ∆ φ and ∆ S obtained using the 1/N expansion, Eq. 4.1. The 1/N expansion results (shown as black crosses in Fig. 2) are consistent with our bound at large N. Moreover, the change in the slope of the boundary line is sharper for large N and it occurs very close to the O(N) vector model points.

Bounds on O(N ) Symmetric Tensors
Bounds on the dimension of the first symmetric tensor scalar operator ∆ T are obtained following a similar procedure as for the singlet bounds. We assume that all the symmetric tensor scalars have dimension greater than ∆ T , all singlet scalars have dimension greater than 1, and all other operator dimensions satisfy unitarity conditions. The allowed values of ∆ T as a function of ∆ φ are shown in Fig. 3 -only the values of ∆ T below the solid line are allowed. Just like in the case of the singlet operators, the boundary lines start at the free field point, ∆ φ = 0.5, ∆ T = 1 and grow monotonically. The change of the slope is more gradual than for the singlet bound. The O(N) vector model points from Table 2 appear to be consistent with the bound to within error bars for N = 2, 3, 4, with some mild tension between the N = 5 bound and the quoted results from [52] and [49].
At large N the symmetric tensor dimension ∆ T approaches 1 in the O(N) model. This is reflected in the fact that our bounds become more and more constraining as N → ∞. Moreover, we find that the order 1/N 2 computation of ∆ T in Eq. 4.1 becomes very close to our bound (and occurs near changes of slope) at large but finite values of N, shown in black crosses at N = 10, 20, ... in Fig. 3.

Bounds on the Central Charge
The central charge c is defined as a coefficient in the two-point correlation function of the canonically normalized stress tensor: where S D = 2π D/2 /Γ(D/2). In our notation, the central charge is related to the OPE coefficient of the stress tensor λ S,3,2 by where c free = D/(D − 1) is the central charge of a free scalar field. We can find an upper bound on this OPE coefficient as follows. Rewrite the sum rule (2.8), separating the contribution of the stress tensor: Applying a functional α such that α(V O ) ≥ 0 for all operators in the spectrum other than unit operator and normalized so that α(V S,3,2 ) = 1, Eq. (4.5) then yields the inequality Finding the functional α that minimizes −α(V unit ) then gives the strongest upper bound on the OPE coefficient. By Eq. (4.4) this implies a lower bound on the central charge.
The most general bound would be obtained by making no assumptions about the operator spectrum, except that they obey unitarity conditions. However, we can obtain a somewhat stronger bound by making additional assumptions about the spectrum. In particular, we can assume there are gaps in the spectrum of singlet and symmetric tensor operators, as long as they are consistent with the results of previous subsections. Here we will assume mild gaps, ∆ S ≥ 1 and ∆ T ≥ 1. This assumption on the operator dimension spectrum is not too stringent; for example, we know from previous determinations that O(N) vector models satisfy these conditions, see Table 2.
The central charge bound as a function of the scalar dimension ∆ φ is shown in Fig. 4. The central charge approximately scales linearly with N (exactly in the non-interacting theory), so we have plotted c scaled to Nc free . At low values of ∆ φ , all of the bounds approach the same asymptote. The slope of the asymptote is −10/3, which is the same curve that one obtains in the N → ∞ limit from Eqs. To obtain stronger bounds on the central charge we can introduce larger gaps in the operator spectrum. In the plots of Fig. 5 we assumed that the gap in the singlet scalar spectrum saturates the bound obtained in subsection 4.1, while the gap in the symmetric tensor scalar spectrum is kept at ∆ T ≥ 1. At low values of ∆ φ , the bounds again approach the same asymptote and in general don't differ too much from the bounds in the Fig. 4. However, here the bounds exhibit a change in the slope at certain value of ∆ φ . At larger values of ∆ φ the bounds are much stronger than the ones in slope is more gradual, but still occurs at ∆ φ close to the known values in the O(N) models.
Just like the operator dimensions in the previous subsections, the central charge values obtained using 1/N expansion in the case of large N lie very close to the boundary line. Encouraged by these facts, we can conjecture that the central charge values will lie on the boundary even for low values of N. Using the values of ∆ φ determined by other methods, we can then make a prediction for the values of c. These are shown, along with our bootstrap predictions for ∆ S and ∆ T , in Table 3. In calculating these values of c we used the bounds of Fig. 5, since these are our strongest bounds for the O(N) vector models.

Discussion
Let us take a moment to reflect on the results of the previous section. First, we have discovered the remarkable fact that operator dimensions in the critical O(N) vector models take on values which saturate general constraints from crossing symmetry and unitarity. This gives us an organizing principle by which we can understand why these theories are special -gaps in the spectrum of operator dimensions are maximized in a way that  is consistent with unitarity. It will be interesting to verify that this trend continues for operators of higher dimension. A promising approach to extracting more of the spectrum in these theories is to consider the locations of the zeros of α(V O ) along the boundary [6,14,16]. We hope to develop this approach in the O(N) models in future work.
As far as we are aware, we have presented the first predictions for the central charge in the O(N) models at small values of N. In doing so we have verified that there is approximately linear growth with N and that c < Nc free for each value of N. It will be interesting if these predictions can be verified in lattice simulations of the O(N) models -this will require a robust lattice construction of the stress-energy tensor, which is a worthwhile task in its own right. One can also easily extend these methods to determine the flavor central charges, In this work we only considered the constraints from crossing symmetry of φ i φ j φ k φ l . It is very interesting to extend this analysis to include constraints from other correlators, such Such extensions are e.g. necessary in order to learn about the Z 2 odd operators in the spectrum of these theories. Studying these correlators may also help to give a sharper criterion that can be used to determine the value of ∆ φ in the O(N) models, going beyond the fact that it appears to take a value near (somewhat smooth) changes in slope of the general bounds.
The recursion representation for the conformal blocks presented in section 3 is a powerful and efficient method of computing conformal blocks in any number of dimensions. This representation can for example be utilized in studies of CFTs that interpolate between 2 < D < 4, in D = 5, 7 where conformal blocks are similarly complicated, or perhaps in constructing an argument (extending [54]) that nontrivial CFTs in large D do not exist. The representation as a sum over poles in ∆ may also be useful for making general analytic arguments (going beyond the large spin arguments of [55,56]) in the context of the conformal bootstrap. The sum over poles may also be particularly interesting from the perspective of Mellin amplitudes [57][58][59][60][61][62][63][64][65].
Finally, we'd like to emphasize that the O(N) models at N = 2, 3 have numerous beautiful realizations in experimental condensed matter systems. E.g., the O(2) model describes the superfluid transition in 4 He and the bicritical point in uniaxial magnets such as GdAlO 3 , while the O(3) model describes the Curie transition in simple isotropic magnets such as Ni, Fe, and EuO. Many more examples can be found in [19]. Thus, the conformal bootstrap in 3D allows one to realize the physicists' dream -it makes quantitative predictions in strongly-interacting systems that can be experimentally tested!

A Improving Rational Approximations
When truncated to a finite number of poles ∆ i , the recursion relation gives a rational approximation for h ∆,ℓ = r −∆ g ∆,ℓ as a function of ∆. The precision of this approximation increases as we include more and more poles ∆ i . However, the degree increases as well, and this can be problematic for computation. Larger degree polynomials slow down semidefinite program solvers.
A useful compromise is to keep n poles ∆ 1 , . . . , ∆ n with the largest residues, and use these as a basis to approximate other poles ∆ j with smaller residues. That is, for poles with small residues, we write where the coefficients a i are chosen to make the approximation as good as possible. In this way, we can approximately include the contribution of the pole at ∆ = ∆ j without increasing the degree of our rational function. Note that the ∆ i lie below the unitarity bound ∆ unitarity , so ∆ itself never approaches a pole when we compute CFT bounds.
How should we choose the coefficients a i ? We need Eq. (A.2) to hold to high accuracy across all ∆ ≥ ∆ unitarity (away from the singularities on both sides). A method that works well in practice is to ensure that Eq. (A.2) and its first n/2 derivatives hold exactly at ∆ = ∆ unitarity and at ∆ = ∞. These conditions give n linear equations which determine the a i .
In practice, including poles with residues less than 10 −2 yields rational approximations to conformal blocks which are correct to within 10 −9 . Including poles with residues less than 10 −10 yields approximations correct to within 10 −22 .

B Implementation in Mathematica and SDPA-GMP
A brief summary of our implementation is as follows. All steps but the last are performed in Mathematica.
1. We compute a rational approximation for derivatives of conformal blocks at the crossing symmetric point ∂ m r ∂ n η g ∆,ℓ ≈ r ∆ P (m,n) ℓ (∆)/Q ℓ (∆), where r = 3 − 2 √ 2. This can be done using either the Gegenbauer expansion and Casimir equation described in [41], or more efficiently using our recursion relation (3.10). The recursion relation can be implemented numerically in the space of vectors of (r, η) derivatives, where multiplication by r k is a matrix on this space. We compute approximations up to order 60 in the r-expansion.