Bootstrapping Heisenberg Magnets and their Cubic Instability

We study the critical $O(3)$ model using the numerical conformal bootstrap. In particular, we use a recently developed cutting-surface algorithm to efficiently map out the allowed space of CFT data from correlators involving the leading $O(3)$ singlet $s$, vector $\phi$, and rank-2 symmetric tensor $t$. We determine their scaling dimensions to be $(\Delta_{s}, \Delta_{\phi}, \Delta_{t}) = (0.518942(51), 1.59489(59), 1.20954(23))$, and also bound various OPE coefficients. We additionally introduce a new"tip-finding"algorithm to compute an upper bound on the leading rank-4 symmetric tensor $t_4$, which we find to be relevant with $\Delta_{t_4}<2.99056$. The conformal bootstrap thus provides a numerical proof that systems described by the critical $O(3)$ model, such as classical Heisenberg ferromagnets at the Curie transition, are unstable to cubic anisotropy.


Introduction
Numerical bootstrap methods [1,2] (see [3,4] for recent reviews) have led to powerful new results in the study of conformal field theories (CFTs). In [5,6] we developed an approach to largescale bootstrap problems which allowed for precise determinations of the CFT data of the 3d critical O(2) model. In this work, we continue the exploration of large-scale bootstrap problems by applying the technology introduced in [5] to the study of the 3d critical O(3) model.
Concretely, we apply these methods to study correlation functions of the lowest-dimension singlet, vector, and rank-2 scalars in the three-dimensional critical O(3) model. Using the "cutting surface" algorithm introduced in [5], we compute the allowed region for the CFT data of these leading scalar operators. Our results, together with comparisons to results from Monte Carlo simulations, are summarized in table 1. We also introduce a new algorithm and software implementation called tiptop, which allows us to efficiently test allowed gaps for other operators across this region. We use it to determine an upper bound on the dimension of the lowest-dimension rank-4 scalar.
The 3d O(3) model is a well-studied renormalization group (RG) fixed point, and its critical exponents have been computed using many methods, both theoretical and experimental. This model describes the critical behavior of isotropic magnets, such as the Curie transition in isotropic ferromagnets, and antiferromagnets at the Néel transition point. Moreover, since disorder corresponds to an irrelevant perturbation, 1 the model also describes isotropic magnets with quenched disorder.

One of the main open questions about the O(3) model is its stability under cubic deformations.
The majority of magnets present in nature are indeed not isotropic: this means that the microscopic Hamiltonian describing the system in the ultraviolet (UV) is not invariant under the full O(3) symmetry group but only under a discrete subgroup, such as the cubic symmetry group. This implies that additional terms will be generated at the microscopic level that are invariant under cubic symmetry but transform in a non-trivial representation of O (3). If any of those deformations turn out to be relevant, the O(3) fixed point would be unstable and could not be reached without further tuning in the UV theory. The attractive, stable, fixed point would instead be the 3d cubic model. Field theory computations and Monte Carlo simulations have shown that these two models have very similar critical exponents: hence, if the cubic perturbation is relevant, it should be very close to marginality and the RG flow connecting the two theories is very short. We will come back to this point in section 1.1.1.
We give a definite answer to the above question: the O(3) model is unstable under cubic deformations. This information is encoded in the dimension of the lowest rank-4 scalar t 4 , which in the O(3) model satisfies ∆ t 4 < 3. As we will discuss, this implies that the O(3) model is also unstable with respect to the biconal fixed point with Z 2 × O(2) symmetry. Relevance of t 4 has been previously suggested by Monte Carlo [9] and perturbative expansions [11], but the proximity to marginality and near degeneracy of the critical exponents between the cubic, biconal, and O(3) fixed points makes this a subtle question ideal for the precision and rigor of the conformal bootstrap.  We denote the leading rank-0, rank-1, rank-2, and rank-4 scalars by s, φ, t, t 4 , respectively. Bold uncertainties correspond to rigorous intervals from bootstrap bounds. Uncertainties marked with a * indicate that the value is estimated non-rigorously by sampling points.

