A Monte Carlo approach to the conformal bootstrap

We introduce an approach to find approximate numerical solutions of truncated bootstrap equations for Conformal Field Theories (CFTs) in arbitrary dimensions. The method is based on a stochastic search via a Metropolis algorithm guided by an action $S$ which is the logarithm of the truncated bootstrap equations for a single scalar field correlator. While numerical conformal bootstrap methods based on semi-definite programming put rigorous exclusion bounds on CFTs, this method looks for approximate solutions, which correspond to local minima of $S$, when present, and can be even far from the extremality region. By this protocol we find that if no constraint on the operator scaling dimensions is imposed, $S$ has a single minimum, corresponding to the Free Theory. If we fix the external operator dimension, however, we encounter minima that can be studied with our approach. Imposing a conserved stress-tensor, a $\mathbf{Z}_2$ symmetry and one relevant scalar, we identify two regions where local minima of $S$ are present. When projected in the $(\Delta_\sigma, \Delta_{\epsilon})$-plane, $\sigma$ and $\epsilon$ being the external and the lightest exchanged operators, one of these regions essentially coincides with the extremality line found in previous bootstrap studies. The other region is along the generalized free theories in $d = 2$ and below that in both $d = 3$ and $d = 4$. We empirically prove that some of the minima found are associated to known theories, including the $2d$ and $3d$ Ising theories and the $2d$ Yang-Lee model.


Introduction
Starting from the pioneering work [1], numerical bootstrap methods have become a widespread work-tool for studying conformal field theories (CFTs) in d > 2 space-time dimensions, see [2] for a review with an extensive list of references. The most established and developed methods make use of semidefinite programming algorithms applied to 4-point functions of local operators to rule in or rule out theories, the state of the art program being currently SDPB 2.0 [3,4]. In this approach one makes assumptions on the CFT spectrum (typically on the low-lying scaling dimension of primary operators) and is able to rigorously establish if that assumption is allowed or not in the space of all possible CFTs. By imposing suitable assumptions to a set of 4-point functions one can restrict the allowed region to a small island in parameter space, so that precise and rigorous properties of theories living at the extremality region can be determined [5][6][7]. Such CFTs at the edge of the extremality bounds can also be studied using the so called extremal functional method, first discussed in [8] and later developed in [9,10]. Although not as rigorous as the previous method, the extremal functional allows a quite accurate and precise determination of the CFT data, see e.g. [11] for an application to the 3d Ising model. These powerful methods generally require that the theory in question lies at the boundary of the allowed region in parameter space to start with. In contrast, the "interior" of the space of allowed CFTs is harder to characterize using the current techniques.
In fact, it could be useful to have some numerical tool to explore the structure of the space of allowed CFTs and, in principle, a way to construct approximate CFT data that are consistent with crossing symmetry of 4-point correlation function of local operators. 1 Locating, even approximately, putative CFTs could be a useful first step to understand which assumptions to impose so that these theories become extremal and hence amenable to, e.g., a fully rigorous bootstrap analysis. The aim of this paper is to propose a novel numerical method to analyze the bootstrap equations and possibly find a subset of its approximate solutions. More specifically, we will try to answer to the following question: Given a finite set of n i operators with spin l i below a given scaling dimension ∆ * , exchanged in some 4-point correlation function, is there a set of CFT data (scaling dimensions and OPE coefficients) for the n i operators which best satisfy the truncated crossing equations? If so, where is the putative solution located in the allowed CFT parameter space?
In one-dimensional systems, where there is no spin, accurate solutions can be constructed [13]. For d > 1 this problem is much harder because one has to look for solutions by choosing among the many possible partitions of the operators in the different spin channels. In addition to that, it is unclear how to choose the initial conditions for any search and which constraints should be imposed to possibly get a finite set of solutions. Finally, it is far from clear that the search would lead to an actual theory. In particular, deterministic methods such as gradient descent, would uniquely give an answer, but most likely this answer won't correspond to an approximate solution to a physical CFT. This is easily understood as follows. Suppose we find a suitable set of constraints which guarantee the existence of a finite set of approximate solutions (we will see below that these are indeed necessary, but easy to identify). If we define a positive definite function F as, say, the modulus square of the conformal bootstrap equations for a single 4-point correlator, we see that in the infinite dimensional space of all possible CFT data, scaling dimensions ∆ i and Operator Product Expansion (OPE) coefficients λ ijk , generally F ≥ 0. Putative physical theories can only be associated to points where F = 0, where crossing symmetry is realized. 2 Suppose now to use a deterministic numerical method that allows us, starting from an arbitrary point in the (truncated) CFT data parameter space where F > 0, to minimize F . The details of how to actually achieve this do not matter for the argument we want to make here. In a truncated numerical approach, of course, F will never be zero, but we will nevertheless uniquely get the minimum value F 0 , given the chosen initial conditions. It is clear that generally F 0 will be an approximation to a local minimum of F (with F min > 0) and not an approximation to one of the global minima of F (where F = 0). Thus, the solution found will not correspond to any CFT; it will be in the "swampland".
In order to alleviate the problem of being trapped at fake CFT minima, we use a stochastic minimization via a Metropolis algorithm. The function S we choose to minimze (the "action") is the logarithm of the truncated bootstrap equations for a single scalar field correlator, S ∼ log F . Configurations with higher actions are from time to time accepted with a probability which depends on a parameter T ("temperature"). In this way, thermal fluctations allows us to "escape" from local minima in the quest for the global ones. 3 The problem of choosing among the many possible partitions of the operators in the different spin channels is solved by adding as many operators as possible to each spin channel and then allowing operators to possibly decouple along the minimization process. The complex dependence of S on the truncated CFT data requires further non-trivial numerical manipulations, which will be discussed in section 2, to finally determine the approximate CFT data of a putative theory.
We consider in this work the simplest bootstrap set-up, namely a single 4-point function of identical scalars in euclidean CFTs and focus on the exploration of the CFT parameter space where the external scalar is close to the unitarity bound. We start in section 2 by presenting the method. In section 2.1 we explain how the bootstrap equations are truncated and discretized, in section 2.2 how the stochastic minimization is defined and performed and finally in section 2.3 the various steps entering in the minimization process are reported and explained.
The results of our numerical explorations, all based on severe truncations of the bootstrap equations with ∆ * 15 and a total number of operators included N Ops 20, are reported in section 3. The first general qualitative feature that emerged from our analysis is discussed in section 3.1 and is the predominance of the free theory scalar theory. If no constraints are imposed -aside unitarity, the presence of an energy momentum tensor and a Z 2 global symmetry -the minimization of S leads only to a single global minimum, associated to the free theory. It is enough to impose a single constraint, such as fixing the scaling dimension of the external operator, to get a finite set of minima. At fixed ∆ σ we first study in section 3.2 the simplest spin partition with one operator per spin. We test the method by rediscovering the free theory and give indications that there are no other CFTs with one operator per spin, at least in their lowest spin channels and lowest lying operators. Our most interesting results are contained in section 3.3, where we study the more realistic set-up of more operators per spin. Now the situation is reversed and we get generally too many minima, unless extra assumptions are imposed. Demanding the presence of only one Z 2 -even relevant scalar drastically reduces the number of minima and makes the analysis feasible. A general qualitative feature that emerged from our analysis is the presence of two regions in the space of CFT data where the bootstrap equations are more easily satisfied for CFTs admitting a conserved energy-momentum tensor, a Z 2 global symmetry and one Z 2even relevant scalar. When projected in the (∆ σ , ∆ )-plane, σ and being the external and 3 Very recently numerical methods based on reinforcement learning have been used to find approximate solutions to the bootstrap equations [14,15]. Though this is an interesting line of research to pursue, it is not clear to us if and how such techniques can alleviate the problem of fake local minima. the lightest exchanged operator, respectively, one of these regions essentially coincides with the extremality line found in previous bootstrap studies. The other region is along the generalized free theories (GFTs) in d = 2 and below the GFT line in both d = 3 and d = 4. The minima found in the (∆ σ , ∆ )-plane are reported in figures 5, 9 and 12 for the d = 2, 3, 4 cases. The extremal minima in the d = 2 case are identified with the generalized minimal models of [16], one extremal minimum in d = 3 corresponds to the 3d Ising model, while all the other minima have not been identified. Unfortunately the numerical accuracy we have does not allow to reach higher values of ∆ * to establish more firmly if (and which of) these minima are numerical artifact or not. CFTs without an energy momentum tensor are considered in section 3.4. The minima now cluster around the GFT line and GFTs themselves are well reproduced, especially in d = 2 where our numerical accuracy is higher. We finally briefly show in section 3.5 how our algorithm can be applied to non-unitary theories by rediscovering the 2d Yang-Lee model.
We discuss in the final outlook the limitations of our algorithm and possible ways to improve it in the future. Three appendices contain some other details of our numerical procedure and further more technical results.

