The S-matrix Bootstrap III: Higher Dimensional Amplitudes

We consider constraints on the S-matrix of any gapped, Lorentz invariant quantum field theory in 3+1 dimensions due to crossing symmetry, analyticity and unitarity. We extremize cubic couplings, quartic couplings and scattering lengths relevant for the elastic scattering amplitude of two identical scalar particles. In the cases where our results can be compared with the older S-matrix literature they are in excellent agreement. We also extremize a cubic coupling in 2+1 dimensions which we can directly compare to a universal bound for a QFT in AdS. This paper generalizes our previous 1+1 dimensional results of arXiv:1607.06109 and arXiv:1607.06110.

In [1] and [2] we initiated a bootstrap analysis of massive quantum field theories. In particular, we obtained bounds on couplings of a quantum field theory compatible with a given spectrum of stable particles.
Physically, one expects such bounds to exist since increasing the interaction strength will typically increase the attraction between particles. As such, we expect to have maximum values for couplings beyond which the masses of bound states must decrease, or new boundstates should emerge from the continuum, or both.
Mathematically, this problem is also very natural once we make the non-trivial assumption that scattering amplitudes are described by functions that are analytic away from the usual physical poles and cuts. The point is that analytic functions always attain their maximum at a boundary of their domain of definition. In the context of scattering amplitudes, these boundaries are the cuts generated by multiparticle intermediate states. For physical kinematics the amplitude along the cut is constrained by the conditions that probabilities add up to one -i.e. by unitarity. For this reason we focus on the two body scattering of the lightest particle in the theory since then all the usual cuts of the amplitude correspond to physical kinematics. In 1 + 1 dimensions where unitarity can be directly applied at the level of the S-matrix (simply, |S(s)| ≤ 1 for s along the cuts) we are faced with a clean problem in the theory of complex functions of a single variable. As we have an analytic function on a domain with a boundary along which it is bounded, so we are able to constrain its values inside this region and in particular the various physical couplings which we define as residues of factorization poles. Section 2 contains a derivation of the two dimensional bound which is a significant refinement of that in [2].
In this paper we move the focus to higher dimensions which contains a plethora of very interesting and difficult elements absent in the simpler 1 + 1 dimensional case. An essential difference is that the most convenient way to formulate unitarity requires introducing partial waves and these are not bounded by unitarity along their entire boundary (only along the so-called "right cut"). Therefore the simple complex analysis argument of 1 + 1 dimensions cannot directly apply. Furthermore, the analyticity and crossing symmetry requirements involve the amplitudes rather than the partial waves, which forces one to use both descriptions of the scattering event. Still, it is possible to overcome these technical obstacles. We shall introduce a kind of uniformization coordinates where the full space of physical kinematics is mapped to (a few) unit circles. This will allow us to Taylor expand the amplitudes in a convergent and manifestly crossing symmetric way in the full physical plane and then to numerically impose unitarity along the physical boundaries.
We start by revisiting the two dimensional results with this new approach in section 2, setup the higher dimensional problem in section 3 and present and analyze the corresponding numerical results in section 4. In section 4.4 we compare our numerical results with the completely orthogonal approach of [1] which is based on QFT in AdS and in particular does not require any analyticity assumptions. We conclude in section 5. A number of appendices are included to complement the main text presentation. Figure 1: Mapping from the cut s-plane to the unit disk given in equation (1). The mapping associates the points z(2 + iy) = z(2 − iy) and maps the half plane Re(s) > 2 to the full unit disk. The grey, dashed curves on the left map to those on the right and are included to help the reader visualize the mapping.