Theoretical approaches to the 3d O(3) model
We start by briefly reviewing past approaches to the 3d O(3) model, including field theory studies, Monte Carlo, and past results obtained by conformal bootstrap techniques. We also describe related models and motivate the calculations in this work.
The simplest continuum field theory in the O(3) universality class is the theory of a scalar field φ transforming in the fundamental representation of O(3) with Lagrangian A large negative mass-squared for the scalar induces spontaneous symmetry breaking and leads to the ordered phase, while a large positive mass-squared leads to the disordered phase. The critical point is achieved by tuning the UV mass so that the infrared (IR) correlation length diverges. The β function of the coupling g has been computed in the ε−expansion and in a fixed-dimension scheme. After a Borel-resummation, both methods predict the existence of an IR stable fixed point. We will review these results in the next sections.
The IR limit of the above field theory captures the same physics as the Heisenberg model.
This model consists of a lattice of classical spins S i , which can take values on a three-dimensional sphere. The Hamiltonian has only nearest-neighbor interactions: where we also introduced an external magnetic field H in the third (z)-direction. When the parameter J is positive, the ground state corresponds to all spins aligned, corresponding to ferromagnets. When J < 0, the energy is minimized when neighboring spins are anti-aligned, corresponding to antiferromagnets.
For small J, the line H = 0 separates a ferromagnetic phase from the paramagnetic one. This line represents a first-order transition and terminates at a value J = J c , where the correlation length of the system diverges, and the transition becomes second order. For J > J c , there is only a disordered phase. At J = J c , the theory in the IR is in the same universality class of the field theory defined in (1). The critical exponents are related to operator dimensions at the fixed point as Here, s ∼ | φ| 2 denotes the lowest-dimension singlet scalar, while t ij ∼ (φ i φ j − trace) denotes the lowest rank-2 scalar. More generically the exponents Y r are associated to the dimensions of the lowest rank-r scalar operator. 2 In the O(N ) model, the dimension of the lowest traceless symmetric operator t describes the instability of the theory against anisotropic perturbations. Because of this, it plays an important role in the description of multicritical phenomena. For instance, the critical behavior near a bicritical point where two critical lines with O(n 1 ) and O(n 2 ) symmetry meet givies rise to a critical theory with enlarged O(n 1 + n 2 ) symmetry.
The Hamiltonian (2) is a simplified model of magnetic interactions, since in a real crystalline solid other interactions are present. For instance, the crystal lattice structure could give rise to magnetic anisotropy. In cubic-symmetric lattices this effect produces an interaction localized at each lattice point i of the form 3 k=1 (S k i ) 4 . This perturbation breaks the O(3) global symmetry of the Heisenberg Hamiltonian, and therefore it cannot be generated by an RG transformation. As such, the IR fixed point of (2) will be described by an O(3) invariant CFT.