Method
We present in this section the numerical method we used to find approximate solutions of the bootstrap equations for a 4-point function of identical scalar primary operators in a generic CFT in d dimensions.

Discretization of the Bootstrap Equations
The crossing equations for 4 identical scalar fields σ with scaling dimension ∆ σ can be written as (notation as in [2]) ∆>0, Here λ σσO ∆ are OPE coefficients, ∆ and are the scaling dimensions and spin of the exchanged primary operator O ∆ (even spin traceless symmetric tensors), z andz are the usual Dolan-Osborn conformal invariant cross ratios [17,18], I(z,z) = |z| 2∆σ − |1 − z| 2∆σ is the identity contribution to the OPE expansion, and finally with g ∆, (z,z) the conformal blocks, normalized as in [1] with an extra 2 factor. As discussed in the introduction, we do not look for rigorous disallowed regions in the space of CFT data, but rather for approximate solutions of the above equations, namely for approximate possibly allowed CFT data. To this aim, we severely truncate (2.1) to the finite set of primary operators with scaling dimension up to a given ∆ * (chosen to be integer), and assign a given number of operators n to each spin channel up to some value max . If unitarity is imposed, ∆ * and max are related. We take so that the scaling dimension of the operators with spin max have a large enough range where they can be varied, compatibly with unitarity: ∆ * − 3 ≤ ∆ max ≤ ∆ * . The total number of exchanged operators O a that we consider is We then choose a finite sample of N z points in the complex z-plane in the neighborhood of z =z = 1/2 and define the N z × N Ops matrix M and the N z vector I i We define an "action" as (spin dependence omitted for simplicity) where ρ a ≡ λ 2 σσOa , by taking a weighted sum of the absolute square of the truncated crossing equations (2.1) at the N z points, with (2.7) The factor (2.7) is introduced to take into account that the points with the best convergent conformal block expansion are around z = 1/2, point where on the other hand the function F vanishes. For any finite truncation, our aim would be then to minimize S in the space of CFT data.

Metropolis Montecarlo and stochastic minimization
Finding the 2N Ops CFT data (ρ a , ∆ a ) that minimize S is a challenging task. First of all, in d > 1 dimensions, where we have spin, one should in principle look for solutions in all possible partitions of the n operators in each spin channel. In addition to that, deterministic methods, such as gradient descent, can be highly inefficient because they converge to a minimum of S depending on the initial conditions. The actual dependence of S on the CFT data is expectedly complicated and depends on the constraints imposed, and can be characterized by the presence of several local minima. One of the main motivations of this work is in fact to shed some light on the "landscape/swampland" space of CFT data. As we will explain in what follows, both problems (avoiding to land on a local minimum and the need to consider all spin partitions) are significantly alleviated by the use of stochastic searches based on Monte Carlo (MC) methods.
We will in particular sample the parameter space via the Metropolis Monte Carlo method, whose structure we now briefly describe.
Given a function f (x) : U → R, where U ⊂ R n , the minima of f in U are sampled by taking an initial point x (0) ∈ U and proposing an update which is accepted with a probability given by The parameter T can be interpreted as a temperature and the exponential factor in (2.9) as a Boltzmann suppression factor. While a gradient descent method would uniquely find the local minimum of f connected to the point x 0 , the Metropolis algorithm, accepting from time to time moves where the function increases rather than decreases, is able to overcome barriers between different minima and sample wider regions of U . The "efficiency" in overcoming barriers depends on T . At fixed δx, the higher T the easier it is. Of course, if T is chosen too large, we would explore the entire region U but we might not detect the minima, because with such a diffusive dynamics the probability density P (x) will tend to be uniform in x.
In our MC implementation, the variables x which are changed at each iteration procedure are the scaling dimensions ∆ a (and possibly the external scaling dimension ∆ σ ), while the function f is identified with the action S in (2.6). The procedure works as follows. Starting from some initial scaling dimensions (∆ a ), we analytically minimize S with respect to the OPE coefficients ρ a , given that they enter quadratically in the action. We get where all the terms in the right hand side of (2.10) are evaluated at (∆ a ). We then plug (2.10) into (2.6) to get an action a ) which depends on operator scaling dimensions only (spin and external operator dependence omitted for simplicity). At this stage, the scaling dimensions (∆ a ) are changed as in (2.8), the OPE coefficients recomputed as above and then the move is accepted or not according to (2.9). This procedure is then iterated N steps times. In this way we explore different the local minima of S, and then study the collective trends of these minima as a function of ∆ * . Of course, for the CFT data (∆ σ , ∆ a , ρ a ) associated to a physical theory we should have lim ∆ * ,N Ops →∞ S(∆ a , ρ a ) = −∞ . (2.11) It is then important to establish that, in a finite truncation, the value of S at the minima decreases as ∆ * increases. This is in a nutshell how our protocol works. In the next subsection we will describe in detail the various steps involved.