Two Dimensions Redux and Unit Circles
In this section we revisit the much simpler two dimensional problem. In two dimensions we can solve things analytically, and so it is a great training ground for developing intuition and testing any new numerical approaches. Nonetheless, for the braver readers eager to learn about the higher dimensional story, this section can be skipped without compromising the logic of the paper.
Most of the mathematical analysis of [2] boils down to minor variations of the following simple problem: Q: Consider all real analytic functions f (z) = [f (z )] with no singularities inside the unit disk apart from a simple pole at z = 0 and which are bounded on the unit circle as |f (e iφ )| ≤ 1. 1 What is the maximum possible residue at z = 0 and which function has that residue?
A: The maximum residue is 1 and the corresponding function is f = 1/z.
Indeed g(z) = f (z)/(1/z) has no singularities inside the disk and obeys |g(z)| ≤ 1 at the boundary of the unit disk. By the so-called maximum modulus principle, it satisfies |g(z)| ≤ 1 everywhere inside the disk. Its value at the origin -which is nothing but the residue of fis therefore at most 1. This maximum value is attained when g is constant everywhere, that is when g(z) = 1 corresponding to f (z) = 1/z.
To see how this simple problem relates to the analysis in [2,3] consider the 2 → 2 Smatrix S(s) for scattering of identical neutral particles of mass m considered in [2]. Assume also that there is a single bound-state showing up in this S-matrix element and for simplicity assume its mass m b > √ 2m. Because of crossing symmetry S(s) = S(4m 2 − s) and we can focus on the region Re(s) > 2m 2 without any loss of generality. In this half plane we have a threshold cut starting at s = 4m 2 , the bound-state pole at s = m 2 b and no other singularities. Consider then the change of variable (1) which maps this half plane into the unit disk, the bound-state pole into the origin of that disk and finally the threshold cut -where unitarity is to be imposed -to the boundary of the disk, see figure 1. In terms of z the S-matrix is therefore exactly constrained by the conditions of the previous point; it has a pole at z = 0 and obeys |S(z)| ≤ 1 at the boundary of the disk. 2 Its maximum residue -which is where we measure the (square of the) coupling to the bound-state -is therefore 1 and the corresponding optimal S-matrix is therefore S(z) = 1/z.
To recover the results of [2] -see e.g. formula (36) therein -we simply need to take into account the Jacobian to go from z to s, the simple kinematical multiplicative factors relating the S-matrix and the T-matrix and a factor of m 4 to render the coupling dimensionless. All other results of [2] for more complicated bound-state spectra can be treated through simple generalizations of this simple example! 3 Although redundant at this point, it is instructive for what will come next in higher dimensions to set up this exactly solvable problem numerically. We define a function S(z) in the unit circle as a pole plus a convergent Taylor expansion which we truncate at some large power z M . Then we simply maximize the residue with the constraint that in a tightly spaced grid of K points on the unit circle unitarity is satisfied. In Mathematica, the simple code below does the job: This nicely yields residue 1 and c n 0 with great numerical accuracy which can be always improved. The reader is encouraged to copy/paste this and try by him/herself. It should take about 2 or 3 seconds to run.
As a last warm-up it is very useful to solve this very same problem in a third way since this last approach is the closest to what we will do in higher dimensions. In this last approach 2 Note that this condition also holds on the lower half of the disk due to real analyticity. 3 Strictly speaking the map to the unit circle is not even needed here. It suffices to assume there is no essential singularity at infinity so that the unitarity cut is the boundary of the region where S(s) takes values. Then S(s)/z(s) −1 is free of singularities in the physical region and obeys |S(s)| ≤ 1 on the cuts which are the boundaries of this region. Hence it can at most be one inside by the maximum modulus principle and the bound on the residue of S follows. This is the argument in [3]. We still found the unit circle discussion to be useful as a warm-up to the higher dimensional case. Figure 2: Mapping from the cut s-plane to the unit disk given in equation (2).
to the problem we start by thinking of the S-matrix as being a function of both s and t as if they were independent variables; they are not since s + t + u = 4m 2 and u = 0 in two dimensions. 4 Then S(s, t) is a function with a cut for s > 4m 2 , another cut for t > 4m 2 as well as poles for single-particle processes in the s-and t-channels. Next we use a very convenient change of variable which maps the full complex plane with those cuts removed into the unit disk. This is the map where s 0 < 4m 2 is a free parameter that we can choose according to convenience. In the present case, it is convenient to choose s 0 = 2m 2 so that ρ s = 0 corresponds to the crossing symmetric point s = t = 2m 2 . A similar map is also very useful in conformal bootstrap studies [4]. It is illustrated in figure 2. The top of the cut maps to the upper boundary of the unit disk and the bottom of the cut maps to the lower boundary of the disk. The interval [0, 4m 2 ] maps to the interval ρ ∈ 2 √ 2 − 3, 1 so this is where we find the poles associated to stable particles.
Apart from the poles corresponding to single particle exchanges, S(ρ s , ρ t ) is analytic for both ρ s and ρ t inside the unit disk and thus we can write Crossing symmetry is guaranteed provided the coefficients of the convergent Taylor expansion are symmetric, c ab = c ba . Since we are going to evaluate the S-matrix on the constraint surface s + t = 4m 2 we can simplify this ansatz further. In terms of ρ s and ρ t this constraint yields This means the representation (3) has a big redundancy. We can always add to it polynomials in the left hand side of the constraint (4). To remove this ambiguity, we can set to zero many of constants c ab (in appendix B we explained in detail which c ab can be set to zero).  Figure 3: Comparison of the exact optimal S-matrix (given by 1/z(s) with z given by (1)) to numerical results using the ansatz (3) with the a, b series truncated at maximum degree N = 5 and m 2 BS = 3. We plot the physical region ρ = e iφ with φ ∈ [0, π). The numerical results (red dashing) are indistinguishable from the exact results.
Numerically, we set a cut-off in the sum (3) and impose unitarity for s > 4 which corresponds to the upper half circle where ρ s = e iφ with φ ∈ [0, π]. We evaluate |S(s, t)| 2 in a uniform grid in the φ interval which gives a set of quadratic constraint equations on the c ab and the residues of the poles. We optimizeĝ 2 in the usual way using FindMaximum for example. The outcome of this third approach is in perfect agreement with our previous analytical and numerical results as illustrated in figure 3.
To summarize: In two dimensions we can find the optimal S-matrix with largest possible residue analytically. 5 We do so by dividing the S-matrix by a clever guess and using the maximum modulus principle to show that this ratio should be one. We recovered the same analytic results numerically in two ways. In the first one we start from a parametrization of the kinematics where we can Taylor expand the S-matrix and then truncate that expansion to obtain a finite algebraic problem which we can put on a computer. The second numerical approach is a small variation where we think of the S-matrix as a function of s and t as if they were independent and then consider a double Taylor expansion in each of them.
What we implicitly used in the last method can be called an analytic extension -note that it is not an analytic continuation as we are increasing the number of variables and not just moving into the complex plane keeping the number of variables fixed. In this extension we promoted the S-matrix to a more general function of two variables which has 5 Notice that if we allow essential singularities at s = ∞ then there is no upper bound onĝ 2 . To see that consider the ansatz where δm 2 = 1 2 min(4m 2 − m 2 b , m 2 b ). For any value ofĝ 2 , we can find a (large) positive integer n such that this ansatz satisfies the unitarity constraint |S(s, 4m 2 − s)| ≤ 1 for s > 4m 2 . We thank Etienne Granet for raising this point. We exclude such essential singularities at s = ∞ because they are incompatible with causality (see for instance appendix D of [5]). no singularities in the cut s and t planes. 6 Equivalently, in terms of the ρ variables, we assumed the existence of an extension into a function S(ρ s , ρ t ) which has no singularities in the polydisk {ρ s , ρ t such that |ρ s | ≤ 1 and |ρ t | ≤ 1} while all we know a priori is that such a regular function exists only in the intersection of the polydisk with the constraint (4). Why do we have the right to assume that such an extension exists at all? For instance, it could happen that such an extension would inevitably introduce new singularities in the full polydisk domain which would then invalidate the convergence of the double expansion (3). Numerically, using this extension method we seem to find perfect agreement with the analytic results so somehow we should be safe. Indeed, the polydisk is a so-called Stein manifold 7 and the constraint (4) is an holomorphic embedding and as such defines a submanifold inside the polydisk which is also Stein. As discussed in greater detail below, there is a rather remarkable mathematical result which states that regular analytic extensions from Stein sub-manifolds inside Stein manifolds to the full Stein manifold do exist! The perfect numerical agreement is thus to be expected.
Of course, in two dimensions this discussion is a clear use of excessive force. On the other hand, in higher dimensions we will also make use of such analytic extensions and there we will not have the luxury of the analytic results to cross-check our numerics. The theorem alluded to above generalizes to that case as well and is key in providing confidence for the higher dimensional numerics.
There is also another more pedestrian explanation of why the double Taylor expansion numerics had to work which we present in appendix A; however, contrary to the discussion above, it makes use of particular features of the two dimensional problem and is not that useful as a warm up to the higher dimensional case.