O(N ) vs. multi-critical models
The O(3) model described above can be generalized to O(N ) by promoting φ to an N component field. We can also consider the closely related cubic model, which describes the continuum limit of the Hamiltonian (2) with the addition of the O(N ) breaking term N k=1 (S k i ) 4 . This interaction is indeed invariant under the symmetries of a hyper-cubic lattice, namely permutations and reflection of the three axes. The field φ i , i = 1, . . . , N , transforms in the fundamental representation of the permutation group S N . Moreover, each component is odd under a reflection of the corresponding axis. The composition of these transformations gives rise to the hypercubic symmetry group Compared to (1), the Lagrangian of the hypercubic model has an additional term in the potential: The computation of the two β-functions β g and β h reveals the existence of four fixed points [12]: the trivial fixed point (g = h = 0), the N decoupled copies of the Ising model (h = 0, g = 0), the O(N ) fixed point (g = 0, h = 0) and the cubic model (g = 0, h = 0). It is straightforward to see that the first two are unstable since the quartic operator parametrized by g is relevant in both theories. 3 Determining which one of the other two fixed points is stable is a more complicated issue, and it turns out to be N dependent.
One way to rephrase the above question is to notice that the additional term in (4) can be rewritten as where t ijkl 4 is the traceless symmetric combination of four fields. The added term in the potential, in O(N ) notation, can be written as a combination of a rank-4 field and a singlet. We know that the singlet is irrelevant at the O(N ) fixed point, by definition. Thus the stability of the O(N ) fixed point or the cubic point is linked to the value of the dimension of the operator t 4 .
In the O(2) model the operator t 4 is irrelevant. A simple proof of this is to notice that for N = 2, as long as h = 0, the cubic Lagrangian can be mapped in the Lagrangian of two decoupled Ising models. This cubic fixed point coincides with the decoupled Ising fixed point, which is unstable. Field theory and Monte Carlo determinations of the dimension of t 4 agree with this argument. This is also consistent with the assumptions made in [5]. On the contrary, at large N , the operator t 4 is relevant, and the cubic fixed point is stable. Thus, it is important to know at which value N = N c > 2 the operator t 4 becomes relevant.
A second closely related model is the multi-critical point with O(n 1 ) × O(n 2 ) symmetry [13]. A field theory description is given in terms of two sets of scalar fields φ 1 and φ 2 , transforming respectively in the fundamental representation of O(n 1 ) and O(n 2 ), with Lagrangian: where we have already set to zero all the mass terms. The analysis of the perturbative β functions shows the existence of six fixed points. Some we already know: the free one (g i = 0, h = 0), the two Wilson Fisher fixed points (g 1 = 0, g 2 = h = 0 and same with 1 ↔ 2), the decoupled fixed point (DFP, g i = 0, h = 0), the symmetry enhanced O(n 1 + n 2 ) Wilson Fisher fixed point, and lastly the biconal fixed point (BFP). The latter one also has all couplings nonvanishing, but the global symmetry is not enhanced.
The problem of understanding the stable fixed point can be again reduced to studying the (ir)relevance of given deformations in the various CFTs. For instance, by inspecting the dimension of the composite operator built out of the lowest dimension scalar singlets in O(n 1 ) and O(n 2 ) theories, one can conclude that the DFP is stable for any N = n 1 + n 2 ≥ 4. It is unstable for N ≤ 3, although the perturbation is close to being marginal. 4 The issue of stability of O(N ) vs. the BFP is again related to the dimension of a certain operator. In the Lagrangian formulation, this is a combination of quartic interactions. At the O(N ) fixed point this term is mapped in a combination of the second-lowest rank-0 (S ), secondlowest rank-2 (t 2 ) and leading rank-4 scalar operator t 4 . If any of these operators is relevant, then the O(N ) fixed point is unstable. While the former two are known to be always irrelevant for any N , the latter is the object of investigations. In particular, if ∆ t 4 < 3 for N = 3, then among the fixed points, the BFP will be the stable one.  We denote the leading rank-0, rank-1, rank-2, and rank-4 scalars by s, φ, t, t 4 , respectively. Another estimate of ∆ t in the fixed-dimensional expansion can be found in [19] in terms of the crossover We do not report it here because the errors depend on the value of ν used.
Both the O(3) model and the cubic model have been extensively studied using different expansion techniques. β-functions for these models are known up to high order in both the -expansion and fixed-dimension expansion, and critical exponents have been computed by Borel resumming the respective series. 5 We report in table 2 the latest results obtained with field theory techniques. 4 The most precise bootstrap determination [8,5] gives ∆ [s Z 2 s O(2) ] = ∆ s Z 2 + ∆ s O(2) = 2.92398(23) < 3. 5 Both approaches are based on a perturbative expansion in the quartic interaction g up to a certain loop order. In the fixed dimension approach one works directly in d = 3 dimension and looks for solutions of the Borel resummed β-function β BR (g * ) = 0. Critical exponents are then computed as ν BR (g * ). In order to remove the divergences one imposes suitable renormalization conditions. Historically the term "fixed dimension" refers to renormalization The question of stability of fixed points has also been discussed in the literature. As we discussed in the previous section, this question can be addressed in two ways: by computing the dimension of the lowest dimension rank-4 scalar in the O(3) model, or by computing the value N c at which the dimension of the second-lowest rank-0 scalar in the cubic model becomes exactly marginal. 6 Results from both methods support the conclusion that O(3) is unstable while the cubic model is stable. The formula for N c in the ε-expansion is [20,21]: and after resummation gives N c = 2.89 (2).
Analysis of the ε-expansion or fixed-dimension perturbative series in the cubic model [22,17,23] shows that the critical exponents of the two models are very close: 7 These differences are much smaller than the typical experimental error (e.g., [26]). This makes distinguishing the two models experimentally very challenging. Curiously, the first few terms in the of the ε-expansion of the critical exponents in (8) are quite different, and it is only after the Borel resummation that the two values appear quite close.
Similarly, also the biconal Z 2 × O(2) model and the O(3) critical exponents are very close, as the flow connecting the two is driven by the same almost marginal operator as in the cubic case: where η BFP and η BFP correspond to the two relevant order parameters charged respectively under Z 2 or O(2).

Monte Carlo results
Using Monte Carlo (MC) techniques, it is possible to obtain precise estimates of the critical exponents for both the O(3) model and the cubic model, as well as information about their stability. Such determinations can also be improved when combined with finite-size scaling (FSS) or hightemperature expansion (HT) methods. A precise determination of the ν and η critical exponents was made using MC and FSS methods in [27]. A more precise analysis combining MC with HT techniques was carried out in [28], while a more precise MC and FSS study was performed in [9]. A very precise MC and FSS analysis of an icosahedral model, as well as improved MC and HT analyses were recently presented in [7]. Several other less precise determinations can be found in [26]. The dimensions relevant to anisotropic perturbations of rank-2,3,4 were computed in [9,29] using MC and FSS methods, and support the conjecture that the O(3) model is unstable under cubic deformations. These results are summarized in table 1.
conditions at zero momentum; the use of a minimal subtraction scheme is instead called a "minimal subtraction scheme without ε-expansion". In the proper ε-expansion approach one works in d = 4 − ε dimensions and solves the condition β(g * ) = 0 order by order in ε. Plugging the solution g * (ε) into the expression for the critical exponent one gets a series in ε that can be Borel resummed. The final critical exponents are then computed as ν ε-BR (ε = 1). 6 The lowest dimension one corresponds to the mass deformation and is always relevant; the second-lowest corresponds to a combination of the two quartic interactions. The orthogonal combination is related to the Lagrangian operator via the equation of motion and is irrelevant. 7 Field theory estimates have also been obtained in [24,25] 1.

