The Effective Bootstrap

We study the numerical bounds obtained using a conformal-bootstrap method - advocated in ref. [1] but never implemented so far - where different points in the plane of conformal cross ratios $z$ and $\bar z$ are sampled. In contrast to the most used method based on derivatives evaluated at the symmetric point $z=\bar z =1/2$, we can consistently"integrate out"higher-dimensional operators and get a reduced simpler, and faster to solve, set of bootstrap equations. We test this"effective"bootstrap by studying the 3D Ising and $O(n)$ vector models and bounds on generic 4D CFTs, for which extensive results are already available in the literature. We also determine the scaling dimensions of certain scalar operators in the $O(n)$ vector models, with $n=2,3,4$, which have not yet been computed using bootstrap techniques.


Introduction
There has recently been a great revival of interest in the conformal bootstrap program [2,3] after ref. [4] observed that its applicability extends to Conformal Field Theories (CFTs) in d > 2 dimensions. Since ref. [4], considerable progress has been achieved in understanding CFTs in d ≥ 2 dimensions, both numerically and analytically. Probably the most striking progress has been made in the numerical study of the 3D Ising model, where amazingly precise operator dimensions and OPE coefficients have been determined [5][6][7].
Essentially all numerical bootstrap studies so far have used the constraints imposed by crossing symmetry on 4-point correlators evaluated at a specific value of the conformal crossratios, u = v = 1/4, or equivalently in z-coordinates at z =z = 1/2 [8]. This is the point of best convergence for the combined conformal block expansions in the s and t channels. Taking higher and higher derivatives of the bootstrap equations evaluated at this point has proven to be very effective and successful in obtaining increasingly better bounds. We will denote this method in the following as the "derivative method". A drawback of the derivative method -both in its linear [4,6,9] or semi-definite [10,11] programming incarnations -is the need to include a large number of operators in the bootstrap equations. This makes any, even limited, analytical understanding of the obtained results quite difficult.
A possible approximation scheme is in fact available: ref. [12] has determined the rate of convergence of the Operator Product Expansion (OPE), on which the bootstrap equations are based. This allows us to extract the maximal error from neglecting operators with dimensions larger than some cutoff ∆ * in the bootstrap equations and thus to consistently truncate them.
These truncated bootstrap equations can then be evaluated at different points in the z-plane. This method, which we denote as the "multipoint method", has been previously advocated by Hogervorst and Rychkov in ref. [1] but has not yet been numerically implemented. The aim of this note is to provide such an implementation and study the resulting bounds. It is important to emphasize that the method of ref. [1] combines what are in principle two independent ideas: i) multipoint bootstrap and ii) truncation of the bootstrap equations. One could study i) without ii), or try to analyze ii) without i). We will not consider these other possibilities here.
We begin in section 2 with a brief review of the results of refs. [1,12,13] on the convergence of the OPE. We use generalized free theories as a toy laboratory to test some of the results obtained in ref. [12]. We then generalize the results of ref. [12] for CFTs with an O(n) global symmetry.
We write the bootstrap equations and set the stage for our numerical computations in section 3. Our results are then presented in section 4. For concreteness, we study bounds on operator dimensions and the central charge in 3D and 4D CFTs, with and without an O(n) global symmetry (with no supersymmetry). For these bounds, extensive results are already available in the literature (see e.g. refs. [5-7, 10, 14-22]). In particular, we focus our attention on the regions where the 3D Ising and O(n) vector models have been identified. We show how the results depend on the number N of points in the z-plane at which we evaluate the bootstrap equations and the cut-off ∆ * on the dimension of operators in the bootstrap equations. Using values for the dimension of the operator φ in O(n) vector models available in the literature and a fit extrapolation procedure, we then determine the dimensions of the second-lowest O(n) singlet and symmetric-traceless operators S and T for n = 2, 3, 4. To our knowledge, these have not been obtained before using bootstrap techniques. Our results are consistent with those from analytical calculations using the -expansion [23,24] with a mild tension with the result of ref. [24] for the dimension of T in the O(2) model. We notice from our results that the "kink" in the bound on the dimension of the lowest scalar (singlet) operator in 3D Ising and O(n) vector models is already visible for relatively small ∆ * , while the minimum in the central-charge bound is very sensitive to ∆ * . For our numerical implementation, we discretize the spectrum and formulate the bootstrap equations as a linear program which we solve using the optimizer CPLEX 1 by IBM. Since we focus on the truncated bootstrap equations with relatively low cutoff ∆ * , double precision as used by CPLEX is sufficient for our purposes. More refined implementations with higher numerical precision, possibly adapting the method and optimizer of refs. [6,9], are certainly possible. More details on the numerical implementation are given in section 5. We conclude in section 6.