Higher Dimensions
We now move on to scattering amplitudes in d + 1 spacetime dimensions. Consider again the elastic scattering process of two identical real scalar particles of mass m. In our conventions the S-matrix element is with normalization such that where E p = m 2 + p 2 . The Mandelstam invariants are given by 6 Of course we still have the poles associated with stable particles but these can be easily treated separately as in (3). Here, we focus on the parametrization of the analytic part of the S-matrix. 7 The unit disk is an open Riemann surface and those are Stein manifolds. Products of Stein manifolds are also Stein so the polydisk is also Stein.
which of course obey s + t + u = 4m 2 , and we henceforth work in units such that m 2 = 1. We often write M (s, t) ≡ M (s, t, 4 − s − t). In the channel under consideration s is the squared center-of-mass energy and the scattering angle is given by Physical values of the Mandelstam invariants are therefore 4 ≤ s and 4 − s ≤ t ≤ 0. We can project onto channels with definite angular momentum by introducing the partial amplitudes: where P (d) (x) is proportional 8 to the Gegenbauer polynomials. In our conventions, with P (x) the usual Legendre polynomials, normalized such that P (1) = 1. We note that S (s) = 1 for odd because Bose symmetry implies invariance under the reflection θ → π −θ.
Although the S-matrix element (6) has all kind of distributional properties, the amplitude M (s, t, u) is a regular function (see e.g. [6, section 4.3]). We will assume that M (s, t, u) obeys three further constraints: • Crossing Symmetry: M (s, t, u) is completely symmetric in its arguments. The symmetry u ↔ t follows from the aforementioned Bose symmetry, but the other generator of the crossing symmetry group can only be found from a more sophisticated analysis and requires the LSZ prescription.
• Analyticity: M (s, t, 4 − s − t) is analytic for arbitrary complex s and t, except for potential bound-state poles at s = m 2 b with 0 < m 2 b < 4, a cut along the real axis starting at s = 4, and the images of these singularities under the crossing symmetry transformations. It further obeys the usual reality condition M (s * , t * 4 − s * − t * ) = M * (s, t, 4−s−t). We note that the analyticity assumption is actually rather optimistic, since this 'maximal' analyticity has not been proven from axiomatic field theory. 9 On 8 In general spacetime dimension, we have Certain analyticity properties are known to be valid very generally, derived either to all orders in perturbation theory or from axiomatic field theory; the latter case sometimes requires the Wightman axioms and other times merely requires the validity of the LSZ prescription and causality. Typically one can prove two-variable analyticity for all s (modulo the known poles and cuts) but only for some finite range of values of t or of x which in particular includes the physical values. A standard result is that the proven analyticity is sufficient to analytically continue the amplitude from the s-channel to the t or u channels, establishing crossing symmetry [7]. We refer to [8,9] and references therein for more extensive discussions. the other hand some a posteriori justification is provided by the remarkable agreement between some of our results and those obtained without maximal analyticity in the older literature. We therefore believe that this assumption is sufficiently mild to generate physically meaningful results. We offer some further comments on this point in section 4.4 and the conclusions section below.
• Unitarity: From S † S = 1 we find that the unitarity constraint for elastic scattering takes the form |S (s)| ≤ 1 for all s ≥ 4 and ∈ {0, 2, 4, . . .}. Generically no other channels are available for a finite window of values of s, starting at 4 and ending at a higher threshold (like s = 9 for three-particle scattering). In such a window the above inequality should in fact be saturated. In this work we will not impose such saturation, but our numerics in principle allows for it.
The aim of the S-matrix bootstrap program (as we envisage it) is to use these general conditions to obtain concrete constraints on the behavior of the function M (s, t, u) or the partial amplitudes S (s) at interesting points. Many results from the previous century can be found in the textbook [10] and the reviews [8,11].
The recent works [12,13] pursue a bootstrap analysis of scattering amplitudes of weakly interacting higher spin theories, where the amplitudes are meromorphic functions of the Mandelstam invariants. Analytically, they beautifully explore the large s and t regime of weakly interacting higher spin scattering amplitudes and observe remarkable universality there. In contrast, our analysis is fully non-perturbative and the only poles of the scattering amplitudes are associated with stable particles (below the 2-particle continuum). Nevertheless it would be very interesting to investigate the same large s and t regime within our numerical approach.

Ansatz
In this subsection we explore the consequences of our analyticity assumption in some detail. As a toy model we can start with a single-variable function f (z) which is analytic in a simple domain D ⊂ C. If we define ρ : D → ∆ as a biholomorphic map between D and the unit disk ∆ = {ρ ∈ C : |ρ| < 1}, then any such f (z) has a Taylor series expansion of the form which converges as long as |ρ(z)| < 1. Our multi-variable problem is unfortunately not so easy, since for M (s, t) the moving cuts imply that the domain of analyticity in one variable, say s, depends on the other variable t. We will remedy this as follows. First we relax the constraint s + t + u = 4 and consider three-variable functions M (s, t, u). Then we transform the variables (s, t, u) → (ρ s , ρ t , ρ u ) using the map (2) which is, with m 2 = 1, In this case, it is convenient to choose s 0 = 4 3 so that ρ s = ρ t = ρ u = 0 corresponds to the crossing symmetric point s = t = u = 4 3 . Now, since the transformation ρ s maps the s-plane minus the right cut starting at s = 4 to the unit disk, we see that in the ρ variables all the cuts lie outside the polydisk ∆ 3 defined by |ρ s | < 1, |ρ t | < 1 and |ρ u | < 1. The only remaining singularities are then the poles and it is natural to write where the triple ρ series converges inside ∆ 3 , and for definiteness we have put in the poles for a single scalar bound state of mass m b . The demands of crossing symmetry are implemented by demanding that the coefficients α abc are totally symmetric in their indices. When restricted to the surface defined by s + t + u = 4 the ansatz (15) obeys the analyticity and crossing symmetry constraints. It is perhaps more surprising that the converse is also true: any function obeying the analyticity constraints on the surface s + t + u = 4 can be extended to a function on ∆ 3 , analytic modulo the poles, and therefore can be written in the form (15). This follows from a mathematical theorem known as Cartan's theorem B, which is a statement about the vanishing of higher cohomologies of coherent analytic sheaves on Stein manifolds (see e.g. [14]) -in the case at hand this implies that there is no obstruction to an extension away from the surface s + t + u = 4. 10 The triple ρ expansion in equation (15) is the starting point for our numerical work. Our approach is to restrict the expansion to a finite sum by imposing and then further restricting to the constraint surface s + t + u = 4 which is given by a polynomial equation and which in practice allows us to eliminate many terms in (15) (in appendix B we explain in detail which terms can be set to zero). The remaining freedom in our ansatz then consists of the finitely many remaining α abc together with the bound state parameters; since this is a finite-dimensional space we can use a computer to numerically explore the space of scattering amplitudes. Of course we want to keep N max as large as possible. As we will see, in fortunate cases the numerical results stabilize already for feasible values of N max , while in other cases we can extrapolate. 11 It will be the job of the computer to impose the unitarity constraints, which are quadratic constraints in the parameters g 2 and α abc . Rather than checking the infinity of constraints for all s and , we impose a cutoff and check that unitarity constraints are obeyed only for ≤ max and along a grid of values for s. Experimentally we observe that our results remain meaningful if max is not much smaller than N max and if the grid is sufficiently refined. In appendix F we discuss the dependence on these parameters in more detail, and outline the numerical implementation.

Results
In this section we present our numerical results for several maximization problems using the S-matrix bootstrap method explained above. For most of this section we restrict our attention to 3+1 dimensional QFTs, i.e. d = 3 in our notation. In the final subsection 4.4, we consider 2 + 1 dimensional QFTs.