The conformal bootstrap
Three dimensional O(3) models have been studied with bootstrap methods in a series of papers [30,10,8], first by considering the correlation function φ i φ j φ k φ l , where φ i is the lowest-dimension scalar transforming in the vector representation of O(N ) , and then by also including correlation functions involving the lowest-dimension singlet scalar s. The most precise determination of the critical exponents was obtained in [8], which isolated a three dimensional region in the space {∆ φ , ∆ s , λ sss /λ φφs } = {0.51928(62), 1.5957(55), 1.205(9)}, under the assumption that φ i and s are the unique relevant scalar operators in their representations. In addition, by scanning over this island, [8] determined the magnitude of the leading OPE coefficient to be λ φφs = 0.5244 (11).
Theories invariant under the cubic symmetry group were also studied using bootstrap methods using single correlators [31,32], and mixed correlators [33,34]. In particular [32] analyzed the bootstrap equations assuming the hypercubic symmetry group C N = Z N S N and observed a series of kinks for various values of N . However, the locations of the kinks in the singlet sector were degenerate with the O(N ) kinks (and hence compatible with (8)), likely reflecting a symmetry enhancement in the extremal bootstrap solutions [35][36][37][38][39], while the bounds in other sectors did not seem to be saturated by the cubic model. The mixed-correlator analysis of [33,34] also did not manage to isolate the cubic model but rather found evidence of a new theory, called the "Platonic CFT," with cubic symmetry and operator dimensions not matching any known CFT.
In this work, we study the O(3) model with numerical bootstrap techniques using a larger system of correlation functions than before: in addition to φ i and s, we incorporate the lowestdimension rank-2 scalar t ij ∼ φ (i φ j) . This setup is similar to the one leading to the successful results obtained in [5] for the O(2) model. Following the strategy detailed in [5], we first scan over the three operator dimensions {∆ φ , ∆ s , ∆ t } and the OPE coefficients {λ sss , λ φφs , λ tts , λ φφt , λ ttt } (or more precisely their ratios) and we determine a three dimensional island in the space of operator dimensions, along with an associated allowed set of OPE coefficient ratios. Next, we compute upper and lower bounds on the magnitude λ φφs , as well as on the current and stress-tensor central charges C J and C T . Finally, we enlarge the parameter space to include one more parameter: the dimension of the lowest rank-4 scalar ∆ t 4 . Using the new tiptop algorithm, which we describe in section 3, we carve out the allowed region in the enlarged four-dimensional space and obtain an upper bound on ∆ t 4 .

Structure of this work
The remainder of this work is structured as follows. In section 2 we describe the crossing equations and relevant O(3) representation theory. In section 3 we describe the new tiptop algorithm that we use in order to bound ∆ t 4 . In section 4 we describe the results of our numerical bootstrap calculations and in section 5 we describe directions for future research. Various appendices describe the code availability, software setup, details about our tensor structures, and give a list of allowed and disallowed points that we have computed.

The O(3) model 2.1 Crossing equations
We begin by describing the representation theory of O(3) = Z 2 × SO(3). We label the irreducible representations q ± of O(3) by the usual SO(3) rank q tensor of dimension 2q + 1 for q ∈ 1 2 Z ≥0 , as well as the Z 2 parity ±. Tensor products of these irreps are given by where if q 1 ± = q 2 ± , then the even/odd q a are in the symmetric/antisymmetric part of the tensor product.
Operators O q ± (x) in the irrep q ± can be written in terms of SO(3) fundamental indices i = 1, 2, 3 as rank-q symmetric traceless tensors O where ∆ ij ≡ ∆ i − ∆ j , the conformal cross ratios u, v are and the operators O that appear both OPEs ϕ 1 × ϕ 2 and ϕ 3 × ϕ 4 have scaling dimension ∆, spin , and transform in an irrep R that appears in both the tensor products R 1 ⊗ R 2 and R 3 ⊗ R 4 . For each R, the SO(3) structure T R can be constructed using the SO(3) Casimir and normalized to give consistent OPE coefficients under crossing using the free theory as described in appendix C. The Z 2 irrep of O follows from trivial multiplication of ± 1 and ± 2 , and so does not require a structure. If ϕ 1 = ϕ 2 (or ϕ 3 = ϕ 4 ), then Bose symmetry requires that O have only even/odd for R in the symmetric/antisymmetric product of We are interested in four-point functions of the lowest dimension scalar operators transforming in the 0 + , 1 − , and 2 + representations, which we will denote following [30,10,8] as s, φ, and t, respectively. 9 These operators are normalized via their two point functions as Correlator s-channel t-channel Eqs Table 3: Four-point function configurations that give independent crossing equations under equating their s-and t-channels, along with whether even or odd spins ± appear for each irrep in each channel, and the number of crossing equations that each configuration yields.
where x 12 ≡ |x 1 − x 2 | and all indices of the same letter should be symmetrized with their trace removed. In table 3 we list the 4-point functions of s, φ, and t that are allowed by O(3) symmetry 10 and whose s and t-channel configurations lead to independent crossing equations, along with the irreps and spins of the operators that appear in the OPE, and the number of crossing equations that they yield. These 4-point functions can be written explicitly as in (11), where the explicit SO (3) structures T R are computed in appendix C. Equating each of these s-channel 4-point functions with their respective t-channels yields the crossing equations where ± denotes which spins appear, and the V 's are 28-dimensional vectors of matrix or scalar crossing equations that are ordered as table 3 and written in terms of The explicit form of the V 's is given in the attached Mathematica notebook. 11