The search protocol in detail
We choose to sample the N z points in the z-plane from a uniform grid with the constraint with λ(z) defined as in (3.11) of [19]. This condition defines a compact and convex region around z =z = 1/2, the point of best convergence of the bootstrap equations. The parameter λ 0 is, together with ∆ * and T , the most relevant "hyperparameter" of our computational pipeline, and has to be chosen carefully in order for the method to work. The number of MC iterations N steps needed to get sensible results can be quite large (∼ 10 8 ). In order to be able to perform exhaustive searches in parameter space we developed a numerical approach which allows estimating efficiently the action S defined in (2.6). Indeed, conformal blocks are computationally-expensive functions of ∆, and z, which would make an extensive sampling impossible with moderate computational resources without an appropriate numerical optimization. We tabulate the blocks for a fixed z−point sample and each , and then use cubic splines [20] to interpolate their values. In this way we get a speedup of O(10 4 ). The interpolation grid is chosen in such a way that the numerical error induced by this procedure is negligible compared to the truncation error (see figure 17 and the discussion in appendix B.2).
It is important here to note that this speed up comes with the tradeoff of limiting our working precision to that of double type in Fortran (i.e. 64−bit). This might come as a surprise to our more technical bootstrap readers, but an important result of this paper is indeed that approximate solutions to the bootstrap equations can be found without resorting to 200 digits of precision. Nevertheless, limited precision is a bottleneck to go to larger values of ∆ * and to improve the quality of our approximate solutions.
We now describe each step of the protocol more thoroughly. Further details can be found in appendix B.
1. Determining the hyper parameters Once a physical system of interest (understood as a global symmetry structure and the relevant bounds on the operators' scaling dimensions) has been chosen, some exploratory work must be done in order to make our search protocol as efficient as possible. In appendix B we list the choices for each search discussed in section 3. We will here concentrate on motivating our choice for T , since we have found empirically that the other hyperparameters are less determinant to the success of the search.
Since we will explore different truncations, the order of magnitude of δS, see eq.(2.9), can be very different from spectrum to spectrum 4 and thus the temperature T must be adjusted accordingly. Although we will use an adaptive step that guarantees that in the stationary state roughly half of the moves will be rejected, choosing a temperature that is too high will make δ∆ too coarse and the MC will "miss" narrow features of the landscape. On the other hand, if T is too small our method will move too slowly around the landscape or might even "freeze" in the basin of a local minimum.
In the following we will empirically choose T in such a way that the projection of the trajectory to each coordinate axis (which in our case correspond to the operators' scaling dimensions) covers the whole allowed range at least 10 times, but still low enough so that the dynamics does not become diffusive, namely the local minima can still be detected as metastable states.
For further detail, we show in figure 18 an example of archetypal trajectories. This figure and further considerations about the choice of T are presented and analyzed in B.3.

Finite temperature MCs
For each ansatz ∆ (0) a 5 we sample the CFT parameter space at the temperature T chosen according to the criteria discussed in the previous step. This is the most important and computationally expensive step, but one must bear in mind that it gives all the necessary information to understand the landscape of possible CFTs at several different truncations.
The variation of the scaling dimensions δ∆ a is taken randomly with uniform distribution over an hyperparallelepiped centered around the origin. This was observed to be more efficient than moving one operator at a time since the landscape shows highly correlated features and thus attempting moves in which all the operators change simultaneously improves the sampling efficiency. The user can choose the length of each side but the overall scale is adaptive and changes at each MC step. This is done in order to guarantee a 50% acceptance rate in the stationary state. In practice, when a move is rejected (resp. accepted) this overall scale is decreased (resp. increased) by a fixed factor. Each operator scaling dimension is confined within a finite region within two boundary values decided by the user. The MC then explores the landscape for ∼ 10 8 steps. Our choice of the temperature guarantees that after this number of iterations the initial conditions are irrelevant. For each of these spin partitions we store each 1000th configuration visited. We will refer to these as frames.

Separation by sectors
In a given MC run not all the operators in the ansatz will have sizable OPE coefficients at each iteration step. In this work we assume that if the OPE coefficient associated to a primary operator is smaller than a given tolerance (we found 10 −8 to be a sensible choice for the max used in this paper) this operator is decoupled and can be neglected. If unitarity is imposed, OPE coefficients that turn out to be negative from (2.10) are set to zero in that iteration. These two effects lead to a classification of the different points visited by the MC in terms of the effective spectrum contributing to the crossing equations. In order to describe these different sectors we introduce the following shorthand notation to indicate each one of them: own. 5 For simplicity of notation, we do not explicitly report ∆ (m 0 ) (m 2 ) · · · (m max ). For example, a point with 3 scalars, 2 spin-2 operators, 2 spin 4 and one spin-6 would belong to spin sector 3 2 2 1. Crucially, m ≤ N Ops . (2.13) This implies that there is no need to consider different partitions of operators in each spin channel. All partitions with m ≤ n will be automatically searched for by MC. Moreover, the same sector can be visited in two different MC runs with different n . Thus, this step effectively recombines all the information from the different spin partitions studied in the previous step.
4. Identification of putative local minima Now that we have sampled the landscape we must identify the points that are more likely to be close to a local minimum of the action. For a low-dimensional parameter space we could visualize the action as a contour plot and identify the points by inspection. It is clear, though, that for more than 3 scaling dimensions, this approach cannot be followed.
Our proposal is thus the following. For each set of frames obtained in the previous step (namely, those corresponding to spin sectors with more than 120k frames and those describing the trajectory of every search) 6 we sort their elements according to the number of nearest neighbors with higher action, in the philosophy of Density Peak clustering [21]. We consider as putative local minima the 20 configurations with highest rank. We have observed empirically that in each sector the number of actual minima is always much smaller than this number, and that many of the smaller rank frames actually converge to other minima, so this is a conservative estimate. To clarify, the point with the highest rank will be the global minimum of that set of frames because every other point will have higher action.

Local minimization
For each putative local minimum we launch Newton -Rhapson (NR) minimizations. This deterministic algorithm is justified because the points found in the previous step are likely contained in the basin of attraction of a local minimum. More details about the implementation of this algorithm can be found in appendix A.
All the endpoints from the NR minimizations are collected and those that violate the bounds imposed on the spectrum are discarded. 7 We then assign to the corresponding sector the minima where one or more operators decouple (in the sense described above). In the case of two or more minimizations converging to the same point (i.e. the relative difference for any scaling dimension is smaller than 20%) we keep only the point with the lowest action. At the end of this step we have a list of approximate solutions to crossing that are local minima of S. 6 It might seem redundant to analyze sector by sector frames that are already contained in the trajectories that result from the MC search described in point 2. Since exhaustivity is paramount and the definition of local minimum depends on the neighboring points, it makes sense to consider the same frame in two different settings in order not to miss any interesting features. 7 Spectra where a whole spin sector is missing, such as 2 1 0 1 1, have been ignored for simplicity.

Convergence and boundary check
Since we usually impose boundaries on the subleading operators in order to better control the type of CFTs that we might find, it is important to rule out the possibility of some of the local minima found in the previous step actually being induced by the presence of this boundary.
To this end, a final NR minimization is performed from all the minima found in the previous step, relaxing all the boundaries except strict unitarity (if imposed). We consider as true minima the ones that reach convergence according to our NR implementation and are contained within the region of interest: this means a minimum can be at the boundary of unitarity, but not at one of the boundaries introduced in step 1.

7.
Identifying akin minima at different truncations Minima belonging to diferent spin sectors but with similar CFT data in their common operators should be identified. In this case the minimum with more primaries should have a value of S lower than the minimum with less primaries. This is a necessary (but by no means sufficient) condition to associate such minima to approximate solutions of the crossing equations of one and the same CFT. In order to find akin minima and quantify their similarity we organize all the minima in a directed graph where the minima are the nodes and A → B if: (i) the sector to which A belongs is contained in the one of B, namely if B has at most maxMismatch more operators than A and (ii) the operators which are in common between the two sectors differ in their scaling dimensions by at most relTol in relative error. The final result of our protocol are the leaves of this directed graph, namely the nodes with no outgoing links. One can also study the branches (i.e. subsets of this directed graph) in order to understand the asymptotic behaviour of certain values, such as the action S or the ∆ of some leading operators.
allowing ∆ T to vary. We focus mostly on unitary theories, though our method does not rely on it. As a proof of concept, we study the non-unitary landscape around the 2d Yang-Lee model.
As already anticipated in the introduction, when ∆ σ is allowed to vary, the landscape of S becomes very simple: the global minimum is at the lower edge of the unitarity region, which in d = 3, 4 coincides with the Free Theory (FT) of a scalar field. Moreover, for generic spectra we find that it is in fact the only reachable minimum. For this reason in most of our analysis we will study "slices" of the landscape by keeping fixed ∆ σ in each MC run.
The section is structured as follows. We start in section 3.1 by showing that the FT is the global (and in many cases the only) minimum of S. From section 3.2 we look for local minima by keeping the external operator fixed. We first consider the simplest possible scenario, namely very sparse CFTs where just one operator per spin appears below ∆ * . We then discuss in section 3.3 more general CFTs with multiple operators per spin. In section 3.4 we discuss non-local theories with a special focus on the GFT. Finally, we briefly discuss the 2d Yang-Lee model as a proof of concept of the validity of the method for non-unitary theories in section 3.5.

The Free Theory is the global minimum of S
Before attempting to identify the local minima of S it is useful to have a qualitative idea of how S behaves as a function of the CFT data. This analysis can be made by launching MC runs with a given temperature and spin partition of operators, and simply looking for the values of the ∆'s corresponding to the lowest value of S explored during the run.
The most important property that comes out from this analysis, common to every scenario studied in this paper, is that the global minimum of S corresponds to the FT (or the point ∆ σ = ∆ = 0 for d = 2). This has been observed in all our searches, and in particular for all the spin partition of operators which we considered. Of course, this does not imply the non-sensical result that the FT is the only consistent CFT. Rather, it implies that among all the physical CFTs allowed within the given assumptions at fixed truncation, within our numerical accuracy and in absence of further constraints, only the FT shows up as an isolated minimum of S.
In order to test the (non-)existence of other minima, we performed a very large number of MC minimizations at T = 10 −4 . This temperature, being of the order of the numerical error of the action, 9 guarantees that the MC can only diffuse towards configurations of lower action, since the temperature is too low to cross any relevant barrier. Independently of the choice of n and of the operator scaling dimensions (∆ a ), all the low temperature MC runs we performed converge to the FT point within 10 8 steps, although longer times might be needed for higher N Ops . 10 We report as an example in figure 1 the low−T evolution for a 3d CFT spectrum with max = 6, N Ops = 8 and an initial spin partition 3 2 2 1. Here the first spin two operator is fixed 9 See appendix B.2. 10 There is an important exception described in detail in section 3.2. In d = 3 for spectra of the kind 1 1 1 and 1 1 1 1 there are actually two local minima, though the FT still is the global one. at ∆ T = 3 and no gap on the number of relevant scalars has been imposed. 11 The starting point of each colored line (indicated with a black bullet) corresponds to a MC with initial values of (∆ (0) σ , ∆ (0) ) at that point. The black arrows on the colored lines indicate the direction of the MC evolution. As can be seen, the final point of all the MC runs (highlighted by red circles) have ∆ σ ≈ 1/2, ∆ ≈ 1. The contour plot shown in the background is indicative of the value of S for the same spin partition but sampled at a much higher temperature (roughly 100 times larger). A two-dimensional representation of S is then obtained by binning the points visited on the (∆ σ , ∆ )-plane and plotting the minimum value of S in each bin. 12 They are reported as a visual aid. As can be seen from the figure, the initial values of (∆ (0) σ , ∆ (0) ) can be taken in the excluded region, in which case they quickly approach the extremality line. The discrepancy between the rigorous bootstrap extremal line and the trajectories of our MC solutions is due to the severity of our truncation. During the MC evolution the OPE coefficients of the extra scalars and higher spin operators decrease in such a way that at the final point we effectively have the FT spectrum with one operator per spin. We illustrate this phenomenon by plotting in figure 2 11 Note that in d = 3 (the case shown in figure 1), if a scalar gap is imposed (for example, demanding that there can only be one exchanged relevant scalar), depending on the initial CFT data it can happen that the FT cannot be reached and a local minimum is induced. 12 Understood as a cell of fixed width in ∆σ and ∆ . dimensions. This is the first important result of our analysis: for most of the spin partitions studied in this work, the landscape of S has a single minimum, which can be reached by trivial gradient descent from (generically) any initial configuration. Note that most trajectories in figure 1 reach the FT point following two broad paths in the (∆ σ , ∆ )-plane: i) the extremality line from above, ii) a curve, almost constant in , from below. The presence of regions in the (∆ σ , ∆ )-plane where the bootstrap equations are more easily satisfied is going to be a common theme in our results.
As we are going to show next, if one fixes at least another ∆ one can find local minima along these slices. In the following subsections we will apply the protocol described in section 2.3 for fixed values of ∆ σ and we will show that some of the local minima that are present correspond to known theories (for example the Ising model). However, it is important to remark that none of the minima that we will describe are stable if ∆ σ is allowed to vary, not even the minimum corresponding to the Ising models in d = 2 and d = 3, in agreement with the results presented in this subsection.