Cubic coupling
For our first result we consider a scattering amplitude with a single pole corresponding to the exchange of a scalar particle of mass m b , exactly as in our ansatz (15), and maximize the value of the residue g 2 as a function of m b . 12 In figure 4 we plot the maximum absolute value of the coupling |g| defined as the residue of the pole, with the different curves corresponding to different values of N max . We have obtained this plot by maximizing |g| for a sequence of values of m b and the indicated curve is an interpolation through our data points. The plot is rather rich; we discuss its key features one by one.
• Convergence with N max . For m b √ 2 we see that |g| max is nearly stationary as we vary N max , whereas for m b √ 2 we observe more significant improvements with N max . We have no explanation for this disparate behaviour (although we suspect it to be related to some subtler higher energy behaviour to which our ansatz is struggling to converge -see also discussion section 5 and appendix G). Numerically we find that we can extrapolate to infinite N max and appear to get a finite answer in either domain. We expect this value to correspond to an upper bound on |g| for any scattering amplitude that obeys the constraints of the previous section. 13 The clear peak is reminiscent of two-dimensional scattering amplitudes, where it was easily explained because in that case the s-and u-channel poles cancel precisely at m b = √ 2 and the number |g| becomes meaningless -so no upper bound can be obtained. 14 In greater than two dimensions the cross-channel poles are smeared into a cut by the projection onto the partial waves. One can easily see from (10) that this cut starts at s = 4 − m 2 b thus we find in the partial amplitudes the s-channel pole starts to overlap with the t-and u-channel cut when m 2 b ≤ 2. While there is a singularity at the branch point of this cut with the correct sign to "screen" the s-channel pole, this singularity is not strong enough to fully cancel the pole as in 1 + 1 dimensions. The singularity is a log(s − 4 + m 2 b ) in 3 + 1 and (s − 4 + m 2 b ) −1/2 in 2 + 1 (see appendix D for the expicit expressions). We thus expect the peak in figure 4 to remain finite as N max → ∞. This is borne out by some crude extrapolations (not shown).
• Behavior near threshold, m b ∼ 2. As explained in appendix E, when m b − 2 is parametrically small we can analytically constrain the behavior of |g| max as a function of m b . This result is plotted in the figure as the dashed red line segment. Figure 5 shows a closer analysis of this limit. We see that it accurately traces our numerical results, with the agreement improving as m b approaches 2.
• Behavior for m b < 1. In this region the scattered particle is no longer the lightest particle in the theory and on physical grounds we expect the two-particle cut in A(s, t, u) to begin at 2m b rather than at 2m. For small enough m b this is corroborated by our numerics since |g| max ∼ 0 so no pole can be present without modifying our ansatz. It would be interesting to understand in more detail the kink near m b ≈ 0.5.
For m b = 1 we can identify the pole with an exchange of the external particle. Reference [15] (see also [10]) discusses an analytic upper bound on |g| for that case which in our conventions takes the value: |g| 16π √ 1.5 · 10 6 ≈ 61562.4 (18) 13 As for any of the results in this paper, it might very well be possible to derive even stronger bounds by including the constraints from other processes involving more particles.
14 In our ansatz (15) this is easily observed by recalling that t = 0 in two dimensions, so also u = 4 − s. There is also only one partial wave with = 0. which is far weaker than our current bounds. 15

Quartic coupling
Our second set of results concerns the scattering amplitudes M (s, t, u) without any bound state poles, as for example would be the case in π 0 scattering. We will constrain the value of the amplitude at the symmetric but unphysical point s = t = u = 4/3 and therefore define: Historically λ was taken to be a measure of the quartic pion interaction strength. In previous works [17] it was constrained both from above and below, in our conventions: These constraints stem only from the use of axiomatially proven analyticity, crossing and unitarity. Another data point is provided by the explicit "amplitudes" constructed by Auberson and Mennessier, one with λ = 2.62 [16] and one with λ = −1.69 [18], both of which obey analyticity, crossing and unitarity. This provides a lower bound for any upper bound and vice versa. It is particularly remarkable that there exists a fairly narrow interval [2, 62, 2.75] in which the best upper bound must reside.  (15) with g = 0. We impose the unitarity constraint (12) for all ≤ max . Convergence requires larger max for higher values of N max . With this ansatz, the maximal quartic coupling continues to increase significantly with N max even for N max = 20. The black line indicates the value 2.262 achieved in the solution of [16], while the red line indicates the rigorous upper bound 2.75 of [17]. For large enough max and N max our curves must eventually form a plateau between these two lines, however the convergence is so poor that this cannot be inferred from the plot. 3 ), now using the ansatz (15) with g = 0, supplemented with the term (21). With this improved ansatz, the maximal quartic coupling effectively saturates for N max 6. A few values of max are shown to demonstrate that the value of the plateau is independent of this cutoff -the data points for various max are indistinguishable until around N max 12 where the plateau is lost for max = 10 (this is just the usual loss of the plateau when N max becomes too large compared to max ).
Let us first discuss the case of the upper bound. Figure 6 shows the largest possible value on λ using the ansatz (15) (with g = 0). One can see that the convergence with N max is quite slow which suggests the presence of a singularity near or on the boundary of the ρ discs. Indeed, as pointed out in [16,17] the amplitude which achieves the upper bound naturally has a singularity of the form (s − 4) −1/2 corresponding to a bound state sitting precisely at threshold. Physically this is intuitive: the positive sign of the amplitude corresponds to an attractive interaction. 16 The situation in which the interaction is as attractive as possible without introducing new bound states occurs just at the point where a resonance is pulled all the way to the threshold. Mathematically it is natural that to make the amplitude as big as possible at the symmetric point it should be made as big as possible at threshold. Figure 7 shows the bound on λ with the threshold bound state included in the anstaz. This amounts to adding to the ansatz (15) where now α is another parameter to be varied. This singularity does not cause a violation of unitarity because it is canceled by the phase-space volume factor in (12). More precisely, we find that the = 0 partial amplitude near threshold behaves like and therefore − 32 The unitarity constraints for the higher spin partial amplitudes do not lead to further restrictions on α.
Once the threshold bound state (21) is included we find that convergence is now quite rapid as indicated by the plateau in figure 7 already seen at modest values of max and N max . The height of the plateau is 2.6613... and since 2.62 < 2.6613... < 2.75. (24) it falls beautifully below the rigorous bound of [17] but above the solution constructed in [16]. Given the flexibility of our anstaz we expect this value to represent the strictest possible bound that derives from unitarity, crossing and analyticity of a single amplitude.
An interesting feature of the optimal solution is what appears to be a tendency toward saturation of unitarity. In right plot in figure 8 one can see that |S 0 | increasingly saturates unitarity for increasing values of N max . A related fact is that we observe numerically α = −32 √ 6π to great accuracy indicating that unitarity is saturated at threshold. Unitarity saturation is also observed in the higher partial waves.
Let us now consider the lower extremum for which our results are shown in figure 9. As in the previous case (without the threshold singularity) the convergence is quite slow in N max . Unfortunately the addition of a threshold bound-state of the form (21)  since we would need α > 0 to lower the value of λ but according to (23) this is not allowed by unitarity of the spin 0 partial amplitude at threshold. Physically this makes sense -if λ < 0 then this indicates a repulsive force which does not favour the creation of bound states nor moving resonances down to the threshold value. Unfortunately we were not able to identify the relevant singularity in this case and thus were not able to improve the slow convergence.
Notwithstanding these convergence issues, we did already significantly improve the lowest possible value of −1.69 that was explicitly constructed in [18]. As the authors of that paper already noted, the discrepancy between their −1.69 and the lower bound −8.2 of [17] means that either the lower bound is quite far from optimal, or that the behaviour of the amplitude which provides this bound is quite "wild" so as to not be contained within the space of functions they explored. Our results indicate that the latter scenario is the correct one since we do seem to be approaching a value in the ball park of the lower bound in (20).