Ward identities
The OPE coefficients of J µ and T µν are constrained by Ward identities in terms of the two-point coefficients C J and C T . In our conventions, we have where C free J,T are the two-point coefficients of J and T in the free O(3) model described in appendix C. Thus, the contribution of these operators to the crossing equation can be parametrized purely in terms of C T and C J , together with the dimensions and charges of the external scalars φ, s, t.

The tiptop algorithm
While our primary search for the O(3) bootstrap island will follow the same methods and software tools used for the O(2) model described in [5], we will also need to compute the maximum value of the scaling dimension ∆ t 4 over this island. This employs a new search strategy and software implementation that we describe in this section.

Software and algorithm
tiptop is a program to assist in finding the maximum value of a coordinate achieved in a region in (N + 1)-dimensional space, where testing whether a point is in the region is computationally expensive. Given a set of inside (allowed), outside (disallowed), and unknown points, tiptop generates successive points to narrow down the boundary of the top of the region. It is meant to be invoked by a driver that takes these points, computes whether they are allowed, and then asks for more points to check. tiptop is freely available (see Appendix A).
The number of dimensions is arbitrary but fixed at compile time. For concreteness and ease of visualization, we assume that N + 1 = 3 for the rest of this discussion, where the dimensions are ∆ φ , ∆ s , and ∆ gap . The algorithm operates unchanged for higher dimensions.
We start with at least one allowed point, a cloud of disallowed points, a cloud of points that are in-progress, and a maximum gap (∆ max gap ). In-progress points are points that the driver already knows about and is working on, but does not yet know if they are allowed. For example, those calculations may have been submitted as calculations to an HPC cluster but not yet completed.
We assume that there are no allowed points with ∆ gap ≥ ∆ max gap . We also assume that islands only shrink at larger gaps. This allows us to assert that if a point is disallowed at one gap, it will continue to be disallowed at larger gaps.
The last assumption is that each N -dimensional island at a fixed value of ∆ gap is convex and simply connected, so each island never becomes a horseshoe or splits into two pieces. We have observed this behavior for a wide variety of theories, as long as the theory is well approximated.
The basic outline of the algorithm for generating points is: • Set ∆ allowed to the largest ∆ gap with an allowed point.
• Explore parameters at ∆ allowed to find the size of the island there. If there are any corners of parameter space left to map out, return one point from there (Section 3.2).
• If the island at ∆ allowed is thoroughly mapped out, generate one point at a higher gap (Section 3.3).
The non-gap dimensions (∆ φ , ∆ s ) are represented as regular floating-point numbers, while the gap dimension (∆ gap ) has been rescaled to a 64 bit integer. This reduces numerical errors where two points at very similar gaps are mistakenly considered to be at the same gap.
tiptop will not return a point in two cases: • The current gap ∆ allowed might be fully explored, but it needs to know the outcome of some in-progress points to be sure.
• There are no valid larger gaps left. For example, consider the case where ∆ allowed = 10000 and tiptop has ruled out any jumps to ∆ gap = 10001. There are no integers between 10000 and 10001, so the algorithm terminates. Figure 1: The max coordinates ∆ max for a collection of allowed, disallowed, and inprogress points.

Rescaling
The islands often have extreme aspect ratios in the 'natural' coordinates. This causes difficulties when exploring an island, so tiptop rescales the coordinates. The first step in rescaling is to get an overall scale for all of the points (allowed, disallowed, and in-progress) from all gaps. We define ∆ max as a scalar equal to the largest absolute coordinate value in all dimensions, as shown in figure  1.
For a given ∆ allowed , we define ∆ previous as the largest gap with allowed points but less than ∆ allowed . This is usually a previous value for ∆ allowed . Figure 2 shows an example of allowed points at ∆ allowed , ∆ previous , and lower gaps.
Using the m allowed points at ∆ previous , we scale the points using a principle component analysis. Specifically, we construct the matrix We then compute the singular value decomposition (SVD) of this matrix where Σ is a rectangular m × N diagonal matrix with non-negative real numbers σ i = Σ ii on the diagonal ranging from the smallest (σ min ) to the largest (σ max ). U is an m × m unitary matrix, and V is an N × N (here 2 × 2) unitary matrix.
We define the N × N matrix Ω as the first N rows of Σ. This is a diagonal matrix with the entries σ i , so the inverse is trivial. Putting this all together, we define the rescaling matrix It may be that there are so few points at ∆ previous that they are not linearly independent. For example, in the beginning, there may not be any points at ∆ previous . If the ratio between the smallest (σ min ) and largest (σ max ) of these singular values is less than a tolerance (we use 10 −8 ), then we only scale by ∆ max where I is the identity.
Everything is scaled by the largest coordinate value ∆ max to guarantee that all points are mapped into a box with extents (-1,1) in every dimension. The factor of 1.75 (about √ 3) is to ensure that all points will fit into the unit box even after rotation.
The transformation has the effect of a rotation and then rescaling of the rotated coordinates, so the allowed region remains convex. However, the allowed points should outline a more circular shape than the extended ellipse we started with, as shown in figure 3.
One concern with this rescaling algorithm is that it weighs dense regions with more points more than equivalently sized regions with fewer points. So it may not produce an optimal transform. In practice, the later steps spread out the points very evenly, so this concern turns out not to be a problem in practice.

