Bootstrapping N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 4 super-Yang-Mills on the conformal manifold

We combine supersymmetric localization results with numerical bootstrap techniques to compute upper bounds on the low-lying CFT data of N\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \mathcal{N} $$\end{document} = 4 super-Yang-Mills theory as a function of the complexified gauge coupling τ. In particular, from the stress tensor multiplet four-point function, we extract the scaling dimension of the lowest-lying unprotected scalar operator and its OPE coefficient. While our method can be applied in principle to any gauge group G, we focus on G = SU(2) and SU(3) for simplicity. At weak coupling, the upper bounds we find are very close to the corresponding four-loop results. We also give preliminary evidence that these upper bounds become small islands under reasonable assumptions.


Introduction and summary
Four-dimensional N = 4 super-Yang-Mills (SYM) theory is one of the most well-studied models of gauge dynamics in theoretical high energy physics. Part of the reason for the continued interest is that this theory has a large amount of symmetry. Indeed, not only does N = 4 SYM have the maximal amount of supersymmetry possible in four dimensions for a non-gravitational theory, but it also possesses conformal symmetry as its beta function vanishes identically. In fact, for any simple gauge group G, this theory possesses a conformal manifold parameterized by the complexified gauge coupling τ ≡ θ 2π + 4πi g 2

YM
, with g YM being the Yang-Mills coupling and θ the theta-angle. The case G = SU(N ) also features prominently in gauge/gravity duality [1][2][3], where in the combined large N and large 't Hooft coupling λ = g 2 YM N limit, N = 4 SYM has a dual description in terms of weakly coupled string theory on AdS 5 × S 5 . N = 4 SYM has been studied using a variety of methods in various limits. Indeed, for any gauge group G, perturbatively at small g YM one can calculate various observables using Feynman diagrams -see, for instance, [4][5][6][7][8] for cutting-edge results at four loops.

