Bootstrapping Mixed Correlators in the 3D Ising Model

We study the conformal bootstrap for systems of correlators involving non-identical operators. The constraints of crossing symmetry and unitarity for such mixed correlators can be phrased in the language of semidefinite programming. We apply this formalism to the simplest system of mixed correlators in 3D CFTs with a $\mathbb{Z}_2$ global symmetry. For the leading $\mathbb{Z}_2$-odd operator $\sigma$ and $\mathbb{Z}_2$-even operator $\epsilon$, we obtain numerical constraints on the allowed dimensions $(\Delta_\sigma, \Delta_\epsilon)$ assuming that $\sigma$ and $\epsilon$ are the only relevant scalars in the theory. These constraints yield a small closed region in $(\Delta_\sigma, \Delta_\epsilon)$ space compatible with the known values in the 3D Ising CFT.


Introduction
The conformal bootstrap [1] in D > 2 dimensions has produced remarkable results, including numerical bounds on operator dimensions and OPE coefficients , analytical constraints [26][27][28][29][30][31][32], and recently a precise conjecture for the Z 2 -even spectrum of the 3D Ising model [22]. However, all previous studies have focused on a single four-point function φφφφ containing identical operators (sometimes in a nontrivial global symmetry representation). It is extremely important to ask how other correlators, such as φφφ 2 φ 2 or more generally φ 1 φ 2 φ 3 φ 4 for different operators φ i , additionally constrain the space of CFTs.
An immediate complication is that the unitarity properties of such mixed correlators are more intricate because coefficients in the conformal block expansion are not necessarily positive. In section 2, we describe how these unitarity properties can be captured by a semidefinite program with a continuously infinite number of constraints (as opposed to the linear programs that arise in the single correlator case). Applying methods introduced in [9,15], we rewrite this as a higher-dimensional semidefinite program with a finite number of constraints, which can be solved on a computer. 1 In section 3, we specialize our discussion to the case of scalar correlators with a Z 2 global symmetry.
To formulate our semidefinite program, we need approximations for conformal blocks as rational functions of the exchanged operator dimension ∆. Such approximations follow from a rapidly convergent expansion for the blocks as a sum over poles 1/(∆ − ∆ i ). We describe this expansion in section 4, generalizing the results of [15] to non-equal external operator dimensions. Our expressions give higher dimensional analogs of a recursion relation for Virasoro conformal blocks developed by Alyosha Zamolodchikov in [33,34].
In section 5 we apply our formalism to numerically study operator dimensions in 3D CFTs with a Z 2 global symmetry, a class of theories that includes the 3D Ising model. We focus on the system of four-point functions { σσσσ , σσǫǫ , ǫǫǫǫ } containing the lowest Z 2 -odd scalar σ and Z 2 -even scalar ǫ, and ask: What are the allowed dimensions (∆ σ , ∆ ǫ ) assuming that σ and ǫ are the only relevant scalars in the CFT?
The existence of only two relevant scalars is an obvious experimental fact about the 3D Ising CFT -it follows from the observation that the phase diagram of water is twodimensional. Nevertheless, despite the mild assumptions of Z 2 symmetry and two relevant scalars, we find a striking result: the dimensions (∆ σ , ∆ ǫ ) are almost uniquely fixed! The allowed region is a tiny sliver around ∆ σ = 0.51820 (14) and ∆ ǫ = 1.4127 (11), in agreement with the most precise Monte-Carlo simulations [35], and also the c-minimization conjecture of [22]. Indeed, our results give strong support for c-minimization -support that can be strengthened with further numerical work.
It is plausible that with arbitrary computational power, the constraints we study are strong enough to uniquely determine the spectrum of the 3D Ising CFT. In particular, our results support the conjecture that the 3D Ising CFT is the only Z 2 -symmetric 3D CFT with exactly two relevant operators, giving a higher-dimensional example of critical universality [36,37].
Because our study includes the correlator σσǫǫ , we can additionally learn about Z 2 -odd operators (appearing in the σ × ǫ OPE), which were not accessible in previous bootstrap studies. In section 5 we also compute an upper bound ∆ σ ′ ≤ 5.41 (1), where σ ′ is the secondlowest dimension Z 2 -odd scalar. A precise determination of the complete Z 2 -odd spectrum of the 3D Ising model is a fascinating problem that we leave for future work. We describe other important future directions in section 6.

Bootstrapping Mixed Correlators and Semidefinite Programming 2.1 What is New About Mixed Correlators?
It was shown in [2] that the bootstrap constraints for a four-point function of identical scalars φφφφ can be transformed into a system of linear inequalities. Studying solutions to these inequalities leads to bounds on CFT data. By contrast, the bootstrap constraints for mixed correlators φ 1 φ 2 φ 3 φ 4 cannot be written in terms of linear inequalities -rather the mixed correlator problem is intrinsically quadratic. 2 In this section, we describe how this quadratic problem arises in a simple case, and how the procedure of [2] can be modified to solve it.

Review of the Bootstrap Argument for Identical Scalars
Let us first recall the original bootstrap argument of [2] for a four-point function of identical real scalars φ(x 1 )φ(x 2 )φ(x 3 )φ(x 4 ) . Using the OPE, we can write the four-point function as a sum over conformal blocks Here, u = are conformal cross-ratios, O runs over real primary operators appearing in the φ × φ OPE, ∆ = dim O, and ℓ = spin O. The OPE coefficients λ φφO are real by unitarity, implying that their squares λ 2 φφO are positive. The four-point function should be independent of how we pair the operators to perform the OPE. Specifically, swapping 1 ↔ 3, we find the crossing equation Grouping terms that multiply λ 2 φφO , we obtain a sum rule with positive coefficients O λ 2 φφO F ∆,ℓ (u, v) = 0, (2.3) Positivity of the coefficients in (2.3) is the key property that leads to universal bounds on CFT data, without having to know the precise details of the operators O entering the OPE. Let us first assume the dimensions ∆ and spins ℓ lie in some specified range. For example, we might assume that all scalars have dimension larger than some ∆ 0 . Consider linear functionals α acting on functions of u, v, and suppose there exists an α satisfying the conditions α(F ∆,ℓ ) ≥ 0 for all ∆, ℓ in the spectrum, and α(F 0,0 ) = 1. (2.5) If such an α exists, the sum rule (2.3) cannot be satisfied with any choice of operators O, and the hypothetical CFT is ruled out.
We can always express α in terms of a basis of functionals, for example derivatives around the crossing-symmetric point, where z, z are defined by u = zz, v = (1 − z)(1 − z). Then, (2.5) becomes a set of linear inequalities (and one affine equality) for the coefficients a mn .