Adaptively meshing the box
While the distribution of points in figure 3 no longer has extreme aspect ratios, the points are still clustered in a small region of the unit box.
Based on the assumption that the allowed island only shrinks as the gap increases, we now only consider three sets of points: allowed at the current ∆ allowed , disallowed at ∆ gap ≤ ∆ allowed , and in-progress. For the rest of this step, we will be treating in-progress and disallowed points identically.
The strategy is to place points in regions that are empty. To quantify this emptiness, we create a hierarchy of meshes covering all the points. Empty regions are then cells that have no points.
The mesh hierarchy starts with a single coarse cell covering the entire unit box. The next level has 2 N cells (4 for our example), the level after that has 2 2N cells, and so on. This subdivision continues up to our predefined limit of 47 levels. This implies a relative minimum cell width of 2 −47 . This is quite small, but still significantly larger than the minimum resolution of an IEEE-754 double-precision number (2 −53 ). This helps reduce errors as points get rotated and scaled.
Next, we compute the coordinate extents of the allowed points in the transformed frame. For example, the extents in the transformed horizontal direction is the difference between the smallest and largest horizontal coordinates for the allowed points. This gives us a list of N numbers. In figure 3, these extents are about 0.2 in both vertical and horizontal directions.
We then define a minimum cell size as the minimum of these N extents multiplied by a fixed fraction f cutoff . We use f cutoff = 2, which is deliberately very coarse. This tends to double the size of the allowed region at each step. Also, if f cutoff is too small, then the algorithm will completely fill in internal regions, even though, by assumption, the internal spaces do not need to be checked. This minimum cell size corresponds to a level l max in the hierarchy of meshes.
Just after a jump to a higher gap, there is only one allowed point at ∆ allowed . In this case, the extents are zero, so we set l max = 47, the finest level.
Restricting our search to cells at level l < l max , we look for any empty cells adjacent to the allowed points. There will, in general, be multiple empty cells at multiple levels. We choose the largest empty cell. We only check diagonals, so points get laid out in a checkerboard pattern as in figure 4.
The new point is not placed at the center of the new cell, but rather simply offset from the existing allowed point. So if the allowed point is in a corner of a cell, the new point will be in the same corner of the empty cell.
In practice, the implementation does not explicitly create the mesh at all levels. Rather, the points are stored in a tree. A node in the tree can have up to 2 N leaves, but leaves are only created if there is a point in that leaf. Adding a point to the tree adds it at all levels, so that it is easy to determine if a cell is occupied at any level.
The observed behavior of this algorithm is that it quickly finds a rough estimate for the boundary between allowed and disallowed, but can spend a lot of effort finding the exact boundaries. In-progress points are treated as disallowed, so too many in-progress points will lead to extra work. In practice, we have up to 16 points in-progress at any one time.

Jumping to a larger gap
If the previous section does not yield a new point, and there are no in-progress points at that gap, then we try to jump to a larger gap.
We start by rescaling the points as in section 3.2.1. We draw a coordinate box around all of the points allowed at ∆ allowed and shrink it by a factor of 2. Then we find the largest gap ∆ ceiling that can accommodate this box without containing any disallowed or in-progress points with ∆ gap ≤ ∆ ceiling . At the beginning, there are no disallowed points at large gaps, so ∆ ceiling = ∆ max gap .
Defining ∆ diff ≡ ∆ ceiling − ∆ allowed , we return the center of the box at ∆ gap = ∆ allowed + ∆ diff /2, thus bisecting the range of allowed gaps. This underscores the need for a good estimate of ∆ max gap . If the estimate is too high, then the algorithm will recommend too many points that are far too large.
One thing to note is that when running with multiple in-progress points, each subsequent point will be at a smaller ∆ gap . So if there are 3 points running concurrently, they will be placed at ∆ allowed + ∆ diff /2, ∆ allowed + ∆ diff /4, and ∆ allowed + ∆ diff /8. If the first point at ∆ allowed + ∆ diff /2 succeeds, then any effort towards verifyinig the latter two points will be wasted.
This partition procedure means that a rough estimate for the upper bound is ∆ ceiling . This assumes both that the island itself is convex and that the allowed region near the tip of the allowed region is convex, as in Figure 7.
Overall, we have found this approach to work reasonably well. More importantly, it is very robust. It is very easy to be too clever, resulting in odd failures.