One operator per spin: d = 3, 4
We now describe the results of a systematic search of local minima at fixed ∆ σ by the protocol described in subsection 2.3. It is instructive to begin by exploring the simplest possible scenario where n = 1, = 0, 2, · · · , max . (3.1) Although this assumption looks very restrictive, some basic features of the shape of the landscape of CFTs is already visible in this set-up.
Step 3 of the protocol is clearly unnecessary when n = 1, since all the configurations in each run will belong to just one spin sector. We focus on d = 3 We searched local minima of S at different fixed values of ∆ σ = 1/2 + n/20, n = 0, 1, 2, 3, with input spin sectors 1 1 1, 1 1 1 1, 1 1 1 1 1 and 1 1 1 1 1 1. The minima that we found are represented in figure 3 (a). In order to better appreciate the value of S at the minima i, the When unitarity is imposed we find two disjoint basins. Allowing for non-unitary configurations gives rise to a single connected basin.
latter are indicated with a circle whose radius r i is given by where S max is the value of S in the minimum with larger action, and α is an offset so that this minimum can still be visible in the figure.
For n = 0, 1 two classes of minima are found, one for each of the two basins of attraction. For n = 2, 3, instead, we find only minima along the lower basin, with ∆ < 1. It should be clear that there is nothing special to the above chosen values of ∆ σ . For each fixed value of ∆ σ we expect to find two local minima if ∆ σ 0.6 and one for higher values.
We then use the protocol described in section 2.3 with relTol= 0.1 and maxMismatch= 2 to identify the different branches, which correspond to sets of minima with different input sectors that can however be considered to represent the same "theory" at different truncations. In this way we find that, among these minima, some are isolated and do not from any branch (orange diamonds in panel (b) of figure 3), while others form branches (green, pink and blue bullets). The only branch in which the action decreases as a function of ∆ * is the one associated to the free theory, the blue branch in panel (b) of figure 3.
In d = 4 a similar high-T search with max = 6 and free ∆ σ gives rise to the contour lines in the background of panel (c) of figure 3. A search of the minima at different fixed values of    As we will show in the following, the main structure of the S landscape with two elongated basins of attractions (valleys) will appear repeatedly in the analysis of more complex situations. Note the striking difference between the d = 3 and d = 4 cases: in the former the basins are two, in the latter the basin is only one. One could argue that this difference is due to the presence of another actual theory in d = 3 close to the FT point, the Ising theory, and none in 4d within the above assumptions. Indeed, our results in sections 3.2 and 3.3 further point to this explanation.
We also find that if the unitarity bounds are relaxed (both on ∆ and on the OPE coefficients squared), the aforementioned disjoint valleys in d = 3 join in a non-unitary point of the parameter space. This is shown in figure 4, where we compare the lowest-lying points in an unitary MC search to those of a non-unitary setting.
The main take-away of this test is that the evolution with respect to ∆ * of the action in branches formed by akin minima seems to be a discriminant factor between spurious minima and physically meaningful ones (the FT in this case); namely, the spurious minima do not create long branches that go down into small values of the action S, while the FT is present for every truncation with very high consistency. We find that in the FT minimum high spin operators are determined quite accurately (∼ 1% error) except for the highest−∆ at each truncation, as can be seen in table 1.
Although it is clear that the minima at ∆ σ = (d − 2)/2 are deeper than those away from it, we performed a MC search at max = 10 without any constraint on ∆ σ besides unitarity. This allows us to understand the precise shape of the landscape around the FT point. After 2 × 10 8 steps of an extensive search (T ∼ 1), we chose the global minimum and performed a T = 10 −4 minimization in order to refine the precise location of this minimum. In table 2 we report these values. For completeness we also show in table 3 the OPE coefficients determined in this case.
We finally note that for higher ∆ * , namely higher N Ops , the FT minimum is much deeper than the other minima. This means that the latter are very poor solutions to the crossing equations. In other words, we did not find other consistent approximate solutions to the crossing equations, besides the FT, where (3.1) is satisfied up to max = 10.

More operators per spin
We now use our protocol to study CFTs with more than one operator per spin with the assumptions discussed at the beginning of the section (presence of an effective Z 2 symmetry and only one relevant Z 2 −even scalar).
We here consider viable solutions to the crossing equations the end-point of branches 13 that arrive up to a given max and include at least 2 different minima with strictly decreasing S. For brevity, we will refer to these end points simply as "end-minima". For example, in the case of only one operator per spin, the "end-minima" are the blue, pink and green circles at the end of the branches in the right panels of figure 3. The rationale behind concentrating on end-minima is that they summarize the information from minima at higher max and thus allow for a cleaner visualization. Since this approach might exclude possible branches that start only at the largest ∆ * here considered, we will also include in our discussion isolated minima at the largest max considered. They will be also called end-minima in a slight abuse of notation.
Our whole search space includes spectra with max = 4, 6, 8, 10. For d = 2 we limit the analysis to max = 8, since at max = 10 the end-minima that we find are so many that they obscure the visualization. Here we only observe that, qualitatively, all the minima that we find at max = 10 differ only by the least relevant operators.
The role of the stress-energy tensor deserves some consideration. For the study of local theories it is necessary to impose the presence of a spin−2 operator with ∆ T = d. In practice, we have also found necessary to impose a gap on the scaling dimension ∆ T of the second spin−2 operator so that it does not come too close to d and create the effect of a non-local theory by "supplanting" the stress energy tensor. In section 3.4 we will study the effect of allowing for non-local theories but hereafter it should be understood that unless specified otherwise we fix ∆ T = d and impose the gap ∆ T ≥ d + 1.
We report the end-minima found in the (∆ σ , ∆ )-plane in figures 5, 9 and 12 for the d = 2, 3, 4 cases, respectively. We discuss these results separately in the sections below, but a few considerations apply to all cases. The contour lines of the lowest-lying configuration sampled by an high−T MC is shown in the background as a visual aid. The orange dashed line corresponds to rigorous bootstrap bounds, while the green continuous line is the line of GFTs. We use the radius of the points to represent the action as before (see (3.2)) and use transparent markers so that coincident minima appear darker. To understand the fact that several minima are superimposed in this representation it is important to keep in mind that there can be up to 10 other dimensions not shown in the plot. Indeed, in many cases there are end-minima which have almost identical values of ∆ σ and ∆ but differ in the scaling dimensions of the other operators.