Exploring scattering lengths
Another set of observables that received interest in days long gone were the scattering lengths a . These are defined as the behavior of the partial waves when s approaches its threshold value 4. We will restrict ourselves to four spacetime dimensions, i.e. d = 3, where it is typically defined as with the limit taken from above in order to make direct contact with experiment. The power of s − 4 in the denominator arises as follows. One assumes that lim s→4 M (s, t) is finite for all t in some neighborhood of zero. Analyticity in t then allows one to write down a Taylor series expansion in t whose radius of convergence remains strictly positive as s → 4. Substituting t = 1 2 (s − 4)(x − 1) and doing the x integral in (10) to project onto the partial waves of spin then gives a finite scattering length for all precisely with the given prefactor (recall that we are considering d = 3). The factor of i is included to make the scattering length real if M (s, t) is real-analytic. In this section we will investigate constraints on the scattering length for amplitudes without bound state poles, so we will be using the ansatz (15) without the pole terms.
Let us begin with the largest possible values of the scattering length. We first recall that, in ordinary quantum mechanics, scattering lengths are known to diverge when a resonance crosses the threshold value s = 4. In the ρ-variables in d = 3 this can be seen by considering scattering amplitudes that locally take the form with the dots denoting subleading terms, which include permutations to make the amplitude crossing symmetric and other terms to make the amplitude unitary for s away from 4. From unitarity near s = 4 we obtain the constraint where we recall that s 0 in our ansatz is equal to 4/3 and we used that [512π 2 (2 +1)] −1 in our conventions. The important observation here is that unitarity bounds µ independently of the value of , whereas the contribution to the spin scattering length is given by (16π) 2 µ (2 + 1) so by sending to zero from above we can get an infinitely large positive scattering length. Notice that < 0 creates a pole on the physical sheet and this is disallowed by our ansatz. 17 The unboundedness from above is borne out by our numerical results. In figure 10 we plot the largest possible values we can obtain for the spin 0, 2 and 4 scattering lengths with our usual ansatz (15), again with g = 0. We observe no convergence to a finite value as we increase N max .
We can also consider the lowest possible values of the scattering lengths. For spin 0 the best known lower bound dates from 1980 and is given by [20] a 0 −1.7 , which slightly improves on a more precise bound obtained five years earlier in [21]: These result were the culmination of a series of works, starting with the observations in [22] which were followed by a series of intermediate improvements in e.g. [10,20,23,24]. 18 Our numerical results are shown in figure 11 and are clearly converging in the neighborhood of the above lower bounds. This shows that the lower bound can more or less be saturated (with an amplitude that falls within our ansatz), which is actually a new result: the best known constructible value was -0.88 [18].
In fact, it may appear that we get dangerously close to the value −1.7 and that further increasing N max may push us over the edge. However for this particular bound the convergence with max is quite slow and the value corresponding to infinite max may in fact increase a little bit. It would be interesting to perform a precision study with larger values of max and N max and to simultaneously re-compute with higher precision the lower bound of −1.7 obtained in [20]. We leave this to the future.
For the higher spin scattering lengths one can use the Froissart-Gribov representation, see e.g. [11], to arrive at the simple lower bound: This is borne out by our numerics but we do not show the results since a plot consisting of nothing but zeroes is not very interesting.

Bonus feature: three spacetime dimensions and QFT in AdS
In our previous work [1] we outlined another method for constraining QFT data, based on putting a QFT in AdS. The main idea is to investigate the boundary correlation functions, which behave exactly like CFT correlation functions (except there is no stress tensor) and are therefore amenable to an ordinary conformal bootstrap analysis. As we explained in [1], the translation between boundary and bulk quantities parallels the standard AdS/CFT dictionary, for example m 2 R 2 = ∆(∆ − d), and furthermore we found precise formulae that dictate how the boundary correlation functions morph into flat-space scattering amplitudes upon sending the AdS curvature to zero. In [1] we numerically tested these equations in 1+1 dimensions and found a quantitative match between the two approaches to the S-matrix bootstrap.
For this paper we set out to repeat this exercise for QFTs in 2+1 dimensions. We focused on the 2+1 dimensional version of the maximal possible coupling that we discussed in section 4.1. This setup was called scenario I in [1]. We discuss the salient points of the methodology before presenting the results.

S-matrix bootstrap approach
For the S-matrix bootstrap, the only difference in the implementation between the 3+1 dimensional analysis of section 4.1 and the present one is that we were no longer able to compute the partial amplitudes (10) analytically. The method explained in appendix D fails because the factor (1 − x 2 ) d− 3 2 in (10) introduces an additional square-root cut in 2 + 1 dimensions (d = 2 in the conventions of this paper). Thus we are forced to evaluate the partial amplitudes by brute force use of Mathematica's NIntegrate. Although slow, this approach is manageable with the use of multiple computing cores. This leads us to the: • First approach: maximal three-point coupling g 2 for any flat-space QFT, obtained by assuming a flat-space scattering amplitude captured by our ansatz (15) and obeying the unitarity condition (12), as a function of m b /m.

QFT in AdS approach
For the QFT in AdS approach we refer to [1] for a detailed exposition of the method, except that presently we consider two-dimensional rather than one-dimensional conformal four-point functions. This implies that there is an extra cross ratio, since z is no longer kinematically equal toz, and conformal blocks are labelled by a pair (∆, ) rather than just the scaling dimension ∆. The combined effect of these modifications is simply that the numerics is computationally much more demanding. 20 Now, in [1] we obtained a precise match in 1+1 dimensions by taking the raw numerical QFT in AdS results and performing a double extrapolation: first to "infinite computational power" and then to infinite ∆ corresponding to the flat-space limit. For our 2+1 dimensional results we unfortunately run into trouble at the first step: our numerical results, obtained for 1 ≤ ∆ ≤ 20 with functionals with up to 136 components, were not amenable to reliable extrapolations. We therefore chose to present the result directly for a QFT in AdS. We chose ∆ = 17 as a representative value. 21 Altogether this gives the: • Second approach: maximal three-point (bulk) coupling g 2 for a QFT in AdS, obtained by assuming boundary correlation functions consistent with unitarity and a spectrum with the natural two-particle gaps, again as a function of m b /m.