Dimension bounds with OPE scans
Next we present our conformal bootstrap island computed using SDPB [41,42], along with its comparison with various Monte Carlo results. Computing the conformal bootstrap island requires scanning over the three operator dimensions {∆ φ , ∆ s , ∆ t } using the Delaunay search algorithm described in [5], and for each point using the "cutting surface" algorithm presented in [5] to decide if there exists an allowed point in the space of OPE coefficient ratios { λsss λ φφs , λtts λ φφs , λ φφt λ φφs , λttt λ φφs }. When computing the island we make the following assumptions about the spectrum unless stated otherwise. We assume that φ, s, and t are the only relevant operators in their respective symmetry representations, so that ∆ φ ,s ,t ≥ 3. In addition, we assume that the leading rank-4 scalar has a dimension satisfying ∆ t 4 ≥ 2. We assume an O(3) current with ∆ J = 2 and stress tensor with ∆ T = 3, with coefficients satisfying the Ward identity constraints. We also impose a twist gap above them, as well as in all other sectors not mentioned above, of size 10 −6 .
In figure 5 we show the conformal bootstrap island we have computed at Λ = 43 using these assumptions, compared to the Monte Carlo results of [7,9]. In figure 6 we show various 2d projections of the bootstrap island. In appendix D we give the full set of allowed and disallowed points we computed at Λ = 43, along with figure 9 which shows the convergence of the allowed points as a function of Λ after performing an affine transformation to make the allowed regions roughly spherical.
In these plots we show our best determination of the allowed region at a given Λ, constructed by computing a Delaunay triangulation of the tested points, choosing triangles that contain both allowed and disallowed points, and plotting the convex hull of the points that are midway between the allowed and disallowed vertices in these triangles. At Λ = 43, this "best-fit" region gives A more rigorous determination can be made by taking the convex hull of the disallowed points in these boundary Delaunay triangles. This region gives the rigorous error bars which we have quoted in table 1.
The allowed points at Λ = 43 are associated with OPE coefficient ratios which live in the ranges 12 These should be viewed as an approximation to the full allowed region of OPE coefficients, which may be slightly larger.

Central charges and λ φφs
Next, we compute upper and lower bounds on the magnitude of the OPE coefficient λ φφs , the central charge C T , and the current central charge C J . We compute these bounds over a small sample of points in our allowed region so the results will be inherently non-rigorous. However, we believe that this treatment gives reasonable estimates for these quantities that are more precise than previous results.
The strategy is similar to the method we employed in [5]. We take seven primal points in the Λ = 43 island, consisting of the scaling dimensions and allowed OPE coefficients. The points are chosen to be sufficiently symmetrized and sparse across the Λ = 43 island we have computed. For each of these points, we extremize C T , C J , and the external OPE norm parameterized by λ φφs , to obtain upper and lower bounds. This calculation was limited to Λ = 35 due to our available computational resources. The data points and SDPB parameters we used are summarized in tables 8 and 5, respectively.
There is an important comment we want to make about the upper bound computation on C T and C J (a similar comment was made in [5]). For computing upper bounds on C T and C J , we have to assume a gap ∆ T /J above the unitarity bound for the next operators in the T or J sectors. Note that this gap was not assumed in our OPE scan, so this extra constraint might turn an allowed point into a disallowed point. If we do not have such a gap, the upper bound is loose and may not give reasonable results. On the other hand, large gaps can make SDPB unable to find a solution.
In table 8, we summarize the gaps ∆ T /J we assume in the upper bound calculations. From spectrum determinations using the extremal functional method (see [43,44]), we have noticed that a gap ∆ T /J = ∆ T /J + 1 above T and J is generally favored. We were able to compute bounds with this gap for three of the points, but for the other four we could not find solutions. For those points, we adopted the weaker assumption ∆ ext,T (J) = ∆ T /J + 0.1.
These results agree with and are more precise than previous determinations of these quantitites (see [30,10,8]).

Upper bound on ∆ t 4
Our last result is the maximum value of the rank-4 scalar dimension ∆ t 4 . In conjunction with the tiptop algorithm described in section 3, we computed points at Λ = 19, 27, and 35. Allowed points at lower values of Λ were used to initiate the search at larger values of Λ. Figure 7 shows a projection of a subset of the 1311 disallowed points and 172 allowed points at Λ = 35, and Figure  8 shows how the island shrinks as we approach the maximum ∆ t 4 .
The largest allowed value of ∆ t 4 was a single point at 2.99052 with the scaling dimensions {∆ φ , ∆ s , ∆ t } = {0.518962, 1.59527, 1.20969}. The limit from ∆ ceiling (see Section 3.3) is This implies that the leading rank-4 tensor in the critical O(3) model is relevant, in agreement with other studies. We have superimposed a convex hull encompassing the allowed points on top, obscuring some of the disallowed points. We can see the behaviour of the tiptop algorithm, exploring the island at one ∆ t 4 before jumping to a larger ∆ t 4 . The jumps become progressively smaller, indicating convergence. We computed 16 points simultaneously, and this calculation took several months during which the tiptop algorithm was being developed. So the points reflect occasional crashes and small inefficiencies in the set of computed points.