Applying the Bootstrap Argument to a Mixed Correlator
Let us try to apply this procedure to where g is a conformal block for scalars with possibly unequal dimensions, ∆ ij ≡ ∆ i − ∆ j , and λ ijO denotes the OPE coefficient of O ∈ φ i × φ j . We will often abbreviate g 0,0 ∆,ℓ (u, v) as g ∆,ℓ (u, v). The left-hand side of (2.7) has manifestly positive coefficients |λ 12O | 2 . However, on the right-hand side there is no a-priori relation between λ 22O and λ 11O , so their product can have either sign. Consequently, we cannot simply apply linear functionals to both sides and derive conclusions about the allowed spectrum. We must modify the bootstrap logic above.
To obtain some kind of positivity condition, we can combine the crossing equation for ∆,ℓ are combinations of conformal blocks analogous to (2.4), and we have suppressed their u, v dependence for brevity. The quantity multiplying λ 11O and λ 22O above is a 2 × 2 matrix of functions of u, v. Formally, it is an element of R 2×2 ⊗ F , where F is the space of functions of u, v in the region where both conformal block expansions converge. 3 The right-hand side of (2.7) contributes to off-diagonal elements of this matrix.
We can now consider linear functionals α acting on functions of u, v α : F → R, (2.9) such that ∈ R 2×2 is a positive semidefinite matrix, (2.10) and α is additionally positive acting on each λ 2 12O term in (2.8). Applying α to both sides, we again have termwise positivity and the bootstrap logic can proceed.
An optimization problem that includes positive semidefiniteness constraints of the form (2.10) is a semidefinite program, as opposed to the linear program (2.5) that appears in the case of identical operators. 4 Semidefinite programming has appeared in the conformal bootstrap before: it was applied to 4D CFTs in [9], and later extended to arbitrary spacetime dimensions in [15]. However, its appearance here is qualitatively different. In [9,15], semidefinite programming was a useful trick for efficiently encoding the infinite number of constraints α(F ∆,ℓ ) ≥ 0 (one for each ∆ and ℓ). This trick is not strictly necessary, and alternative methods have also been successful, for example the discretization of ∆ in [2] and the modified simplex algorithm in [22]. By contrast, the appearance of semidefinite programming in (2.10) is unavoidable, stemming from the intrinsically quadratic nature of the crossing constraints. 5 In this work, we will combine the semidefinite programming trick of [9,15] with the novel appearance of semidefinite programming in (2.10). Thus, semidefinite programming appears in two ways: one optional, the other obligatory. An alternative approach that may be fruitful would be to adapt the modified simplex algorithm of [22] to work with semidefiniteness constraints. 6 3 The convergence region F includes a finite open neighborhood of the point z = z = 1 2 [38]. 4 Because of the infinite number of constraints (one for each ∆, ℓ), one technically has a semi-infinite program in the identical operators case. We will not bother with this distinction. We do not know the correct terminology for a semidefinite program with an infinite number of constraints, as in the case of mixed operators. 5 It is possible to approximate the semidefiniteness constraint with a finite number of linear constraints by approximating the cone of semidefinite matrices as a polytope, as explored in [39]. 6 Semidefinite programs enjoy a duality similar to the linear programming duality underlying the primal simplex method in [22].

Spin and Global Symmetry Representations
Before describing our construction in detail, let us note that the appearance of semidefinite programming is generic in the conformal bootstrap, and previously considered problems are special cases where it can be avoided. Semidefinite programming appears whenever a conformal block expansion has coefficients with indefinite sign.
An important example is a four-point function of operators with spin, where multiple structures can appear in the OPE. For example, let J µ (x) be a conserved current in a 3D CFT. A primary operator O µ 1 ···µ ℓ with even spin ℓ can appear in the J × J OPE with two different parity-even tensor structures [40][41][42], The tensors t 1 (x) and t 2 (x) are fixed by conservation, symmetry under exchanging the J's, and conformal invariance. Each tensor structure has an independent OPE coefficient λ 1 JJO , λ 2 JJO , and thus the conformal block expansion of J µ J ν J ρ J σ contains terms proportional to λ 1 JJO λ 2 JJO (which can have either sign). Another example is a four-point function of operators in representations r i of a global symmetry group G, where the tensor product r i ⊗ r j contains irreducible representations with nontrivial multiplicity. The fact that bootstrap conditions become semidefiniteness constraints in this case was first noticed in [7].
where φ r is a scalar operator transforming in the representation r, and φ † r transforms in the dual representation r. Suppose O s,ℓ ∈ φ r ×φ r is a spin-ℓ operator transforming in the representation s of G. The three-point function φ r φ r O † s,ℓ must be proportional to an invariant tensor t O of G. Specifically, where Sym 2 and ∧ 2 denote symmetric and antisymmetric tensor squares, and (·) G denotes the G-invariant subspace. The space of such invariant tensors may be multidimensional: its dimension counts the multiplicity of s in the decomposition of Sym 2 r and ∧ 2 r into irreducibles. If so, we can expand t O in a basis of invariant tensors t i , each with an independent coefficient The conformal block expansion of our four-point function in the (12) → (34) channel can then contain products λ i φφO λ j * φφO , which are not necessarily positive. Previous bootstrap studies of CFTs with global symmetries have focused on small representations: either G = SO(n) with r the vector representation, or G = SU(n) with r the fundamental representation. In each of these cases, the spaces Sym 2 r, ∧ 2 r, and r ⊗ r (in the other channel) decompose into irreducibles with multiplicity at most 1.
However, it is easy to find examples with higher multiplicities. For instance, take G = SU(n), n ≥ 3, and let r to be the largest irreducible representation in Sym 2 Ad G , of dimension 1 4 n 2 (n − 1)(n + 3). In this case, Sym 2 r = 2 r ⊕ . . . (2.14) so that there are two independent OPE coefficients for an operator O r,ℓ with even spin. The three-point structures for this example are written explicitly in appendix A.