Convergence of the OPE
We begin with a brief review of the results of refs. [12,13] (see also ref. [1]) about the convergence of the OPE in a euclidean, reflection positive, CFT in any number of dimensions. 2 For more details see the original references. Consider the 4-point function of a scalar primary operator φ with scaling dimension ∆ φ : are the conformally-invariant cross-ratios (x ij ≡ x i − x j ). Applying the OPE to the operator pairs φ(x 1 )φ(x 2 ) and φ(x 3 )φ(x 4 ) in the 4-point function, one can write where u = zz, v = (1 − z)(1 −z) and the sum runs over all primary operators O that appear in the φ × φ OPE with ∆ and l being respectively their dimension and spin. For each primary, the sum over all its descendants is encoded in the conformal block function g ∆,l (z,z). In a euclidean CFT,z = z * and the conformal blocks are regular everywhere in the complex z-plane, with the exception of a branch-cut along the real line [1, +∞). 3 Thanks to reflection positivity, the OPE coefficients λ O are real and thus λ 2 O > 0. Crucial for our considerations will be a bound on the remainder of the sum in eq. (2.3) when it is truncated at some primary operator of dimension ∆ = ∆ * . To determine this bound, one first uses that as follows e.g. from a representation of the conformal blocks in terms of Gegenbauer polynomials [1]. It is therefore sufficient to estimate the remainder for real z =z. As was found in ref. [12], the most stringent bound is obtained by using the coordinate Bounds on the OPE convergence are obtained in an alternative way using crossing symmetry in ref. [25].
Interestingly, ref. [25] sets bounds which are also valid for finite values of ∆ * at z =z = 1/2, though they are relative and not absolute bounds. It would be interesting to explore the approach followed in this paper further. We thank Slava Rychkov for having pointed out this reference to us. 3 The branch-cut is best seen in Lorentzian signature, where z andz are two independent variables. At fixedz (z), g ∆,l (z,z) is a true analytic function in z (z) with a branch-cut along the line [1, +∞).
The z-plane is mapped to the unit disk in ρ and the branch-cut is mapped to the boundary of the disk. The conformal blocks in ρ are then defined for |ρ| < 1. In the manifestly reflection positive configuration withρ = ρ = r, the function g(u, v) in eq. (2.3) can be written as 4 where c n (∆, l) are positive coefficients whose explicit form is not important here and the sum over n takes into account the contributions from the descendants of each primary. It is convenient to rewrite g(r) as Here β ≡ − log r, k runs over all operators (primaries and their descendants) which are exchanged in the OPE and f (∆) is a spectral density with positive coefficients ρ k . Again, their explicit form is not relevant for our considerations. The behaviour of g(β) in the limit β → 0 (corresponding to the OPE limit x 3 → x 2 , in which case z → r → 1 and 1 − z → β 2 /4 → 0) is dominated by the exchange of the identity operator and one finds: Here a ∼ b means that a/b → 1 in the considered limit. The key observation of ref. [12] is that since the coefficients ρ k are all positive, this asymptotic behaviour determines the leading, large-∆ behaviour of the integrated spectral density by means of the Hardy-Littlewood tauberian theorem (see e.g. [26]): 6 .
The remainder (2.4) can then be bounded as follows: We first note that 4 For simplicity, we use the same symbol to denote the functions g(u, v) andg(r) = g(u(r), v(r)) etc. here and below. 5 This is true in general only in d > 2 dimensions. In d = 2, one has to be careful since scalar operators can have arbitrarily small dimensions. See also the discussion after eq. (2.23). 6 It is in fact sufficient that the coefficients are all positive for operators with dimension larger than some fixed value ∆0.
since the r.h.s. contains contributions from all operators with dimension larger than ∆ * , whereas on the l.h.s. only primaries with dimension larger than ∆ * and their descendents contribute. Using eq. (2.11), the r.h.s. can in turn be bounded as where Γ(a, b) is the incomplete Gamma function. Clearly, this bound applies for parametrically large values of ∆ * , where eq. (2.11) holds. Using eq. (2.5), we finally get the bound on the remainder This is valid in any number d > 2 of dimensions for 4-point functions with identical scalars. It was pointed out in ref. [13] that the conditions for the applicability of the Hardy-Littlewood tauberian theorem in both 3 and 4 dimensions are also fulfilled for the rescaled conformal blocks with γ = 1. Repeating the derivation reviewed above for a remainder involving the rescaled conformal blocks, it is straightforward to get the alternative bound For −∆ * log |ρ(z)| 1, eq. (2.17) can be approximated as We see that for |ρ(z)| not too close to 1 and ∆ * 8∆ φ , the bound is more stringent for γ = 1 than for γ = 0. It was furthermore shown in ref. [13] that in d = 3 dimensions, γ = 1 is the maximal allowed value such that the Hardy-Littlewood tauberian theorem remains applicable, whereas it was conjectured without proof that the maximal allowed value in d = 4 dimensions is γ = 3/2. Correspondingly we use eq. (2.17) with γ = 1 for the remainder both in 3 and 4 dimensions in our numerical implementation. 7 The above derivations were based on the existence of a configuration for which the function g(u, v) turns into a positive definite function of a single variable. The remainder is then estimated using the Hardy-Littlewood tauberian theorem. One cannot naively apply these arguments to arbitrary derivatives of g(u, v) w.r.t. u and v, unless the resulting functions remain positive definite and derivatives can be brought inside the absolute value in the l.h.s. of eq. (2.16). See the appendix of ref. [27] for a recent discussion on how to estimate the remainder on derivatives of g(u, v). It would be interesting to verify if this allows us to also study truncated bootstrap equations with the derivative method.