d = 2
One can wonder if our method can work in 2d given that the landscape of CFTs is notoriously very rich. On the other hand, 2d is the ideal playground to test the method, since entire classes of CFTs are exactly soluble.
We have considered various MCs for 8 different fixed values of ∆ σ : When ∆ σ is too small the numerical analysis becomes quite noisy. For this reason, we did not explore values of ∆ σ < 1/8. 7 of the values in (3.3) have been chosen to be "small" rational numbers, while one was chosen irrational, e.g. numerically equivalent to a "large" rational number. Rational scaling dimensions can be associated to exactly known rational conformal field theories (RCFT). As we will see, in order to limit the number of allowed theories and make the analysis feasible, the assumption of having only one relevant Z 2 -even scalar turns out to be particularly important in 2d. Before presenting our results, it is useful to report some (by no means exhaustive) instances of known CFTs admitting scalars with the above scaling dimensions.
A glimpse at the 2d CFT Landscape The landscape of the known CFTs in 2d is amazingly rich. In contrast to the situation for d > 2, in particular, 2d CFTs can have a continuum spectrum (e.g. Liouville field theory). It is clear that we cannot detect such theories. 14 Exactly marginal deformations are ubiquitous in 2d. These theories will not give rise to minima of S, but to a valley of minima, each corresponding to a CFT as we vary the marginal couplings. Since we look for isolated minima in our search, such theories will also be difficult to detect. The CFTs we can hope to find should have a discrete, relatively sparse, spectrum and should be isolated. Their central charge c is expected to be of order one. Prototypical theories of this kind are RCFTs. Even if we focus on these theories only, we need to further restrict our search space to avoid a proliferation of theories. We will do so by assuming, as mentioned, that only one relevant scalar can appear in the σ ×σ OPE. The importance of this assumption can be appreciated by focusing on specific exactly soluble CFTs. We will briefly consider in what follows the space of diagonal unitary minimal models (c < 1), the S 1 /Z 2 orbifold (c = 1) and SU (2) k Wess-Zumino-Witten (WZW) models (c > 1).
Recall that the scaling dimensions of the (Virasoro) primary fields φ r,s in diagonal unitary minimal models M(m + 1, m) is (conventions as in [24]): 15 ∆ r,s = h r,s +h r,s , where h r,s =h r,s = ((m + 1)r − ms) 2 − 1 4m(m + 1) The range in s in (3.5) can be restricted to 1 ≤ s ≤ r using the field equivalence Minimal models feature a faithful discrete Z 2 symmetry for any m (see e.g. [25] for a pedagogical exposition). 16 Using the equivalence class (3.6), for any r and s we can choose the representative field φ r,s with r + s an even integer. On these representatives, Z 2 -even fields are those with r and s odd, while Z 2 -odd fields have r and s even (the identity operator is φ 1,1 ). It is a straightforward exercise to show that there are one or more unitary minimal models which contain a Z 2 -odd field σ with scaling dimension equal to each of the first 7 values of ∆ σ reported in (3.3). 17 For instance, for ∆ σ = 3/20 within minimal models with m ≤ 100 we find two which contain a Virasoro primary with that dimension: φ 8,8 ⊂ M (15,14) and φ 14,14 ⊂ M (26,25). The fusion rules σ × σ contain respectively 42 and 132 Virasoro primaries. Aside from the identity operator, 11 and 21 of these Z 2 -even fields have ∆ < 2, respectively. The smallest Virasono primaries in the two cases are φ 5,5 and φ 3,3 , with scaling dimensions ∆ = 2/105, 2/325, respectively. Models in such kind of field correlators cannot be detected by our method, because i) their OPE contain too many light operators and ii) operators with a too small scaling dimension makes the numerical analysis unfeasible. Imposing that we have only one relevant Z 2 -even scalar removes most of these field identifications. The only identification left, for any m ≥ 3, is the one where 15 In the 2d CFT literature, Virasoro and SL(2, C) primaries are usually denoted primaries and quasi-primaries (respectively). In numerical bootstrap analysis we usually care about SL(2, C) primaries, which are the analogues of the primary fields in d > 2 CFTs. Most of the considerations that follow can be made at the level of Virasoro primaries with no need to decompose them in SL(2, C) primaries. 16 The presence of such a Z2 symmetry is clear from the equivalence of minimal models with Wess-Zumino-Witten cosets of the form [SU (2)m−2 × SU (2)1]/SU (2)m−1, where Z2 is the combination of the two Z2 center symmetries of SU (2)m−2 and SU (2)1 which is not gauged by the center of SU (2)m−1.
17 Of course, minimal models associated to the irrational value of ∆σ reported in (3.3) can also be found provided we take m large enough and approximate √ 5/10 with an appropriately chosen rational number. σ = φ m−1,m−1 (∼ φ 1,2 ) and = φ 1,3 , already considered in [26], which have and the simple fusion rules σ × σ = 1 + . However, as far as we consider the four-point function of σ only, given the fusion rules (3.8), we can effectively assign an unfaithful Z 2 charge to σ such that Z 2 (σ) = −σ and identify such Z 2 as the global discrete symmetry of the four-point correlator.
Another notable example which shows the importance of restricting to one relevant scalar is given by the theory of a free compact scalar on a S 1 /Z 2 orbifold. This CFT has two Virasoro primaries with ∆ σ = 1/8, which are the twisted ground states of the orbifold. They are odd under the quantum Z 2 global symmetry that emerges after the gauging and hence we can consider any of them as our field σ. As well-known, the S 1 /Z 2 orbifold is a c = 1 theory admitting an exactly marginal deformation, given by the radius R of S 1 . At R 2 = 1 the theory is equivalent to a decoupled pair of Ising theories. 18 Deforming the radius gives rise to a family of theories related to the Ashkin-Teller model, namely the theory obtained by coupling two Ising models [27]. The OPE of two twist fields give rise to an infinite number of Kaluza-Klein and winding vertex operators (Virasoro primaries) with scaling dimensions [28] ∆ m,n = m 2 R 2 + n 2 R 2 4 , m ∈ Z , n ∈ 2Z . (3.9) At R 2 = 1, the two Ising theories are decoupled and this theory is indistinguishable from a single Ising copy. When R = 1 the two couples to each other and we get a continuum of c = 1 CFTs parametrized by R. For any R such theories contain more than 1 relevant Z 2 -even scalar, as we see from (3.9) by setting either m or n to zero. In this case, restricting to only one relevant Z 2 -even scalar is crucial to avoid a continuum of theories, which would be a challenge for our numerical methods. The external operator σ could also be part of a CFT with a continuous global symmetry G. CFTs with continuous global symmetries should of course be analyzed in a covariant way exploiting the symmetry, see e.g. [30]. However, if G ⊃ Z 2 the field σ could be considered some (real) component of a multiplet in a representation of G for which Z 2 (σ) = −σ. For example, σ could be the real component of a spin j multiplet, with half-integer j, of a SU (2) k WZW model.  [29] assuming only one Z 2 -even relevant scalar. All the end minima correspond to local theories ∆ T = 2, ∆ T ≥ 3.
More specifically, we could have where φ j,m (z) andφ j,m (z) are the holomorphic and anti-holomorphic highest-weight state components of the spin j field. The symmetry Z 2 could be identified with the center of either the holomorphic or of the anti-holomorphic SU (2) gauge group. With the identification (3.10), we have This quick detour on some well-known 2d CFTs shows the importance of reducing the space of viable theories by demanding only one relevant scalar. Needless to say, the space of CFTs is vastly larger than the three classes we discussed, so it is likely that other CFTs with σ operators with scaling dimensions as in (3.3) and only one relevant Z 2 -even scalar exist.
Results in d = 2 We report in figures 5 and 6 the location of the end-minima found with max ≤ 8 in the (∆ σ , ∆ ) and (∆ σ , c)-planes, respectively. The end-minima belong to different sectors, all of them contained in 4 4 3 2 1.
Most of the end-minima are aligned along two approximately straight trajectories in the (∆ σ , ∆ )-plane. The blue end-minima are along the upper extremal bootstrap bounds, while the green ones are aligned along the GFT line, despite the constraint ∆ T = 2.
Let us first discuss the blue end-minima. Two less significant blue end-minima appear above the extremal line in the disallowed region, indication of being fake end-minima not associated to CFTs. The rest are in good agreement with the bootstrap bounds.
It is reassuring that, among such end-minima, we reproduce the first 4 minimal models (red crosses in figure 5), with the field identification as in (3.7) discussed before. The remaining four points along the line correspond instead to the so called generalized minimal models [16], i.e. Comparison between the OPE squared coefficient λ 2 σσ numerically determined from the blue end-minima in figure 5 and its analytic expression for generalized minimal models (red line) [16], as a function of ∆ . non-unitary minimal models obtained by analytically continuing the CFT data appearing in the 4-point correlation function from integer values of m (minimal models) to real values of m (generalized minimal models). Such theories are visible in the conformal bootstrap, despite the enforcement of unitarity, because the OPE coefficients squared of the exchanged quasi-primaries have been conjectured in [16] (see also [31]) and then proved in [29], to be positive for any m.
Given (3.7), the blue end-minima should be aligned along the line (3.14) In terms of the central charge, we have We compare in figure 6 the central charge (3.15) with the one extracted from the OPE coefficient λ σσT by means of the relation valid for general d dimensions. We see that the agreement is very good, despite the severity of our truncation, and confirm also the identification of such end-minima with generalized minimal models.
As a further check we compare in figure 7 the OPE coefficient squared λ 2 σσ with the one analytically found for such models in appendix B of [16]. The good agreement found is yet another indication of the accuracy of our determination and of the nature of the blue endminima. The only exception is the minimum at ∆ σ = 7/40, which falls slightly above the SDPB bounds in figure 5 and its values of c and λ 2 σσ in figures 6 and 7 are slightly misaligned with those of the generalized Minimal Models. We expect this to be due to some numerical artifact and that a more thorough study will find a lower action minimum that better agrees with the theoretical values.
Let us now discuss the green end-minima of figure 5. In contrast to the blue ones, we do not have a convincing explanation for the nature of these end-minima as local CFTs, assuming they are not numerical artifacts. Most of them lie on the GFT line ∆ = 2∆ σ . A notable class of local CFTs that have in their spectrum operators with scaling dimensions related in this way are N = 2 CFTs. It would be interesting to understand if the green points along the GFT line (or some of them) can be identified as N = 2 Supersymmetric Minimal Models, provided an appropriate identification of σ is made. This analysis would probably require also to understand whether generalized N = 2 minimal models, in the spirit of [16], exist and can be constructed.
Despite we enforced the presence of an energy momentum tensor, the green end-minima along the GFT line could simply be an approximation of GFTs themselves. This observation is supported by the fact that if we remove the gap assumption ∆ T ≥ 3 most of the green minima are unstable. In particular, λ 2 σσT ∼ 0 and ∆ T ≈ 2∆ σ + 2 effectively replaces T as the first spin-2 operator. Similarly, if we do not assume the existence of an energy momentum tensor, most minima turn out to be GFTs. As we will see in the next subsections, the analogue of the green points in figure 5 will not occur in d > 2 when we enforce the presence of an energy-momentum tensor.
It is useful to also report the spectrum of SL(2, C)-primary operators which appear at the end-minimum at a given truncation. We report as an example in figure 8 the full spectrum of some end-minima found at ∆ σ = 1/8, 1/5 (two in each case). The operators associated to the same end-minimum are connected with lines. In the Ising case, where the identification of the end-minima with the Ising model is quite convincing, we see how the agreement of the whole low-lying spectrum is quite accurate, with small deviations occurring for higher-spin operators in the first Regge trajectory (the spin-8 orange diamond) or for operators in the second Regge trajectory (the second spin-0 green diamond, the second spin-4 orange diamond).
In the GFT-like case, as already mentioned, we do not have a convincing explanation for the nature of such theories. Despite having enforced the presence of a spin 2 operator at ∆ = 2, we can see from the right panel of figure 8 that the spectrum resembles roughly the one of GFTs. The presence of more operators with respect to the Ising case makes the accuracy of the numerics less accurate. In particular, it is not possible to disentangle if the significant deviations from the GFT spectrum at higher Regge trajectories are, in addition to numerical noise, due trivially to the fact that we have imposed the presence of a non-existent operator (the stresstensor), or because the theory under question is in fact distinct from a GFT. See appendix C for some further details, in particular for examples of selected branches, the endpoints of which correspond to the blue and green end-minima discussed above.
It is well possible that other theories with a relatively sparse spectrum satisfying our assumptions (or effectively doing it, like the generalized minimal models) might have not been found by our search protocol, given also the limitations imposed by machine precision computations and the small values of ∆ * and N Ops we sampled. Possible theories of this kind are for instance the SU (2) k WZW models (3.13) with the field identification (3.10).
(3.17) We plot in figure 9 the end-minima found in the (∆ σ , ∆ )-plane, with ∆ T = 3, ∆ T ≥ 4, ∆ ≥ 3. All the end-minima lie at the extremality bound, like in the d = 2 case. In contrast to the latter, no end-minima along the GFT line appear, while we find a couple of faint end-minima with ∆ < 1, below the GFT line. 20 Reassuringly, we get end-minima in correspondence of the 3d Ising model, which sits at the kink of the dashed orange lines. End-minima are found for all the ∆ σ in (3.17) with the exception of ∆ σ = 0.505. The lack of an end-minimum at such value is due to the imposed gap ∆ ≥ 3. Indeed, by repeating the MC analysis of section 3.1 with ∆ σ unconstrained assuming ∆ ≥ 3, the "flow" towards the FT passing through the valley connecting the Ising point with the FT in figure 1 gets interrupted around ∆ σ = 0.510. This is an indication that possible extremal CFTs with 0.5 ≤ ∆ σ ≤ 0.510 do have more than one Z 2 -even scalar. 21 We plot in figure 10 the central charge of the end-minima, computed using (3.16), as a Ising reference values taken from [11]. Right: Full spectrum of the minima found at ∆ σ = 0.53, ∆ T = 3 that have ∆ below the GFT line. The colors and the lines connect operators within the same end-minimum.
function of ∆ σ . The dashed orange lines delineate the rigorous lower bound of [31]. Most of the end-minima in figures 9 and 10 lie in the forbidden region, slightly in the former case and more evidently in the latter. This difference should not come as a surprise, since it is known that the bounds on the central charge are more sensitive than those on ∆ to higher dimensional operators. Compare e.g. figures 3 and 4 of [19] to appreciate the different convergent properties of the two bounds as the degree of truncation is varied. The spectra associated to the branches with ∆ σ = 0.5181489, 0.53 is reported in figure 11. In the left panel we also report as reference values the Ising spectrum found in [11] using the Extremal Functional Method (purple lines). The operators associated to the same end-minimum are connected with lines. In the Ising case we see how the accuracy is significantly worse with respect to the 2d case. We have a reasonably good agreement of the low-lying spectrum in the first Regge trajectory, but the remaining operators have a significant indeterminacy. Similarly, in the right panel of figure 11 the operators in the first Regge trajectory agree among the two orange and green branches, but significantly differ at higher levels.
In d = 4 we run our protocol for the following values of ∆ σ : Due to the limited numerical accuracy of our protocol in d = 4 we just report in figure 12 the end-minima found in the (∆ σ , ∆ )-plane with ∆ T = 4, ∆ T ≥ 5, ∆ ≥ 4, with max ≤ 10. Once again we have end-minima along the extremality bound but it is interesting to see how our minima satisfy the older bounds (with less derivatives) of [1] but are slightly excluded by the bounds in [32]. Interestingly enough, for each value of ∆ σ sampled we also find end-minima below the GFT line.
∆T = 2∆σ + 2 (the GFT value), as expected. 21 Note that this region is devoid of not only end-minima, but of minima at all.