General Semidefinite Programs for the Bootstrap
In the previous subsection, we saw examples of different ways that semidefinite programming can arise in the bootstrap. Now let us generalize these examples and show how the generic statement of crossing symmetry and unitarity can be phrased as a semidefinite program. The discussion here is somewhat abstract. In section 3, we will specialize to the case of interest for the remainder of this work.
Consider a CFT whose symmetry group H is a product of the conformal group and global symmetry groups. 7 Primary operators O i transform in unitary irreducible representations R i of H. Let H 0 be the isotropy subgroup of H (the group of symmetries that fix a spacetime point), which is generated by Lorentz transformations, dilatations, special conformal transformations, and global symmetry transformations. The representations R i are induced from finite dimensional representations R i0 of H 0 . Hence, each O i carries an index a i for R i0 . For example, an uncharged spin-1 operator J µ has a i = µ, a Lorentz index.
The four-point functions of the theory are given by For each set of representations R i , R j , R k , R l , the four-point function can be expanded in a finite basis of four-point structures I The number of four-point structures depends on the representations R i , R j , R k , R l . For example, when all the operators are scalars there is a single four-point structure. Its definition is ambiguous up to multiplication by conformal cross ratios u, v. We choose Crossing symmetry is the statement that (for bosonic operators) swapping (i, a i , x 1 ) ↔ (k, a k , x 3 ) leaves the four-point function unchanged, 8 G a i a j a k a l ijkl For instance, for scalar operators, using (2.17) we have We should think of G f ijkl (u, v) as a vector with indices given by ijkl, f, u, and v. The matrix S acts on the f and u, v-indices. For convenience, let us additionally define the operator T : (2.21) As usual in the bootstrap, to constrain solutions to crossing symmetry, we can look at linear functionals acting on the crossing equation. A linear functional acting on G has the form where each α ijkl f acts on functions of two variables, as in (2.6). The dual of the crossing equation is the statement that for all α.
The final ingredient is the conformal block expansion for the functions G f ijkl . For an operator O p appearing in the OPE O i ×O j , there are in general several three-point structures t m that can appear, and each structure has an associated OPE coefficient λ (m) , Consequently, conformal blocks in the ij → kl channel are labeled by pairs of three-point structures (m, n), The conformal blocks g depend on the representations (dimensions, spins, and global symmetry charges) of the external and exchanged operator, together with the threepoint structures (m, n) and four-point structure f .
Plugging the conformal block expansion (2.25) into the dualized crossing equation (2.23) gives ) . (2.26) All repeated indices are summed, but we have chosen to indicate some of the sums explicitly. Let us think of λ ijp as a vector λ p with indices ij, m. For each p, the quantity in parentheses is then a matrix A p with left indices ij, m and right indices kl, n, and the above equation can be written From (2.27), it is clear that the crossing equation can be studied using semidefinite programming. Specifically, we can make an assumption about the CFT spectrum, and then search for α such that 0, for all p which can appear in the spectrum. (2.28) If such an α exists, the assumed spectrum is ruled out. As usual in the bootstrap, the advantage of studying this dual formulation of crossing symmetry is that one can make progress by restricting to a finite-dimensional space of α's. By contrast, truncating the operator spectrum itself can introduce approximations that are difficult to control.