Comparison with Generalized Free Theories and Asymptotics for z → 1
The results reviewed in the previous subsection are based on eq. (2.11) which holds in the limit ∆ * → ∞. Of course, for any practical use, we need to know the value of ∆ * beyond which we can trust eq. (2.11) and thus the bound eq. (2.16). It is difficult to determine this value for a generic CFT. But we can get useful insights by considering exactly calculable CFTs, like generalized free theories (sometimes called mean field theories) for which the CFT data are known and the function g(u, v) in eq. (2.1) in any number of dimensions reads ( In fig. 1, we show η as a function of ∆ * evaluated at the symmetric point z =z = 1/2. Notice that at the point of best convergence the actual remainder is always significantly smaller than R, and that the ratio gets bigger and bigger as ∆ * increases for large ∆ * . In particular, η is greater In order to make the bound more stringent, one could then alternatively use the series representation in ref. [1] which includes contributions from primary operators and their descendants separately. Using this series truncated at contributions corresponding to dimension ∆ * instead of the full conformal blocks g ∆,l would make the r.h.s. of the inequality (2.12) the actual remainder to be bounded. This would thus make eq. (2.16) with γ = 0 more stringent. Here, however, we choose not to follow this approach. The reason is that the representations for the full conformal blocks g ∆,l can be considerably faster calculated than (our implementation of) the truncated series representation of ref. [1]. When z → 1, both the numerator and the denominator of η in eq. (2.20) blow up, since the OPE is not convergent at z =z = 1. Operators with high scaling dimension are no longer suppressed and the remainder completely dominates the OPE. 8 More precisely, we have independently of γ. Notice that this limit is universal for any CFT that includes in its spectrum a scalar operator with dimension ∆ φ , because z =z → 1 selects the universal identity contribution in the t-channel. This class of CFTs always includes a GFT for the operator φ itself. In this case the universal nature of the limit is trivially checked using eq. (2.19): where in the last equality we have used that |1 − z| → (log |ρ(z)|) 2 /4 in the limit. It was found in refs. [28,29] that the spectrum of any Lorentzian CFT resembles that of a GFT for parametrically large spin operators. In particular, in ref. [28] this has been established by analyzing crossing symmetry in the limit z → 0 andz fixed for d > 2, where large twist operators are suppressed. The two-dimensional case is more subtle, because there is no longer a gap between the identity (which has the minimum twist zero) and the other operators. Indeed, the results of refs. [28,29] and those of ref. [12] in the euclidean do not straightforwardly apply for d = 2.
In the euclidean, operators of any twist should be considered. However, given the results of refs. [28,29], it is natural to expect that the leading behaviour (2.22) is expected to come from operators with parametrically high dimension and high spin for any CFT, asymptotically approaching the GFT spectrum in this regime. It would be interesting to understand within euclidean CFTs, where the twist does not play an obvious role, why this is so.

Remainder for CFTs with O(n) Symmetry
The generalization of the OPE convergence estimate to CFTs with O(n) global symmetry is straightforward. For concreteness, let us consider scalars φ i in the fundamental representation of O(n). The only non-trivial point is to identify a proper linear combination of 4-point functions that leads to a positive definite series expansion, otherwise the Hardy-Littlewood tauberian theorem does not apply. A possible choice is where for simplicity we have omitted the x-dependence of the fields. The parameter η can in general take an arbitrary complex value, but it is enough for our purposes to consider η = ±1. Forρ = ρ = r and any η, this correlator is manifestly positive definite, because it corresponds to the norm of the state φ 1 |φ 1 + ηφ 2 |φ 2 . (2.26) The leading term in a η (u, v) for x 2 → x 3 is given by the exchange of the identity operator in the first two correlators and hence is independent of η. On the other hand, expanding in conformal blocks in the (12)-(34) channel, we have [19] A η = 1 where S and T denote operators in the singlet and rank-two symmetric representations of O(n), respectively. Both sums run over even spins. We can now repeat essentially verbatim the derivation below eq. (2.6). For η = −1, this gives rise to the bound where R is given in eq. (2.17). The factor 1/2 with respect to the non-symmetric case arises because the identity operator is exchanged in two correlators but a factor 4 is present in the second term in the r.h.s. of eq. (2.27). For η = 1 we similarly get Another positive definite linear combination of correlators is corresponding to the norm of the state Again, we consider η = ±1. In the (12)-(34) channel the correlator B η can be written as 9 where A stands for operators in the rank-two antisymmetric representation of O(n). The first sum runs over even spins, whereas for the second one they are odd. As before, the leading term in b η (u, v) for x 2 → x 3 is given by the exchange of the identity operator in the first two correlators and is independent of η. For η = 1, eq. (2.32) gives rise to the same bound given in eq. (2.28), while for η = −1 we have It is straightforward to see that the bounds (2.28), (2.29) and (2.33) are the best that can be obtained. Indeed, in the free-theory limit one has λ 2 S = λ 2 /n, λ 2 T = λ 2 A = λ 2 /2 with λ 2 being the OPE coefficients for a single free field (see e.g. eq. (5.11) in ref. [20]). The above three bounds then reduce to eq. (2.16) which is known to give the best bound on the r.h.s. of eq. (2.12) (see however footnote 7) [12]. Any potentially better bound for O(n) theories should in particular apply to the free theory, but would then be in contradiction with the results of ref. [12].
The above bounds will be used in the next section to bound the remainder of the bootstrap equations in CFTs with an O(n) global symmetry.

Bootstrapping with Multiple Points
The bootstrap equation for a 4-point function with identical scalars φ with scaling dimension ∆ φ in any number of dimensions is given by the sum rule (see refs. [30,31] for pedagogical reviews) Splitting the sum into two parts, for dimensions smaller and larger than a cutoff ∆ * , we can write Using eq. (2.16), the remainder E of the sum rule is bounded by where we have omitted the dependence on ∆ * , ∆ φ and γ. The truncated sum rule (3.2) still involves a generally unknown spectrum of operators up to dimension ∆ * . In order to make it amenable to numerical analysis, we discretize the spectrum and make the ansatz 10 for the quantum numbers (spin,dimension) of the operators that can appear in the truncated sum rule. For each spin l, the dimension runs in steps of size ∆ step from the unitarity bound ∆ d,l min ≡ l + (d − 2)/(1 + δ l0 ) to the cutoff ∆ * (or a value close to that, depending on ∆ step ). Accordingly, l max is the largest spin for which the unitarity bound is still below the cutoff, ∆ d,lmax min < ∆ * . In practice, we vary the step size ∆ step somewhat depending on the spin and dimension. This is discussed in more detail in sec. 5. We find that the bounds converge when going to smaller ∆ step , meaning that the discretization does not introduce any artifacts into our calculation.
We similarly choose a finite number of points z i in the z-plane where the sum rule is evaluated. The details of our choice for this distribution of points are discussed in sec. 3.1. Together with the discretization of operator dimensions, this turns eq. (3.2) into the matrix equation The elements of the matrix M are the functions F ∆ φ ,∆,l (z,z) evaluated for the different quantum numbers in eq. (3.4) along the rows and for the different points z i along the columns. Furthermore, the vector ρ consists of the squared OPE coefficients λ 2 O of the operators corresponding to the quantum numbers in eq. (3.4) and Using the bound (3.3), we then obtain the matrix inequality where max is defined as but with E replaced by E max . This is the starting point for our numerical calculations. In order to determine bounds on OPE coefficients, we search for vectors ρ which satisfy eq. (3.7) and extremize the entry corresponding to that OPE coefficient. For bounds on the dimension of the lowest-lying scalar operator, on the other hand, we make an assumption on this dimension and drop all scalar operators with smaller dimension from our ansatz (3.4). This gap then allows for a consistent CFT only if there exists a vector ρ which satisfies eq. (3.7) with the reduced ansatz. By trying different assumptions, we can determine the maximal allowed gap. Both problems are linear programs which can be solved using fast numerical routines. An advantage of solving eq. (3.7) is that the vector ρ gives us the spectrum of operators and their OPE coefficients of a potential CFT living at the boundary of the allowed region. This has been used before in ref. [6]. 11 We also consider CFTs with an O(n) global symmetry. For an external scalar operator in the fundamental representation of O(n), the sum rule reads [19] and we have suppressed the arguments of the functions F and H. Splitting the sums in eq. (3.8) into two parts, for dimensions smaller and larger than a cutoff ∆ * , we can write Using eqs. (2.28), (2.29) and (2.33), we obtain the bounds on the remainders with E max defined as in eq. (3.3). Discretizing the space of operator dimensions as in eq. (3.4) and evaluating the sum rule at a finite set of points z i , we again obtain a matrix inequality of the form (3.7). This is the starting point for our numerical calculations for CFTs with O(n) global symmetry.

Choice of Points
An important choice for the multipoint method is the distribution of points in the z-plane at which the bootstrap equations are evaluated. Using the symmetries z ↔z and z ↔ (1 − z), z ↔ 1 −z of the bootstrap equations, we can restrict these points to the region Re(z) ≥ 1/2 and Im(z) ≥ 0 of the z-plane. The remainder of the truncated sum rule is controlled by |ρ(z)| and |ρ(1 − z)| (cf. eqs. (2.18) and (3.3)). Guided by this, we introduce the measure and consider points with λ(z) ≤ λ c for some constant λ c . It is desirable to choose λ c and the distribution of points within that region in such a way that the obtained bounds are as stringent as possible. We have performed extensive scans over different values for λ c and distributions with different density profiles and have found that a flat profile leads to as good or better bounds than more complicated profiles. We therefore choose the former and put points on a grid centered at z = 1/2. The grid spacing is chosen such that the desired number of points is within the region λ(z) ≤ λ c , Re(z) ≥ 1/2 and Im(z) ≥ 0. We have then found that λ c = 0.6 (3.12) gives the best bounds for all cases that we have studied. 12 In fig. 2, we show the corresponding region in the z-plane and a sample distribution of 100 points. In order to test the influence of the choice of measure on the bounds, we have performed further scans with λ(z) ≡ max(|ρ(z)|, |ρ(1 − z)|) proposed in ref. [1] and λ(z) ≡ |z − 1/2| (for the latter we have removed points at or close to the branch-cuts). We have found that, once the optimal λ c is chosen, the bounds obtained with these measures are indistinguishable from those obtained with eq. (3.11). This indicates that the precise form of the region within which points are sampled has only a marginal effect on the quality of the bounds.

Results
We now present the results of our numerical analysis. In subsection 4.1, we study bounds on the dimension of the lowest-dimensional scalar operator in the OPE and bounds on the central charge in 3D CFTs, focusing in particular on the regions where the 3D Ising and O(n) models have been identified. In subsection 4.2 we then study the same bounds for generic 4D CFTs. We analyze in particular how our results depend on the number N of points chosen in the z-plane, and on the cutoff ∆ * . In subsection 4.3 we give a closer look at the spectrum of the 3D O(n) models and determine the operator dimensions of the first two scalar operators in the singlet and rank-two symmetric representation of O(n).
Before presenting our results, it is important to emphasize an important difference between the multipoint and the derivative bootstrap methods. As mentioned in the introduction, in the latter we do not have a reliable way of truncating the OPE series defining the bootstrap equations at some intermediate dimension ∆ * , because we do not have a reliable estimate of the resulting error. We are therefore forced to have ∆ * as large as possible to minimize this error and can only check a posteriori if the chosen ∆ * was sufficient. 13 More than ∆ * (or its analogue), the key 12 In more detail, we have considered bounds on the central charge and the dimension of the lowest-dimensional scalar operator, in 3D and 4D, with O(n) and without symmetry, and with different choices for the number of points N and the cutoff ∆ * . It is remarkable that λc = 0.6 (within ±0.02, the resolution of our scan) comes out as the optimal choice for such a variety of cases. 13 We are a bit sloppy here in order to keep the discussion simple and get to the point. For instance, in numerical methods based on semi-definite programming one is able to include all operator dimensions continuously up to infinity. The rough analogue of our ∆ * in that case is the maximum spin of the primary operators entering the parameter that controls the accuracy of the method is given by the total number of derivatives N D that are applied to the bootstrap equations. Of course, the larger N D is, the better are the bounds. The accuracy is then limited by the largest N D that allows the calculation to be performed within an acceptable amount of time with the available computing resources.
In the multipoint method, on the other hand, we can reliably vary ∆ * due to the bound on the remainder of the truncation discussed in sec. 2. In addition, we can also vary the number N of points in the z-plane which is the analogue of N D in the derivative method. The parameter region for the multipoint method corresponding to the typical bootstrap analysis with the derivative method is then very large ∆ * and N as large as possible given the available computing resources. In this paper, on the other hand, we are mostly interested in the regime where ∆ * is not very large, with values O(10)-O (20). We find that for this range of ∆ * , the results converge for N ∼ O(100) and do not improve further if N is increased. This corresponds to the fact that the rank of the matrix M in the discretized bootstrap equation (3.5) is then O(100). Note that since CPLEX is limited to double precision, we also cannot take ∆ * arbitrarily large. Due to the excellent speed of CPLEX, on the other hand, we have found that taking N large enough so that the bounds converge is no limiting factor.
OPE which are taken into account for the numerical implementation.

3D Ising and O(n) Models
The most remarkable numerical results from the conformal bootstrap have been obtained in 3D CFTs. One interesting bound to study is on the dimension of the lowest-dimensional scalar operator appearing in the OPE. We denote this operator by and the operator that is used to derive the bootstrap equations by σ. It was noted in ref. [5] that the 3D Ising model sits at a special point, a kink, at the boundary of the allowed region of ∆ as a function of ∆ σ . The Ising model is similarly special with respect to the bound on the central charge c as a function of ∆ σ , sitting again at the boundary of the excluded region, at the point where c is minimized [5,6]. Note, however, that the theory minimizing c does not actually correspond to the 3D Ising model, but rather to some exotic theory with ∆ < 1. Most likely this theory is unphysical (though we are not aware of a solid argument to dismiss it). In practice this theory is removed by assuming a gap in the operator spectrum such that ∆ > 1. Independently of the nature of this theory, the condition ∆ > 1 is satisfied by the Ising model and can be legitimately imposed if we are interested in this particular 3D CFT. In fig. 3, we show the bound on ∆ as a function of ∆ σ for N = 100 points and different values of ∆ * . Notice how the kink shows up already for ∆ * = 13 and converges quite quickly as ∆ * increases. In fig. 4, we show the bound on the central charge c (normalized to the central charge c free of a free scalar theory) as a function of ∆ σ for N = 100 points and different values of ∆ * . The gap ∆ > 1.1 is assumed in the operator spectrum. A lower bound on c is obtained even for ∆ * = 10, but the convergence when going to larger ∆ * is now much slower than for the bound on ∆ . A minimum is visible starting from ∆ * = 16 but even at ∆ * = 22 it is a bit shifted to the right with respect to its actual value. We have still not reached the asymptotic value for ∆ * . Unfortunately, we cannot get reliable results for much higher ∆ * because the numerical accuracy of CPLEX is limited to double precision. Nevertheless, it is clear from comparing figs. 3 and 4 that the lower bound on c is more "UV sensitive" than the bound on ∆ . In both figures, the crosses mark the location of the 3D Ising model, as determined in ref. [6].
In order to quantify the dependence of our results on the number N of points, we show in figs. 5 and 6 the bounds on respectively ∆ and c as a function of ∆ σ for different values of N at fixed ∆ * = 16. We see that in both cases the convergence in N is quite fast, with N = 40 for ∆ and N = 60 for c being already an excellent approximation. Notice that for increasing N , the bound on ∆ converges faster than the bound on c, similar to the dependence on ∆ * . We have studied the dependence on N also for different values of ∆ * and have found as expected that the value N * beyond which no significant improvement in the bounds is observed increases with ∆ * . The dependence is however very mild for the central charge c and barely observable for ∆ . This is still a reflection of the different "UV sensitivities" of the two quantities. In all cases, N * O(100) up to ∆ * = 24.
Let us now turn to 3D CFTs with O(n) global symmetry. We consider a primary operator φ in the fundamental representation and denote the lowest-dimensional scalar singlet operator in the φ × φ OPE by S. It was found in refs. [14,16] that these CFTs have kinks in the bound on ∆ S as a function of ∆ φ similar to that found for the Ising model. Moreover, the kinks coincide, for all values of n that have been studied, with the values of ∆ φ and ∆ S associated with the 3D O(n) models. On the other hand, a minimum in c no longer occurs for generic O(n) models and the lower bound on c instead monotonically decreases for n > 3 (see ref. [14] for details).
In figs. 7 and 8, we show respectively the bound on ∆ S and c (the latter normalized to the central charge nc free of n free scalars) as a function of ∆ φ for different O(n) symmetries, at fixed N = 80 and ∆ * = 16. For the central charge, gaps ∆ S > 1 and ∆ T > 1 in the spectrum of respectively singlet operators S and rank-two symmetric-traceless operators T are assumed as in ref. [14]. This assumption is satisfied for the O(n) models and leads to more stringent bounds. The dashed line corresponds to the leading large-n prediction. All the qualitative behaviours found in ref. [14] are reproduced, though with milder bounds, as expected. 14 In particular, the kinks in the (∆ φ -∆ S ) plane are not well visible at ∆ * = 16. In figs. 9 and 10, we show the same bounds on ∆ S and c as a function of ∆ φ at fixed N and n, for different values of ∆ * . We see the  same qualitative behaviours regarding the "UV sensitivities" found for 3D CFTs with no global symmetry (the Ising model). In particular, in fig. 9 we see how the kink in the bound becomes well visible at ∆ * = 18 and does not significantly improve for ∆ * = 20. Its location is in very good agreement with that found in ref. [14]. On the other hand, the central-charge bound in fig. 10 is still monotonically decreasing for ∆ * = 18 and a minimum appears only for ∆ * = 20.
There are no signs of convergence comparing the bounds at ∆ * = 18 and 20, indicating the need to go to larger ∆ * to approach the optimal bound.

4D CFTs
All the above considerations can be repeated for 4D CFTs. There are no known non-supersymmetric CFTs at benchmarks points but it is still interesting to study general bounds on operator dimensions and OPE coefficients. See e.g. refs. [4,10,[17][18][19][20][21][22]33], where bounds of this kind (and others) have been determined with the derivative method using both linear and semidefinite programming. In figs. 11 and 12, we show bounds respectively on the dimension ∆ φ 2 of the lowest-dimensional The analysis of 4D CFTs with O(n) global symmetry also closely resembles its 3D counterpart. We again take the external field φ to transform in the fundamental representation of O(n) and denote by S the lowest-dimensional singlet scalar operator that appears in the φ × φ OPE. For illustration, we report in fig. 13 the bound on ∆ S as a function of ∆ φ for CFTs with O(4) symmetry, at fixed N and for different values of ∆ * . By comparing figs.11 and 13 we notice that the convergence in ∆ * of the operator-dimension bound in 4D CFTs with O(4) symmetry is slower than its analogue with no global symmetry.

A Closer Look at the Spectrum of 3D O(n) Models
In the last subsections, we have shown how previously determined bounds are reproduced using the multipoint method. Here we present some new results for the spectrum of O(n) models. To this end we assume, as previous analyses indicate, that the 3D O(n) models sit precisely at the kink on the boundary of the excluded region in the (∆ φ -∆ S ) plane (∆ S -maximization). The vector ρ that we obtain from solving the linear program (3.7) then gives us the spectrum and OPE coefficients of the operators that are exchanged in the φφφφ correlator of the O(n) models. Here we report the scaling dimensions of the first two operators in respectively the singlet and rank-two representation of O(n), S, S and T , T , for n = 2, 3, 4. Scalar operators with larger scaling dimensions are physically uninteresting, whereas S and T are important in determining the stability of the fixed points of the O(n) models (being marginal operators in the underlying UV 4D Landau-Ginzburg theory) [24]. 15 Actually, one additional operator should be considered, denoted as P 4,4 in ref. [24], but it transforms in the rank-four representation of O(n) and hence cannot appear in the OPE of two scalar operators φ in the fundamental representation. Its dimension might be bounded (or computed) by considering a correlator involving, e.g., four T 's. As far as we know, the scaling dimensions of S and T have not been previously determined using the conformal bootstrap. The best determinations of these parameters have been made   using a five-loop computation in the -expansion in refs. [23] and [24]. 16 In table 1, we report the values of ∆ φ , ∆ S , ∆ S , ∆ T , ∆ T determined in the literature, for n = 2, 3, 4. They should be compared with the values in table 2 which have been determined in this paper as follows: We take the values of ∆ φ for O(n) models with n = 2, 3, 4 calculated in refs. [34][35][36] as input and determine the scaling dimensions ∆ S , ∆ S , ∆ T and ∆ T using ∆ Smaximization. We repeat this procedure for the lower, central and upper value of ∆ φ given in these references and for different values of the cutoff ∆ * ∈ [18,23] and the number of points N ∈ [60, 120]. 17 At fixed N and ∆ * , we then take the average over the scaling dimensions obtained with the different input values of ∆ φ . Sometimes the same operator appears twice in the spectrum, at two different but close values of the scaling dimension. In this case we take the average of these values, weighted by the size of the corresponding OPE coefficient. Let us denote the resulting scaling dimensions by ∆ O (N, ∆ * ) for O = S, S , T, T . Each of these values is associated with an error, resulting from the averaging. The stepsize ∆ step of our discretization has been set to 10 −4 in the region where the operators were expected to be found (the resulting uncertainty in the scaling dimensions is typically negligible compared to the other errors).
At Carlo simulations. On the other hand, since ∆ T has been determined only using the -expansion, we have decided to omit the other results for ∆ S . The interested reader can find them, e.g., in table I of ref. [24], where the coefficients y4,0 and y4,2 give ∆ S = 3 − y4,0 and ∆ T = 3 − y4,2. For completeness, we also report the relations defining ∆S and ∆T in the notation of ref. [24]: ∆S = 3 − 1/ν, ∆T = 3 − y2,2.
17 Our numerical precision does not allow us to take higher values of ∆ * and N without having issues with numerical stability.   We have then extrapolated to N = ∞ using a linear fit in 1/N which seems to well describe the behaviour of ∆ O (N ) as a function of 1/N . An example of this extrapolation fit is shown in fig.15. We denote the resulting scaling dimensions as ∆ O ≡ ∆ O (∞). 18 We do not have an analytic understanding of why the results should scale as 1/N for parametrically large ∆ * . We simply take it as a working hypothesis. We expect that possible deviations from the linear behaviour should be contained within the errors of our determination (cf. fig.15). Note that having N as large as possible is clearly important for high precision. However, at fixed ∆ * the bounds saturate for sufficiently high N and there is no gain in taking N larger.
We have noticed that, at least for n = 2, 3, 4, ∆ O (N, ∆ * ) decreases as N and/or ∆ * increase (this is obvious for S, but not for the other operators). If we assume that this is true for any N and ∆ * , we may then set rigorous upper bounds without using any fit extrapolation. These bounds are reported in table 3. Comparing them with the results in table 2 gives an idea of the impact of the fit extrapolation on the final results. As can be seen, all the scaling dimensions that we have determined are compatible with previous results in the literature. The only exception is ∆ T for the O(2) model for which our result has an approximate 3σ tension with that of ref. [24]. Our accuracy in the determinations of ∆ S and ∆ T is comparable with that achieved in ref. [14], though it should be emphasized that the results there do not rely on extrapolations. Furthermore, our accuracy in the determinations of ∆ S and ∆ T is comparable with that achieved using the five-loop -expansion. This is an indication that a slightly more refined bootstrap analysis will be able to improve the determinations of these scaling dimensions.
As we mentioned at the beginning of this subsection, ∆ S -maximization also allows us to determine the OPE coefficients λ φφO . We have not performed a detailed analysis with fit extrapolations as above to determine the asymptotic values of λ φφO as ∆ * , N → ∞. Instead we just report λ φφS as determined with the highest values ∆ * = 22, 23 and N = 110, 120 used in We have not determined the error associated with these results and have instead rounded them to the last shown digit. The results for O(2) and O(3) are in agreement with the recent determination in ref. [7], whereas the result for O(4) is new as far as we know.

Details of the Implementation
For the conformal blocks in d = 4 dimensions, we use the closed-form expression from ref. [8], normalized as in ref. [19]. For d = 3 dimensions, on the other hand, we use the recursion relation for the conformal blocks found in ref. [14]. 19 To this end, we iterate the recursion relation up to some cutoff ∆ rec . We choose this cutoff large enough such that the resulting error in the conformal blocks is smaller than the error from neglecting contributions of operators with dimensions larger than the truncation cutoff ∆ * . In practice, we find that ∆ rec = ∆ * + few is sufficient to ensure this. For the ansatz (3.4) of discretized operator dimensions, we closely follow ref. [5]. We generate the discrete spectra T1 to T4 (the latter only for sufficiently large ∆ * ) in their table 2, where we rescale the stepsizes δ by the factor ∆ step /(2 · 10 −5 ). We then remove duplicates from the combined spectrum and restrict to operator dimensions less than or equal to ∆ * . We have performed extensive scans using different stepsizes ∆ step and have found that the bounds converge for sufficiently small ∆ step . This is in particular satisfied for ∆ step = 2 · 10 −3 which we choose for all the plots in this paper. For the determination of the spectra in sec. 4.3 we add additional operators with stepsize ∆ step = 10 −4 around the previously determined scaling dimensions for the operators S, S , T , T in the O(n) models. Furthermore, for bounds on operator dimensions for which the plots extend to bounds ∆ φ 2 > 3 (the largest dimension of T1 of ref. [5]), we have included additional operators in the scalar sector so that the smallest stepsize ∆ step is used up to the largest bound on ∆ φ 2 shown in that plot. We have also performed scans using different parametrizations for the ansatz (3.4) and have found that the bounds become indistinguishable from the bounds obtained with the ansatz discussed above for sufficiently small ∆ step . This gives us confidence that the discretization does not introduce any artifacts into our calculations.
We use Mathematica to evaluate the conformal blocks for the different operators that appear in the ansatz (3.4) and for the set of points in the z-plane. The linear progam (3.7) is then set up by a program written in Python and is subsequently solved with the optimizer CPLEX by IBM using the primal simplex algorithm. Since this optimizer is limited to double precision, it is important to reduce the spread in size of the numerical values in the problem. To this end, note that we can rescale each row of the inequality (3.7) separately by a positive number. Denoting a given row by R, we rescale its elements by Similarly, we can rescale each column of the matrix M separately by a positive number if we 19 Alternatively, we can use the recursion relation also in d = 4 dimensions by setting d = 4 + (to avoid double poles that appear at d = 4). However, Mathematica evaluates the closed-form expression faster than (our implementation of) the recursion relation and we therefore choose the former.
redefine the corresponding (squared) OPE coefficient in the vector ρ. We again choose and correspondingly for ρ. This procedure is iterated three times in our Python code, using precision arithmetric with 120 digits to ensure that no significant rounding errors are introduced in the process (the conformal blocks have been calculated with the same precision). Since we perform our own rescaling, we switch off this option in CPLEX. We find that the above rescaling typically reduces the orders of magnitude in the ratio between the largest and smallest numerical value in eq. (3.7) by about half. Nevertheless, precision is a limiting factor and does not allow us to go to cutoffs ∆ * much larger than 20. The fact that double precision is sufficent for smaller cutoffs, on the other hand, makes our calculations (combined with the excellent speed of CPLEX) very fast.

Conclusions
We have implemented the method proposed in ref. [1] to numerically study the bootstrap equations away from the symmetric point z =z = 1/2. Using this method, we have qualitatively reproduced various results that have been determined in the bootstrap literature using the more common method of taking derivatives at the symmetric point. The main aim of our work was to show that bootstrapping with multipoints works and is a valid alternative to the standard derivative method. In particular, it can be useful at a preliminary stage when one wants to qualitative bound or approximately compute some quantities using the bootstrap. By choosing a sufficiently low cutoff ∆ * , one can get qualitatively good results within seconds of CPU time with a standard laptop! Since the optimizer CPLEX that we use is limited to double precision, we can not achieve the high precision of refined bootstrap codes such as Juliboots [9] or SDPB [11]. Nevertheless we have shown how, using ∆-maximization, relatively precise results can be obtained for the scaling dimensions of operators (though we relied on an extrapolation procedure). In particular, for O(n) models with n = 2, 3, 4 we have determined the scaling dimensions of the second-lowest-dimensional operators S and T in the singlet and symmetric-traceless representation, respectively. To our knowledge, these have not been determined before using bootstrap techniques. We believe that it should not be difficult to go to arbitrary precision and get rid of the discretization (and the extrapolation procedure) by, for instance, adapting the algorithm developed in refs. [6,9] to multipoints. We do not exclude that bootstrapping with multipoints might then turn out to be comparable to (or better than) the derivative method for high-precision computations. From a conceptual point of view, the multipoint method is more rigorous, since the crossing equations are not truncated but bounded by an error. 20 We have also discussed how the multipoint method is useful in understanding to which extent a given numerical result depends sensitively on the high-dimensional operators. In particular, we have noticed that bounds on operator dimensions are less sensitive in this respect than bounds on the central charge.
Ideally, one might want to push the multipoint method to the extreme "IR limit", by choosing a cutoff ∆ * so low that an analytic approach may become possible. This is certainly a very interesting direction that should be explored. Among other things, it requires to improve on the estimate of the OPE convergence given in ref. [12] that applies in the opposite regime, for parametrically large ∆ * . Perhaps the results of ref. [25] might be useful in this respect. 21 An important line of development in the numerical bootstrap is the analysis of mixed correlators which so far are numerically accessible only using semi-definite programming [15]. It would be very interesting to implement mixed correlators in the multipoint bootstrap, either by adapting the semi-definite programming techniques or by extending the linear programming techniques.