Non-local theories
In this subsection we analyze the structure of the minima found when relaxing the condition ∆ T = d. The most notable and exactly calculable non-local theory is the GFT. At a given ∆ * the number of operators in the GFT spectrum is larger than the maximal truncation we considered in this paper. Although in principle this rules out the possibility of finding exactly the GFT as end-minima in our approach, the limitation is only theoretical, considering that our numerical accuracy does not anyhow allow to precisely determine operators beyond the first Regge trajectories, as we have seen. In fact, we will find end-minima which seem quite compatible with the leading operators of the GFT at different values of ∆ σ . The landscape of end-minima is now much richer and cover most of the allowed region. Interestingly enough, along the ∆ T direction we find that there is a clear tendency of the minima to cluster around ∆ T ≈ 2∆ σ + 2, see figure 13 (a) for d = 2 (a similar phenomenon occurs in d = 3 and d = 4), where we report the distribution of the end-minima in the (∆ σ , ∆ T ) plane. As can be seen, we have many more end-minima with respect to those found when ∆ T = d and ∆ T > d+1 were imposed. Interestingly enough, most end-minima have d ≤ ∆ T ≤ d+2∆ σ . Note that only for a limited range in ∆ σ , end-minima with ∆ T = d are found. This is particularly evident in d = 2, where an end-minimum at ∆ T = 2 is found only at ∆ σ = 1/8. This signals the fact that the extremal end-minima associated to the (generalized) minimal models with m > 3 found in subsection 3.3.1 become unstable when we relax ∆ T and move towards the GFT region. In panel (a) of figure 13 two specific end-minima points have been singled out at ∆ σ = 1/8 (green circle) and at ∆ σ = 3/10 (red circle). They have been selected as those resembling most closely the GFT conformal data using (3.19) as criterion. The truncated spectrum at those points is reported in panel (b) of figure 13. Note the good agreement beyond the first Regge trajectory, in particular in the scalar sector where the first three scalars are well reproduced in both theories.
We have also explored the effect of fixing the scaling dimension of the first spin 2 operator at the exact GFT value, still assuming ∆ T ≥ d + 1. This is useful to understand whether minima at fixed ∆ T are "harder" to reach than similar points when ∆ T is allowed to vary. We show these results in figure 14 for the d = 3 case, where the blue circles are obtained at ∆ T = 3 and are the same as in figure 9, while the red ones are those arising at ∆ T = 2∆ σ + 2. The red end-minima align along the GFT line, as expected, but they spread a wider area, which however does not include the extremality region. Similar considerations apply also in d = 2 and d = 4.
We saw in figure 13 (a) that some of the minima found by our protocol at ∆ σ = 1/8 do resemble local theories, although they are clearly disfavoured in front of the GFT-like ones. We would like to determine more precisely how well we reproduce local theories when no constraint on ∆ T is imposed. This can be attained by taking the 2d Ising theory and the GFT as reference points and defining a distance measure between our minima and these exactly known CFTs. We will here take the scaling dimensions of the first scalar and first spin-2 operators and their associated OPEs coefficients, and compute the average relative error with respect to those of the known theory: (3.19) Figure 14: d = 3. Same as figure 9 but we also show with red circles the minima found by fixing ∆ T to the GFT value. Note that we show here in blue all the end-minima that in figure 9 were reported in blue and green.
where v = (∆ , ∆ T , λ 2 σσ , λ 2 σσT ). In figure 15 we show how the blue circles reach configurations that are not too far away from the 2d Ising theory (roughly 10% mean relative error) but also cover configurations which reproduce the leading conformal data of the GFT even better than the MC which had ∆ T = 2 + 2∆ σ to begin with. Finally, figure 15 supports the interpretation of the green minima in figure 5 to be "images" of the GFT, since we can appreciate that they correspond to solutions which reproduce the Conformal Data of this theory with an error of less than 10%.