SDP 2
The semidefinite programs described in previous sections have an infinite number of positivity constraints (2.10) -one for each ∆ and ℓ (and more generally for each conformal representation) that can appear in the spectrum. To encode these constraints on a computer, we write them in terms of a finite amount of data using a trick from [9,15].
Let us briefly review this trick in the case of identical operators. The positivity constraints (2.5) are linear inequalities of the form where ∆ min ℓ is an ℓ-dependent lower bound on the dimension (e.g., the unitarity bound). We take α to be of the form (2.6): a sum of derivatives with respect to z, z around the crossing symmetric point z = z = 1/2.
The trick proceeds by first rewriting our inequalities as positivity of polynomials. This is possible because there exists a systematic approximation The left hand side is now a polynomial in x. The degree of this polynomial grows logarithmically with the accuracy of the approximation (2.30).
A classic theorem [43] states that a polynomial P (x) is positive on the positive real axis, if and only if where a(x) and b(x) are sums of squares of polynomials. That is, a(x) = t p t (x) 2 and similarly for b(x).
It is easy to show that any degree-2d sum of squares is of the form Tr(AQ d (x)), where A 0 is a positive semidefinite matrix. Consequently, (2.33) is equivalent to (2.34) Here and " " means "positive semidefinite." Using (2.34), we can write the constraints (2.31) in terms of positive semidefinite matrices A ℓ , B ℓ with linear relations between the elements of A ℓ , B ℓ and the coefficients a mn . This efficiently encodes an infinite number of constraints in terms of finite matrices. Now let us consider the case of distinct operators, with constraints of the form (2.10), Here, F ∆,ℓ are combinations of conformal blocks, and α is again a sum of derivatives with respect to z, z. Similarly to (2.30), a systematic positive-times-polynomial approximation for each entry of the above matrix exists, Crucially, the positive function χ ℓ (∆) is independent of the matrix indices i, j. Dividing by Analogously to the case of positive polynomials, an N × N matrix polynomial that is positive on the positive real axis, such as the one in (2.37), can be written as [44] a(x) + xb(x), where a(x) and b(x) are sums of squares of matrix polynomials, i.e. a(x) = t p T t (x)p t (x), and similarly for b(x). 10 Let the highest degree of the polynomials p t be d 1 , and write the matrix polynomials p t (x) as where p t;ρ are M × N matrices. Any matrix polynomial a(x) can be written as acts only on the first tensor factor R d 1 +1 and the trace is over this factor, leaving an N × N matrix that depends on x. In terms of components, the elements of a are given by If a(x) is a sum of matrix squares, then A is positive semidefinite on the space R d 1 +1 ⊗R N . To see this, note that for an arbitrary N( The converse also holds: if A is positive semidefinite, then a(x) is a sum of matrix squares, as can be seen by writing A as a sum of outer-products. Therefore, the matrix polynomials of (2.37) can be written in a way analogous to (2.34): where A ij ρσ and B ij ρσ are respectively N(d 1 + 1) ×N(d 1 + 1) and N(d 2 + 1) ×N(d 2 + 1) positive semidefinite matrices. 11 Once again, the infinite number of (matrix) constraints is encoded in terms of finite matrices.

Specializing to Scalars
The analysis of section 2.3 is completely general and abstract. To simplify the discussion, in this section we will focus on four-point functions of scalar operators. Here we do not assume any additional global symmetries. Let φ i = (φ 1 , φ 2 , φ 3 , . . .) be a collection of scalar primary fields with scaling dimensions ∆ i . We consider the correlation function of four scalars, The expression (2.16) for the four-point function simplifies because there are no representation indices a in (3.1) and there is only one four-point structure (2.17) for four scalars. Therefore, G ijkl can be written as The correlation function must be invariant under the exchange (1, i) ↔ (3, k), which gives the crossing equation Following (2.25), the function G ijkl can be written in terms of conformal blocks g as 11 We thank João Penedones for discussions of this idea.
where compared to (2.25) there is no sum over different three-point structures, as only one structure appears in three-point functions containing two scalar operators. The conformal blocks also depend on the representations of the external operators φ i , but since they are scalars only the dependence on the dimensions ∆ i remains. Moreover, the conformal blocks depend only on the differences ∆ ij and ∆ kl , as shown for instance in [45,46] and reviewed in section 4. For real scalar external operators, the OPE coefficients λ ijO are real [2]. The crossing equation is then given by The dual form of the crossing equation is given by (2.26), which for scalar operators becomes and can be analyzed using semidefinite programming, as explained in section 2.4.
We now introduce notation that follows more closely the notation used in the analysis of the single correlator crossing equation. Let us define functions which are respectively symmetric and antisymmetric under the exchange u ↔ v. They are also invariant under the simultaneous exchange of i ↔ k and j ↔ l. In terms of these functions, the crossing equation becomes If all operators are equal, the upper sign case reduces to the single correlator crossing equation (2.3).

Simplest System with a Z 2 Symmetry
Before switching to the dual form of the crossing equation, we will further simplify the system of crossing equations under consideration by assuming that the system has a Z 2 symmetry. Under this symmetry, all of the operators can be classified as even or odd. Let σ and ǫ be the lowest dimension Z 2 -odd and Z 2 -even scalars, respectively. The OPE structure of these operators can be written schematically as Here O + runs over Z 2 -even operators of even spin and O − runs over Z 2 -odd operators of any spin. An important example of a system described by OPEs (3.9) is the critical Ising model, where σ and ǫ can be thought of as the spin and energy density operators.
We now write the crossing equations (3.5) for four-point functions containing σ and ǫ. Due to the Z 2 symmetry, some of the OPE coefficients vanish and we end up with five independent constraints: where in the last equation we used the identity λ σǫO = (−1) ℓ λ ǫσO . We can write this in vector notation as where V −,∆,ℓ is a 5-vector and V +,∆,ℓ is a 5-vector of 2 × 2 matrices: Let α = (α 1 , . . . , α 5 ) be a 5-vector of functionals, α i : F → R. Acting on (3.11) with this functional gives the dual form of crossing equation (3.14) The sum over operators in (3.11) must contain the unit operator. For convenience we will write the unit operator contribution separately. We also assume that the operators are normalized so that λ σσ1 = λ ǫǫ1 = 1. The crossing equation can then be written as

Bounds from Semidefinite Programming
Now (3.15) is in the right form for a semidefinite programming analysis as described at the end of section 2.3. The logic is to make an assumption about the CFT spectrum and then search for a functional α such that If we manage to find such a functional, then crossing symmetry can not be satisfied and we conclude that the initial assumption about the spectrum is wrong. Note that (3.16) represents a sufficient condition for a functional to exclude the spectrum, but not a necessary one. In particular, it does not take into account the symmetry of the OPE coefficients, λ σσǫ = λ σǫσ . Using this identity we could weaken the conditions on a functional that would still violate the crossing symmetry constraints (strengthening the resulting bounds).
The space of all functionals is infinite-dimensional, so for implementation on a computer we must select a finite subspace of functionals as well as a convenient set of basis vectors. We follow the traditional approach, which has given good results in previous studies, and use linear combinations of z and z derivatives at the point z = z = 1/2 as the basis for functionals. The functional α i is given in that basis by where we limit the number of derivatives with the parameter n max by requiring n + m ≤ 2n max − 1. 12 The sum in (3.17) only contains terms with m ≥ n because the functional acts on functions of u, v, which are symmetric under the exchange z ↔ z. Moreover, the functional α 5 acts on F + , which is symmetric under the exchange u ↔ v, so only the terms with m + n even are non-zero. Other components α i act on functions F − antisymmetric under u ↔ v, so for them only the terms with n + m odd are non-vanishing. The total number of independent components of the functional α is then N = 5n max (n max + 1)/2. This is the dimension of the vector space in which we search for a functional α with the properties (3.16).
To turn the constraints (3.16) into a semidefinite program suitable for numerical analysis on a computer, we employ the method described in section 2.4. We use the rational approximation of conformal blocks described explicitly in section 4 to approximate the conformal blocks and their derivatives as where p ∆ 12 ,∆ 34 ;mn ℓ (∆) are polynomials in ∆ and χ ℓ (∆) are functions that are positive for all values of ∆ in the CFT spectrum. Then we can similarly write where P ijkl;mn ±,ℓ (∆) are polynomials related to p ∆ ij ,∆ kl ;mn ℓ . Using the last expression, we can rewrite the conditions (3.16) in the following form: 0, for all Z 2 -even operators with even spin, Y ℓ (∆) ≥ 0, for all Z 2 -odd operators in the spectrum. (3.20) Here Y ℓ (∆) are polynomials and Z ℓ (∆) are matrix polynomials in ∆ defined as  The constraint Y ℓ (∆) ≥ 0 can be written in terms of positive semidefinite matrices as where A ℓ , B ℓ 0 are positive semidefinite, and d 1 , d 2 are chosen appropriately for the degree of Y ℓ . Similarly, the constraint Z ℓ (∆) 0 can be written in terms of positive semidefinite matrices as where C ℓ , D ℓ 0 are positive semidefinite acting on R d 1,2 +1 ⊗ R 2 , and d 1 , d 2 are chosen appropriately for the degree of Z ℓ . Matching terms of equal degree in x on both sides, (3.22) and (3.23) become linear equations relating the variables a i mn , A ℓ , B ℓ , C ℓ , D ℓ . Written in terms of these variables, our optimization problem is now in a form that can be fed to a semidefinite program solver. We give details of our implementation in the solver SDPA-GMP [47,48] in appendix B.

Additional Constraints
There are a few additional constraints on systems of multiple correlators that we have not yet incorporated into our numerical analysis. One is that the coefficient of the stress tensor conformal block should be consistent with Ward identities. In the OPE φ × φ, Ward identities imply that the stress tensor appears with coefficient proportional to ∆ φ . Thus, in the conformal block decomposition for φ 1 φ 1 φ 2 φ 2 in the φ 1 φ 1 → φ 2 φ 2 channel, the stress tensor block appears with a coefficient of λ 11T λ 22T ∝ ∆ 1 ∆ 2 c (up to a theory-independent pre-factor that depends on our definition of the conformal block), where c ∝ T µν T ρσ is the coefficient of the stress tensor two-point function. 13 In the 3D Ising model, this implies that the coefficients of the stress tensor block in σσσσ , σσǫǫ , ǫǫǫǫ must appear in the fixed ratios ∆ 2 σ : ∆ σ ∆ ǫ : ∆ 2 ǫ . We might add this additional information to obtain stronger bounds. However, this condition is useless without an additional assumption of a gap in the spin-2 spectrum, ∆ T ′ ≥ 3 + δ with δ > 0, where T ′ is the second-lowest spin-2 operator in the theory. In the absence of a gap, spin-2 operators with dimension 3 + ε with ε ≪ 1 can mimic the contribution of the stress tensor, obliterating any information about ratios of stress-tensor coefficients. Concretely, the condition Z 2 (3 + ε) 0 in (3.20) for all ε > 0 also implies Z 2 (3) 0, which is strictly stronger than necessary given the Ward identities.
In the 3D Ising model, there is indeed a large gap in the spin-2 spectrum, ∆ T ′ 5.5 [22]. We have not incorporated the existence of this gap and the accompanying Ward identity constraints in this work. It will be very interesting to do so in the future.
Another constraint on the 3D Ising model is that there is precisely one operator ǫ with dimension ∆ ǫ . 14 The condition Z 0 (∆ ǫ ) 0 in (3.20) actually allows for many operators of dimension ∆ ǫ with different OPE coefficients. In particular, it assumes only that the sum is a positive semidefinite matrix (which is a consequence of unitarity). However, because there is exactly one choice for O above, the sum is not a generic positive semidefinite matrix: it has rank one. Let us suppose the vector (λ σσǫ , λ ǫǫǫ ) is proportional to e θ ≡ (cos θ, sin θ) for some angle θ ∈ [0, π). We can then replace the condition Z 0 (∆ ǫ ) 0 with the weaker condition Running our semidefinite program subject to (3.25) will yield some allowed region A θ in the space of CFT data. Since we do not know the actual value of θ in the 3D Ising model, we 13 c is commonly referred to as the "central charge," although it has no (known) relation to symmetry algebras in greater than 2 spacetime dimensions. 14 We thank Slava Rychkov for explaining how to exploit this constraint. must then scan over θ, computing the final allowed region The region A * could in principle be smaller than what one gets by imposing the naïve condition Z 0 (∆ ǫ ) 0. The idea of scanning over θ to exploit the fact that (3.24) has rank one was initially explored in [39]. Unfortunately, performing this scan is infeasible with our current methods, given the time it takes to compute each allowed region A θ . These constraints will be important to study in the future.
Another constraint is the symmetry of three-point coefficients, in particular the relation λ σǫσ = λ σσǫ . This constraint is straightforward to impose within our formalism, and we are currently exploring its consequences.

Rational Representations for Conformal Blocks
In order to study the semidefinite program described in the previous sections, we require systematic approximations for the derivatives of the functions F ijkl ±,∆,ℓ in terms of positive functions times polynomials in ∆. Such approximations directly follow from a representation for the conformal blocks as a sum over poles in ∆. Representations of this type were first developed for 2D (Virasoro) conformal blocks by Alyosha Zamolodchikov [33,34], and a generalization for global conformal blocks with identical external scalars in general spacetime dimension D was developed in [15]. 15 In this section we generalize the formula obtained in [15] to non-identical external scalars.
Poles in the conformal block occur at special (non-unitary) dimensions ∆ = ∆ * where some descendant P n |O of the state created by the primary operator O becomes null. This null state and all of its (also null) descendants together form a conformal sub-representation, and hence the residue of the pole is proportional to a conformal block: 16 Since poles in ∆ determine g ∆ 12 ,∆ 34 ∆,ℓ up to a function that is analytic on the entire complex plane, we can write

2)
15 Other recent studies of global conformal blocks can be found in [11,41,42,45,46,[49][50][51][52][53][54][55], with connections to Mellin amplitudes in [56][57][58][59][60][61][62][63][64][65]. Older work includes [66][67][68][69][70]. Superconformal extensions have been studied in [5,21,[71][72][73][74][75][76][77][78]. 16 At this stage it is not obvious that all such poles must be simple poles. Indeed, when the spacetime dimension D = 2n is an even integer, double poles occur, while outside of even dimensions only simple poles occur. Our formulas assume D = 2n, but reproduce the correct even-D conformal blocks in the limit D → 2n.  whereg ∆ 12 ,∆ 34 ℓ (∆, r, η) is an entire function of ∆ and we describe the conformal cross ratios using radial coordinates [52]. In Euclidean signature, where z = z * , these are defined by The block g ∆ 12 ,∆ 34 ∆,ℓ (r, η) has an essential singularity of the form r ∆ as ∆ → ∞. Stripping this off, we have We see that the entire functionh ∆ 12 ,∆ 34 ℓ (r, η) = lim ∆→∞ h ∆ 12 ,∆ 34 ∆,ℓ (r, η) no longer depends on ∆, since there are no singularities as ∆ → ∞. Now, the functionh ∆ 12 ,∆ 34 ℓ can be easily computed by solving the conformal Casimir equation [46] to leading order in ∆, giving the result 17 where ν = D−2 2 and C ν ℓ (η) is a Gegenbauer polynomial. The locations of the poles ∆ i in (4.2) are the same as in the case of equal external dimensions (though in that case some have vanishing coefficients) since they only depend on the representation theory of the exchanged operator. In both cases we find three sequences of poles, reproduced in table 1. However, the coefficients c ∆ 12 ,∆ 34 i depend on the external dimensions and can be found by solving the conformal Casimir equation. In practice we do this order by order in the r expansion, following the procedure described in [52]. We compute coefficients up to high order in the r-expansion, guess a formula, and check the formula to even higher orders. The resulting coefficients are given by where (a) n = Γ(a + n)/Γ(a) denotes the Pochhammer symbol. It should be possible to analytically derive these coefficients using conformal representation theory. A derivation could shed light on their generalization to other conformal blocks (e.g., for operators with spin) and superconformal blocks.
Using the recursion relation (4.5), it is straightforward to compute derivatives of the conformal blocks around the crossing-symmetric point r = r * = 3 − 2 √ 2 ≈ 0.17, η = 1. 18 These have the form where q ∆ 12 ,∆ 34 ;mn ℓ (∆) is a polynomial in ∆ and a ∆ 12 ,∆ 34 ;mn ℓi are numerical coefficients. Poles corresponding to larger values of n i in table 1 are suppressed by higher powers of r * ≈ 0.17. Thus, we can get a good approximation by truncating to a finite number of poles n i ≤ ν max , with the result where p ∆ 12 ,∆ 34 ;mn ℓ (∆) is a polynomial obtained by combining the poles in the partial fraction expansion (4.8). 19 This has the form (3.18), with χ ℓ (∆) = r ∆ * i (∆ − ∆ i ) −1 . The accuracy of this approximation depends on ν max . Bounds involving more derivatives (higher n max ) require more precise expressions for the blocks, and consequently higher ν max . We have verified that the bounds computed in this work are essentially unchanged if ν max is increased.

General Bound on ∆ ǫ
Previous numerical bootstrap studies considered only a single four-point function σσσσ . One of the results of the single correlator bootstrap is a rigorous upper bound on the dimension of the lowest-dimension scalar ǫ appearing in the σ × σ OPE. The bound on ∆ ǫ as a function of ∆ σ was obtained for D = 3 in [11]; we reproduce it here and plot it in several figures, as explained below.
In the case of multiple correlators, both σ and ǫ appear as external operators, resulting in the system of equations (3.11). It is clear that the bound resulting from this system will be at least as strong as the single-correlator bound because we can set α 2,3,4,5 = 0 in (3.15) to reduce it to the single-correlator problem. A priori, it is possible for the complete system (3.11) to give an even stronger bound on ∆ ǫ . However, after proceeding with the computations described in section 3.3, assuming nothing else about the spectrum except that σ and ǫ are respectively the lowest dimension Z 2 -odd and -even scalars, 20 we find that the multi-correlator and single-correlator bounds agree (at n max = 6).
While it is not obvious what the allowed region should be, it is clear that Z 2 -symmetry combined with multiple correlator constraints should not fix the dimensions (∆ σ , ∆ ǫ ) without additional assumptions. For example, the O(N) vector models admit a Z 2 symmetry, so they should lie in the allowed region, but below the 3D Ising model. 21

Bound on ∆ σ ′
A new feature of the multiple correlator bootstrap for the Ising model is the access to the Z 2 -odd spectrum. In particular, this allows us to place an upper bound on the dimension of the second Z 2 -odd scalar σ ′ . We assume that all Z 2 -odd scalars (other than σ, σ ′ ) have dimensions greater than ∆ σ ′ and try to find a contradiction with the crossing equation (3.11). In the notation of section 3.3, we must find a functional satisfying the constraints: 20 Additionally, we always include the usual unitarity bounds on operator dimensions. 21 To apply our bounds to the O(N ) models, we may take σ = φ 1 and ǫ = S 11 , where φ i is a vector under O(N ) and S ij is the lowest-dimension traceless symmetric tensor of O(N ) in φ i × φ j . Of course, much stronger bounds can be obtained by using the full information of O(N ) symmetry, as in [15].
where ∆ min ℓ is the unitarity bound for a spin-ℓ operator. For given ∆ σ and ∆ ǫ , we find the minimal value of ∆ σ ′ for which the spectrum is excluded. This value, depending on both ∆ σ and ∆ ǫ , is an upper bound on the dimension of σ ′ .
Instead of making a three-dimensional plot of ∆ σ ′ vs. (∆ σ , ∆ ǫ ), we choose a curve in the (∆ σ , ∆ ǫ ) plane and plot the ∆ σ ′ bound along it. In particular, we choose (∆ σ , ∆ ǫ ) to lie on the (single-correlator) upper bound on ∆ ǫ computed at n max = 10 -the black dotted line in figures 3 and 4. This choice of curve is somewhat arbitrary. An advantage is that it should pass near the 3D Ising point. Our primary goal is to get a general picture of the constraints on ∆ σ ′ .
The bound on ∆ σ ′ as a function of ∆ σ (with ∆ ǫ following the n max = 10 single correlator bound) is shown in figure 1 at n max = 6, corresponding to a linear functional with N = 105 components. 22 It has an almost rectangular peak centered near the Ising value of ∆ σ . The sides of the peak are close to vertical and the top is relatively flat. The width of the peak is partially an artifact of our choice of (∆ σ , ∆ ǫ ) curve, as will be clear in the next subsection.
This peak provides another example of interesting behavior in bootstrap bounds near the Ising point. Away from the Ising point, for ∆ σ < 0.517 or ∆ σ > 0.52, the bound on ∆ σ ′ is quite strong, implying ∆ σ ′ < 3 for the range of ∆ σ plotted in figure 1. Just at the Ising point, the bound ∆ σ ′ 6.5 is relatively weak compared to the expected value of ∆ σ ′ = 3 + ω A in the 3D Ising CFT, which has been estimated to be ∆ σ ′ 4.5 by resumming the ǫ-expansion [79], ∆ σ ′ ≈ 5.4 using the scaling field approach [80], and ∆ σ ′ ≈ 4.7 − 5.2 using exact RGE methods [81].
Enlarging the search space of functionals by icreasing n max , the peak gets narrower and smaller. At n max = 10, we computed a bound on ∆ σ ′ to precision 0.01 at the three points (∆ σ , ∆ ǫ ) = (0.5181, 1.41206), (0.51815, 1.41267), (0.5182, 1.41312) near the expected values in the 3D Ising model. At each of these points we find ∆ σ ′ ≤ 5.41(1), closer to the estimates of other methods, but perhaps still somewhat larger. The difference in the n max = 6, 10 results indicates that this bound has yet to converge to its optimal value and could likely be improved with further numerical work.

Bounds on (∆ σ , ∆ ǫ ) with Gaps in the Operator Spectrum
The ∆ σ ′ bound in figure 1 indicates that away from the Ising point the spectrum must contain a Z 2 -odd scalar of dimension ∆ σ ′ < 3. Conversely, imposing ∆ σ ′ ≥ 3 will exclude values of ∆ σ sufficiently far from the Ising model. The assumption ∆ σ ′ ≥ 3 has physical meaning: it implies that there is only one relevant Z 2 -odd operator, σ. This assumption is known to hold for the critical Ising model, where the only relevant operators are σ and ǫ. Using this as input, we can additionally assume a gap ∆ ǫ ≥ 3 in the Z 2 -even spectrum to obtain even stronger bounds on scaling dimensions of the Ising model.
With the assumption of a gap in the Z 2 -odd spectrum, we find a strong constraint on the values of ∆ σ and ∆ ǫ . The allowed region in the (∆ σ , ∆ ǫ ) plane is shown shaded light  3 and 4). The sharp spike in the bound occurs when the values of ∆ σ , ∆ ǫ lie within the allowed region in figure 3, with no assumption about Z 2 -even gaps (medium-blue shaded region). Only when ∆ σ , ∆ ǫ take these values, is it possible to have ∆ σ ′ ≥ 3, and indeed the upper bound on ∆ σ ′ becomes extremely weak. This bound is computed at n max = 6, ν max = 8. We expect that as n max increases, it becomes more sharply peaked, while the top moves down to the correct value of ∆ σ ′ in the 3D Ising model. At n max = 10, ν max = 14 we have computed the stronger bound ∆ σ ′ ≤ 5.41(1) (the dashed line) at the points (∆ σ , ∆ ǫ ) = (0.5181, 1.41206), (0.51815, 1.41267), (0.5182, 1.41312). blue in figure 2 for n max = 6. As expected from the ∆ σ ′ bound, it consists of a small closed region around the Ising point and another big region at ∆ σ 0.54. The dashed line is the single correlator bound. Note that since we expect ∆ σ ′ ≈ 4.5 in the Ising model, we could have assumed a bigger gap in the Z 2 -odd sector to get an even smaller allowed region in the (∆ σ , ∆ ǫ ) plane. However, due to the fact that the sides of the peak in the ∆ σ ′ bound in figure 1 are almost vertical, we expect that the allowed region around the Ising point is not significantly affected by the exact value of the Z 2 -odd gap, as long as it is greater than 3.
The effect of the gap in the Z 2 -even sector is shown in figure 3. Without the Z 2 -odd gap, we get exactly the same bound that was obtained for the single correlator bootstrap in [11]. The allowed region for that case is shaded light blue in figure 3. Assuming gaps in both the Z 2 -odd and -even parts of the spectrum, we find the allowed region around the Ising point (shaded dark blue) to be of similar shape, but somewhat smaller size than the allowed region when assuming only the Z 2 -odd gap (shaded medium blue). A zoomed version of this region is shown in figure 4. : Allowed region of (∆ σ , ∆ ǫ ) in a Z 2 -symmetric CFT 3 where ∆ σ ′ ≥ 3 (only one Z 2 -odd scalar is relevant). This bound uses crossing symmetry and unitarity for σσσσ , σσǫǫ , and ǫǫǫǫ , with n max = 6 (105-dimensional functional), ν max = 8. The 3D Ising point is indicated with black crosshairs. The gap in the Z 2 -odd sector is responsible for creating a small closed region around the Ising point.
The allowed region around the Ising point shrinks further when we increase the value of n max . Finding the allowed region at n max = 10 (N = 275) is computationally intensive, so we tested only the grid of 700 points shown in figure 5. The disallowed points in the figure were excluded by assuming both ∆ σ ′ ≥ 3 and ∆ ǫ ′ ≥ 3. On the same plot, we also show the n max = 14 single-correlator bound on ∆ ǫ computed in [22] using a very different optimization algorithm. The final allowed region is the intersection of the region below the n max = 14 curve and the region indicated by our allowed multiple correlator points.
Since the point corresponding to the 3D Ising model must lie somewhere in the allowed region, we can think of the allowed region as a rigorous prediction of the Ising model dimensions, giving ∆ σ = 1/2 + η/2 = 0.51820 (14) and ∆ ǫ = 3 − 1/ν = 1.4127 (11). In figure 6 we compare our rigorous bound with the best-to-date predictions using Monte Carlo simulations [35] and the c-minimization conjecture [22]. Although our result has uncertainties greater than c-minimization by a factor of ∼10 and Monte-Carlo determinations by a factor of ∼3, they still determine ∆ σ and ∆ ǫ with 0.03% and 0.08% relative uncertainty, respectively. Increasing n max further could potentially lead to even better determinations of ∆ σ and ∆ ǫ . Indeed, the single correlator bound at n max = 14 passing through the allowed region in figure 5 indicates that the n max = 10 allowed region is not yet optimal. At this point, it is not even clear whether continually increasing n max might lead to a finite allowed allowed region with various gaps in ∆ ǫ ′ , ∆ σ ′ (n max = 6) ∆ σ Single Correlator ∆ ǫ ′ ≥ 3 Figure 3: Allowed regions in a Z 2 -symmetric CFT 3 , assuming various gaps in the scalar spectrum. The dashed line is an upper bound on ∆ ǫ using crossing symmetry and unitarity of σσσσ , with no assumptions about gaps, at n max = 6. The black dotted line is the same bound with n max = 10. The light blue shaded region assumes a gap ∆ ǫ ′ ≥ 3 in the Z 2 -even sector. The medium blue shaded region assumes a gap ∆ σ ′ ≥ 3 in the Z 2 -odd sector, and uses crossing symmetry for the system of correlators σσσσ , σσǫǫ , ǫǫǫǫ (same as figure 2). The dark blue region assumes both ∆ σ ′ , ∆ ǫ ′ ≥ 3, and uses the system of multiple correlators. All bounds other than the black dotted line are computed with n max = 6, ν max = 8 (21 components for single correlator bounds, 105 components for multiple correlator bounds). The 3D Ising point is indicated with black crosshairs.
region or a single isolated point.
We note that in our determinations we did not assume the c-minimization conjecture or anything similar. The only assumption besides unitarity and conformal symmetry was the existence of a Z 2 symmetry and the assumption that σ and ǫ are the only relevant scalars. It is therefore encouraging that the two methods are in such good agreement.

Discussion
In this work we have elucidated the power of mixed correlators in the context of the 3D conformal bootstrap. While the simplest upper bound on the leading Z 2 -even operator dimension ∆ ǫ does not differ from the single correlator bootstrap, mild assumptions about the number of relevant operators give rise to very tight constraints on the allowed values of allowed region with various gaps in ∆ ǫ ′ , ∆ σ ′ , zoom (n max = 6) ∆ σ ∆ σ and ∆ ǫ , almost uniquely determining their values. Our results support the conjecture that the 3D Ising CFT is the only Z 2 -symmetric CFT in 3 dimensions with exactly two relevant operators. No other such CFT has been found experimentally, and it appears that using bootstrap techniques a numerical "proof" may be forthcoming.
Moreover by considering the mixed correlator bootstrap we are also able to gain information about the Z 2 -odd spectrum, finding a general upper bound on ∆ σ ′ . We fully anticipate that further studies of the mixed correlator bootstrap will yield an accurate picture of the complete low-lying spectrum of the 3D Ising CFT.
There are several directions for future research. First, given the vital role that semidefinite constraints play in general formulations of the conformal bootstrap, it is important to find and implement improved algorithms for high-precision solutions of semidefinite programs of the type encountered in this work. Such improvements will make it much easier to perform broad explorations of the space of conformal field theories in general dimensions.
More concretely, it would be interesting to perform similar studies of the simplest multiple correlator constraints in D = 3, as well as in CFTs with different global symmetry groups. For example, in 2D one could understand what assumptions are needed in order to isolate the minimal model solutions, in 3D one could perform similar studies of the O(N) vector models, and in 4D one could try to better understand the space of CFTs with a small number of relevant operators which may have phenomenological interest. Moreover, the time is ripe to begin including constraints from 4-point functions of operators with spin allowed region with ∆ ǫ ′ , ∆ σ ′ ≥ 3 (n max = 10)  Figure 5: Allowed and disallowed (∆ σ , ∆ ǫ ) points in a Z 2 -symmetric CFT 3 with only one relevant Z 2 -odd and Z 2 -even scalar, using the constraints of crossing symmetry and unitarity for σσσσ , σσǫǫ , ǫǫǫǫ at n max = 10 (275 components), ν max = 14. The light grey points are ruled out, while the dark blue points are allowed. The light blue shaded region shows the region allowed by crossing symmetry and unitarity of the single correlator σσσσ at n max = 14, computed in [22]. The final allowed region is the intersection of this shaded region with the region indicated by the dark blue points (see figure 6). 23 -such studies will likely use techniques similar to what we have developed in this work.
It is also interesting to study mixed correlators in theories with supersymmetry. In particular, previous numerical bootstrap studies have focused on 4-point functions containing the lowest component of a given supersymmetry multiplet, while mixed correlators could allow one to incorporate the full constraints of supersymmetry on the external operators. 24 Such studies may help to clarify the origin of the "kink" observed in previous studies of the 4D N = 1 superconformal bootstrap [9] and may also reveal rich new structure in theories with N = 2, 4 supersymmetry, extending the results of [14,18,25]. Finally, there is significant room for incorporating mixed correlators into general analytical studies of the bootstrap, both in the context of large N theories [26,27] and in constraining the spectrum at large spin [28,29]. 23 The computed points in figure 5 lie on a grid, where each row has constant ∆ ǫ − ∆ σ , because ∆ ǫ − ∆ σ is the quantity entering the conformal blocks g ∆12,∆34 ∆,ℓ (u, v). Restricting it to a small number of values means we have fewer tables of blocks to compute. We thank Slava Rychkov for this idea. 24 One can think about this in two ways: in components, we have four-point functions of different operators in the same SUSY multiplet; in manifestly supersymmetric notation, multiple superconformally covariant structures can appear in a three-point function. : Allowed values of (∆ σ , ∆ ǫ ) in a Z 2 -symmetric CFT 3 containing only two relevant scalars. The blue region is our rigorous bound from figure 5, computed at n max = 10, 14. The dark grey rectangle is the Monte-Carlo prediction of [35]. The light grey rectangle is the prediction of the c-minimization conjecture [22], using single-correlator results at n max = 21. There may be additional disconnected regions for ∆ σ 0.54, as in figure 2, but we have not computed them here.
While the conformal bootstrap involving identical external operators has already shown itself to be surprisingly constraining, our results demonstrate that the larger system of mixed constraints, combined with mild assumptions about gaps, may be sufficiently powerful to uniquely locate isolated CFTs. Indeed, if one previously did not know about the 3D Ising CFT, one would have discovered it following the general logic of this paper! There may be many more isolated CFTs waiting to be discovered, perhaps theories without Lagrangian descriptions or supersymmetry. The space of such theories can be mapped out in a systematic way using the conformal bootstrap, by inputting gaps and searching for small closed allowed regions in the space of operator dimensions. There is much exploration to be done!

A Multiple SU(n) Three-Point Structures
In this appendix, we give a concrete example of multiple structures appearing in a threepoint function of operators charged under a global symmetry. Let G = SU(n) and let r be the largest irreducible representation in Sym 2 Ad G , of dimension 1 4 n 2 (n − 1)(n + 3). The symmetric tensor square of r decomposes into irreducibles as Consequently, there are two independent three-point structures (and hence OPE coefficients) in the three-point function φ r φ r O r when O r has even spin. Let us write these structures explicitly in the case where O = φ.
In terms of SU(n)-indices, φ kl ij has two symmetric upper and two symmetric lower indices, and satisfies the tracelessness condition φ il ij = 0. To write correlators of φ, it is convenient to use index-free notation, where we contract φ with auxiliary bosonic vectors U, V in the fundamental and dual representations, respectively, The operator with explicit indices can be recovered by differentiating with respect to U, V and subtracting traces of the form δ k i , δ l j , δ k j , δ l i , Any expression of the form V ·U does not contribute after subtracting traces. Hence, φ(U, V ) should only be defined modulo the ideal of functions proportional to V · U. Quotienting by this ideal is equivalent to restricting φ(U, V ) to the locus V · U = 0, so we will henceforth impose this condition. To apply (A.3), we can choose an arbitrary extension of φ(U, V ) away from V · U = 0 and then differentiate. Similar index-free techniques were used for classifying correlators of operators with spin in [40,41,50].
A correlator of φ(U m , V m , x m )'s must be a function of the SU(n)-invariants Q mn ≡ V m ·U n (with n = m since V m · U m = 0) which is quadratic in each of the V m , U m . For a three-point function, there are two such structures consistent with permutation symmetry, Each structure comes with its own OPE coefficient λ 1 , λ 2 . The explicit SU(n)-indices for this three-point function can be recovered by applying (A.3) for each operator.

B Implementation in SDPA-GMP
In this appendix, we follow the notation of the SDPA manual [82]. In section 3.3, we expressed our semidefinite program in terms of the variables a i mn and the positive semidefinite matrices A ℓ , B ℓ , C ℓ , D ℓ . These are subject to linear constraints (3.22, 3.23), where we equate coefficients of each power of x on both sides. The a i mn are unconstrained. We can write them in terms of positive variables by introducing a "slack variable" s ≥ 0, and writing a i mn = b i mn − s, where b i mn ≥ 0. All unknowns can now be grouped into one block-diagonal positive semidefinite matrix Here, m, n, i run over the 5n max (n max + 1)/2 components of the functional α, and ℓ runs over spins up to some large maximum value. In this work, we take ℓ = 0, 1, . . . , 25, 26, 49, 50. Derivatives of the conformal blocks converge rapidly as ℓ → ∞, and in practice ℓ max ≈ 50 is enough to ensure appropriate positivity conditions for all ℓ (one can check this by plotting functionals at high ℓ once they are determined).
Y is the matrix of unknowns in the "dual" formulation of a semidefinite program defined in [82]. The equality conditions (3.22, 3.23) can be expressed in the form where F i are symmetric matrices with the same block-structure as Y , i runs over spins ℓ and powers of x entering which can also be written in the form (B.2) with c i = 1, where only the entries in F i corresponding to b i mn , s are nonzero. Since we are only interested in determining whether a feasible solution exists for Y , and not in optimizing a particular function, we take the dual objective function F 0 to be identically zero.
Sometimes we have isolated operators in the OPE (for instance, in the Ising model, ǫ is isolated from the remaining Z 2 -even scalars which have ∆ ≥ 3). Demanding that α be positive on the contribution of these operators gives additional semidefiniteness constraints on the variables b i mn , s. To accommodate these, we extend the matrix Y with semidefinite matrices E k , The semidefiniteness constraints now become equalities relating b i mn , s, and E k , which can again be written in the form (B.2). This suffices to write our semidefinite program in the form required by SDPA-GMP.
If our semidefinite program is feasible (i.e., if a functional α exists satisfying our constraints), then it should be possible to reduce the primal objective function to zero (the   dual objective function is identically zero). By decreasing the parameter epsilonDash, we can force SDPA-GMP to make the primal objective smaller and smaller. Failure to decrease the primal objective below a given finite value means the SDP is infeasible. In practice, we set epsilonDash very small (∼ 10 −30 ), and use the final value of the primal objective as a measure of whether the problem is feasible or infeasible. Our SDPA-GMP parameters are summarized in table 2. Our condition for feasibility is |primalObjective| ≤ 10 −13 .
Our computational setup is as follows. A Mathematica program computes tables of derivatives of conformal blocks using the recursion relation described in section 4. A separate Mathematica program reads these tables and writes a semidefinite program to a file in sparse SDPA format. This file is read and solved by SDPA-GMP itself. We have modified SDPA-GMP to allow checkpointing: it periodically saves its state to a file so it can be started and stopped at will. 25 Thus, computations taking several days can be interrupted safely without having to start over again from the beginning. The n max = 10 computations in this work are quite time-intensive. Writing the SDP to a file takes about 30 minutes, and solving it takes approximately 2 weeks. The n max = 6 computations are much less time-intensive, taking about 8 hours each. It is extremely useful to run several computations in parallel on a computing cluster. Checkpointing allows us to set small time-limits on each individual process, continually freeing up the cluster for jobs from other users. 26 25 The checkpointed version of SDPA-GMP is available at https://bitbucket.org/davidsd/sdpa-gmp-checkpointed/overview. 26 Our cluster management software uses Cloud Haskell [83,84] and MongoDB [85].