JHEP01(2023)038
As mentioned above, the G = SU(N ) theory can be studied in the combined large N and large λ limit using the holographic duality. Certain observables at large N and finite λ have also been computed using the technique of integrability (see, for instance, [9,10] for reviews and references). Lastly, exact results at finite N and finite λ so far have been restricted to properties of supersymmetry-protected operators that can be deduced either from anomaly considerations, the 2d chiral algebra subsector [11], or using the technique of supersymmetric localization (see for instance [12][13][14][15], as well as [16] for a collection and reviews and references). In addition to these exact results, at finite N and finite λ, one can also obtain bounds on the low-lying local CFT data using the technique of conformal bootstrap [17][18][19][20][21]. In these studies, however, it was not known how to make the bootstrap equations sensitive to the coupling τ , so these bounds necessarily apply to the entire conformal manifold.
In this paper, we will combine supersymmetric localization with the numerical conformal bootstrap to put bounds on local CFT data for any τ . 1 We study explicitly the G = SU (2) and SU(3) theories, but our analysis generalizes straightforwardly to any G. In [23,24], exact relations were derived between certain integrated stress tensor multiplet four-point functions and derivatives ∂ τ ∂τ ∂ 2 m F | m=0 and ∂ 4 m F | m=0 of the sphere free energy F (m, τ,τ ) of the N = 2 * theory, 2 which was computed using localization in terms of a rank(G)-dimensional integral that depends nontrivially on the full complex τ [12] (see also [25][26][27]). These two integrated constraints were used in [23,24,[28][29][30] to constrain the large N expansion of the correlator. In particular, by combining these constraints with analyticity in Mellin space and crossing symmetry, refs. [23,24,[28][29][30] determined the first few 1/N 2 corrections at finite τ , which correspond to protected higher-derivative corrections to the supergravity action on AdS 5 , and matched the type IIB S-matrix in the flat space limit.
Here, we explain how to combine these two integrated constraints with the infinitely many constraints from crossing symmetry and unitarity. When expanded in conformal blocks, the integrated constraints and the crossing symmetry constraints receive contributions from blocks of large ∆ that differ exponentially in ∆. This leads to difficulties in combining these constraints. In the typical numerical formulation of the bootstrap problem, one would approximate the crossing symmetry constraints by a polynomial in ∆ times a positive factor. If one were to do the same for the integrated constraints, the difference in the large ∆ behavior would require approximating an exponential by a polynomial. This approach was in fact used in an analogous 2D study [31], but in our case we found that such an approximation made the bootstrap insensitive to the integrated constraints. Instead, we employed the approach to the conformal bootstrap based on linear programming [32], in which the values of ∆ are discretized up to a cutoff. In this way, both the crossing and integrated constraints become a set of linear constraints for which feasibility can be evaluated using standard linear programming methods.

JHEP01(2023)038
In practice, because the localization input used in the integrated constraints for the SU(N ) gauge theory is an (N − 1)-dimensional integral, we only applied this approach to N = 4 SYM with gauge group SU(2) and SU (3). We computed upper bounds on the scaling dimension and OPE coefficient of the lowest dimension scalar unprotected operator (at weak coupling, this is the Konishi operator). At weak coupling, we found that our upper bounds are nearly saturated by the four-loop results. When scanning over all couplings, we found that the maximal value of the scaling dimension bound occurs at the self-dual point τ = e iπ/3 with enhanced Z 3 symmetry, and that this maximal value was strictly lower than the bound obtained without using the integrated constraints [17][18][19][20][21]. Furthermore, after imposing a gap above the lowest dimension operator, we found lower bounds that converge to the upper bounds as the gap is increased.
The rest of this paper is organized as follows. In section 2 we review the conformal block decomposition of the stress tensor four-point correlator, the definitions of the integrated constraints, and weak coupling results for local CFT data. In section 3, we show how to combine the integrated constraints with the usual crossing symmetry constraints to numerically bootstrap N = 4 SYM as a function of τ for any gauge group. We carry out this procedure for the SU(2) and SU(3) theories. Finally, in section 4 we end with a discussion of our results and of future directions. Various technical details are discussed in the appendices.

Stress tensor four-point function
The main object of study in this paper is the stress tensor multiplet four-point function. We begin by discussing general constraints coming from N = 4 superconformal invariance, including the expansion in superblocks. We then review the two exact constraints on the correlator that relate certain integrals of the correlator to derivatives of the mass deformed sphere free energy, which can be efficiently computed using supersymmetric localization. Finally, we review analytic weak coupling predictions for CFT data that appears in the correlator, which we will compare to our numerical finite N and τ results in the next section.

Setup
Let us denote the bottom component of the stress tensor multiplet by S. This operator is a dimension 2 scalar in the 20 irrep of the su(4) R ∼ = so(6) R R-symmetry algebra, and can thus be represented as a rank-two traceless symmetric tensor S IJ ( x), with indices I, J = 1, . . . , 6. For simplicity we will contract these indices with null polarization vectors Y I , where Y · Y = 0. We are interested in studying the four-point function SSSS , which is fixed by conformal and su(4) symmetry to take the form

JHEP01(2023)038
The constraints of superconformal symmetry are given by the Ward identity in [33], whose solution can be formally written in two different ways. The first expression is where S free (U, V ; σ, τ ) is a free theory correlator in a theory of 4c = N 2 − 1 scalar fields. In particular, in such a free theory, where a = 1, . . . , 4c, and using Wick contractions with the propagator X a Thus, all information about the interacting N = 4 SYM theory is encoded in the function T (U, V ). Another way of solving the superconformal Ward identity is [18,33] where Θ is the same as in (2.3) and we define (2.6) If we then equate (2.5) to (2.3), we find that the free theory correlator fixes the holomorphic functions f j to be while the remaining free part contributes to G(U, V ), which is then related to T (U, V ) as We can expand G(U, V ), 3 in terms of long and short multiplets as short (U, V ) , (2.9) 3 We could also in principle expand T (U, V ) in blocks, but then the OPE coefficients squared would not necessarily be positive, because the free theory correlator in the definition of T (U, V ) in (2.3) contributes to both long and short multiplets.

JHEP01(2023)038
where G ∆, (U, V ) with scaling dimension ∆ and spin are 4d conformal blocks (2.10) The long multiplets that appear in the S × S OPE must have even and satisfy the unitarity bound ∆ ≥ + 2. The short multiplet OPE coefficients do not depend on the coupling and so can be computed from the free theory to give the exact expressions F (0) short and F (1) (2.11) All non-trivial interacting information in the block expansion (2.9) of the correlator is then given by ∆ and for the long multiplets. This data is constrained by the crossing equation derived by swapping 1 ↔ 3 in (2.1), which in terms of the block expansion in (2.9) leads to =0,2,... ∆≥ +2 where we defined (2.13)

Integrated constraints from supersymmetric localization
Another source of nonperturbative constraints on the CFT data comes from the two integrated constraints discussed in the introduction, which relate certain integrals on S 4 of 4 For c < 3 4 , unitarity requires that certain short multiplets have modified OPE coefficients, and that in particular higher spin currents appear, so all such theories must be free theories [18]. We will not consider c < 3 4 in this paper. For c ≥ 3 4 , these formulae apply to both free and interacting theories, where for the free theory there are long multiplets that saturate the unitarity bound with OPE coefficients given in (2.22), at which point these long multiplets become conserved current multiplets as expected for a free theory.

JHEP01(2023)038
SSSS to derivatives of the mass-deformed sphere free energy F (m, τ,τ ). These constraints can be written in terms of the function T (U, V ) defined in (2.3) as [23,24] 1 8c (2.14) Here, 5 whereD 1,1,1,1 (U, V ) is a scalar 1-loop box integral that can be written explicitly as Using (2.8) and (2.9), the function T appearing in (2.14) takes the form (2.17) The left-hand sides of the constraints (2.14) are written in terms of the mass-deformed free energy F (m, τ,τ ). The partition function Z(m, τ,τ ) ≡ exp(−F (m, τ,τ )) was computed using supersymmetric localization in [12] for N = 4 SYM with gauge group G in terms of a rank(G)-dimensional integral. 6 We are mostly interested in the G = SU(N ) expression, where a ij ≡ a i − a j . Here, the integration is over N real variables a i , i = 1, . . . , N , subject to the constraint i a i = 0 imposed by the delta function, and H(m) ≡ e −(1+γ) where G is the Barnes G-function and γ is the Euler-Mascheroni constant. The Nekrasov partition function Z inst encodes the contribution from instantons localized at the north pole of S 4 , and can be conveniently expanded as inst (m, a ij ) representing the contribution of the k-instanton sector and normalized such that Z  in [34,35], 7 and the resulting integral of this expansion converges rapidly for any τ in the SL(2, Z) fundamental domain (2.20) In appendix A we provide more details about this calculation and explicit expressions for the SU(2) and SU(3) cases written as one-and two-dimensional integrals, respectively, that can be easily computed numerically for any τ in (2.20). We plot the localization inputs in figure 1 as a function of τ for SU(2) and SU(3). Figures 2 and 3 additionally show cross sections of these inputs along the imaginary axis and on the arc from e πi/3 to e 2πi/3 .

Weak coupling expansion
In the next section, we will find it useful to compare our results to the small g YM perturbative expansion, where instantons do not contribute. At g YM = 0, we have the free theory correlator (2.4), which, when expanded in superconformal blocks, gives long multiplets of twist t = 2, 4, 6, . . . with OPE coefficients [33]:
Later in [36,37], the SU(N ) instanton partition function was obtained by factoring out U(1) contributions Z inst U(1) motivated by the AGT correspondence [36]. Since Z inst U(1) is holomorphic in τ and independent of the ai, this contribution is killed by the m and τ derivatives that we consider, so we will simply use the results for U(N ) instanton partition functions in the following. 8 The leading instanton correction has also been computed in [38][39][40], and takes the form γ1-insnt =  The localization inputs for SU (2) for θ = 0 as a function of g 2 YM (above), and as a function of τ along the arc from e πi/3 to e 2πi/3 (below). We include terms with up to 10 instantons in (2.19), which as shown in appendix A is more than sufficient for convergence inside the fundamental domain. Note that the slope of these inputs goes to zero at the Z 2 self-dual point τ = i and at the Z 3 self-dual point τ = e πi/3 . where the 't Hooft coupling is λ ≡ g 2 YM N . At twist 4, there are generically four degenerate operators, but when N = 2 there are only two due to trace relations. At one loop, the scaling dimensions of these twist 4 operators whose scaling dimension are written in terms of a quantity ω which obeys a quartic equation: [41,42]

Numerical bootstrap with integrated constraints
In this section we will describe how to combine the numerical bootstrap of the stress tensor correlator with the two integrated constraints. We will begin by applying the functionals  . The localization inputs for SU (3) at θ = 0 as a function of g 2 YM (above), and as a function of τ along the arc from e πi/3 to e 2πi/3 (below). We include terms with up to 2 instantons in (2.19), which as shown in appendix A is sufficient for convergence inside the fundamental domain. Note that the slope of these inputs goes to zero at the Z 2 self-dual point τ = i and at the Z 3 self-dual point τ = e πi/3 . in (2.15) to each block G ∆, (U, V ) in the expansion of the correlator for any spin and dimension ∆. We will then describe the abstract bootstrap algorithms that can be used to bound the scaling dimensions and OPE coefficients of the unprotected operators that show up in the S × S OPE. We then discuss how to implement these algorithms in practice, and in particular why the conventional semidefinite programming approach [43] cannot efficiently impose both integrated constraints and the crossing equations, so instead we use an updated version of the original linear programming approach of [32]. We conclude with non-perturbative bounds on low-lying CFT data coming from these constraints for finite N and τ , which we compare to the weak coupling predictions.

Integrated constraints for block expansion
The integrals in (2.15) are over R 2 (written in polar coordinates (R, θ)), while the block expansion of a four-point function has a finite radius of convergence. Thus, we divide R 2 into domains that are permuted under the S 3 crossing symmetry of the four-point function, such that the conformal block expansion converges within each of these domains. This division is more clearly seen in the (z,z) coordinates used in the conformal block expansion, which are related to the polar coordinates (R, θ) used in the integrals in (2.15) by The S 3 crossing symmetry is generated by z → 1 − z (coming from the interchange of the first and third operators in the four-point function) and by z → 1/z (coming from the interchange of the first and fourth operators in the four-point function), along with their complex conjugates. 9 We define the three regions that are interchanged by the S 3 symmetry, and that are each marked with a different color in figure 4. 10 Out of the three regions, the s-channel block expansion that we considered in (2.9) converges in D 1 [44]. For any crossing-invariant quantity, we can restrict the integration range in (2.15) to one of these domains, or indeed any union of two of the fundamental domains under crossing symmetry, and then multiply by 3 to compensate. We will use a domain of integration denoted D (b), to be defined shortly. To motivate the definition, we first work out how to efficiently integrate the conformal blocks. 9 Defining by s1 the map z → 1 − z, by s2 the map z → 1/z, and by e the identity transformation, it is clear that s 2 1 = 1 and s 2 2 = 1. One can also check (s1s2) 3 = 1, which is the defining property of the permutation group S3. 10 Each of these regions is the union of two fundamental domains on which the S3 permutation symmetry acts.

JHEP01(2023)038
It is useful to expand the blocks such that one integral can be done analytically, leaving just one numerical integration. 11 The block expansion is most easily expressed in terms of the radial coordinates r, η defined in [46] as We can then use the recursion relation in appendix B 12 to expand the conformal blocks in a small r expansion as where U s (η) are Chebyshev polynomials, B n,s are numerical coefficients, and the block is normalized such that B 0,0 = 1. The integrals (2.15) that act on the blocks U −2 G ∆+4, (U, V ) that appear in the integrated constraints (2.14) can then be written as the r, η integrals whereD 1,1,1,1 can also be written in terms of r, η using its explicit expression (2.16) and the changes of variables from (z,z) to (r, η) obtained by combining (2.6) and (3.3). For a given ∆, , we can expand the blocks, the integration measures, andD 1,1,1,1 , 13 all at small r to some order p, perform the integrals in r exactly, and then the remaining integral in η numerically. The error from this expansion goes like r p max , where r max is the maximum value of r in our integration region D . So long as r max < 1, we see that the method converges quickly. For instance, in the region D 1 we have r max = 2 − √ 3 ≈ 0.268. However, this maximum value occurs at η = 0 where U s (0) = −(−1) s/2 , so if we integrate over D 1 then the integrated constraints will oscillate with spin. This presents a problem for the bootstrap, because for large ∆ the integrated constraints grow as (4r max ) ∆ while the constraints from crossing symmetry grow as (4(3 − 2 √ 2)) ∆ , where 3 − 2 √ 2 is the value of r at the crossing-symmetric point. For D 1 , we have r max > 3 − 2 √ 2, and so for large ∆ the integrated constraints dominate. The sign oscillations in spin would then prevent us from including the integrated constraints in any positive functional, and so we could not use the integrated constraints to improve the rigorous bootstrap bounds. 14 JHEP01(2023)038 We will thus follow the approach of [31] to modify the integration region such that the maximum value of r occurs at η = ±1 instead of at 0, thus removing the sign oscillations. We first extend D 1 to form a region D (b) defined by where α and β are defined by requiring the boundary to intersect the points (r, η) is then achieved at η = ±1 as desired. To correct for the expansion of the integration region, we should subtract off the image of D (b)\D 1 under 1 ↔ 3 crossing, which is bounded from below by the solution to and from above by the boundary of D 1 , namely (3.6) with α = 0 and β = 1. We define D (b) to be D (b) with this crossing image of D (b)\D 1 removed, as depicted in figure 5.
In the (r, η) coordinates, we can define this region by where r 0 (η) = |η| + 2 − η 2 + 4|η| + 3 (3.9) and r 1 (η, b) and r 2 (η, b) come from solving (3.6) and (3.7), respectively. Since the four-point function we are integrating is crossing-invariant, the choice of integration region should not affect the bounds we can achieve in the infinite-precision limit. However, as we have just mentioned, different regions can behave very differently for the finite truncations needed for the numerical bootstrap; this is why we choose a non-oscillating region like D (b) instead of D 1 . Furthermore, although constraints coming from two different integration regions are redundant in the infinite-precision limit, for our numerical bootstrap approach we find that including constraints from two different regions significantly improves the bounds. We thus use both D (1/3) and D (7/25) in what follows. Both of these values of b exceed 2 − √ 3, and thus do not exhibit the sign oscillations of D 1 , so we can use them together to build positive functionals. The values of the integrated constraints computed using these regions are shown in figure 6.

JHEP01(2023)038
In addition to the integrals of the conformal blocks, we need the integrals of the short contributions. We numerically integrate the closed-form expression for T short given in (2.17) for both of our integration regions, giving the following results:

Numerical bootstrap algorithms
Let us now explain the numerical bootstrap algorithm that incorporates both the integrated constraints and the crossing equations. For a given function f of the cross ratios (U, V ), we (3.11) The first four components of v f are the functionals I 2 and I 4 defined in (2.15) that appear in the integrated correlators, for both of the two non-oscillating integration regions we employ, acting on the function f . The remaining components correspond to the function We can denote by V ∞ the vector space consisting of the vectors v f of the form (3.11). The constraints (2.12) and (2.14) then impose that an infinite sum of vectors with positive coefficients λ 2 ∆, plus the contribution from the protected terms must vanish: where v ∆, is a shorthand notation for v f with f (U, V ) = F ∆, (U, V ). To find instances in which the equation (3.12) cannot be obeyed, consider a functional (written as an infinitedimensional row vector) where α 2,i and α 4,i are numbers, and α ∞ is a functional acting on functions of (U, V ).

JHEP01(2023)038
The functional α can be used to bound scaling dimensions ∆ by the following algorithm: Scaling dimension bound: which contradicts eq. (3.12), and so our assumptions ∆ ≥∆ must be false. If we cannot find such an α, then we conclude nothing.
For instance, we can set the lower bounds∆ to their unitarity values∆ = + 2 for all except a certain , then by varying∆ , c, and τ , this algorithm can be used to find an upper bound on∆ , i.e. on the scaling dimension of the lowest dimension operator with spin , as a function of c and τ . Without further assumptions we cannot do better than finding upper bonds on scaling dimensions. If we furthermore assume that only a few operators have dimensions below a certain value, then the allowed region as a function of the scaling dimensions of these operators could be more complicated. For instance, consider adding the additional constraint on α to the algorithm above: for some spin˜ and dimension∆˜ . One can then fix∆˜ and vary∆˜ . For instance, at weak coupling we saw in section 2.3 that there are only two relevant scalar unprotected operators, while at large N and large 't Hooft coupling λ there is only one such operator. Thus, we may assume that at any value of N and τ there are at most two unprotected scalar operators. We can impose this in our bootstrap study by setting a gap∆ 0 = 4, inserting two operators∆ 0 and∆ 0 below, and varying these latter scaling dimensions for given c and τ . We can also get upper bounds on a certain λ 2 ∆ , , and also lower bounds if we assume a gap, by a slight variation of the scaling dimension algorithm: JHEP01(2023)038 OPE coefficient bound: where s = ±1 for upper/lower bounds.
2. Assume that the scaling dimensions ∆ of all unprotected operators other than the one with spin and dimension ∆ obey lower bounds ∆ ≥∆ (e.g. unitarity).
If we just set∆ to the unitarity bounds, then we can get upper bounds but not lower bounds, because the potential existence of operators with ∆ arbitrarily close to ∆ makes step 1 inconsistent with step 3. To avoid this we must set a gap above unitarity for∆ and below this insert the operator with ∆ , as in (3.16). Note that this algorithm involved maximizing α[w], whereas the scaling dimension algorithm only involved the existence of α for fixed α[w].

Linear versus semidefinite programming
To implement the above algorithms numerically, we must perform three truncations / discretizations. The first truncation is the number of spins, which we can achieve by imposing a simple cap max on the range of spins we consider. We then check afterwards that the functionals we find are positive for all spins. The second truncation is on the space of functions of (U, V ), or equivalently on the size of the functionals used to probe this space. As is customary in conformal bootstrap studies, instead of functions of (U, V ), we consider a finite number of derivatives at the crossing symmetric point z =z = 1 2 . Thus, the vector space V ∞ is replaced with a finite-dimensional vector space V Λ defined by replacing a general vector with m + n = 1, 3, 5, . . . , Λ and m ≤ n.

JHEP01(2023)038
Here, we restrict to odd derivatives with m ≤ n because F ∆, (U, V ) (2.13) is symmetric in z ↔z and odd under crossing. 16 For sufficiently high max , the numerical bounds will monotonically improve with increasing Λ, which is what makes the numerical bootstrap rigorous. Lastly, one should discretize the continuum of ∆ ≥∆ in the space of allowed operators.
Without integrated constraints, the most efficient way of dealing with the continuum in ∆ is to use the semidefinite programming approach of [43]. The strategy is to use an approximation to the conformal blocks that makes the crossing equations polynomial in ∆.
The positivity constraints α[v ∆, ] ≥ 0 then reduce to the positivity of certain polynomials, which can be implemented using semidefinite programming software such as SDPB [50]. To write (2.12) as a polynomial in ∆, one uses the small r expansion of the blocks in (3.4) (truncated at some orderp) to write the where we used the fact that r = 3 − 2 √ 2 at the crossing symmetric point and the B n,s (∆, ) are meromorphic functions of ∆. We can then multiply the crossing equations by (4(3 − 2 √ 2)) −∆ as well as by the finite product of all positive factors that appear in the denominators of B n,s (∆, ) to get a polynomial in ∆.
The problem with using a similar approach here is that, as discussed in section 3.1, the integrated constraints scale at large ∆ as because the blocks scale as (4r) ∆ U (η) and b is the maximum value of r in the integration region D (b). Since the asymptotic scaling of the integrated blocks and the crossing equations differ by an exponential factor, writing both as polynomials in ∆ would require approximating an exponential by a polynomial. 17 While approximating an exponential by a polynomial seemed to give reasonable results in [31], in our 4d case we found that such an approximation made the bootstrap insensitive to the integrated constraints. We will thus avoid this problem with semi-definite programming by returning to the original linear programming method of [32]. In this case, we truncate each ∆ to some upper cutoff ∆ max , and then discretize ∆ with spacing ∆ sp , so the bootstrap algorithms become a finite set of linear constraints. We can then solve these linear constraints with linear programming software. In practice, we will still use SDPB for a linear programming problem. Since we have cut off ∆ max , the exponential scaling and oscillation of the integrated constraints relative to the crossing constraints do not matter as long as we use enough 16 We can efficiently compute ∂ m z ∂ n z F ∆, (U, V )| z=z= 1 2 using the scalar_blocks code, available online from the bootstrap collaboration [47][48][49]. Note that their blocks differ from our conventions as G ours ∆, (U, V ) = ( + 1)G theirs ∆, (U, V ). 17 One cannot make the crossing equations scale the same as the integrated constraints by expanding the former around a different point, because then the crossing equations at large ∆ would oscillate with spin as (−1) /2 , and cause a similar oscillation problem. We thank David Simmons-Duffin for discussion on this topic.

JHEP01(2023)038
numerical precision so that we are sensitive to all the constraints. Also, for sufficiently small ∆ sp and large ∆ max , the numerical bootstrap will become insensitive to the precise values of ∆ sp and ∆ max we use, similar to the cutoff on spin. See appendix C for more details on our numerical implementation.

Results
Using the method outlined above, we can obtain an upper bound on the scaling dimension ∆ 0 of the lowest scalar as a function of complex coupling τ . Figure 7 shows, for SU (2) and SU (3), the maximal upper bound on ∆ 0 that we obtain at any value of θ, as a function of g YM . For this plot and all others, we work at a bootstrap resolution of Λ = 39, as defined below (3.19).
For g 2 YM 1, our scaling dimension bounds are approximately saturated by the weak coupling result (2.23). In figure 7 we plot (2.23) to three and four loops, and a Padé approximant of order (2,2) to the four-loop expression. All of these expressions nearly coincide with our upper bound at very small coupling, and the Padé approximant tracks our bound over a somewhat larger range of coupling values.
In figure 8 the upper bound on ∆ 0 is plotted across the standard fundamental domain (2.20) in terms of the complex coupling τ , and extended via SL(2, Z) duality to the remainder of the upper half-plane. Within the standard fundamental domain, the dependence on θ is very weak. The bounds vary by less than 1% between θ = 0 and θ = π. The bottom panels of figure 8 show the bounds on the bottom boundary of the standard fundamental domain, namely along the arc of the unit circle between e πi/3 and e 2πi/3 .
As described above, we can also obtain bounds on the squared OPE coefficient of the lowest scalar operator, denoted by λ 2 0,∆ 0 . Figure 9 shows the upper bounds we obtain on this OPE coefficient as a function of g 2 YM for θ = 0. These bounds also appear close to being saturated by weak-coupling estimates. The agreement is less precise in this case because we have to input the scaling dimension of the lowest operator as estimated in figure 7.
We plot the bounds on the OPE coefficient as a function of τ in figure 10. Like in figure 8, our results are extended by SL(2, Z) duality to the remainder of the upper half-plane. We see again that the dependence on θ in the fundamental domain is weak. The bottom panels of figure 10 show how the OPE coefficient bound varies on the arc between e πi/3 and e 2πi/3 .
Without imposing a gap, we are not yet able to obtain lower bounds on any of these pieces of CFT data. This is because, at the highest bootstrap resolution we have attempted (Λ = 39), theories with a single relevant scalar seem to still be allowed. Let ∆ 0 (τ ) denote an allowed scaling dimension for this scalar. If we then insert two relevant operators, the second operator could have any scaling dimension so long as the first operator has scaling dimension ∆ 0 (τ ). In particular, the second operator could have a dimension as low as the unitarity bound, thus preventing us from obtaining lower bounds.
This behavior is shown in figure 11 in the SU(2) case for θ = 0. For three values of g 2 YM , we plot the allowed region in the space of the lowest two scalar operator dimensions, denoted ∆ 0 and ∆ 0 , respectively. The gap on all other scalar operators in these plots is set to 4, so that the presence of allowed points with ∆ 0 = 4 indicates the feasibility of having a  The bounds are approximately saturated by the weak-coupling estimate (2.23) in the regime where it is valid. We also show the upper bounds that would be obtained at the same bootstrap resolution without using the integrated constraints. Note that, like the localization inputs, the slope of these upper bounds go to zero at the Z 2 self-dual point where g 2 YM = 4π.  The dependence of the bounds on θ is weak within the fundamental domain, varying by less than 1% between θ = 0 and θ = π. Note that, like the localization inputs, the slope of these upper bounds go to zero at the Z 2 self-dual point τ = i and the Z 3 self-dual point τ = e πi/3 . single relevant scalar. As a result of this fact, we see that there are also allowed points with ∆ 0 as low as 2, and so we cannot establish a lower bound. Nevertheless, if we knew that ∆ 0 is always above some threshold, then we could obtain a lower bound on ∆ 0 . For instance, when g 2 YM 1, we can use (2.24) to estimate the second-lowest scalar. It will be below 4, but barely so, and thus by imposing a gap slightly below 4 we should obtain lower bounds that include the weak-coupling result (2.23). Figure 12 shows the lower bounds on ∆ 0 we obtain for various gap assumptions for the SU(2) theory at θ = 0. If the second-lowest scalar is greater than one of these gaps for all coupling values, then the bounds corresponding to that gap are rigorous; otherwise, the bound applies only where the second-lowest scalar is above the gap value. In the  YM at θ = 0, obtained with bootstrap resolution Λ = 39. We plot from g 2 YM = 0 to the self-dual point at g 2 YM = 4π. The weak-coupling estimates for these coefficients are included, showing that our bounds are nearly saturated by these values at small g 2 YM . We also show the upper bounds that would be obtained at the same bootstrap resolution without using the integrated constraints, and without using bounds on the scaling dimensions obtained with integrated constraints. Note that, like the localization inputs, the slope of these upper bounds go to zero at the Z 2 self-dual point g 2 YM = 1.

JHEP01(2023)038
SU ( bottom panel of figure 12, we zoom in on the weak-coupling region to show that the Padé approximant that appears in figure 7 is indeed included within the lower bounds obtained by this method.

Discussion
In this paper we have shown how to combine two integrated constraints from supersymmetric localization with the infinitely many constraints from crossing symmetry and unitarity to numerically bootstrap the stress tensor multiplet correlator SSSS of N = 4 SYM as a function of the complexified coupling constant τ . Since the localization inputs are easier to compute for small rank gauge groups, we restricted our analysis to the cases G = SU(2) and SU (3). For these cases, we determined upper bounds on the scaling dimension ∆ 0 and OPE coefficient of the lowest dimension scalar unprotected singlet operator. At weak coupling, we found that our upper bounds are close to being saturated by the four-loop weak coupling estimates [4][5][6][7][8] for a remarkably large range of couplings. Interestingly, we found that the maximal values of our scaling dimension bounds occurs at the Z 3 point τ = e iπ/3 , and that these maximal values are appreciably smaller than the conformal bootstrap bounds obtained without inputting the integrated constraints. Lastly, we also found preliminary evidence for the entire range of τ that after imposing a gap ∆ 0 above ∆ 0 , a lower bound appears that meets the upper bound as ∆ 0 is increased.
Looking ahead, our first goal is to upgrade these upper bounds to islands for every G and τ . As was the case for the O(N ) model bootstrap [48,[50][51][52][53][54], we expect that putting a gap to the spacetime dimension (4 in our case) and inserting the number of relevant operators below should give precise islands for the scaling dimensions of these operators. For the N = 4 SYM theory, weak coupling and large N estimates for G = SU(N ) suggest that there are at most two relevant operators for any N and τ . In order to get islands, the bootstrap must be able to accurately restrict the range of ∆ 0 to an interval below 4 for all τ , which we have not yet achieved. It is possible (but not likely) that we might accomplish this by simply increasing the parameter Λ that measures the size of our functionals. Alternatively, we could consider a mixed correlator between S and the dimension 3 half-BPS primary S 3 , as was considered without integrated constraints in [20,21]. Not only will mixed correlators improve the bootstrap accuracy by allowing us to impose the uniqueness of the relevant operators as they appear in multiple correlators, which proved very effective for the O(N ) bootstrap [52], but this setup will also give us access to the integrated constraint on SSS 3 S 3 as derived in [23], where the localization input could be computed for any G = SU(N ) and τ following [35,[55][56][57]. In fact, there are infinite possible integrated constraints that we could access by considering SSS p S p for each integer p, where S p is the dimension p half-BPS operator considered in [23].

JHEP01(2023)038
In order to efficiently apply our algorithm to these more complicated correlators, we will need to improve the numerical implementation. Instead of discretizing ∆ up to some ∆ max and applying linear programming, one could directly apply semidefinite programming to the continuous interval of ∆ between the unitarity bound and ∆ max , as was previously explored in [58]. 18 This will allow us to use the full power of SDPB, instead of treating it as a linear programming software in our current approach.
We would also like to apply our algorithm to any gauge group G, not just for G = SU(2) and SU(3). Since brute force evaluation of the matrix model integral for the localization input involves rank(G) integrals, we will need more compact expressions. The ∂ τ ∂τ ∂ 2 m F | m=0 input was already computed for G = SU(N ) at finite N and τ in [59], but we will also need a finite N, τ expression for the other input ∂ 4 m F | m=0 . It would also be nice to generalize these finite N, τ expressions to general gauge group along the lines of [21], and to the localization inputs that apply to SSS p S p in [23].
Finally, we would like to emphasize the broad applicability of our new method of combining numerical bootstrap with integrated constraints from localization to study SCFTs. This method should apply to any SCFT where the sphere free energy is computable from supersymmetric localization, which includes 3d N = 2 [60,61], 4d N = 2 [12,13], and 5d N = 1 [62,63] CFTs. In particular, 4d N = 2 SQCD also has a conformal manifold with one complex parameter, and we expect to be able to derive integrated constraints on correlation functions of some of the half-BPS operators. Another theory that would be interesting to study is 3d N = 6 ABJ(M) [64,65] with gauge group U(N ) k × U(N + M ) −k and Chern-Simons level k, which has dual descriptions in either M-theory, string theory, or higher spin gravity depending on the values of N , M , and k. In [66], we combined the numerical conformal bootstrap with two protected OPE coefficients, which could be computed from supersymmetric localization [67,68], to fix a 2-parameter subspace of {k, N, M }. To fully fix the theory, we would need a third constraint, which could come from the integrated constraints derived in [69,70]. This would allow us to compute small precise islands for each N, M, k, as was already done for the N = 8 case with k, M = 1, 2 in [71]. For the N = 8 case itself, adding this integrated constraint as well as another possible constraint coming from the squashed sphere free energy [72,73] could make the islands so precise that we could read off large N corrections to the correlator, as initiated in [74,75], which would allow us to derive the M-theory S-matrix following [76]. We have truly only begun to explore the power of combining the conformal bootstrap with supersymmetric localization.

A Supersymmetric localization
In this appendix, we will compute the mass derivatives of the SU(N ) N = 2 * sphere free energy F (m, τ,τ ) that appear in the integrated constraints (2.14) as defined in (2.21). We start by considering the massless partition function, where the instanton term cancels so that we simply get where the second line was evaluated using orthogonal polynomials in [28] and G(x) is the Barnes-G function. We can then write F 2 as an (N − 1)-dimensional integral: where K(z) ≡ − H (z) H(z) and the expectation value with respect to the unit-normalized massless partition function is We can similarly write F 4 as To compute the mass derivatives of the Nekrasov partition function Z inst (m, τ, a ij ), we use the explicit expression for the Z (k) inst (m, a ij ) in (2.19), which as shown in [34,35] is given by a sum over N -tuples of Young diagrams Y = (Y 1 , Y 2 , . . . , Y N ) with k boxes in total:

JHEP01(2023)038
Here, s labels a box (α, β) (α-th column and β-th row) in a given Young diagram, while h i (s) and v i (s) denote the arm-length and leg-length, respectively, of the box s in the diagram Y i . Each individual Young diagram Y consists of columns of non-increasing heights λ 1 ≥ λ 2 ≥ · · · ≥ λ l with λ l ≥ 1. The transpose (conjugate) diagram Y T has columns of heights λ T 1 ≥ λ T 2 ≥ · · · ≥ λ T m with λ T m ≥ 1. Then the arm-length h and leg-length v of the box s in Y are given by Note that the definitions of h and v extend beyond boxes in Y to the entire quadrant (α, β) ∈ Z 2 + in the obvious way, so they can be negative (e.g. when Y is empty). The integrals in these expectation values in (A.2) and (A.4) can then be computed numerically for low N , where recall that there are N − 1 integrals because we set a N = − N −1 i=1 a i as required by the SU(N ) measure. For numerical stability, it is useful to change integration variables so that a i − a j terms do not appear in the Vandermonde determinant, and also to rescale the variables by g YM so that it does not appear in the Gaussian factor. We will now show some explicit results for the SU(2) and SU(3) cases that we consider in the main text.

A.1 SU(2)
For SU(2), after setting a 2 = −a 1 = g YM x, we find that inst (m, a ij )| m=0 for the lowest couple k are: inst (m, a ij )| m=0 we get k = 1 : 12 In both cases it is easy to compute higher values of k, but we found that the numerical integrals for τ in the fundamental domain (2.20) converged quickly even with just a few terms. For instance, at the Z 3 self-dual value τ = e iπ/3 we get including k = 0 : including k ≤ 1 : including k ≤ 2 : including k ≤ 3 :

B Conformal block expansion
In this appendix we describe how to efficiently compute the small r expansion of 4d conformal blocks G ∆, (U, V ). We start by defining the variables which are simply related to z andz as z = 4ρ (1 + ρ) 2 ,z = 4ρ (1 +ρ) 2 . (B.2) Expanding in small ρ,ρ thus coresponds to expanding at small r.
Recall that the 4d conformal blocks are written in terms of the functions k β (z) defined in terms of hypergeometrics in (2.10). From the standard series representation of the hypergeometric, we can define a recursion relation for k β (z) in terms of ρ as where u ∞ (ρ) = 1 1 − ρ 2 to efficiently expand each term in (2.10) at small r to get the block expansion (3.4).

C Numerical bootstrap details
As described in section 3.2, our bootstrap calculations reduce to solving linear programs, following the setup in [32]. We search through a space of functionals α parametrized by The parameters α m,n , along with α 2,1 , α 2,2 , α 4,1 , and α 4,2 which define the action of the functional on the integrated constraints, are the variables of the linear program. The constraints of the linear program implement the positivity conditions (3.14). For each spin , we wish to include constraints corresponding to a continuum of possible operators from some gap∆ all the way to infinity. To make the problem finite, we consider operators with dimensions spaced by some separation ∆ sp , and up to some cutoff ∆ max . We also include one additional constraint for each spin at ∆ = 500, which approximates the asymptotic large-∆ behavior.
In practice, once we make ∆ sp sufficiently small and ∆ max sufficiently large, our set of constraints approximates a continuum well enough that the results are no longer sensitive to the values of these parameters. In all the plots shown in this paper, ∆ sp = 0.1 is small enough to see this behavior, and ∆ max = 100 is large enough.
In figure 13, we give an example of a positive functional calculated at Λ = 39, giving an upper bound of 2.8688 on the lowest scaling dimension for the SU (2) theory at the Z 3 symmetric point. This corresponds to the rightmost point of the SU (2) plot in figure 7, and is the maximum upper bound we find throughout the conformal manifold. The bound we obtain at this point is insensitive to increases in ∆ max or the maximum spin, or decreases in ∆ sp .
We used two linear program solvers, Gurobi and SDPB [50,77]. Gurobi is heavily optimized for solving linear programs, and is very efficient for obtaining results at low Λ. However, it can only work at machine precision. For high enough Λ, the large dynamic range of conformal block derivatives and integrated constraints leads to numerical instability. This is most straightforwardly corrected by working at higher precision. When this is required, we use SDPB with a precision of 256 bits, and all other parameters set to their default values. The parameters used are summarized in table 1.
Empirically, we find that Gurobi is only reliable at Λ ≤ 29. We use Gurobi to estimate the allowed region of scaling dimensions using lower values of Λ, and then use SDPB to obtain the exact allowed region at Λ = 39. All the bounds in this paper are obtained via binary search with a resolution of 2 × 10 −4 , except for plots on the bottom arc of the fundamental SL(2, Z) domain for which we increase the resolution to 2 × 10 −5 . Positivity is imposed only in the yellow shaded region, but we find that the functional is positive for all ∆ ≥∆ for each spin . (b): the functional above evaluated only for spins 0, 2, and 4, showing in more detail the low-lying extremal spectrum for these spins. Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited. SCOAP 3 supports the goals of the International Year of Basic Sciences for Sustainable Development.