A taste of non-unitarity: The 2d Yang-Lee Model
The MC search protocol of section 2.3 does not rely on unitarity. Although we mostly focus on unitary models in this paper, we would like to briefly show here a proof of concept of its use in the non-unitary realm, by checking that the d = 2 Yang-Lee (YL) model can easily be found using our method. Non-unitary theories have received less attention from the bootstrap community due to their non tractability with linear or positive semidefinite programming methods. The best approach so far to solve non-unitary theories with the conformal bootstrap is given by [33,34], which also severely truncates the crossing equations (as we do). As well-known, the d = 2 YL model can be described by the non-unitary minimal model M(5, 2), which contain only two Figure 15: Distance of different end-minima to the d = 2 Ising and the GFT with ∆ σ = 1/8. The quantity Q is computed according to (3.19) for end-minima with ∆ T = 2, 2 + 2∆ σ and those in which ∆ T is determined dynamically.
Virasoro primaries related by the following fusion rule: (3.20) This means in practice that the external dimension ∆ φ will also be the scaling dimension of the first exchanged scalar and that we must not impose positivity of the OPE coefficients. We launched a Metropolis Monte Carlo search with T = 0.5 and a spectrum 3 1 2 1 (∆ * = 9) to sample a parameter space close to the one containing the exact CFT data of the 2d YL model. We find that the global minimum is compatible with the exact values, as can be seen in figure  16 for ∆ φ and ∆ 4 (the first spin 4 operator). At the global minimum the other operators show also similar agreements with the exact values.
We note in passing that in d = 3, if we relax the positivity constraint on the OPE coefficients, a search with a 1 1 1 1 spectrum shows another basin of attraction besides those discussed in section 3.2. The basin is in a region well below the unitarity bounds for ∆ σ and lies roughly along the line ∆ σ = ∆ . Besides being an interesting feature of the landscape in itself, we consider this to be an indication that the d = 3 YL model should also appear as a minimum if investigated with our method.

Outlook
Motivated by the pressing need of a tool to explore the inner part of the allowed space of CFTs, we have introduced in this paper a method based on a stochastic minimization via a Metropolis algorithm to find approximate numerical solutions with a sufficiently sparse spectrum of operators to truncated bootstrap equations for CFTs in arbitrary dimensions. Given an initial value of scaling dimensions for a reduced set of primary operators up to some scaling dimension ∆ * , the action (2.6) is minimized and an approximate set of CFT data found.
While for all physical CFTs we should eventually have F = e S → 0 in the limit when all operators are taken into account, at fixed operator truncation the rate in which this limit is reached does matter. For instance, we have seen that for arbitrary severe truncations of the operator spectrum, the free theory minimum dominates over all and no other theory can be found, see e.g. figure 1. In order to find other theories, extra conditions, such as fixing the external operator dimension, have to be imposed. It is clear that, at fixed truncation, not all CFTs are equally accessible, even assuming the same sparsity of spectrum. We expect that richer set of solutions can be found by generalizing our method to multiple correlators and by imposing further constraints that would halt the MC evolution towards specific attractors (e.g. the extremality line or the GFT line in d = 2). Though we mostly focused on unitary theories, the method does not rely on it and can efficiently be used in the non-unitary realm. This is particularly important considering that for non-unitary theories we still do not have rigorous numerical protocols. It would be interesting to see how our numerical method performs with respect to Gliozzi's method [33,34] and its subsequent versions [35,36]. The aim of this work is very similar to that of [14,15], where approximate solutions to crossing have been achieved using a similar logic but a different technique (reinforcement learning instead of MC techniques). It would be interesting to systematically compare the two approaches.
In this work the bootstrap equations have been evaluated in a given sample of points in cross-ratio space. This choice has been primarily dictated by the flexibility, when developing the algorithm, of changing the number of points and the possibility (used in earlier versions of the code) of sampling stochastically these points. The actual performance of the code sensitively depends on how points in cross-ratio space are chosen. It would be interesting to further improve on this aspect. In particular, an improvement could result by choosing derivatives evaluated at the crossing symmetric point z =z = 1/2, as currently done in functional bootstrap studies.
The FORTRAN-code used in this work is limited to machine-like precision. While this allowed us to perform extensive MC searches at little cost, limited precision does not allow us to systematically increase the value of ∆ * and N Ops . For example we do not have enough accuracy to understand the nature of the green end-minima we found in 2d (figure 5), whether they are unidentified local CFTs of some kind or some sort of numerical artifacts resembling GFT theories. For d > 2 the accuracy is even less and the nature of the found end-minima is yet to be understood. Despite the limited accuracy of the algorithm, we think the results of this paper have shown that the method works. Upgrading the current code to arbitrary precision will be the key for a fully working program and is the most important thing to do in the near future. We believe this is feasible, given the ample room for improvement of the code performance in terms of speed.
In summary, we hope that the non-rigorous approach initiated in this work, properly extended and improved, might be a useful tool to navigate in the space of allowed CFT and guide subsequent analysis performed using more rigorous numerical algorithms.
In short, the method uses the hessian matrix to approximate the function to be minimized up to quadratic order and solves iteratively for the minimum of this hyper paraboloid. In this work a modified version of this method is used (known as the Levenberg-Marquardt algorithm [37,38]), in which a damping factor β interpolates between gradient descent and the bare Newton-Rhapson method, avoiding the overshooting of the minimum when very close to it. This is necessary because due to numerical noise some eigenvalues of the hessian might not be positive, which violates the conditions for the convergence Newton-Rhapson method.
More concretely, let f (x) : R n → R be a twice-differentiable function. Given an initial guess x 0 we find the next point by using the following formula: At each step, if f (x 1 ) < f (x 0 ) then β is decreased by a factor of 3. Otherwise it is increased by a factor of 2. This is known in the literature as "Delayed gratification" [39] and gives stable results in the minimization problems explored in this work. The factor β also gives the termination condition for the algorithm: when β > 10 6 it is clear that a stationary point has been reached.