Future directions
In this work we have applied the methods developed in [5,6]  3.01. Because this term is marginally irrelevant with δ = ∆ X − 3 0.01, if we want to reach the cubic fixed point by a Monte Carlo simulation, the size of the lattice has to be around the order of 2 1/δ , which is impractical to implement. An alternative way to estimate the cubic CFT data is using conformal perturbation theory. We start with the perturbed action S = S O(3) + g d 3 xX. Using the formalism in [45], one finds the beta function to be The dimension of an operator O at the cubic fixed point is then given at linear order in δ by CFT. Specifically one obtains ∆ X = ∆ 0 + 2δ, which justifies the estimate ∆ X 3.01.
The OPE coefficient λ XXX is proportional to λ t 4 t 4 t 4 . Unfortunately, using the setup of the present paper, we don't have access to λ t 4 t 4 t 4 . To access λ t 4 t 4 t 4 , one needs to bootstrap all fourpoint functions involving {φ, s, t, t 4 }, which is a concrete task for future research. Here we can estimate that the correction to ∆ t in the cubic fixed-point is of order δ = 0.01. On the other hand, the corrections to ∆ φ , ∆ s start at order δ 2 0.0001 since λ φφt 4 = λ sst 4 = 0. Note that in this work the error bar for ∆ t is much smaller than δ. Therefore a careful study of the {φ, s, t, t 4 } system should yield a solid prediction for the correction to ∆ t in the cubic model. Of course, it will also be interesting to understand how to isolate the cubic fixed point more directly using the conformal bootstrap, perhaps using a larger system of correlators than was considered in [31,33,34]. One can also straightforwardly apply the large-scale bootstrap techniques we have developed to other O(N ) models, as well as to 3d CFTs with fermions (using the newly developed software [46]) or to study conserved currents [47][48][49]. Using these methods one can also continue exploring larger systems of correlators that may help us to isolate CFTs containing gauge fields, such as 3d QED [50,51] and 4d QCD. Now that we have precisely isolated the O(3) model, we are also in position to do a more detailed study of its low-twist trajectories of operators as a function of spin, which can be compared to analytical calculations using the Lorentzian Inversion formula [52,53], following the approach of [54,6,55]. Such analytical techniques can also be used to estimate the leading Regge intercepts and related Lorentzian data of the O(3) model. In future work it will also be important to understand how to better incorporate insights from the analytical bootstrap, such as our precise understanding of the large spin asymptotics, into making large-scale numerical methods even more powerful.
The SDPB parameters used in our computations are given in tables 4 and 5.

C Tensor structures
In this appendix we compute the SO(3) tensor structures T R that appear in the conformal block expansions (11) for the 4-point functions listed in table 3. We start by defining a basis of tensors for each configuration in table 3: where the indices for each of the four operators are labeled as i, j, k, l respectively, all indices with the same letter should be symmetrized with the trace removed, and for simplicity we suppress the indices on the left-hand side. For the first four configurations with non-trivial bases, we can find the tensor structure using the rank-2 SO(3) Casimir C acting on a basis B with n, m number of i, j indices, respectively, as: where G are the usual SO(3) generators. The K eigenvectors (T K ) J of M IJ are eigenvectors of C: where the eigenvalue c for a rank q SO(3) irrep is q(q + 1), which allows us to identify each T K with an irrep. Up to an overall normalization, these T K are then the desired tensor structures. For the last 8 configurations there is only one basis element, so the tensor structure is simply that element also up to an overall normalization. The final list of tensor structures is then φφφφ : tφtφ , φttφ : ssss : T 0 + 0 + 0 + 0 + 0 + = B ssss , φsφs , sφφs : T 1 − 1 − 0 + 1 − 2 + = B φsφs , tsts , stts : T 2 + 2 + 0 + 2 + 0 + = B tsts , ttss : T 0 + 2 + 2 + 0 + 0 + = B ttss , φφss : sttt , ttst : The overall normalization of these tensor structures has been chosen so that the OPE coefficients λ ϕ 1 ϕ 2 O and λ ϕ 3 ϕ 4 O in (11) are consistent under permutation of their subscripts. This can be checked using the free theory, where we have the operators which have been normalized consistent with the 2-point function normalization in (13). We can then compute all the 4-point functions in  in (11) using the tensor structures in (32) to verify this consistency. 13

D Computed points
In In figure 9 we show a plot of the allowed regions at Λ = 19, 27, 35, 43 after performing an affine transformation which makes the Λ = 19 region roughly spherical. The precise affine transformation is given by:  Figure 9: The convex hulls of the allowed points in the affine space (35) Table 8: Allowed points in the Λ = 43 island used to obtain bounds on λ φφs , C T , and C J , along with the gaps above ∆ T and ∆ J that were assumed.