Results
The resulting bounds are shown in figure 12. Notice the logarithmic scale. 22 It is clear that the upper bound obtained from QFT in AdS is way larger than the largest value obtained from the S-matrix bootstrap, but the AdS results have not converged yet and one may hope that the numerical upper bound can decrease much further. The good news, however, is the remarkably similar shape of the two curves, both having a somewhat asymmetric peak slightly above m 2 b = 2. In this sense we see a repetition of the results in 1+1 dimensions, namely that we can obtain similar bounds on the residue of a pole in a scattering amplitudes using two drastically different methods.
Physically, it is important to realize that our QFT in AdS approach is completely devoid of any assumptions about the analyticity of the flat-space scattering amplitude. If one agrees that the result in figure 12 provides evidence of the equivalence between the two approaches, then either our S-matrix bounds on the coupling do not require the amount of analyticity 20 The introduction of spin does lead to one new subtlety, namely the magnitude of the two-particle gap for spinning particles. If there is a single scalar particle corresponding to a boundary operator with dimension ∆ then we chose to set the two-particle gap at 2∆ + as in free field theory. Notice that the flat-space limit merely dictates that the gap tends to 2∆ for very large ∆, but there is freedom in choosing the subleading terms. 21 For ∆ = 17 we find that m 2 R 2 = ∆(∆ − 2) = 255 so the reduced compton wavelength of the particle is about 16 times the AdS radius of curvature in our setup -in this sense space is already quite flat. 22 On a regular scale the shape of the peak is very similar to the one shown in figure 4. that we have imposed or the analyticity (at least of the extremal scattering amplitudes) is a property that we may hope to derive from the QFT in AdS construction. Either option would be very interesting and should be investigated further. 23 Although ∆ = 17 was the largest value for which we had a full set of results, let us briefly discuss the result for 0 < ∆ ≤ 20. In line with the results in [1], the absolute value of the numerical bounds decreases quickly upon decreasing ∆. For ∆ 4 the curve always has a peak hovering around m 2 b = 2, which broadens a bit upon decreasing ∆. For 0 < ∆ 4 the peak moves more or less linearly towards m 2 b = 4 as ∆ → 0. In the future it would be interesting to invest more computational resources and explore in more detail both this behavior and the general convergence of the bounds.