B Implementation details
The protocol detailed in section 2.3 was implemented in a series of automatized bash scripts that used Julia for the data processing and FORTRAN90 (compiled with ifort v2021.2) for the numeric evaluation of the action and the different minimization algorithms. LAPACK was used for the linear algebra manipulations. They were run on Intel Xeon TM processors in an institutional cluster. The CPU time to run 10 8 steps was of roughly 12 hours for representative spin partitions. The parameters used in the searches are shown in table 4.
In order to provide the technical details we will discuss in the following subsections the different elements of our pipeline and their implementation.

B.1 Numerical evaluation of S
We start by discussing the approximation and numerical evaluation of (2.6). We describe each of the approximations used in order to render computationally feasible sampling O(10 8 ) points in ∼ 10 hours on a single CPU.
i) Discretization in z The grid is parametrized as z = x 0 + an + y 0 + ibm with m, n ∈ N and x 0 = 1/2 + 1/1000, y 0 = 1/1000. In this work we used a = 1/100 and b = 1/50 but in the computation of S only a randomly chosen sample of nZ points was considered. This was done by ordering them with respect to their value of λ and then sampling this list. ii) Truncation of the spectrum This was implemented by passing to the FORTRAN routine the list of the spins of each operator in the truncation under study. It is in this step that we input important parameters such as the spin partition, ∆ * and max .
iii) Evaluation of the conformal blocks The conformal blocks were tabulated using 160 digits of internal precision in Mathematica and stored as plain text files for each . In d = 2, 4 we used the closed expressions in terms of hypergeometric functions [17]. For d = 3 the recursion relations implemented as in [40] were used. For each point in z we computed the CBs at 5000 points evenly spaced in [∆ unit ( ) − 99/100, ∆ * ] with ∆ * = d + 15. We note in passing that this gives a lower bound on the error in our determination of the scaling dimensions (∆) of the exchanged operators. In the FORTRAN code, the interpolation is done with cubic splines.
iv) Determination of the OPE coefficients Once the matrix of the conformal blocks for each z and operator (∆, ) has been computed, the determination of the OPE coefficients is then a weighted least squares problem which we solve using LAPACK (see 2.7). If positiveOPEs is true, this is not the end of the story because the OPE coefficients that solve the quadratic problem might as well be negative. To enforce unitarity, we iteratively decouple the offending operators until every OPE coefficient is positive. It is important to underline that this is done at each point evaluated in our algorithm, and thus operators do not decouple for the whole test but only for those particular configurations. v) Computation of S As described in section 2, we obtain S by taking the logarithm of the sum of the squared residuals. If any operator has a ∆ outside its allowed boundaries we impose a quadratic penalty that is then added to the action. This guarantees the analyticity of our potential and thus avoids noxious boundary effects.
More concretely, suppose an operator bound to be in the interval [∆ 0 , ∆ 1 ] has ∆ < ∆ 0 . In that case, the effective action will be where the constant wall is set to 10 4 .

B.2 Numerical error
After all the approximations discussed above, reasonable doubt could remain about the accuracy of S computed with our method. Thus, we consider it necessary to quantify the magnitude of the numerical errors in our implementation. We accomplish this by computing S at representative points with our framework and by comparing them to the values of S obtained in Mathematica with arbritrary precision.
From this analysis we conclude that the action computed with Mathematica is systematically lower than the approximated one using our FORTRAN implementation. However, the discrepancy is such that at the values of ∆ * studied in this paper we can still trust the results. Naturally, the errors become larger for bigger N Ops and ∆ * , but even in the cases with the worst agreement the relative errors are smaller than 10 −3 . Moreover, the largest discrepancy sets a natural lower bound for T , since we must allow for fluctuations of this height in order not to wrongly identify artifacts as minima. For example, in d = 3 the minimum temperature for max = 6 is T ∼ 10 −4 , whereas for max = 8 is T ∼ 10 −3 . We note that this resolution is likely enough for any physical minimum, given that in the ∆ * → ∞ limit any physical theory should have S → ∞.
As an example, we show in figure 17 the relative error computed for 10 4 gaussian perturbations of a local minimum in sector 2 2 1 1 1 for d = 3, λ 0 = 0.42 and nZ = 200. The fact that the error is mostly negative can be easily understood from the delicate numerics of the bootstrap equations. It is known that in order to solve (2.1) to a high degree of precision, very fine-tuned cancellations between the terms must occur. When using our double-precision interpolation of the blocks, there will be points where these errors will spoil the cancellations and thus the value of S can only be bigger than the "exact" one.
Another point that deserves discussion is the determination of the OPE coefficients (ρ in eq. (2.6)). Solving for ρ amounts to finding the minimum of a quadratic form. We find that the curvature of this form is very anisotropic, with hierarchies of as much as 10 orders of magnitude. This is in fact one of the main bottlenecks to adding too many operators: the curvature of the subleading operators at some points becomes too small to be resolved with double precision and the OPE coefficient of these operators cannot be reliably determined.

B.3 Ideal Temperature
As mentioned earlier in section 2, the key step of our protocol is the wide search performed using the Metropolis Monte Carlo algorithm. In order for this step to be efficient the temperature must be chosen wisely. To this end we performed test searches for different values of d, N Ops and ∆ * , where several different temperatures were used. Then, we determined by inspection the ideal In figure 18 we show a concrete example of how to identify the ideal temperature. We show "cold", "hot" and just right temperatures in blue, red and green (respectively). We can see clearly that the blue line "freezes" into the first local minimum found by the MC and then stays there for the whole test, whereas the red trajectory travels back and forth randomly. We consider the green trajectory to be representative of an efficient search because it oscillates around local minima for some time before making a transition into another one. While this could mislead the reader into thinking that the time spent oscillating around a local minimum is in some sense "wasted" one must bear in mind that during those steps, the MC effectively samples the rest of the scaling dimensions, thus refining the solution to crossing in that neighborhood.

C Minima and End-Minima
In this appendix we illustrate the relation between minima and end-minima. As discussed in the main text, end-minima are defined as those special minima which are the end-points of branches. The latter are a set of akin minima, supposedly associated to the same CFT.
We report in the left panel of figure 19 the individual minima found in d = 2 for some selected branches at ∆ σ = 1/8 (Ising) and at ∆ σ = 1/5 (GFT-like). Note the difference with respect to figure 5. There we reported only the end-minima associated to the end-points of all the branches. Here, we pick up a given branch and show all the minima contained in it. Since most of these minima overlap in figure 19 and similarly several branches overlap in figure 5 (detectable from the darker color of the circles), the difference between the two figures can be overlooked at first glance. In the right panel of figure 19 we report the value of the action S for each of these minima making the branches manifest. We also note in passing that the spectrum at the end-minimum for each value of ∆ σ is shown in figure 8 in orange (the green spectrum corresponds to a similar spectrum not belonging to the same branch).
In order to have an estimate of how the action is supposed to decrease with ∆ * , we also report the branch obtained by starting low T -MCs from several truncations of the exact Ising scaling dimensions (purple squares in the right panel of figure 19). This implies that there are several minima in the vicinity of the exact Ising values and that our protocol is able to find at least a subset of them. It is also reassuring to see that the decrease in S of the branches found and of the Ising benchmark branch are fully compatible. The good behaviour of the GFT-like branches is the main reason why we believe such theories are not numerical artifacts.
As further example we show in the left panel of figure 20 the individual minima for three selected end-minima in d = 3. The blue one corresponds to the Ising model, the green one to the end-minimum below the GFT green line and the olive one to one of the end-minima at ∆ σ = 0.55. In the right panel we show the value of the action S for each of these minima (the color code indicates the different branches).
We point out in passing that a consistent branch should not necessarily be monotonically decreasing as ∆ * increases, because it can happen that minima with less operators have a slightly higher ∆ * than a related minima with a denser spectrum. For instance, the green branch in figure  20 (b) would seem to "turn around" at the beginning, but the important fact is that we make sure that minima with more operators have strictly smaller S.