Discussion
Here we continued our exploration of the space of S-matrices of gapped quantum field theories initiated in [1,2]. We present a fresh approach to an old question of constraining S-matrix elements based on unitarity, crossing and analyticity. The former two properties are firmly established properties of the S-matrix whose meaning requires no clarification. By analyticity we mean the rather simplistic (but perhaps most natural) assumption that M (s, t, u) is an analytic function of each of its variables with no singularities in their respective cut planes. We make no assumption about the properties of the S-matrix outside of this union of cut planes -i.e. off the physical sheet.
Of course there are many open questions in S-matrix theory pertaining to analyticity. Are all singularities in the complex Mandelstam variables s, t, u associated to Landau diagrams (as expected based on perturbation theory) or should we be open to more exotic possibilities especially in strongly coupled theories? What is the most general possible large energy behaviour of scattering amplitudes? Finally, if we bravely cross the gates and delve into the various Riemann sheets of non-perturbative scattering amplitudes by crossing its various cuts in the physical sheet, what kind of scary Chimeras await us down there?
We tried to be optimistic -by assuming the minimal expected singularities in the physical sheet -and cautious at the same time -by assuming as little as possible about the uncontrollable world of the other unphysical sheets or the large energy behaviour of scattering amplitudes. In short we mapped the physical sheet into a few unit disks and assumed little about the behaviour of amplitudes on the boundary of those disks which is where both the large energy behaviour as well as the various physical thresholds lie. Inside these disks we assumed that the only singularities were poles associated to stable bound states.
In the future, it would be interesting to develop new numerical investigations relying on more rigorous analyticity assumptions. Perhaps our results are not too sensitive to this distinction, or perhaps we will encounter exotic S-matrices which make use of the allowed non-analyticity to allow for a wider range of values. Both would be very interesting! To this end, it is worth noting that in the case of the quartic pion coupling and the lower bound on the spin zero scattering length we can say with confidence that we are in the former scenario -our results approach the bounds obtained in [17,18] and in [20] which are based on rigorously proven analyticity properties. More evidence for the first scenario is the at least qualitative match between our maximal coupling and the upper bound on the same observable for a QFT in AdS, since the latter computation relied on no analyticity properties whatsoever. Finally we can point to the consistency of our approach with a Mandelstam representation expansion discussed in appendix C.
As for the behaviour at the boundary of the disks the idea here is that we can be agnostic about it and let regular Taylor expansions in the bulk converge towards whatever they want to. Of course, without inputing the correct singularities at the boundary of the disk, the numerics should still work but their convergence will suffer considerably. We encountered two examples of this already in the main text. The first is the quartic coupling numerics whose convergence increased substantially once we allowed for a bound state singularity at threshold. Another example is probably the four dimensional bound state coupling numerics when the bound-state mass is less than √ 2 times the mass of the lowest particle. The numerics are converging much slower for that range as clearly seen in the left curves in figure  4. We suspect in this case it is rather related to a non-trivial large energy behaviour of the S-matrix which the ansatz has a hard time reproducing. 24 It would be interesting to investigate this further.
It is also at the boundary of these disks where we read physical amplitudes with any s > 4m 2 and negative t. Multi-particle production will show up as further cuts at larger s such as 9m 2 , 16m 2 , etc and infinitely many others like (m + m ) 2 , etc if there are other stable particles. We saw no signs of these singularities in our numerics. As we for example show in figure 8, our optimal S-matrices do not seem to open multi-particle production cuts in any significant way. A priori this sounds very strange. How could we have no particle production of four particles from two particles if -by crossing one particle to the past -that amplitude is related to a 3 → 3 process which obviously must exist? 25 Indeed, it is known [26] that particle production is mandatory. It can not be strictly zero or it would lead to important contradictions. Unfortunately, the same work [26] -or any other work as far as we knowdoes not put a lower bound on how much particle production one must have and as such we could not reach a sharp contradiction with the numerics which by definition can never rule out an arbitrarily low particle production. 26 Nonetheless, absence of particle production is unphysical in spacetime dimension greater than 2. We would like to describe more realistic theories where particle production naturally arises. One way of forcing such particle production in a natural way is to study multiple S-matrix elements where we consider a system of scattering elements involving not only the lightest particle but also the next-to-lightest etc. We are currently working on this and finding some very encouraging preliminary results in two dimensions where the bounds are often significantly improved and the corresponding S-matrices do exhibit particle production and thus must correspond to genuinely non-integrable theories in contrast to our previous work [2].
The analyticity properties of scattering amplitudes of several particles of different mass are more intricate than what we considered here. The optimistic scenario is that all singularities on the physical sheet follow from Landau diagrams describing propagation of on-shell particles. This Landau analyticity is far from being rigorously established but it is a reasonable physical conjecture to start from. Even with this assumption, we will have to deal with anomalous thresholds (singularities that arise from Landau diagrams that are not on a line). A simple example is the scattering amplitude of particles of mass greater than √ 2 times the mass lightest particle. We plan to analyse this issue in the future, starting in 1+1 dimensions.
in part by the Government of Canada through NSERC and by the Province of Ontario through MRI. This research received funding from the grant CERN/FIS-NUC/0045/2015. This work was additionally supported by a grant from the Simons Foundation (JP: #488649, BvR: #488659, PV: #488661) JP is supported by the National Centre of Competence in Research SwissMAP funded by the Swiss National Science Foundation.
A x(s) vs ρ s , ρ t in 1 + 1 dimensions Consider the map which maps the full s-plane minus the cuts s > 4 and s < 0 into the unit disc |x(s)| ≤ 1 and the map s which maps the full s-plane minus a single cut s > 4 into the unit disc |ρ s | ≤ 1. An analytic function in the s-plane minus the cuts s > 4 and s < 0 -such as the S-matrix once we subtract out its known poles -can be written as Now, we have which admits a convergent expansion in powers of ρ s and ρ t provided they are both inside the unit list (and hence so is their product in the denominator). Hence the function f (s) can also be cast as As such, our 1 + 1 numerics had to work.

C Mandelstam Representation
The double dispersion representation proposed by Mandelstam [28] implies that the amplitude can be written as follows where B(s, t) = ds dt C(s , t ) (s − s)(t − t) .
If there are no stable particles below threshold, the double discontinuity C(s, t) has support inside the region s > 4m 2 and t > 4m 2 . In practice, this form of the double dispersion relation is not valid and one needs to include subtractions. A simple trick to derive the form of the dispersion relation with n subtractions is to use the identity in equation (41) for both factors in the denominator. This leads to In general the integrals (44) do not converge. The subtracted dispersion relation is (43) considering c k (t) and c k,l as independent functions from the double discontinuity C(s, t).
Stable particles correspond to delta-function pieces in the single discontinuities c k (s). 27 Besides these delta-functions, the support of c k (s) is s ≥ 4m 2 . Therefore, the analytic properties of equation (43) imply that with a convergent double ρ series in the product of two unit disks. This is a more restricted form of formula (15) where we set to zero all coefficients α abc with a > 0, b > 0 and c > 0.
In order to test the validity of Mandelstam representation, we reconsidered the problem discussed in section 4.2 using the more restricted ansatz In figure 14, we show the maximal value of the quartic coupling λ obtained with this ansatz. The maximal value λ ≈ 2.6613... is obtained for N max 6. This result suggests that in the limit of large N max both ansatze cover the same space of functions. Figure 14: Comparison of upper bound on pion coupling using ansatz (15) with g = 0 and the threshold singularity (21)

D Partial Wave Integrals D.1 Pole contributions
Here we will consider the contribution to partial waves coming from poles of the scattering amplitude. Consider It is easy to compute the partial wave decomposition of this expression. For d = 3 we get with x b = x(s, t = m 2 b ) and Q (z) the Legendre function of the second kind with branch cut along z ∈ (−1, 1). For d = 2 we instead get Now consider the contribution to the amplitude from a threshold bound state. The pole part is If we focus on the case d = 3, we must compute integrals of the form: with t(x) = − 1 2 (s − 4)(1 − x). Introducing the generating function for the Legendre polynomials +∞ n=0 z n P n (x) = 1 it is not difficult to obtain Adding up contributions from s, t, u the partial amplitudes are Here we will show how to obtain the contribution to the partial amplitudes from terms of the form ρ a s ρ b t ρ c u analytically in d = 3. While the calculation is somewhat tedious, the underlying concept is simple: the integral that we want to do has only one cut (of square-root type) in the integrand and thus with a simple trigonometric change of variables the integrand can be converted to a rational function and computed by partial fractions (or some more clever method).
The non-trivial integrals to perform take the form with, as in (2) with m = 1, In applications we typically set s 0 = 4/3. We next introduce our first inspired change of variables from x to φ which is given by In these variables we get: where we also introduced We should now do the usual change of variables, This gives We have and the integral runs from y i to y f with The trick now is to rewrite the integration region using the discontinuity of a logarithm, where (y 1 , y 2 ) is a clockwise contour wrapping the line segment from y 1 to y 2 . In our case f (y) is a rational function, therefore we can pull the contour to infinity so that it picks up the poles of f (y) to obtain exact expressions.

D.3 Large energy
Let us consider the large energy limit s → ∞ of our ansatz. Since unitarity is imposed for each spin , we are interested in the limit s → ∞ with fixed scattering angle θ. In this limit, we find The contribution from the pole terms in our ansatz are real and of order 1/s in this limit and therefore can be neglected. The leading term in (65) only contributes to the spin 0 partial wave. The large s expansion of S 0 (s) is given by Unitarity implies that (for d > 4 the inequality must be saturated) If d > 2 then unitarity also implies that For d = 2, the correct condition is For > 0 (even) we find where 28 Therefore, unitarity implies (for d > 5 the inequality must be saturated) where we used that I < I 2 = √ 4−s 0 10π for > 2. Where applicable, we have verified the above constraints a posteriori for our numerical solutions and found them satisfied to very good numerical accuracy.
As a final comment, we remark that the unitarity constraints dictate that lim s→∞ S (s) = 1 for any amplitude within our ansatz with finite N max . 29 This property is likely to be too restrictive, and it is therefore worthwhile to try to improve our ansatz with more singular terms compatible with unitarity and analyticity. As a first attempt we added an extra term of the form (ρ s + 1)(ρ t + 1) −1 plus s, t, u permutations, which allows lim s→∞ (S (s) − 1) to be non-zero -this modification however did not significantly change any of the results displayed above. In the future we plan to add other more singular terms and investigate their effect in more detail. 30 Finally, the restricted behavior at large s might also be a source of slow convergence when N max → ∞ we have observed in some cases. This idea is also corroborated by the two dimensional analysis in appendix G.

D.4 Large spin
The partial waves can also be written in terms of an hypergeometric function, x = 0. Generically, this will come from the poles associated with stable particles. More precisely, where In fact, the pole −g 2 which decays exponentially with l. Notice that this gives a purely imaginary contribution to S (s) (see equation (10)), which by itself would violate unitarity. However, unitarity can be restored with a small real contribution of the order of the square of (84). At large l, this requires that we match the exponential behaviour In other words, unitarity can be restored with another particle or particles of total invariant mass squared m 2 2 ≥ 4m 2 1 . This is what happens in perturbation theory where the box diagram restores unitarity of the tree-level exchanges.
Let us now study the contribution from the polynomial terms ρ a s ρ b t ρ c u in our ansatz. The discontinuity of M for x > 1 comes from where and we only kept the leading behaviour of the discontinuity near its lower end x (s). Similarly, we can approximate and find At threshold the poles become constants and are irrelevant. This is not so for threshold poles which are discussed separately below. Define w := √ s − 4. Then for s → 4 + above the cut we have Recall that in our conventions the partial waves take the form: The leading contribution for the spin partial wave corresponds to the k = term in the above, leading to Absence of particle production would imply reality of δ (s), and hence a measure of the inelasticity of the amplitude at the threshold is We see that for positive spin we generically get a divergent result in the threshold limit. This means that our ansatz does not automatically give an amplitude which becomes purely elastic as we approach threshold, unlike what we would expect on physical grounds. In order for purely elastic scattering to hold, we would have had to impose order linear constraints on the coefficients of the threshold expansion of the spin partial wave. We did not impose these in our numerical computations. However, experimentally we do find that as the number of parameters in our ansatz is increased, the coefficients in the threshold expansion seem to decrease.

E Non-Relativistic Limit
Consider a scalar φ of mass m interacting with itself via the exchange of a second heavy scalar Φ with mass m b = 2m − with small . We can think of Φ as a loosely bound state of two φ particles with binding energy . The two body amplitude for φ + φ scattering contains a pole at s = m 2 b due to virtual production of a Φ which is just below the the two-particle threshold at s = (2m) 2 . The residue of this pole g 2 is the square of the φφΦ coupling. Now consider low energy φ + φ scattering and write s = (2m + E) 2 where E is the centre of mass energy after subtraction of the rest mass. The s-channel pole of the amplitude is given by 31 where we have assumed small E and . The l = 0 phase shift inherits this pole through the relation Plugging (102) into (103) and zooming in on the pole of the phase shift at s = (2m − ) 2 we have We write the pole of the phase shift as g 2 N R /(E/ + 1) where g 2 N R is the residue in units of the binding energy . We then have We will show below that there is a bound on the non-relativistic coupling g 2 N R ≤ 2 2 . Note that this correctly predicts the behaviour in 1 + 1 dimensions [2]. Moreover, this limit has been studied extensively in 3 + 1 dimensions (d = 3) [29,30]. These authors find (adding a factor of 2 to their results to account for identical particles) and thus we find perfect agreement with (105).
Let us now derive the bound on g 2 N R quoted above. Recall that we are considering a very weakly bound state with binding energy . We wish to obtain the behaviour of g 2 max ( /m) for small /m. Thus we concentrate on "slow" physics at energies E ∼ (recall E is the centre of mass energy after removal of the rest mass). Formally, in the phase shift we consider s →s 2 and consider finites as → 0. Any singularites of the phase shift that are a finite 31 The factor m 5−d is to make the coupling g 2 dimensionless.
distance (in s) from the two-particle threshold -e.g. the left cut and inelastic thresholdswill be infinitely far away ins and thus only contribute through positive powers of . We can thus neglect these singularities to obtain the leading behaviour of g 2 max ( /m) and consider a non-relativistic phase shift S NR (Ē) with only a right-hand cut starting atĒ = 0 and a single bound-state pole atĒ = −1, whereĒ = E/ . Since this phase shift is bounded by unitarity along the cut and cannot grow faster than a constant at infinity then the residue of the pole can easily be bounded by maximum modulus theorem. Perhaps the cleanest way to derive the precise value of the bound is to consider the change of coordinates which maps the E-plane minus the positive real axis to the unit disk and maps the bound state pole to the origin Now note that the function f (x) = x S NR (x) is analytic throughout the unit disk and obeys |f | ≤ 1 on the boundary due to unitarity. Thus maximum modulus theorem implies 1 ≥ f (0) = g 2 N R /4 which is the desired bound.
The eigenvalues of this matrix are precisely As befits a Hermitian matrix, they are always real since U ≤ 1 by construction. It is now clear that and the unitarity constraints are therefore precisely those of a semidefinite program.
We need to choose a grid of values of s and a finite set of spins for which to test the unitarity constraints. We found it sufficient to take approximately 200 values of s, interspersed uniformly along the upper half of the unit circle in the ρ s variable defined in the main text. We observed no significant change in the results by taking a more refined s grid, or by distributing the points differently along the unit circle. The maximal value of the spin max is indicated in the various plots. Notice that max needs to be sufficiently big since otherwise the extremal value completely destabilizes -see for example the data points in figure 6 with max = 10 for large N max . In practice we observed convergence by taking max at least as large as N max , and for the scattering length computations we needed at least N max + 4. Increasing max beyond these values did not affect our results.
In our numerical computations we did find it necessary to retain very high precision, generally at least 1000 binary digits. This appears to stem from the approximate redundancy that remains even after imposing the polynomial constraint B. To illustrate this we can for example compute a derivative like In a typical solution we find that this derivative is rather modest in magnitude, of order 10 2 or so, whereas the individual coefficients can be very large, of order 10 24 in some solutions. These kind of cancellations require high precision.
We have performed all the numerical computations in section 4 with sdpb [31]. Details of the computations like parameter settings are available from the authors upon request.

G Slow convergence on a simple 2D example
In this appendix we revisit once more the two dimensional problem considered in section 2 but this time done in the language of the M amplitude rather than S. In two dimensions the two are simply related by S(s, t) − 1 = 1 2 √ st × M (s, t) , s + t = 4m 2 . Great convergence Poor convergence Great convergence Figure 15: Plot of |(S num − S analytic )/S analytic |, that is of the relative mismatch in the numerical solution in this two dimensional example where the analytic solution is available. In all these plots we use Λ = 20 and check unitarity in a small grid of 40 points. With these parameters, mathematica's built-in FindMaximum suffices and produces an outcome in about two or three seconds. We see on the left that for m b > √ 2m the agreement is spectacular with the most naive ansatz (120) while in the middle we see that the agreement is much worse (a few percent off) with the same ansatz when m b < √ 2m. On the right we see that this is neatly fixed -leading again to a perfect convergence -by simply adopting an improved ansatz as in (121). and unitarity then reads This discussion will provide us with a simple example of numerics which work yet converge very slowly until we slightly improve our ansatz and thus completely solve this convergence issue.
To be concrete we consider here the case where there is a single bound-state with mass m b whose coupling we maximize. The S-matrix with the largest coupling and such bound-state is given by [2] S max g = sign(m At high energies the S-matrix approaches +1 for m b > √ 2m and −1 for m b < √ 2m and this leads to a very different behavior when translated to the amplitude M . In particular, for a light bound state m b < √ 2m we see that the amplitude M in (117) must diverge at high energies so that the right hand side approaches −2. This is hard for an ansatz a la (3) to achieve, that is it would require that the sum in to develop a divergence as s = 4m 2 − t → ∞ which corresponds to ρ s , ρ t → −1. Such non-analytic behavior at the boundary of the unit disc can be achieved but a numerically sufficiently accurate approximation requires very large Λ.
In this case there is however a very obvious improvement which is to simply allow for a divergence at large energies which is after all allowed by unitarity and write down instead an ansatz of the form c ab ρ a s ρ b t + β √ ρ s + 1 √ ρ t + 1 +β (ρ s + 1)(ρ t + 1) (121) This immediately allows for a more general high energy behavior and thus an extreme improvement in convergence as illustrated in figure 15.
The moral of this story seems to be that we better allow for flexible ansatze which can easily capture various analytic properties of scattering amplitudes if we want to achieve optimal convergence. In this simple two dimensional example, allowing for an ansatz with a more flexible high energy behavior led to a drastic improvement in the numerics.