An Analytical Toolkit for the S-matrix Bootstrap

We revisit analytical methods for constraining the nonperturbative $S$-matrix of unitary, relativistic, gapped theories in $d \geq 3$ spacetime dimensions. We assume extended analyticity of the two-to-two scattering amplitude and use it together with elastic unitarity to develop two natural expansions of the amplitude. One is the threshold (non-relativistic) expansion and the other is the large spin expansion. The two are related by the Froissart-Gribov inversion formula. When combined with crossing and a local bound on the discontinuity of the amplitude, this allows us to constrain scattering at finite energy and spin in terms of the low-energy parameters measured in the experiment. Finally, we discuss the modern numerical approach to the $S$-matrix bootstrap and how it can be improved based on the results of our analysis.


Introduction
The idea of bootstrapping the S-matrix of a unitary, relativistic, gapped theory in d ≥ 3 was actively pursued in the 60's. While many interesting results have been derived [1][2][3][4], no nonperturbative physical S-matrices have been computed. The main reason being that no solid, nonperturbative calculation scheme was ever put forward without relying on some unreliable approximations. In addition to that, analytic properties of the multi-point amplitudes were never fully understood. Moreover, even at the level of the two-to-two scattering amplitude, the region of analyticity that is usually assumed in the bootstrap analysis has not been rigorously established.
While the problem of analytic properties of multi-point scattering amplitudes is still widely open, the question of finding a good calculation scheme has recently acquired an interesting twist with the development of the conformal bootstrap [5,6]. In this context, by exploring bounds on the OPE data in the space of solutions to the CFT bootstrap equations, it was found that sometimes the physical theories of interest saturate the bootstrap bounds [7][8][9] and are in this sense solvable. 1 Remarkably, a similar phenomenon was observed for 2d S-matrices [10][11][12], where some previously known two-dimensional scattering amplitudes of physical theories were found to saturate the bootstrap bounds. We can loosely call such theories bootstrap-solvable. A fundamental, open question in the S-matrix theory program therefore is: Are there bootstrap-solvable S-matrices in d ≥ 3?
A distinguishing feature characteristic to scattering in d ≥ 3 dimensions is a remarkable connection between scattering and particle production. Here, by scattering we mean a non-trivial 2 → 2 amplitude and by particle production, we mean a non-zero 2 → n amplitude with n > 2. While scattering without particle production is a commonplace in d = 2 [13], in d ≥ 3 it is widely believed that scattering implies production [14]. 2 The underlying reason is an elegant interplay between analyticity, elastic unitarity and crossing symmetry.
The simplest non-trivial S-matrix element is a 2 → 2 connected scattering amplitude T (s, t) as a function of Mandelstam invariants s and t. In a gapped theory, this amplitude is subject to an exact non-linear equation called elastic unitarity. This equation originates from the fact that when we scatter the lightest particles in the theory at energies below the first multi-particle threshold (s 0 ), only two particles can be produced in the final state. 3 As a result, unitarity of the S-matrix becomes a non-linear equation satisfied by T (s, t) for 4m 2 < s < s 0 . It is the purpose of the present paper to revisit the implications of elastic unitarity, when combined with crossing and analyticity, on the structure of the nonperturbative amplitude T (s, t).
One motivation for our analysis is recent numerical investigations of higher-dimensional scattering in [15,16]. In these works elastic unitarity was not imposed and it was observed that various bootstrap bounds tend to be saturated by purely elastic functions. It is therefore an interesting, open question how to efficiently implement elastic unitarity and particle production in the current S-matrix bootstrap program. An obvious way to tackle the problem is to include higher-point amplitudes in the bootstrap analysis explicitly. However, due to unknown and complicated analytic properties of higher-dimensional amplitudes it is not clear if it is feasible in d ≥ 3. Another possible way to make progress, which we will follow in the present paper, is to focus on how to implement structures of the amplitude that are dictated by elastic unitarity into the current nu-merical approach of [15,16]. In this way we hope to be able to zoom in closer on the physical higher-dimensional S-matrices and, if we are lucky, maybe eventually solve them. Here, we discuss various possibilities of doing that and will report the numerical results in [17]. We start by laying out our assumptions that serve as the basis for the further analytic study.

Assumptions
We assume that in the far past and in the far future states of the system are described by a set of free particles. The Hilbert space therefore is taken to be the Fock space of free particles. 4 For simplicity we assume that the spectrum contains a single scalar particle of mass m together with its multi-particle states. Under these conditions it has been rigorously proven that the amplitude satisfies: Our extra assumptions for the connected two-to-two scattering amplitude T (s, t), which have not been rigorously established, are the following: 3. Extended analyticity: T (s, t) is an analytic function for complex s and t in some region D, except for potential poles in 0 < s < 4m 2 and a cut starting at s = 4m 2 , as well as images of these singularities under crossing. We will often assume that the region of analyticity D extends to the full complex s-plane (t-plane) for some finite region in t-plane (s-plane), but many of our arguments can be adopted to the situation when D is bounded in both variables simultaneously. 5 As usual in this paper T (s, t) stands for the analytic continuation of the amplitude from the physical regime to the principle sheet -without going through the multi-particle cuts. 4. Polynomial Boundedness: for fixed t |T (s, t)| < |s| J 0 (t) , |s| → ∞ , (s, t) ∈ D . (1. 3) The formula above assumes that D includes s = ∞ for fixed t. More generally, we will assume that the amplitude is polynomially bounded on the principal sheet (away from the bound states poles and multi-particle thresholds).
We also have two extra technical assumptions:

5.
Continuity: partial waves f J (s) are real analytic functions in the elastic region 4m 2 < s < s 0 . 6 6. Z 2 symmetry and no bound states: for simplicity we assume no bound states in the spectrum. We also assume that we have a single stable particle of mass m which is odd under Z 2 symmetry. Therefore only an even number of particles can be produced in the scattering of two particles.

Plan of the Paper
In section 2 we review basics of the two-to-two scattering amplitude. We consider scattering of identical scalar particles and briefly review known analyticity results. We discuss unitarity and elastic unitarity in terms of T (s, t). We introduce partial wave expansion and the Froissart-Gribov formula. We briefly discuss continuity properties of T (s, t).
In section 3 we consider analytic continuation of elastic unitarity. We first review the derivation of the Mandelstam equation for the double spectral density ρ(s, t), or equivalently double discontinuity of T (s, t), based on analytic continuation of elastic unitarity in one of the Mandelstam variables (s or t). We then discuss its relation to analytic continuation of partial waves in spin J via the Froissart-Gribov formula. We also consider analytic continuation of elastic unitarity in energy s. Finally, we exhibit the structure of the Landau-Karplus curves along which the double spectral density develops a nontrivial support in the elastic region.
We then analyze various implications of elastic unitarity: • In section 4 we discuss positivity properties of double spectral density ρ(s, t). We review the argument that scattering implies production in d ≥ 3. This result follows from the combination of crossing symmetry and positivity of ρ(s, t).
• In section 5 we introduce the notion of the threshold expansion for partial waves, as well as for the first and second discontinuities of the scattering amplitude. The threshold expansion of partial waves and the discontinuity of the scattering amplitude is a consequence of elastic unitarity and it describes low-energy or non-relativistic scattering. The Mandelstam equation then maps it to the expansion of double spectral density close to the boundary of its nontrivial support, the so-called Kaprlus-Landau curve.
• In section 6 we map the threshold expansion in the t-channel to the large J expansion of partial wave coefficients in the s-channel. It comes from "inversion" of the threshold expansion via the Froissart-Gribov formula. Turning to the 1 J corrections we find that remarkably the computations can be sometimes done exactly in 1 J .
• In section 7 we turn the large J results into finite J predictions plus an error estimate. This error estimate is deducted from a local bound on the discontinuity of the amplitude in the region s, t > 4m 2 . No such rigorous bound is known. We discuss natural error estimates and perform the finite spin, finite energy computations in a simple toy model.
In section 8 we consider the modern numerical approach to the S-matrix bootstrap following [15]. We discuss why and what should be improved in the existing approach. We suggest several ways in which this approach can be improved. In particular, we discuss various ways to implement elastic unitarity numerically. In section 9 we briefly comment on the relation between the present analysis and similar ideas in the conformal bootstrap. Finally, in section 10 we conclude and present some future directions. Several appendices contain technical details that should be helpful in understanding the details of our arguments and calculations.

Amplitude Basics
In this section we briefly review the basic kinematics of the two-to-two scattering of identical, scalar particles to set the conventions for the further analysis. As usual we write the S-matrix aŝ S ≡1 1 + iT , (2.1) whereT is zero in the theory of a free massive scalar. We are interested in the matrix elements that describe two-to-two scattering S 2,2 (p 3 , p 4 |p 1 , p 2 ) ≡ p 3 , p 4 |Ŝ|p 2 , p 1 , where the initial and final states are characterized by the on-shell momenta As in (2.1) we can separate the contribution of the disconnected and connected parts of the S-matrix S 2,2 (p 3 , p 4 |p 1 , p 2 ) = S 1,1 (p 3 |p 1 )S 1,1 (p 4 |p 2 ) + S 1,1 (p 4 |p 1 )S 1,1 (p 4 |p 2 ) + S c 2,2 (p 3 , p 4 |p 1 , p 2 ). (2.4) The disconnected part is given by an overlap of the one particle states which is uniquely fixed by Lorentz symmetry

Analyticity
In the discussion above scattering matrix elements were defined for physical momenta that correspond to actual scattering. The starting point of the S-matrix considerations is the statement that the physical matrix element T (s, t) is a boundary value of an analytic function of the relativistic invariants regarded as complex variables where we assumed s and t to be real and −(s − 4m 2 ) < t < 0, 4m 2 < s. This corresponds to scattering in the s-channel 1, 2 → 3, 4. Using the basic principles of QFT correlation functions, reduction formulas that relate them to the S-matrix elements, and techniques of analytic completion one can establish various analytic properties of the scattering amplitudes as functions of complex s and t, see e.g. [4,20,21] for a pedagogical exposition.
The first result of this type is subtracted dispersion relations in s for fixed −t 0 < t ≤ 0 in the physical s-channel region [21,22]. In this case one has to understand analytic properties of T (s, t) as a function of complex s. For the case of π 0 π 0 → π 0 π 0 scattering, for which our treatment applies directly, one can show that T (s, t) is an analytic function of s with two cuts: s > 4m 2 and u > 4m 2 , with t 0 = 28m 2 . See Table 1 of [21] for the processes for which a fixed-t dispersion relation has been proven and the corresponding values of t 0 . 7 Another well-known property, originally due to Lehmann [24], concerns analytic properties of T (s, t) as a function of t for fixed physical s > 4m 2 . Lehmann showed that T (s, cos θ) is analytic inside an ellipse, the so-called Lehmann ellipse, in the cos θ complex plane with foci at cos θ = ±1 and semi-major axis cos θ sL > 1 which depends on the details of the theory, energy and masses of particles, as well as the scattering process. Lehmann also showed that the absorptive part or discontinuity of the amplitude Disc s T (s, cos θ) is analytic in a larger ellipse, the so-called large Lehmann ellipse, with a semi-major axis cos θ LL = 2 cos 2 θ sL − 1.
The third class of results concerns analyticity of T (s, t) when both s and t are complex. Bros, Epstein and Glaser [25] showed that any point (s, cos θ) in the physical region is surrounded by an analyticity neighborhood whose precise form is not known in general, see e.g. [4] for details. An explicit domain of simultaneous analyticity in both variables was derived by Lehmann [26] for elastic processes that obey a fixed-t dispersion relation by continuing the Lehmann ellipse to complex s. For cases with single variable dispersion relations in all three channels (as in ππ → ππ scattering), Mandelstam [27] derived domains of the form |s t| < b for any complex s and t outside the single variable dispersion relation cuts. For pion scattering the largest domain occurs for b = 256m 4 π . These domains have the drawback that when s → ∞ the Lehmann domain shrinks to the line −t 0 < t < 0 and the Mandelstam domain shrinks to the point t = 0.
Based on the results above, the analyticity domain was further enlarged using unitarity by A. Martin [28]. The final result is that the amplitude is analytic within |t| < R and s cut-plane. For scattering of identical particles that we consider in the present paper R = 4m 2 . From this result the standard bounds on high energy behavior of the amplitudes follow. It also follows that for |t| < 4m 2 the scattering amplitude admits fixed t dispersion relations with at most two subtractions.
Finally, let us introduce the notion of maximal analyticity (or Mandelstam analyticity) which states that the scattering amplitude is analytic in the (s, t) complex planes with only singularities on the principal sheet being unitarity cuts s, t, u > 4m 2 and bound state poles for 0 < s, t, u < 4m 2 . This property is consistent, at least for the scattering of lightest particles in the theory, with expectations from unitarity and analysis of perturbation theory. However it is important to keep in mind that it has not been proven. The original attempts to prove maximal analyticity in perturbation theory [29,30] were later found to have a loophole [31] which, to the best of our knowledge, has never been closed even for the scattering of lightest particles. In this paper we freely assume analyticity beyond what has been rigorously proven but we do not assume maximal analyticity.
There are two other properties that we will use. One is crossing 12) which states that in particular that scattering in different channels is described by different boundary values of a single analytic function. For the two-to-two scattering it has been proven in [23].
Beyond the two-to-two scattering only a partial progress has been achieved [32]. Another property is real analyticity This was established within axiomatic quantum field theory in [33]. For the S-matrix argument see e.g. [1].
Finally, let us also mention in a related context the result by A. Martin [34] who showed that maximal analyticity in the form of the Mandelstam representation together with knowledge of the double spectral density in the elastic region fix the scattering amplitude completely. Similarly, knowledge of the amplitude at fixed energy in the elastic region as a function of the scattering angle is believed to fix it almost completely, see e.g. [35].

Unitarity and Elastic Unitarity
In terms of the T -matrix, the unitarity of the S-matrixŜ ·Ŝ † =1 1 reads where in the second line we have inserted a complete basis of asymptotic states and the Lorentz invariant measure is (2.15) Due to the momentum conservation, a 2n-particle intermediate state can only contribute if s > (2nm) 2 . Otherwise, there is no enough energy to create an on-shell state with 2n particles. In particular, for 4m 2 < s < 16m 2 , only two particle states are possible. Hence, in that regime we can replace the sum in (2.14) by the first term and the equation closes on the 2 → 2 transitions only. This is the elastic unitarity regime. Explicitly, we have 2T s (s, t) = 1 2 where we have introduced the notations In writing the above we used real analyticity of the amplitude (2.13). Finally, the overall factor of one half is a symmetry factor for identical bosons.
We will now reduce the integral to an integration over the two scattering angles by performing all the kinematical integrations explicitly. For that aim, we first go to the center of mass frame where q = − q ≡ p n, where n is a unit d − 1 vector and p = √ s − 4m 2 /2. In these variables the elastic unitarity constraint (2.16) becomes s−4m 2 is the Jacobian coming from the energy conservation delta-function. The integrand only depends on the two scatttering angles in terms of which we can write the measure as is the cosine of the external scattering angle. We find that (see appendix A for more details) Using that E k = √ s/2 , we can write (2.18) covariantly as where t(x) ≡ −(s − 4m 2 )(1 − x)/2, not to be confused with the external momentum transfer t, which is held fixed. The step/delta function in the phase space integration kernel P d (z, z , z ) has a simple geometrical origin, see figure 1.
For s > 16m 2 and general t the unitarity constraint involves scattering elements with more than two particles. To get a constraint on the two particle amplitude, we note thatT ·T † on the right hand side of (2.14) is a positive semi-definite matrix. Hence, for any state Ψ we have that 23) Figure 1: In the elastic strip 4m 2 < s < 16m 2 the discontinuity of the amplitude comes from two intermediate particle exchange only. This result in the exact elastic unitarity equation (2.22). The corresponding phase space integration kernel P d (cos θ, cos θ , cos θ ) in (2.21) is proportional to a step/delta function which has a simple geometrical origin. In the center of mass frame we have three (d − 1)-dimensional vectors, p 1 = − p 2 , p 3 = − p 4 , and q = − q . The geometrical angles between these three vectors, {θ, θ , θ } are therefore restricted to the range θ 1 + θ 2 ≥ θ 3 , where θ 1,2,3 are any permutation of {θ, θ , θ }.
and hence, if we drop all the contributions with more than two particles in (2.14) we get an inequality for the 2 → 2 scattering matrix For example, if we pick a wave function that consists of two particles with a specific momenta then we have the the amplitude in the forward limit where p 3 = p 1 and p 4 = p 2 . For this choice of wave function, the unitarity constraint (2.24) becomes 8 where we used that P d (1, z , z ) ∝ δ(z − z ), see (B.7) for the precise formula. We also used that

Partial Wave Expansion
Unitarity of the S-matrix implies the non-linear integral relations that the 2 → 2 T -matrix has to satisfy, (2.22) and (2.24). To simplify these complicated constraints we choose a wave function 8 The unitarity relation in the forward limit is nothing but the optical theorem. In d dimensions it takes the form Ψ that diagonalizes the T -matrix and therefore also the integral kernel in (2.22), (2.24). This can be done using the Lorentz symmetry of the problem. Namely, we decompose the amplitude T (s, t) in a complete basis of intermediate states which transform in irreducible representations of the SO(1, d − 1) symmetry. These representations are characterised by their energy and the little group SO(d − 1) angular momentum in the center of mass frame, E and J. For two particle states the SO(1, d − 1) quantum numbers are enough to characterize the states and we have J, m are the d-dimensional spherical harmonics, and the energies dependant prefactor will not be relevant for us. 9 We can now insert a complete basis to these states to decompose the S-matrix element p 3 , p 4 |T |p 1 , p 2 in all possible spins. Since the operatorT is both, translation and SO(1, d − 1) invariant, due to the Wigner-Eckart theorem we have that where the convention-dependent proportionality factor is independent of the energy and and the angular momentum m. These functions are the so-called partial wave coefficients, in terms of which the amplitude takes the form where the sum runs over all (even) spins and n (d) J are convention-dependent normalization factors. Here, P (d) J (cos θ) are the partial waves. They represents the angular dependence of the amplitude due to the exchange of all the states with spin J. A simple way of determining these functions is to go to the center of mass frame and act with the SO(d − 1) quadratic Casimir on the two outgoing particle, while holding the momentum of the two incoming particles fixed. This equation takes the form where z = cos θ is cosine of the scattering angle (2.20). This second order differential equation has two independent solutions. Spin J unitary representations are composed of states with angular momentum in the plane of scattering ranging between −J and J. Hence, the corresponding solution of (2.30) is a degree J polynomial of cos θ that is given by The partial wave coefficients can be extracted from the amplitude using the orthogonality relation of these polynomials (2.32) 9 For d = 4 the factor is , see [36] for details.
Here we have chosen the convention for which the unitarity constraint presented below takes a simple form. In this convention we have Because the S-matrix is diagonal in the spin basis, so does the unitary constraint. We consider first the elastic regime 4m 2 < s < 16m 2 where this constraints takes the form (2.22). Using (2.34), we project both sides to a fixed spin J. On the left hand side we find the discontinuity of the partial wave coefficient. Real analyticity (2.13) of T (s, t) leads to real analyticity of f J Hence, the discontinuity of the partial wave is equal to the imgionary part 1 2i (f J (s + i ) − f J (s − i )) = Imf J (s). On the right hand side, it is useful to first represent the kernel as a sum over partial waves of z 1 , z 2 and z. Because this kernel represents the angular integration in (2.18), its partial wave decomposition must also be diagonal in spin. It takes the form (see appendix C) Using (2.34) for the three integrals and real analyticity (2.35), we arrive at the elastic unitarity constraint 37) or equivalently Here 1 can be traced to back to1 1 in (2.1). In this way the trivial unitary S-matrixŜ =1 1 becomes S J = 1 in the partial wave basis.
The solution to this is with δ J (s) being real for 4m 2 < s < 16m 2 and is called the scattering phase.
Similarly to the above, for s > 16m 2 we chose ψ(p 1 , p 2 ) = p 1 , p 2 |p, J, m in (2.24). In that way we arrive at the same equation, but with an inequality instead of an equality We close this section with a short discussion on the range of convergence of the partial wave sum (2.29) for fixed physical s as a function of cos θ. 10 It is a well-known fact that the amplitude T (s, cos θ) is analytic inside the small Lehmann-Martin ellipse and its absorptive part, T s (s, cos θ) is analytic inside the large Lehmann-Martin ellipse. These ellipses have foci at cos θ = ±1 and semi-major axis z small and z large . Correspondingly, inside these ellipses the sum (2.29) and its discontinuity converge.
In the case of scattering of identical lightest particles which is our main interest we have (2.42) In section 4.1 we will see that extended analyticity, elastic unitarity and crossing imply that the partial wave expansion converges in a larger region.

Froissart-Gribov Formula
The Froissart-Gribov formula is a representation of the partial wave coefficients in terms of the discontinuity of the amplitude. It has multiple applications and, in particular, it allows us to analytically continue partial wave coefficients in spin. Correspondingly, in section 3.3 we will use the Froissart-Gribov formula to analytically continue elastic unitarity in spin. It will also allow us to better understand the analytic structure of the amplitude in the Mandelstam invariants and relate the threshold expansion, see section 5, to the large spin expansion, see section 6.
Let us introduce the Gegenbauer Q-functions. These are given by the second linearly independent solution of the second order Casimir equation (2.30). They are uniquely fixed by their asymptotic behavior J is a normalization constant. The corresponding Q-function is Our convention is .
(2.45) 10 The convergence of the partial wave expansion can be seen using Neumann's argument [37]. Consider a function f (z) analytic inside some region C which includes the [−1, 1] interval. We can then write where γ ∈ C is some contour that wraps the interval [−1, 1] counterclockwise and contains z inside the integration contour. To exchange the summation and integration, we also used that given z, 1 J (t) converges uniformly in t as long as t is outside the ellipse with foci at −1 and 1 that passes through z. We also used the relation between Q  J (z) which will be explained below, see (2.46). Therefore, the partial wave expansion (2.41) converges when z is inside the ellipse with foci at −1 and 1 and ∈ C. Note therefore that the size of the domain of convergence of the partial wave expansion can be less that the analyticity domain C. We partially open up the contour. Sometimes this representation for partial waves is called the truncated Froissart-Gribov formula. The advantage of this representation is that we only use a finite amount of extended analyticity that has not been rigorously proven. c. We open the contour all the way to infinity and arrive at the usual Froissart-Gribov formula (2.49) with two integrations of the discontinuity of the amplitude along the t-channel and u-channel cuts (in black).
The Q-function has a cut running between z = −1 and z = 1. The fact that there are only two independent solutions to the Casimir equation means that the discontinuity of Q can be expressed in terms of Q and P . The precise relation takes the form 46) or equivalently (for integer J) We can then plug (2.46) into the partial wave coefficient (2.34) as where the integral is counterclockwise around the interval z ∈ [−1, 1]. By blowing up the contour, we get two integrals along the t-and the u-channel cuts, see figure 2 and we have assumed that s > 4m 2 , so that the t channel cut runs from z 1 = z 1 > 1 to infinity.
Here we have dropped the contributions of the arcs at infinity. This is justified for large enough spin J > J 0 (s) using (2.43), where J 0 (s) is the Regge intercept lim |t|→∞ |T (s, t)| < |t| J 0 (s) . (2.51) We can now use crossing to simplify (2.49). We change the integration variable for the u-channel integral from z to −z. Crossing symmetry implies that T u (s, u(z)) = −T t (s, t(−z)), where we have used that z(u) = −z. Under this change of variables We get that f J = 0 for odd J. For even J we get As opposed to (2.34), the Froissart-Gribov representation of the partial waves (2.53) is suitable for analytic continuation in J. It follows from the Carlson theorem that this analytic continuation is the unique continuation that does not grow too fast at large J. The Froissart-Gribov integral (2.53) converges as long as ReJ > J 0 (s) thanks to (2.51) and (2.43).
This integral is written for s > 4m 2 . As s approaches the threshold from above s − 4m 2 → 0 + , the lower end of the integral is pushed to infinity, z 1 → ∞. To analyze f J (s) in this limit, it is useful to use (2.43) and to switch back to an integral over t. In that way one finds This integral should be understood as follows. The large t contribution is finite because |Im t T (4m 4 , t)| < |T (4m 2 , t)| < t J 0 (4m 2 ) and J > J 0 (4m 2 ) by assumption. If the integrand diverges at some finite t, and in particular as t − 4m 2 → 0 + , then we should step back and write it as a contour integral of T (4m 2 , t) around the cut, which is manifestly finite.

Further Continuity Assumption
Based on the standard QFT axioms, scattering amplitudes, cross sections and partial waves are distributions rather than continuous functions. This leads to various subtleties, sometimes known as pathologies a-la Martin [19]. For example, we can imagine total cross-section having singularities which are local in energy variable s, that are not detectable by the finite resolution experiments. As such it is hard to exclude them based on physical grounds. One way to produce such a singularity is to consider an infinite number of resonances that accumulate on the real axis, see [19].
To the best of our knowledge there is no known, first principle argument that can exclude these possibilities. One way to eliminate them is to simply assume that various cross sections are continuous functions of energy s. We adopt this practical approach add this to the list of our assumptions. More precisely, we assume that boundary values of T (s, t) (this includes both single and double discontinuity) are continuous functions. It is common in the literature to impose the condition that boundary values of T (s, t) are uniformly continuous or Hölder continuous, see e.g. [38], but we will not use it in the present analysis.
Another related common assumption is regarding finiteness of scattering lengths which are commonly measured in the experiments or using the lattice. They are defined as follows, see (2.54), (2.55) We will assume that the scattering lengths are finite for J ≥ 2. Through the Froissart-Gribov formula these are related to the assumption of finiteness of the discontinuity T t (4m 2 , t) at s = 4m 2 as well as convergence of the J = 2 Froissart-Gribov integral (2.54).
There is an interesting connection between the continuity of the amplitude and macrocausality [39,40]. Macrocausality is a set of statements about scattering amplitudes when particles grouped according to space and time of interactions and then moved away from each other by large translation. The notion relevant for the continuity of scattering amplitudes is what is called strong asymptotic causality in [39] and it has not been proved within the field theory.

Analytic Continuation of Elastic Unitarity
The elastic unitarity relations (2.21) was derived for energies in the elastic region, above the two particle threshold 4m 2 < s < 16m 2 , and for physical kinematics 4m 2 − s < t < 0. In this section we analytically continue this relation in t, outside of the regime of real scattering angles. We also consider the double discontinuity of the amplitude and the closely related analytic continuation of elastic unitarity in spin.

Mandelstam Kernel
The dependence on the scattering angle z = cos θ enters the right-hand side of the elastic unitarity relation (2.22) through the kernels (2.21). These kernels contain a delta or a step functions and are thus not suitable for analytic continuation. To overcome this difficulty, we use the analyticity of T (±) (s, t(z)) inside the small Lehmann-Martin ellipse to express them as a counterclockwise Cauchy integral around [−1, 1] We can now exchange the order on integrations in (2.22) and perform the z and z integrals explicitly. In this way we arrive at where the new kernel is .
These integrals are evaluated in appendix A. For |η |, |η | > 1 the result is The Mandelstam kernel for |η | < 1 or |η | < 1 is obtained from (3.4) by analytic continuation. Note that in (3.3) P d (z, z , z ) is not analytic in z, see (2.21). On the other hand, the Mandelstam kernel K d (z, η , η ) is analytic in z and therefore is suitable for analytic continuation.
Similarly, the representation of T s (s, z) in (3.2) is now suitable for analytic continuation in t. That is because t only enters through the Mandelstam kernel that is manifestly analytic in z.
As for the kernel P . We then note that the z and z integration has the effect of converting the partial waves P     (2.47). In that way we arrive at

The Double Spectral Density and Crossing
The combination of elastic unitary with crossing symmetry is very restrictive. In the next sections we will explore some of its consequences at length. With this aim in mind, we now represent elastic unitarity in a form that is more suitable for imposing crossing. By taking another discontinuity of (3.2) with respect to t, the left-hand side becomes the double discontinuity of the amplitude 11 = Disc t Disc s T (s, t) = Disc s Disc t T (s, t) = ρ(t, s) . This crossing-symmetric function is known as the double spectral density. 11 Stated in words, the double spectral density ρ(s, t) is defined as a certain combination of boundary values of the analytic function T (s, t) unambiguously specified by i in the definition above. . As s or t approaches the Landau curve from above, the integration region shrinks to zero. As a result, the double spectral density vanishes below the Landau curve z = 2z 2 1 − 1.
By taking the t discontinuity of (3.2) we arrive at where 4m 2 < s < 16m 2 . Outside of the elastic region there are additional non-elastic contributions to the double discontinuity that will be considered elsewhere. Next, we deform the η and η integration to wrap the t and u channel cuts of T (±) (s, η). The discontinuity of the Mandelstam kernel in (3.8) is analytic in η and η along the deformation. In that way, we end with real η and η that are positive on the t channel cut and are negative on the u channel cut. For η η < −1 the integral in the Mandelstam kernel (3.4) starts at η + (η , η ) < −1 and can be chosen to run along the negative real η-axis. This choice make it manifest that the discontinuity of the kernel for z > 1 and η + (η , η ) < −1 is zero. We remain with 12 where we have mapped the u-channel cut to the t-channel cut using T (s, −η) and K(z, −η , −η ) = K(z, η , η ). Here, the lower limit of integration is the point where the t channel cut starts (2.50). The discontinuity of the kernel, for η η > 0 and z > 1 is given by (3.10) 12 The result for z < −1 is obtained by analytic continuation and takes the same form as (3.9) with DisczK(−z, η , η ) instead of DisczK(z, η , η ). This formula was first derived by Mandelstam in d = 4 [41]. In appendix A we present the derivation for any dimension d ≥ 3.
We now discuss the region of support of the double discontinuity ρ(s, t). Because the first discontinuity has only support for energies above the two particle threshold, so does the double discontinuity. Looking at (3.9), it is clear that the situation is more interesting. The discontinuities of the amplitude only start at η , η > z 1 . At the same time, the Mandelstam kernel has a nonzero support only for z > η + (η , η ). This constraint is an upper bound on the η , η limits of integration, see figure 3. Hence, the double spectral density in the elastic region is non-zero only for (3.11) In terms of s and t, the boundary of this region is known as the Landau or Karplus curve and is given by (red in figure 4) As t is increased above t 1 (s), the range of the η and η integration opens up. For example, for any value of t > t 1 (s), the integral over η is bounded in the range Hence, as we increase t inside the elastic region 4m 2 < s < 16m 2 (where (3.9) is valid), more and more channels of T (+) t (s, η ) and T (−) t (s, η ) kick in. Their corresponding contributions to ρ(s, t) start at other Landau curves in the (s, t) plane that are above the elastic one (3.12), see figure 4. In section 3.5 we present a minimal and complete set of Landau curves in the elastic region 4m 2 < s < 16m 2 that appears in the physical S-matrices.
In general, the positions of the Landau curves look somewhat technical. However, they are all kinematic and therefore have a geometrical origin. For example, the constraint (3.11) takes a simple form in terms of the integral scattering angles η = cosh θ , η = cosh θ , and is given by (3.14) In physical kinematics, this constraint follows from a simple geometrical consideration that is described in figure 1. The analytic continuation to the non-physical kinetatical regime of positive s and t effect the range of the angles, but leaves this geometrical constraint unchanged. At the technical level, this is because the kinematical constraint only involves the cosines of the scattering angles.
We end this section with a comment regarding the region below the Landau curve, where the double discontinuity vanishes (the gray region in figure 4). For the two-to-two scattering, existence of this region is a direct consequence of elastic unitarity continued to s, t > 4m 2 . The precise shape of the region depends on the details of the unitarity kernel K d , as well as on the analytic structure of the amplitude that enters into the elastic unitarity relation. A similar phenomenon occurs in the higher-point amplitudes as well. In this case one considers double discontinuity in the so-called overlapping channels, see e.g. [42] for a detailed definition. This time it is possible to consider double discontinuity for physical kinematics directly, as opposed to the two-to-two case which requires continuation in one of the Mandelstam invariants. It then follows that for physical kinematics the double discontinuity vanishes for the overlapping channels. These are known as the Steinmann relations [43,44], and it is again a direct consequence of the multi-particle unitarity. Steinmann relations are useful for constraining the analytic structure of amplitudes with six and more particles, see [45] for a recent discussion. Based on the two-to-two case one can try to analytically continue the relevant unitarity relations to the unphysical values of the relevant kinematical invariants, and find the extended region where the double discontinuity vanishes. It would be interesting to understand the precise shape of this region for the multi-particle case. This would require analytic continuation of multiparticle unitarity kernels, as well as relevant scattering amplitudes, analogous to the one made above.

Analytic Continuation in Spin
We will now argue that the elastic unitarity relation (2.37) holds in the complex J plane, provided that the partial wave coefficients are analytically continued using the Froissart-Gribov representation (2.53). For that aim, we first rewrite the elastic unitarity condition (2.37) with the i prescription explicitly In our previous discussions f J (s) ≡ f J (s+i ). For integer J and real s in that range, real analyticity (2.35) leads to (2.37). The form (3.15) is however more suitable for analytic continuation.
Originally, (3.15) was derived for J being integer and even, however using the Froissart-Gribov representation we can continue partial waves in spin J This representation can in principle be continued to the whole complex J plane, going beyond the Re[J] > Re[J 0 (s)] region. However, the continued partial waves are not guaranteed to coincide with the physical ones for J < J 0 . Moreover, it is clear from (2.13) and the reality properties of Q (d) 16), that real analyticity of partial waves (2.35) is only guaranteed to hold for real J > J 0 .
We can then use (3.16) to separately analytically continue the left and right hand sides of (3.15). It follows from the Carlson theorem that these two analytic continuations have to agree in the whole complex J plane. Namely, the two analytic continuations agree for real integer J > J 0 and do not grow too fast at large |J|, as can be seen from the large J exponential decay of f J (s), see section 6.1 for more details, Let us relate the analytic continuation in J of (3.15) to the Mandelstam equation (3.9). To get the former from the latter we integrate (3.9) with (2.47). We then use the following identity (for derivation see appendix B) which is valid for Re[J] > −1, d ≥ 3 and |η 1 |, |η 2 | > 1. In practice, we will use this identity in (3.9), where η and η are real and positive.
Let us also mention that analytically continued elastic unitarity constrain the possible Regge limit behavior of the amplitude. In particular, the leading Regge singularity of the amplitude in the elastic region 4m 2 < s < 16m 2 cannot we a pole located at some real Regge spinJ 0 (s). We review the derivation of this and a slightly more general result, known as Gribov's theorem, in appendix D.

Analytic Continuation in s
As s is increased above 16m 2 there are new multi-particle cuts contributions to Imf J (s) that are not captured by (2.37). Still, if we denote by f J (s) the partial wave that was analytically continues on the second sheet through the elastic cut 4m 4 < s < 16m 2 , then the equation holds true when analytically continued away from the elastic strip in the full multi-sheet complex s plane. In that sense, elastic unitarity can be analytically continued in s.
Note that there is a difference in the nature of the two-particle cut between odd and even dimensions. Due to the power of (s − 4m 2 ) in (2.38) we get that Correspondingly, for even dimension the continued elastic unitarity condition (3.20) becomes This relation implies that the elastic cut 4m 2 < s < 16m 2 describes a two-sheeted Riemann surface in even d. Indeed, continuing (3.22) around the elastic cut again we conclude that S 2× J (s) = S J (s). In odd d, because of the sign difference in (3.21), we do not have (3.22) and the elastic cut is not two-sheeted. Similar conclusion holds for T (s, t). As we will see below, the nature of the elastic unitarity cut is the one of a simple square-root for even d and is infinitely-sheeted for odd d. It immediately follows from (3.22) that a pole of S J (s) (resonance) corresponds to a zero of S J (s) in

Elastic Landau Curves
In this section we further study the kinematical properties of the double spectral density in the elastic strip. In general, discontinuities of the amplitude result from the exchange of intermediate on-shell states. Correspondingly, the double discontinuity of the amplitude receives its support from intermediate on-shell particles exchanged in both channels. As we have seen in section 3.2, the first double discontinuity in the elastic strip starts at the leading Landau curve (3.12). It comes from the exchange of two particle in the s-channel and four particles in the t-channel. As we increase s or t, more and more processes become accessible. The curves in the s − t plane where they start to contribute are higher Landau curves. We now derive an infinite set of higher Landau curves in the elastic strip that are required by elastic unitarity.
To each of these Landau curves one can associate a Landau diagram. In the elastic strip, all these diagrams have a simple iterative structure that is plotted in figure 5. Correspondingly, the shape of the curve has a relatively simple geometrical origin, (for the case of the first Landau curve, see figure 1 and the discussion around (3.14)). Here instead, we follow a shortcut and derive their shape directly from the Froissart-Gribov representation and the elastic unitarity constraint.
The starting point is the kinematical structure of the discontinuity and the double discontinuity of the amplitude as function of one of the Mandelstam invariants (3.24) Figure 5: a. The Landau diagrams that contribute to the double spectral density in the elastic strip 4m 2 < s < 16m 2 and 4m 2 < t. In this kinematical regime there can be only two particles in the s-channel.
Hence, the corresponding Landau diagrams have a simple structure of iterative two-particle exchange in the s-channel and in between, any number of particles exchange in the t-channel. b. Analogously to figure 1, the corresponding Landau curves originate from a simple geometric constraint on physical kinematics, see discussion after (3.14).
where the functions t (i) (s, t) are the Landau curves we are after. For that aim, we first have to analytically continue the discontinuities of the amplitudes, T 2→2n t (s, t), to the regime of positive s and t. As we do so, they develop new thresholds. Because T t (s, t) is a real analytic function, these new thresholds must coincide with the Landau curves t (i) (s). What is important for us here is that the multi-particle thresholds that are manifest in (3.23) are also present in the non-physical kinematical regime.
To derive the functional shape of the Landau curves in the elastic strip we impose the consistency of (3.24) and (analytically continued) (3.23) with elastic unitarity. This can be done by either plugging them into the Mandelstam equation or by imposing elastic unitarity at the level of the partial waves. Here we follow the latter strategy and in appendix E we present the former, both leading to the same result.
To impose consistency of (3.23) and (3.24) with the partial waves elastic unitarity, we plug them into the Froissart-Gribov projection (2.53), that we quote here for convenience Next, we take the large J limit of (3.25). Using the large J decay of the Q functions that can be schematically written as, see section 6.1, as well as the two-particle step function in (3.23), we conclude that the leading large J behavior of Similarly to (3.27), the n-particle threshold of T t (s, t) and the i'th Landau curve threshold of ρ(s, t) in (3.23) result in contributions to Re f J (s) and Im f J (s) that start at large J as where z n ≡ 1 + 8n 2 m 2 s−4m 2 . Importantly, (3.26) receives only 1/J-power corrections and no nonperturbative exponential corrections, see appendix C. Hence, Im f J (s) does not receive any exponential large J behavior other than the ones that result from Landau curves thresholds (3.28). Similarly, the exponential large J behavior of Ref J (s) can only come from thresholds of T t (s, t) in (3.25), but not from nonperturbative terms in the large J expansion of Q (d) J (z). We now plug these set of large J exponential behaviors into elastic unitarity and derive the minimal set of Landau curves, t (i) (s), that are required to close the equation, together with new thresholds of T t (s, t) that are not present in physical kinematics.
Elastic unitarity (2.37) can be written in the following schematic form where we omitted the pre-factor because it is irrelevant for the present discussion. First, we see that the multi-particle threshold exponents of Ref J (s) in (3.28) result in the following exponents in Im f J (s) Using (3.28), these then lead to the Landau curves For example, the leading Landau curve (3.12) corresponds to the case where m = n = 1. This behavior must then also exist in Ref J (s). To see that, recall that T t (s, t) is a real analytic function of t for fixed s > 4m 2 . Hence, a threshold in the imaginary part of T t (s, t), namely in ρ(s, t), must be accompanied by a corresponding threshold in Re T t (s, t).
Having established the presence of the quadratic in λ(z n ) terms in the large J expansion of Ref J (s) and Im f J (s), we now come back to elastic unitarity (3.29). The presence of these terms together with the linear one (3.28), now induces higher powers of λ(z n ) in the large J expansion of Im f J (s) and, via real analyticity, of Ref J (s) as well. These take the general form This large J structure implies existence of an infinite set of Landau curves labeled by a set of integers {N 1 , N 2 , N 3 , . . . } and given by the following equation (3.34) In terms of the scattering angle, for which z = cosh θ and λ = e θ , these Landau curves take the form with θ n = arccos(z n ). This form is the generalization of (3.14) for all the elastic Landau curves.
In a direct analogy with figure 1, it also has a simple geometrical origin that is discussed in figure  5.b. Using that z = 1 + 2t s−4m 2 , we can translate (3.34) into polynomial equations in s and t with real coefficients. 13 For example, we have 13 In this form they are familiar in the study of perturbative Feynman integrals, see e.g. [1] for more details. Figure 6: The regime of positivity of the double spectral density ρ(s, t) in the elastic strip 4m 2 < t < 16m 2 , (dashed region). This is the region above the first Karplus curve s 1 (t) (in blue) and below the curve There is an identical positive region in the crossed strip of 4m 2 < s < 16m 2 , see (4.4).
The Landau curves have the following asymptotic behavior (3.37) These asymptotes are precisely the t-channel normal thresholds (3.23). This is evident from the Landau diagram interpretation of figure 5 and it would be interesting to see if it can be established rigorously. 14

Positivity of ρ and Multi-Particle Production
Drawing an analogy with the conformal bootstrap, one can expect that combining crossing symmetry of ρ(s, t) with some sort of positivity property leads to nontrivial constraints on scattering amplitudes.
In this section we show that there is indeed a finite region in the elastic strips of the (s, t)-plane, where ρ(s, t) is positive, as noted in [50]. In subsection 4.2 we discuss a direct consequence of this positivity and crossing -the necessity of multi-particle production T 2→2n for any n.

Positivity of The Double Spectral Density
To establish positivity, we assume that the integral in the Mandelstam equation for the double spectral density (3.9) converges and study under what conditions the integrand is positive.
Consider first the discontinuity of the Mandelstam kernel Disc z K d (z, η , η ) given in (3.10). This function is strictly positive for z > 1, which is the case for any s, t > 4m 2 . Next we turn our attention to T (±) t (s, η ). It is clear that T t (s, η ) is non-negative in the region where the partial wave expansion converges 15 where t is not to be confused with the external t in (3.9). Indeed, unitarity implies that Imf J (t ) ≥ 0 and P is positive for s, t > 4m 2 . Consider the function g 1 + 2s t −4m 2 ≡ T t (s, t ) for fixed t . Provided that g(z) is analytic inside an ellipse in the complex z-plane with foci at z = ±1, it then follows from Neumann's argument, see footnote 10, that the partial wave expansion (4.1) converges inside that ellipse. Under the assumption of extended analyticity and for fixed t , the first singularity is where T t (s, t ) develops a discontinuity with respect to s. Namely, the partial wave expansion (4.1) converges below the first Landau curve, red in figure 4, 16 This region of convergence is called the large s-channel Martin-Lehmann ellipse [28]. Note that the standard argument of [28] refers to the convergence of (3.19) in the crossed region of s < s 1 (t ), 4m 2 < t < 16m 2 and does not require the extended analyticity assumption. Here, we are saying that under this assumption, the partial wave expansion (4.1) converges in the union of the t-and s-channel large Martin-Lehmann ellipses.
The maximal value of the η integration is given in (3.13). Hence, we conclude that the double spectral density is non-negative for By solving this condition for t we conclude that where the lower bound comes from (3.12), and the upper bound from (4.3). This region of positivity was first derived in [50] and is plotted in figure 6. 17 Note that at the upper bound of that region ρ(s, t) is strictly positive. Hence, the region of positivity can always be extended from above. A question to which we do not know the answer is whether the region of positivity of the double spectral density can be extended to arbitrary large s and t. Here we only list some known implications of such a scenario. In [51] A. Martin argued that positivity of the double spectral density for all s and t implies that the total cross-section at high energies satisfies the following bound This immediately implies that for theories that saturate the Martin-Froissart bound σ tot (s) ∼ log 2 s [28,52], the double spectral density is not positive-definite.
Similarly, in [53] A. W. Martin explored implications of the positive-definite double spectral density for scattering at finite angles. He showed that amplitudes with positive double spectral density do not exhibit diffraction peak in the near-forward scattering, which is where b(t) is a slowly-varying function. Such a diffraction peak is observed, for example, in the scattering of pions. Therefore, the double spectral density for pion scattering cannot be positive for arbitrary s and t.
It would be interesting to understand positivity properties of ρ(s, t) in physical theories more systematically.

The Aks Theorem and Necessity of Particle Production
We now review an elegant argument for scalar particles by Aks [14]. 18 It states that scattering implies particle production. Namely, provided that T 2→2 = 0, also T 2→n = 0 with n > 2. The theorem applies to any crossing symmetric scalar scattering amplitude in d ≥ 3 that satisfies extended analyticity in a finite region above the leading Landau curve.
To derive Aks's result, let us therefore assume that we have a nontrivial scattering amplitude T (s, t), but T 2→n are identically zero for n > 2. This implies that elastic unitarity in the form of Mandelstam (3.9) holds for any s ≥ 4m 2 , whereas in the the theories with particle production elastic unitarity only holds below the first threshold s 0 > s ≥ 4m 2 . It then follow that ρ(s, t) = 0 below the first s-channel Landau curve 4m 2 < t < t 1 (s) (3.12) for any s, see red curve in figure  4. This region however includes the crossed Mahoux-Martin region of positive ρ(s, t) > 0 (4.4), s < 4m 2 (3t+4m 2 ) 2 (t−4m 2 ) 2 , see gray region in figure 6. 19 Therefore, we conclude that our assumptions of crossing symmetry and scattering without production are not consistent.
Note that in the case of an infinite J = 0 scattering length, the derivation of the positivity of ρ(s, t) in the Mahoux-Martin region does not always apply. That is because in that case, the integral in (3.9) fails to converge for d ≥ 5. However, when the J = 0 scattering length is infinite, ρ(s, t) is already non-zero at the leading Landau curve and we reach the same conclusion that there must be particle production.
A way to relax the assumption of extended analyticity was explained in [55] by exploiting polynomial boundedness and continuity assumptions. The idea is to use elastic unitarity to extend the region of analyticity of the amplitude. One starts with fixed t < 4m 2 dispersion relations and then use elastic unitarity to first extend the analyticity domain of T s and then use dispersion relations to continue the scattering amplitude itself. The key point being that assuming no production we can use elastic unitarity to continue the discontinuity at arbitrary energy, which is necessary if we want to use this inside the dispersion relations.
Note that the argument above also implies that we must have four-particle production. That is because the crossed region of positivity starts at s 1 (16m 2 ) = 64 3 m 2 < 36m 2 , it is enough to assume that T 2→4 = 0 to reach a contradiction. One can wonder if having T 2→4 is enough to fix the problem, or T 2→2n with n > 2 are also necessary? To address this question, we can then proceed via crossing. 20 By unitarity of the 4 → 4 amplitude, ImT 4→4 ∼ |T 2→4 | 2 , non-vanishing T 2→4 implies that we have a non-vanishing T 4→4 amplitude. Applying crossing this becomes T 2→6 . Continuing this recursion we conclude that all T 2→2n should be non-zero. Therefore, not only scattering implies production but it requires all possible production (here we assumed Z 2 symmetry so that only an even number of particles is present in the final state). An alternative argument that does not use crossing symmetry of higher-point amplitudes was put forward in [56]. This argument relies on an unproven assumption that Landau curves can only asymptote to normal thresholds for the scattering of lightest particles, see section 3.5. This assumption is consistent with the perturbative analysis of [29][30][31] and it would be interesting to establish it rigorously.

Bounding Inelasticity
There is another, more intuitive way to think about the result of Aks and necessity of particle production in higher dimensions. Ideally, one would like to take a discontinuity of the 2 → 4 amplitude that is given by a product of two 2 → 2 amplitudes. For physical kinematics however, such a discontinuity only exist for the 3 → 3 setup. Instead, let us discuss impact parameter scattering, which is only possible in d ≥ 3.
As reviewed for example in appendix E of [57], the effect of going to the impact parameter space is the same as continuing the conjugate momentum invariant to the unphysical kinematics. It follows that inelasticity in a gapped theory cannot be exactly zero at very large impact parameters. To see this, decompose the four particles in the final state into a pair of dipoles. Then consider a scattering in which the two dipoles in the final state, as well as the pair of incoming particles in the initial state, are separated by a finite distance b in the transverse space, see figure 7.a. Unitarity becomes, in the impact parameter space, the expansion in Yukawa potential suppressed . At large separation, this expansion is dominated by the one-particle exchange e −bm , while the multi-particle corrections are further exponentially suppressed. 21 In that way, a non-trivial (analytically continued) four-point amplitude imply a non-trivial 2 → 4 amplitude. We would then like to bound the 2 → 4 amplitude from below.
There is a convenient way of bounding the integrated discontinuity of T 2→4 in the kinematical regime of figure 7.b from below. It is based on the discussion in the previous sections, see in particular sections 3.2 and 4.1, 4.2. It also highlights how the crossing of ρ(s, t) discussed in the section above works microscopically. Let us start with the the square of T 2→4 that appears in the discontinuity of the 2 → 2 amplitude, T s (s, t) (4.7) By construction, the unitarity integral in the right-hand side of (4.7) depends only on s and t. For physical scattering we consider s > 16m 2 and t < 0.
We would like next to analytically continue (4.7) to the unphysical Martin-Mahoux region discussed above We also would like to consider 16m 2 < s < 36m 2 to focus on the T 2→4 amplitude. This condition together with (4.8) imply 36m 2 5 < t. For the 2 → 2 case we reviewed the procedure in the sections above. For the multi-particle case it was discussed in [46][47][48][49].
After taking a discontinuity in t and using crossing symmetry of the double spectral density, we arrive at the following schematic form where in the formula above we switched toz = 1+ 2s t−4m 2 andz 1 = 1+ 8m 2 t−4m 2 . Here in the right-hand side each DiscT 2→4 contains a delta-function that puts the exchanged particle in figure 7 on-shell. The phase space integral dLIPS 4 should be understood in terms of the analytic continuation a-la Mandelstam, and we will discuss it in more detail elsewhere [58].
Note that since we are in the Mahoux-Martin region, the partial wave expansion of Disc s T (±) 2→2 (t, η) converges and both the Mandelstam kernel and the Legendre polynomials that enter into partial waves are positive. Therefore we can write (4.10) By plugging (4.10) into (4.9) and using positivity of the Mandelstam kernel DisczK d (z, η , η ) we arrive at the lower bound on the integrated discontinuity Disc t T 2→4 .
While this argument bounds the discontinuity of T 2→4 in the unphysical kinematics of figure  7.b, for 36m 2 5 ≤ t < 16m 2 from below, it does not tell us anything about T 2→4 in the physical kinematics t < 0. Let us however comment why one expects to have T 2→4 of the same order also in the physical kinematics. Schematically, we can write the following representation of the 2 → 4 impact parameter amplitude We can rewrite this as follows where c MP encodes the contribution of the multi-particle exchanges.
Without extra fine-tuning we expect that in a strongly coupled theory T 2→2 , c MP ∼ O(1) and therefore the one-particle exchange to dominate for b 1 m . This would make T 2→4 ∼ O(1). On the other hand, making T 2→4 1 given T 2→2 ∼ O(1) would require a very fine-tuned cancellation between the exchange of one-particle state and the multi-particle state as well as c MP 1. It seems quite possible that such a scenario is not consistent with multi-particle unitarity, but it is very hard to show it explicitly due to the complexity of the latter.
One scenario in which the regime of single particle dominance can be delayed to arbitrary large impact parameters is if the theory contains extended objects, such as strings. In such case we can choose l string 1 m , and the one-particle exchange is expected to dominate only for b l string . In this case however, the spectrum is expected to contain particles of mass 1 l string m which contradicts our assumption about m being the lightest particle.
In some numerical applications that we discuss in more detail below, see e.g. [15], it was observed that small J partial wave converge very well. Therefore we can set in the formula (4.10) J 0 = 2 and then use it in (4.9). We can then evaluate the integral above to rigorously bound from below DiscT 2→4 in the Mahoux-Martin kinematical region. In the numerical analysis when maximizing some couplings, one usually finds Imf 0 (s) ∼ 1 therefore the formulas (4.9), (4.10) will produce DiscT 2→4 ∼ 1. It is however still an open problem to translate this fact into a rigorous statement about T 2→4 for physical kinematics.
Independently of the discussion above, in section 7 we bound inelasticity using an additional input about the discontinuity of the 2 → 2 amplitude.

Threshold Expansion
To further analytically constrain a nonperturbative amplitude a small parameter is needed. In this section, as well as the next one, we will study the expansions of the amplitude in two kinematical small parameters, as well as the relation between the two expansions.
One kinematical small parameter that always exists in a gapped theory is the energy distance from the two-particle thereshold in units of the mass gap, where in the last equality we evaluated σ s in the center-of-mass frame (2.7).
The threshold expansion is, thus, an expansion in powers of σ s . The small σ s expansion is known in nuclear physics as the effective range expansion [59]. As we take σ s 1 or, equivalently, | p| m scattering becomes non-relativistic and we can characterize it by some effective potential whose properties are captured by the threshold expansion parameters [60]. In this section we use elastic unitarity, extended analyticity and crossing to argue that both the discontinuity and the double discontinuity of the scattering amplitude admit a natural threshold expansion. Importantly, the parameters that enter the threshold expansion are controlled by the low energy physics which is well-known from the experiments or lattice simulations, see e.g. [61].
Whenever a theory has a conserved charge, there is a natural small parameter -one over the charge of states that are being exchanged. In the case at hand, the only symmetry we assume is Lorentz symmetry. Correspondingly, the only available conserved charge is the spin. Indeed, in the next section we will show that the partial wave coefficients admit a systematic large J expansion. Moreover, by combining the Froissart-Gribov formula with crossing we will relate the threshold expansion to the large spin expansion.
The logical structure and steps of this and following sections is summarized in figure 8. We start by constructing a threshold expansion of the partial wave coefficients that is consistent with elastic unitarity.

Threshold Expansion of f J (s)
To expand the the partial wave coefficients close to the threshold in a way that is consistent with elastic unitarity, it is first convenient to solve the latter in a different fashion than (2.39). After dividing by |f J (s)| 2 , the elastic unitarity condition (2.37) takes the form Hence, 1/f J has a branch cut at threshold for even d and a logarithmic cut for odd d. The general solution to (5.2) combined with real analaticity (2.35) takes the form where b J (σ s ) is a real analytic and single-valued function in some finite neighborhood around the origin, except for potentially isolated singularity at σ s = 0. Note that in writing (5.3) we used the continuity assumption for f J (s), see assumption 5 at the introduction and section 2.5.
As we have learned from the Froissart-Gribov representation (see (2.54)), as σ s → 0 + from above, f J (s) → a J σ J s (2.55). In terms of the function b J (σ), this behavior is where we have assumed that J 0 (4m 2 ) < 2. For J = 0 we do not expect the Froissart-Gribov representation to hold close to σ s = 0. Therefore we do not have a similar prediction. After factoring out the leading threshold singularity of b J≥2 (σ), it has a regular expansion whereJ 0 is an unconstrained integer at this point and a J , b J,i are real coefficients. Reality of scattering lengths a J and effective ranges b i,J follows from real analyticity property of the partial waves coefficients (2.35). For actual values of these parameters in various QCD processes see e.g. [62].
Similarly, using real analaticity of b J (s), we have the following equations for the imaginary part of the partial wave in the elastic unitarity region To summarize, by plugging (5.5) into (5.3) and (5.6) we have constructed a threshold expansion of f J (s) and Imf J (s) that automatically solves elastic unitarity for any real coefficients. Note that the most singular possible threshold behavior is completely fixed by elastic unitarity, with no free coefficient. It comes from J = 0 and corresponds to setting b 0 (σ s ) = 0 in (5.3), so that the partial wave is dominated by the universal term In particular, this implies that the spin zero scattering length is infinite. In d = 3 the threshold behavior (5.7) is realized in massive φ 4 theory [63]. In d = 4 the asymptotic behavior (5.7) appeared in the study of coupling maximization in [15] and corresponds to having a bound state at s = 4m 2 . 22 For J ≥ 2 ≥ J 0 (s) the singular behavior (5.7) does not occur due to (2.54).
Below we will also be interested in the case of finite Imf 0 (s). In this case, the spin zero scattering length is finite, namely b 0 (s) = m d−4 a 0 , and we get the following leading behavior of Imf 0 (s) As a final remark, note that the analysis of this section can also be generalized to non-integer spin, see appendix F.

Threshold Expansion of Discontinuity
For our purposes we will be interested in a closely related expansion, that of the discontinuity of the amplitude T s (s, t) for t > 4m 2 and σ s → 0. As we will see, the Froissart-Gribov formula (2.53) relates this expansion to the large J behavior of the partial wave coefficients.
We proceed by considering the s-channel partial wave expansion for the discontinuity of the amplitude For s close to 4m 2 , the sum over J above converges inside the large s-channel Lehmann ellipse which is for t < 16m 2 s s−4m 2 . Given some fixed t, and considering the limit s → 4m 2 we stay within the convergence region.
We can now plug (5.6) into (5.9) and perform the threshold expansion under the sum. We get where we can systematically expand T s (s, t) using the threshold expansion of partial waves (5.5) which is based on elastic unitarity. Note that the expansion parameter in even dimensions is simply σ s , whereas in odd dimensions we have in addition powers of log σ s , as well as inverse powers 1 log σs , which we do not write here explicitly.
Let us emphasize that in the argument above it was absolutely crucial to consider the discnontinuity of the amplitude T s (s, t) and not the amplitude T (s, t) itself. Indeed, if we were to try repeating the argument above for T (s, t) itself, we would run into the following problem. In (5.9) we have partial waves f J (s) ∼ (s − 4m 2 ) J and Legendre polynomials behaving as P Therefore the expansion for s → 4m 2 requires re-summation of partial waves of all spins which is beyond our control.
The conclusion is that T s (s, t) admits a systematic threshold expansion in terms of the solution to elastic unitarity in the s-channel. Moreover, the contribution of higher spin partial waves are suppressed by an additional factor of σ J s . We expect the threshold expansion converges as long as we stay below the leading Landau curve, namely for t < 16m 2 s s−4m 2 . Luckily for us, it is the discontinuity of the amplitude and not the amplitude itself that enters the Froissart-Gribov inversion formula. Therefore the results of this section can be readily put to use.

Threshold Expansion for Double Spectral Density
We now turn to the study the double spectral density (3.7) close to the leading Landau curve. We consider the Mandelstam equation (3.9) that is repeated here for convenience where 4m 2 < s < 16m 2 . The crucial observation is that when we are very close to the leading Landau curve t = 16m 2 s s−4m 2 , the integral in (5.11) is again controlled by the threshold expansion of T (+) t (s, η ) and T (−) t (s, η ). The reason for this is that the region of integration starts at the threshold, η , η ≥ z 1 , and ends at the boundary of support of the kernel, η + (η , η ) ≤ z. As t approach the leading Landau curve, this region shrinks close to the threshold for both η and η integrals, see figure 3.
A convenient small parameter in the problem is the dimensionless distance from the Landau curve, see (3.11) We can now plug the threshold expansion of the discontinuity (5.5) into (5.11) and expand the result in powers of δz. At any given order in that expansion only finite number of terms in the threshold expansion of the discontinuity contribute. For even spacetime dimension, the relevant integrals that appears in the expansion arẽ To leading order in δz this integral is given by 23 (1 + O(δz)) .

(5.15)
For n 1 , n 2 ≤ d−5 2 the integral (5.13) diverges, but this apparent divergence is not physical. Namely, the integral in (5.15) originates from the contour integral that wraps around η , η = z 1 . Therefore, the integral can be safely deformed to a keyhole contour around the dangerous region, see figure 9. The result of the keyhole integration is equivalent to analytically continuing (5.15) in m, n (as long as the final result is finite). 23 To derive (5.15) it is useful to do the following change of integration variables close to δz = 0 where the integration in (5.11) for δz 1 is restricted to 0 ≤ x and α ≤ 1. We see that indeed, higher powers in the threshold expansion result in higher powers of δz. For odd spacetime dimension, due to the presence of 1 log 2 σs+π 2 in (5.6), the relevant integral is more complicated, but the power suppression in δz is of course the same. We discuss it more in appendix H.2.
The most singular possible behavior close to the Landau curve comes from the corresponding universal threshold behavior of Im f 0 (t) in (5.7) and corresponds to n 1 = n 2 = 0 in (5.13). In that case we find

Radius of Convergence
Let us briefly discuss the radius of convergence of the threshold expansion introduced above. As usual the radius of convergence is controlled by the analytic properties of the function at hand.
Consider first f J (s). It has a normal threshold cut starting at s = 4m 2 and multi-particle cuts at s ≥ 16m 2 , as well as the u-channel cut for s ≤ 0. In addition to that we can have bound states and resonances (which correspond to poles on the second sheet). In a given theory, we expect the singularity that is the closest to s = 4m 2 to control the convergence radius of the threshold expansion of f J (s). In this paper we assume for simplicity that there are no bound states, however it should be easy to include them in the analysis.
Consider next T s (s, t). In this case, the partial wave expansion (5.9) converges below the first Landau curve. For 4m 2 < s < 16m 2 and t > 16m 2 , that is the regime of σ s < 16m 2 t−16m 2 . Moreover, if we are to first expand each Imf J (s) close to the threshold under the sum, we will have to argue that the sum over J and the threshold expansion commute. It is not clear to us how to do it. Instead below we adopt a different approach. We separate a few low spin partial waves in (5.9) and apply the threshold expansion to them. We then bound the sum over spins without doing the threshold expansion under the sum over spins.

Large J Expansion
In this section we map the threshold expansion developed in the previous section to the large J expansion of the partial waves coefficients f J (s). 24 The basic idea is very simple. The Froissart-Gribov formula (2.53) directly maps the threshold expansion of the discontinuity of the amplitude (5.10) to the large J expansion of the partial wave coefficients in the crossed channel. 25 Interestingly, once combined with crossing, the Froissart-Gribov formula allows us to map the low energy threshold expansion, that is dominated by the low spins, to the large spin behavior of the partial wave coefficients for any energy and in the same channel. In particular, it automatically predicts the amount of inelasticity at large J.
In this section we will use the results for the threshold expansion of the double spectral density and the discontinuity of the amplitude in the crossed t-channel as opposed to the s-channel threshold expansion of the previous section. We hope this will not cause any confusion. In sections 6.1 and 6.2 we work out the leading large J expressions for f J (s) and Imf J (s) correspondingly. In section 6.3 we compute the exact contribution of a given term in the threshold expansion of the amplitude to the partial waves. This includes infinitely many 1 J corrections to the results of 6.1 and 6.2. Note that in this section we only concern ourselves with the threshold expansion of the amplitude close to the two-particle threshold. Of course, in addition there are contributions from multi-particle thresholds. These lead to further nonperturbative in 1 J corrections to the results of this section and are potentially important if we would like to discuss partial waves at finite J which is the subject of the next section.

Large Spin Expansion of f J (s)
We start by analyzing the Froissart-Gribov integral (2.53) in the large spin limit J 1. For convinience, we quoting the integral here The large J behavior of Q 3) 24 In d = 4 this was pioneered by Dragt [64], see also [65], and further developed by A. W. Martin in [66]. Here we consider a slightly more general form of the amplitude, as well as generalize the analysis to any spacetime dimension. 25 This is analogous to the analytic bootstrap in CFTs [67][68][69]. 26 An efficient way to systematically expand Q (d) J (z) to an arbitrary order in the large J, fixed z expansion is to start from its representation in terms of the hypergeometric function [70] Q (d)  J (z) decays very fast for fixed z > 1 at large J. The integral is therefore dominated by the region close to the threshold z 1 . Explicitly, for J 1 we have and therefore the integral in (6.1) is controlled by the region of the size ∼ 1 J close to z * = z 1 . In that way, the large J behavior is controlled by the threshold expansion of T t (s, t(z)) established in the previous section in a manifest way.
Next, we explicitly plug the threshold expansion into (6.1) and compute the leading large J behavior of the partial wave coefficients. For that aim, it is convenient to express the integral as follows where we have 27 Note that in deriving (6.6) we took into account the integration measure (z 2 − 1) d−4 2 and the Jacobian from switching the integration variable to δz.
After performing the crossing transformation s ↔ t, we can now plug the leading threshold expansion expressions from the previous section (5.10) to get the J 1 limit of f J (s). Let us start with the universal and most singular case (5.7). For this case we get If, on the other hand, we consider the situation with a finite spin zero scattering length (5.8), we getf It is clear that we can systematically include 1 J correction in this expansion. Note that in all cases, the leading large J result for f J (s) is exponentially small and purely real. In the next section we compute the leading J contribution to Imf J (s). 27 Note that δz in (6.7) is different from the one used in the previous section.

Large Spin Expansion of Imf J (s)
We now repeat the analogues large J expansion for Imf J (s) in terms of the threshold expansion of the double spectral density. In this case, the Froissart-Gribov formula can be written as It is instructive to separate s into three regions, 4m 2 < s < 16m 2 , 64 3 m 2 > s > 16m 2 and s > 64 3 m 2 .
The Elastic Region 4m 2 < s < 16m 2 In this elastic strip the double spectral density is zero for t < 16m 2 s s−4m 2 or, equivalently, for z < 2z 2 1 −1, see figure 4. The large J limit of Imf J (s) is thus dominated by the expansion of double spectral density close to the leading Landau curve which we worked out in section 5.3. Hence, we can simply use (6.5) with z * = 2z 2 1 − 1, otherwise the consideration is identical to the one in the previous subsection. Explicitly, we have (1 + O(1/J)) . (6.14) An alternative way to arrive at (6.13) and (6.14) is to start from the large J expansion of f J (s) discussed in the previous section (6.8) and use elastic unitarity. Indeed, one can check that (6.8) and (6.13), (6.14) are consistent with elastic unitarity log 2Imf J (s) The Inelastic Region 16m 2 < s < 64 3 m 2 This regime is in the multi-particle part of figure 4. To analyze it one should first identify the leading Landau curve in that segment, which is beyond the scope of this paper.
The Inelastic Region s > 64 3 m 2 As s increases passed 64m 2 3 , we know from crossing that the leading Landau curve is the one of the dual channel (the blue curve in figure 4). A remarkable consequence of this fact is that crossing allows us to measure the non-zero inelasticity at large spin J.
As before the leading large J behavior of Imf J (s) is controlled by the leading Landau curve in the relevant kinematics. For s > 64 3 m 2 the leading Landau curve is at t = t 1 (s) = 4m 2 s s−16m 2 , or equivalently at See blue curve on figure 4.
In this case we get that the Froissart-Gribov integral takes the same form as in (6.10) with where t(δz, J) is a definition of δz that we will use below.

(6.18)
In odd d the situation is more complicated due to the logarithmic nature of the threshold singularity. While the leading power dependence on J can be easily computed, we do not give an explicit form for the dependence on log J that multiplies the leading power. 29 For the case of finite spin zero scattering length (5.8), we get in any d We can use the results above for the leading large J behavior of Imf J (s) for s > 64 3 m 2 to estimate the amount of inelasticity that exists at given s and J. Recall that if the scattering were purely elastic partial waves would satisfy 2Imf J (s) = (s−4m 2 ) d− 3 2 √ s |f J (s)| 2 , see (6.15). 28 By taking s → ∞ with x = m 3 s 3/2 J kept fixed in (6.11), (6.17) one gets the d-dimensional version of the Haan-Mütter scaling law [71,72]. 29 It should be possible to use the integral form for ρ(s, t) given in appendix H.2 to efficiently evaluate the integral numerically for given J but we do not pursue it here.
Using the explicit results of this section we conclude that elastic unitarity together with crossing lead to the universal inelasticity ratio at large J In figure 10 we plotted r J (s) as a function of s/m 2 for fixed J. We see that r J (s) approaches 0 for s = 20m 2 and s = ∞. It acquires its maximal value of about 0.47J at s = 40m 2 .

Threshold Expansion in the J-space
In this section we evaluate exactly the contribution of a given term in the threshold expansion to the partial waves. This generalizes the results above which correspond to the large J limit of the exact formulas that we present in the current section. For simplicity let us first focus on d even, when the threshold behavior is of the square-root type.
Consider the threshold expansion of the discontinuity of the amplitude. In the regime where it converges close to the threshold we can write In practice we truncate this expansion at any desired order and plug it into the Froissart-Gribov formula. Note that the Froissart-Gribov integral goes all the way to infinite z where the expansion (6.21) is no longer valid. However, that regime of integration is exponentially suppressed in J and therefore only affects the nonperturbative corrections at large J. It turns out that the relevant integral can be done exactly. It is given by we see that higher orders in the sum are more suppressed in J. Hence, one can use this representation to explicitly and systematically obtain all the coefficients in the large J expansion of f J (s), the leading 1/J result of course coincides with the previous analysis.
Similarly, we can write down the threshold expansion formula for ρ(s, t). Based on the previous discussion we have ρ(s, t(z)) = 1 The problem of finding d m (s) becomes completely algebraic after we note that (6.25) implies where we again performed the Froissart-Gribov integral exactly.
Elastic unitarity (6.15) relates the two expansions, (6.23) and (6.26). It then allows one to express the d m (s)'s in terms of the c n (s)'s. This is possible because the d m (s)'s do not depend on J. Equivalently, we can use (5.11) to directly map the threshold expansion of T t (s, t) to that of ρ(s, t) near the Landau curve. The details of this expansion are presented in appendix G.

Finite J and Finite s
In practice, one is interested in making statements about partial waves at some finite (but potentially large) spin J and finite energy s. Our ability to make such statements crucially depends on our ability to estimate an error in the Froissart-Gribov integral produced by the approximation to the amplitude. To that extent we can think of the discontinuity of the amplitude as follows is an approximation to T t (s, t) based on our knowledge about it. It can involve expansion coefficients close to various normal thresholds, information about resonances, or about the Regge limit. The more information is available to us, the less we can make T error t (s, t) and therefore the better is our knowledge of f J (s).
Assuming continuity of the amplitude we can try to bound the amplitude as follows. Consider a given s and let us assume that T t (s, t) ∼ t J 0 (s) at large t. Let us consider an integer N > J 0 (s)+ d−3 Continuity of the amplitude then implies that there exists c N (s) such that Indeed, this bound matches the neglected term close to the threshold and is also consistent with the Regge behavior since N > J 0 (s) + d−3 2 . The minimal value of c N (s) depends on the behavior of the amplitude at intermediate s. In particular if at some fixed t = t 0 the discontinuity develops a "bump", see figure 11, then c N (s) should be made large enough for (7.3) to hold. A familiar example of such a bump is a resonance and it is due to a singularity on the second sheet. More generally, we can call such a bump an outlier. Correspondingly, we call the problem of bounding c N (s) "the problem of outliers." 30 Note that for the integral on (7.3) to converge at large z, we need to take J > N + 1. Hence, a better choice of variables to encode our knowledge about the threshold behavior of T t (s, t) are ones that do not grow at large z. For example, we can replace (z 1 − 1) → (z − 1) in (7.2). Doing so is necessary if one wants to go to higher orders in the threshold expansion, but we will not pursue this in the present paper. The advantage of using just (z − z 1 )/(z 1 − 1) is that the Froissart-Gribov integral is known explicitly, see (6.22).
In a related manner, when deriving the Froissart-Gribov formula we do not have to close the contour all the way to infinity. We can instead keep arcs in the complex plane at some finite energy, see figure 2.b. This is sometimes called the truncated Froissart-Gribov formula, see e.g. [72,74]. The problem of deriving finite s and J formulae then requires bounding the contribution of the arcs. The advantage of this approach is that we do not have to assume extended analyticity all the way to t = ∞.
Let us proceed with c N (s) in (7.3) being our phenomenological parameter and see how various 30 For a similar discussion in the context of the CFT bootstrap see [73]. quantities depend on it. First of all, we can immediately bound the error in f J (s). We have n,J (z) is the integral (6.22).
Next, we consider Imf J (s). The relevant regime to study is s > 64m 2 3 , where we can estimate inelasticity using our knowledge about the threshold expansion in the crossed channel. In this case we can start from J (z)ρ(s, t(z)) , (7.5) and split the z-integral in two regions, 4m 2 < t < 16m 2 and t ≥ 16m 2 .
In the region 4m 2 < t < 16m 2 we can use the Mandelstam equation to compute ρ(s, t). It will involve the terms in the threshold expansion that we worked out above together with an error term that is controlled by (7.2) and (7.3). For t ≥ 16m 2 , since Disc s T approx t (s, t) = 0, we have Together with (7.2) and (7.3), this can be used to estimate the bound on Imf J (s).
If one is not willing to make an assumption of the form (7.3) about T error t (s, t) then one can still derive a bound on inelasticity [74], albeit very weak and only asymptotic when s → ∞.

A Toy Model Example
We now study a toy model example that is motivated by the numerical analysis of [15]. We consider four-dimensional spacetime (d = 4) and take T approx t (s, t) to be given by the first universal term in the threshold expansion (5.7) We assume that the Regge limit is bounded by J 0 (s) ≤ 1 2 for real s in some region. Hence, for this case we can take N = 1 in (7.3) and bound the error as where for convinience we have introduced the notation δc(s) = c 1 (s)/(32π), (7.3).
We will now go through the steps in figure 8 and apply them to the toy model (7.7). We will use the bound on the error (7.8) to have a finite J and finite s bound on the error at each step.
We start with the error on the partial waves coefficients. Using the Froissart-Gribov formula, it is given by (7.4) Note that due to our assumption on the Regge trajectory, the Froissart-Gribov formula is applicable all the way down to J ≥ 2 for any s ≥ 4m 2 .
To quantify the error in (7.9) we take the ratio between the error integral I The apprximation for the partial waves is good when this ratio is small. For example, we can focus on the point s = 40m 2 , where the large J inelasticity ratio (6.20) is maximal. In figure 12 we have plotted this ratio for s = 40m 2 as a function of spin. For example, we have err J=4 (40m 2 ) 0.51 δc(40m 2 ) , err J=20 (40m 2 ) 0.08 δc(40m 2 ) . (7.11) Next, we use the Mandelstam equation (3.9) to evaluate the double spectral density in the elastic strip and bound the error there. We get (see appendix G.1 for more details) n 1 ,n 2 was defined in (5.13), and recall that s cross 1,0 + δc 2 (s)Ĩ In estimating the error we used the fact that both the Mandelstam kernel Disc z K(z, η 1 , η 2 ), and T approx t are non-negative.

(7.18)
Let us comment on the origin of the various terms. The error in f J (s) is given by (7.9). In the expression for λ 2J 1 2Imf J=20 (s = 40m 2 ) the term linear in δc el comes from the error in the elastic region, and similarly the term linear in δc(40m 2 ) comes from the integral over t > 16m 2 . The δc 2 el term comes from the error in ρ(s, t) in the elastic region.
To summarize, given a bound on our ignorance about T t (s, t), specified by δc(s) for the case at hand, we can explicitly derive a lower bound on the amount of inelasticity that is present in the toy model. This is of course assuming that the underlying amplitude satisfies elastic unitarity. It would be also interesting if the analysis above can be improved using more refined error estimates. For example, a better approach would be to bound the error distributionally, see section 8.3.1 below. 31 Note that for s = 40m 2 this region is inside the Mahoux-Martin positivity region (4.4).

Elastic Unitarity and Coupling Maximization
In [15] the numerical procedure outlined below in section 8.1 has been put forward and carried out for the problem of maximizing T ( 4m 2 3 , 4m 2 3 ). It was observed that the low spin partial waves f J (s), with J = 0 and J = 2, converge very well as a function of N max which characterizes the number of terms used to approximate the amplitude, see section 8.1 for details. Moreover, it was found that that f J=0,2 (s) tend to saturate elastic unitarity above s ≥ 4m 2 .
It is interesting to apply the finite J and s analysis described above to this case. We find that the function produced by the numerics fits the toy model approximation we considered above in section 7.1. In particular, the relevant bounds on the the discontinuity of the amplitude T t (s, t) come from the threshold behavior of the amplitude The convergence properties of δc(s) at general s > 4m 2 are less clear as they probe the unphysical region. An interesting possibility to avoid the question of their convergence is to add the bounds on (7.19) as extra conditions to the bootstrap algorithm. Note, that in theories that satisfy elastic unitarity δc(s) does not depend on s since in this case only f 0 (t) contributes. Therefore dependence of δc(s) can itself be used as a probe of elastic unitarity.
Nevertheless, we proceed and consider numerics for N max = 11 and J max = 36 with the leading threshold behavior fixed to the universal behavior dictated by elastic unitarity. In this case we find from the numerics [17] δc el = 0.67 , δc(40m 2 ) = 1.17 . (7.20) In figure 13 we plot the results for the inelasticity ratio r J (40m 2 ) (6.20) obtained in the toy model of the previous section with the parameters taken from the numerics (7.20).
As we discuss below, there are several ways to improve the numerical procedure so that the plots would agree at finite N max . We leave the exploration of these possibilities for future work [17].

Analytical Methods and Numerical Bootstrap
In this section 32 we discuss how some of the analytical methods described in the paper can be implemented in the numerical approach to the S-matrix bootstrap that has been put forward in [15]. We will report the actual results of the numerical explorations elsewhere [17].

Review of the Numerical Framework of [15]
Let us briefly review the setup of [15]. The basic idea is to write an ansatz for the expansion of the scattering amplitude which is linear in unknown real parameters α abc maps the complex s-plane minus the s-channel cut to the unit circle and the point s 0 to the origin. Here, the extra terms may be added to make some particular properties of the amplitude manifest. Their presence or absence depends on the particular problem at hand. Crossing symmetry is imposed by demanding that the coefficients α abc are permutation-invariant. Finally, the relation s + t + u = 4m 2 leads to a redundancy in the basis of coefficients that can be addressed systematically.
To approximate an amplitude using the ansatz (8.1), the sum is truncated such that Given a finite N max , unitary in the form is imposed over a finite grid of points and for spins that are truncated by some maximal value J ≤ J max (N max ). As shown in [15], remarkably unitarity in the form of (8.3) can be restated as a semidefiniteness condition as follows. We write for physical J and s wheref J (s) are kinematical objects and all the dynamical information is in the coefficients α. The condition (8.3) can be then rewritten as a semi-definitedness condition for the matrix 32 We thank Madalena Lemos for collaboration on this section.
At this point one can maximize numerically some quantity linear in the α-parameters by imposing unitarity in the form (8.5) over the chosen grip in s and for J ≤ J max (N max ). For example, in [15] the "coupling", T ( 4m 2 3 , 4m 2 3 ), is maximized. If a certain maximization task reliably saturates as a function of N max we stop the process and trivially extrapolate to N max = ∞ to get the actual bound on the space of physical S-matrices.

Why and What Should be Improved
The setup of [15] has several very important and desirable properties: • It is simple and practical. It is not too hard to implement, manipulate and obtain bounds.
• The space of functions (8.1) is complete inside the ρ unit circle. Hence, any function analytic inside the circle (except, maybe a finite number of isolated poles which can be added explicitly) can be expanded in that way.
• Crossing is trivialized and is satisfied exactly.
At the same time, there is some tension and potential issues in applying the procedure above for exploring physical amplitudes that we list below and comment on how one may improve on them: • Physical amplitudes exhibit a set of normal multi-particle thresholds at s > 4m 2 . Hence, the ρ-expansion of such a function is not expected to converge point by point for real s > 4m 2 (|ρ s | = 1). On the other hand, in the procedure above unitarity is imposed point-wise in exactly that regime. Thus, it is unclear if the space of functions that are probed in this procedure includes among them physical amplitudes, with multi-particle thresholds. Below we discuss two possible ways of addressing this point. One is to improve the ansatz and the other is to change the way in which unitarity is imposed.
• The ansatz (8.1) has restrictive behavior at s = ∞. In particular, the partial waves satisfy lim s→∞ S J (s) = 1. There is no reason to expect this to be a correct property of physical amplitudes. Hence, similar to multi-particle thresholds, the ρ-expansion is expected to be a poor approximation in this kinematical regime. One way to fix it is to explicitly add extra terms to the ansatz. Another way, which we discuss in more detail below, is to change the way in which unitarity is imposed.
• Elastic unitarity is not satisfied. This is manifest for any finite N max -elastic unitarity implies that ρ(s, t) = 0 in the elastic strips, below the first Landau curve. If one tries to impose it exactly on the truncated ansatz (8.1) then clearly the only solution is α abc = 0. One may still hope that elastic unitarity will emerge at large N max . However, without extra constraints, we see no reason for that to happen. In practice, there is no conclusive evidence that the functions that emerge at the boundary of the allowed space satisfy elastic unitarity within the numerical error.
Imposing elastic unitarity is hard for the simple reason that this condition is nonlinear in the unknown parameters α abc . Therefore, within the approach of [15], we can only hope to impose elastic unitarity-type constraints. Namely, constraints that go beyond (8.3), but still include physical amplitudes in the space of functions that satisfy them. This is important if eventually we want to explore the space of physical amplitude that in particular do satisfy (8.6).
Below, we elaborate on several ways in which elastic unitarity can be pursued within the numerical approach of [15]. Importantly, each of them is still linear in α abc and therefore possible to implement using the standard solvers.

How the Numerical Framework Can Be Improved
We now suggest a few ways of addressing the issues identified above within the framework of [15].

Distributional Unitarity
Here we address some of the issues identified above, namely the lack of point-wise convergence for real s, or |ρ s | = 1, and too restrictive behavior at s = ∞ by suggesting a different way in which unitarity can be imposed numerically.
The basic idea is that even though the truncated expansion of a function with multi-patricle thresholds does not converge at |ρ| = 1, it can still converge to it as a distribution. Indeed, this is what is expected for the functions that are polynomially bounded as we approach |ρ| = 1 or real Mandelstam variables s and t. A relevant mathematical result that addresses this question is known as Vladimirov's theorem [75]. It was recently discussed in detail in the context of the conformal bootstrap in [76] to which we refer the reader for technical details.
Let us consider negative t needed to compute S J (s) starting from the amplitude. For such t, it is known that the amplitude is polynomially bounded as |s| → ∞. It is also polynomially bounded for real s, namely as Ims → 0, see e.g. [77]. These facts imply that we can apply Vladimirov's theorem to partial waves S J (s) and we expect that partial waves of physical amplitudes computed using the ansatz (8.1) will converge for real s distributionally, namely after integrating S J (s) against some smooth test function g(s). Note that this includes also amplitudes that grow with s or t, which can be modeled by (1+ρt) n (1+ρs) n type terms. Expanding such terms around ρ t = ρ s = 0 leads to a series that converges distributionally for |ρ s | = 1 as explained in [76].
Therefore, if we treat the amplitude (and the partial wave coefficients) as a distribution, we can still approximate it by the ansatz (8.1). Given N max , however, this approximation is expected to be good only for test functions that do not resolve the local in s features of the amplitude/partial waves, or have support for very large s. In other words, given N max we can only hope to have a reliable approximation on average over big enough intervals of s, with the intervals becoming smaller as we increase N max .
Similarly, considering test functions localized at larger s requires taking larger N max . This is a consequence of the accumulation of the multi-particle thresholds at large s, as well as due to the potential growth of the amplitude as s or t become large. Both effects lead to poor distributional convergence of the ansatz (8.1) at large s, or, equivalently, close to ρ s = −1. By restricting the support of the test functions away from that region we can use the ansatz (8.1) to probe functions that have multi-particle thresholds and grow in the Regge limit.
Given a real, non-negative function g(s), unitarity on average takes the form This can be still restated as a semi-definitedness condition and therefore can be implemented numerically using the same methods. In practice, given N max and expected properties of the physical amplitude there is a set of test functions g(s, N max ) that we can use.
To understand what are the reasonable functions to be used let us note that on the boundary of the circle, truncation of the maximal power in ρ n is the same as truncating the Fourier harmonics of ρ = e iφ . Therefore if we truncate n ≤ N max we cannot hope to resolve the amplitudes on scale δφ < 2π Nmax . Through the map ρ(s), this translates into a statement about the s-plane. We leave the detailed discussion of the test functions g(s, N max ) for the future.
Distributional convergence puts the numerical S-matrix bootstrap algorithm of [15] on a much more solid mathematical ground. It justifies imposing distributional unitarity (8.7) for the truncated ansatz (8.1).
In practice, for certain types of problems it can happen that a more naive point-wise analysis of unitarity still leads to correct results. Based on the reasons explained above this is expected to work if the underlying amplitude does not exhibit multi-particle thresholds and does not grow in the Regge limit. This expectation agrees with the numerical results observed in [15].

Extended Basis
Another way of addresing the issue of multi-particle thresholds is to extend the ansatz (8.1). This ansatz makes the structure of the two-particle normal threshold manifest. An obvious extension of the basis which makes the structure of the multi-particle cuts manifest is to add to it any power of the functions Such an extension while making the structure of the multi-particle normal thresholds manifest still has the property that the double spectral density misses the regions where ρ(s, t) = 0 carved out in the (s, t)-plane by the Landau curves.
Ideally, one would like to write down an ansatz which does not only make maximal analyticity manifest but also has a correct structure of the Landau curves. Such functions are naturally generated in perturbation theory. We can use them to write down functions that have expected behavior in the elastic strip. Let us take φ 4 in d = 3 and consider the following diagram Box(s, t, u) = Box(s, t) + Box(s, u) + Box(t, u) , (8.9) .

(8.10)
This function has the property that it is crossing symmetric and has the zero double discontinuity in the expected region.
We can consider for example the following ansatz By choosing n > 2 we make sure that the normal threshold coming from ρ (n) s,t,u starts after 16m 2 and the correct analytic structure inside the elastic strip comes from the sum of the φ 4 diagrams. Hence, this ansatz has an advantage of having the built-in analytic structure consistent with elastic unitarity.
Similarly, other Landau curves can be manifestly incorporated by choosing different perturbative diagrams and dressing them by an appropriate ρ-ansatz. In this way we can hope to have an ansatz which has more structure of the actual scattering built in. At the same time, linearity in α abc as well as crossing are still manifest.
That said, without imposing extra constraints that result from elastic unitarity, there is no a priori reason for the numerics to turn on the α abc 's that are associated to sub-leading Landau curves. We turn to such conditions next.

Elastic-type Unitarity at Non-integer J
Another way of injecting constraints from elastic unitarity into the numerical bootstrap is by imposing unitarity for non-integer spins. As we explained in section 3.3, the elastic unitarity condition (8.6) can be analytically continued to complex J as long as ReJ > J 0 (s), where J 0 (s) is the leading Regge trajectory.
Let us consider real J > J 0 (s). Numerically, we can then impose the following elastic unitaritytype condition  While physical amplitudes will saturate this bound the condition above still includes them as a part of the solution however goes beyond the integer spin unitarity conditions. Importantly, it is still linear in α abc and therefore can be easily implemented numerically.
When imposed in that way, the elastic unitarity-type conditions (8.12) are very similar to the original conditions (8.3). Note that the Froissart-Gribov projection probes the regime of arbitrary large t where the numerics convergence is slower. On the other hand, for integer spin one normally uses the partial wave projection (2.34) that only probes physical t's.

Lower Bound on Inelasticity
As we discussed in section 7 an additional knowledge of the behavior of the discontinuity of the amplitude T t (s, t) leads to a more refined prediction about the amount of inelasticity at finite J and finite s. For example, we can put in some expectation about the Regge behavior and low energy data on the scattering length, as well as structure of bumps or resonances to get a relatively accurate estimate of T t (s, t). Similarly, it can happen that while running the numerical algorithm one observes that the gross features of T t (s, t) saturate quickly as one increases N max . In this way one can get an accurate estimate of c N (s) in (7.3). This in turn allows us to put a lower bound on the amount of inelasticity r J (s) (6.20) for finite J and finite s.
The minimal amount of inelasticity can be easily implemented numerically. We can replace (8.5) by a set of modified matrices The corresponding modified positive semi-definetedness condition (8.13) is equivalent to the inequality where in the original problem a J (s) = 1 and more generally a J (s) specifies the minimal amount of inelasticity at given J and s.
From our discussion in section (7) and the large J analysis we know that a J (s) 1 for large enough J and s > 64m 2 3 , see (6.20). Moreover, given a local bound on the discontinuity of the scattering amplitude of the type (7.3), we can derive a set of improved elastic unitarity-type bounds (8.13) at finite J and finite s, see section 7. Luckily, a bound of the type (7.3) is again linear in α abc . Hence, the improved bound can be implemented numerically.

Mahoux-Martin Positivity
As we reviewed in section 4 elastic unitarity leads to the positivity property of the double spectral density in the so-called Mahoux-Martin region, see 4.4.
Positivity or more generally non-negativity of the double spectral density is obviously linear in α abc and therefore is straightforward to implement numerically. We can think of this either using the improved basis described in section 8.3.2 which directly implements the Steinmann shadow region where ρ(s, t) = 0. Alternatively, we can consider the original ansatz (8.1) and impose the Mahoux-Martin type positivity constraints The condition (8.15) still includes physical amplitudes, which in the elastic strip 4m 2 < s < 16m 2 satisfy the more restrictive conditions ρ(s, 4m 2 ≤ t < 16m 2 s s−4m 2 ) = 0 and ρ(s, 16m 2 s s−4m 2 < t ≤ 4m 2 (3s+4m 2 ) 2 (s−4m 2 ) 2 ) > 0.

Comments on CFTs
In this paper we assumed extended analyticity and we studied the structure of the amplitude for s, t > 0, where the amplitude develops crossing-symmetric double spectral density ρ(s, t). It is interesting to understand what are the analogous statements in CFTs.
The double discontinuity dDiscG and the quadruple discontinuity qDiscG were introduced in [69]. It was shown in [69] that crossing symmetry of qDiscG, see also [80,81], readily implies the presence of multi-twist operators in the OPE. This is the CFT analog of the Aks theorem reviewed above.
Elastic unitarity is a consistency condition of the two-particle sector of the S-matrix. In generic CFTs, the analog of two-particle states in the dual AdS space are double-twist operators [67,68] that are defined at large spin J. In large N CFTs the two-particle states in AdS correspond to double trace operators. Similarly, in AdS QFTs, see e.g. [82], we expect a natural set of operators corresponding to two-particle states to be present in the spectrum. However unitarity, as formulated in the CFT language through the OPE, does not admit a truncation analogous to elastic unitarity that emerges as we take the flat space limit of the theory.
Imposing that the twist spectrum structure of a CFT is the one coming from the light-cone bootstrap leads to the so-called Polyakov conditions, see e.g. [83] for a recent discussion in the nonperturbative context. It is interesting to understand to what extent the consequences of imposing the Polyakov conditions in AdS are analogous to elastic unitarity in flat space. 33 Using this analogy, the exploration of the present paper suggests that an interplay between the Polyakov conditions and crossing symmetry of the quadruple discontinuity can lead to interesting results. It will be interesting to investigate it further.
Let us next comment on extended analyticity. In the context of amplitudes it implies in particular that the analytic structure of the discontinuity T s (s, t) is similar to the one of the scattering amplitudes, modulo interchanging normal thresholds to Landau curves. In the context of CFTs it would require understanding analytic properties of dDiscG. In general this is a complicated problem since it requires going to the region of the (u, v) space in which no OPE channel converges, see e.g. [84] for a detailed, recent discussion. It is however possible to make progress in 2d CFTs. Indeed, in this case thanks to the Virasoro symmetry [85], the OPE converges on an arbitrary sheet [86]. One finds that a statement analogous to Mandelstam analyticity indeed holds, namely the only singularities on an arbitrary sheet of G(u, v) are branch points at u, v = 0, ∞. One does not expect an analogous statement in higher dimensional CFTs due to a more complicated structure of Lorentzian singularities of the correlator. However, it would be very interesting to investigate this analytic structure in more detail.
Finally, let us comment on the validity of the Mandelstam representation in CFTs. Recall, that to obtain Mandelstam representation in flat space we start with the usual dispersion relation (we ignore subtractions for simplicity) (9.1) 33 We thank S. Caron-Huot and J. Penedones for discussion on this point.
We then write the dispersion relation for the discontinuity of T s (s , t) (9.2) and plug in the formula above to get An important ingredient in this argument, apart from maximal analyticity, is polynomial boundedness of T s (s, t) for arbitrary s. 34 Let us now see what is the analogous situation in CFTs. Let us first consider 2d CFTs where maximal analyticity follows from the Virasoro symmetry as described above. Let us recall what were the main ingredients in the derivation of CFT dispersion relations in [79]. There it was shown that given a single-valued G(u, v) analytic in the cut-plane and bounded in the Regge (and Euclidean OPE) limit one can write a dispersion relation. What happens if as above we try to write dispersion relations for dDiscG(u, v)? Using the OPE, as described in [86], one can clearly bound any limit of the correlator or dDiscG(u, v) on any sheet.
Single-valuedness of dDiscG(u, v) is however not obvious using the Virasoro block OPE [85] and in general we do not expect it to hold. It is easy to check explicitly what happens in the case of minimal models. For the critical Ising model it is easy to see that dDiscG(u, v) is single-valued. It therefore satisfies all the necessary properties to write dispersion relations [79]. In this sense 2d Ising model correlators (somewhat trivially) admit Mandelstam representation. On the other hand, already in the tricritical Ising case single-valuedness does not hold, so we cannot apply the dispersion relations of [79].
In higher dimensions the situation is much more complicated due to absence of Virasoro symmetry. Here again there is no reason to expect single-valuedness of dDiscG(u, v). On the other hand, single-valuedness is a true property of the double discontinuity in free field theories, which therefore also admit the CFT analog of Mandelstam representation. One can wonder if this property continues to hold for the theories with slightly broken higher spin symmetry, e.g. Chern-Simons vector models in d = 3.
It would be also interesting to understand if there exists some other, more sophisticated way to think about writing an analog of Mandelstam representation in CFTs. As a different direction, thinking about some other versions of dispersion relations, see e.g. [88], that do not rely on singlevaluedness of the underlying correlator might very well be useful in certain applications.

Conclusions
One of the challenges of the modern conformal bootstrap is to efficiently combine analytical insights with the numerical methods to corner and solve physical theories [80]. Analogously, in this paper we revisited analytical techniques for the nonperturbative S-matrix bootstrap. A natural next step for the S-matrix bootstrap program is to combine them with the existing [15] or future numerical methods to compute physical amplitudes.
Concretely, in this paper we studied the implications of elastic unitarity and extended analyticity for the relativistic, unitary, gapped S-matrix in d ≥ 3. Our goal was to develop the analytical methods to constrain the nonperturbative scattering amplitude, which can be further used in the numerical bootstrap approaches. 35 The analytic bootstrap was the subject of active investigation in the 60's. Most of our ideas and results, when restricted to d = 4, are contained in some form in the old works of Dragt [64] and A. W. Martin [66], as well as more recent work of Roy and A. Martin [74]. We believe however that there is some value in "re-discovering" these methods from the modern perspective and in pushing forward the current incarnation of the S-matrix bootstrap.
As usual, if one wants to do analytic computations in a nonperturbative setting, a small parameter is needed. For a nonperturbative S-matrix, there are two expansions in two small kinematical parameters -the threshold expansion in s−4m 2 4m 2 , and the large spin expansion in 1/J. These two are related via the Froissart-Gribov formula (2.53). Ones combined with elastic unitarity and crossing symmetry, these two expansions lead to the bootstrap scheme outline in figure 8. The upshot of this analysis is that one can start with the low energy, low spin data (the threshold expansion), and use it to bootstrap the amplitude away from this regime. We, however, do not restrict the low energy data. In this sense, the scheme is analogous to the analytic CFT bootstrap [67][68][69].
While the analytic bootstrap methods reveal important structural properties of the amplitude, by themselves, they are not strong enough to "solve" the problem. Correspondingly, the low energy data that enters the threshold expansion and the bound on Regge are taken here as an unconstrained input for the analytic bootstrap scheme. Currently, the only known way of constraining these parameters systematically is using the numerical bootstrap techniques [15] or experiment [61]. As we discussed in the present paper, the numerical methods should be improved by implementing the structure that originates from elastic unitarity and extended analyticity. Indeed, it was observed in the numerical studies that the putative amplitude functions that saturate bounds tend to saturate unitarity. That seems in tension with the Aks theorem of section 4. The problem with the latter is that it does not provide us with a finite energy lower bound on particle production that can be implemented numerically. Provided a local bound on the discontinuity of the amplitude however, one does get a finite lower bound on particle production. Moreover, it can be implemented numerically as an extra constraint. Hence, it is instructive to consider the S-matrix bootstrap in a given class of discontinuity bounded amplitudes. We discussed this in more detail in section 7, where an explicit example is also given. In section 8 we suggest various ways in which inelasticity and other constraints that emerge from elastic unitarity can be implemented.
A related important question is to what extent physics in the elastic regions 4m 2 < (s or t) < 16m 2 studied in the present paper dominates the dynamics of the amplitude? In other words, under which conditions our ignorance of the multi-particle kinematics at s, t > 16m 2 leads to a small controllable error? As we have seen in section 7, when considering the toy model, in the case where the low energy interaction is strong (infinite scattering length) and the Regge behavior is relatively soft, the elastic region strongly constrains the behavior of partial waves at finite s and J. We can easily imagine a different situation, e.g. relevant for pion scattering, when the low energy interaction is weak. Based on our analysis in this case we do not have a reason to expect the physics of the partial waves to be dominated by the elastic region (unless the spin is very large). Correspondingly, in this case the dynamics in the multi-particle region is expected to be important. Bootstrapping such an S-matrix would then potentially require a more detailed understanding of 35 Many remarkable structures were recently unraveled in the study of perturbative scattering amplitudes of both massless [89] and massive [90] particles. In this paper we have focused on nonperturbative aspects of the two-to-two massive particle scattering. It would be interesting to see if any of these new insights can be put to use in the nonperturbative setting.
the analytic constraints that result from the physics in the multi-particle region.
Let us briefly discuss a few future directions. First, it would be very interesting to extend the current numerical approaches to the S-matrix bootstrap by implementing structures that originate from elastic unitarity in one of the ways suggested in this paper. We will report on this in [17]. Second, most of the explorations in this paper were bounded to the elastic strips, where one of the Mandelstam invariants is between the two-and the four-particle threshold energies, see figure 4. This region is particularly manageable because in one of the channels it is controlled by two-totwo amplitude only. It is an interesting and important task to explore the multi-particle region, where the energy is above the four-particle threshold in two of the channels. Finally, it would be interesting to explore the landscape S-matrices, other than d = 4 QCD. Ideally, one would like to find an S-matrix in d ≥ 3 that may play the analogous rule to the one played by the Ising model in the conformal bootstrap. Whether such "bootstrap-solvable" S-matrices in d ≥ 3 exist or not is yet to be shown. If it exists, we expect its solution to teach us a lot about nonperturbative QFT in general. Implementing efficiently the structure of the amplitude that we discuss in the present paper would be an important step towards constructing an example. A natural candidate theory to explore is φ 4 theory in d = 3.
There are also a few technical avenues along which our work can be extended. One is relaxing the Z 2 symmetry we assumed, which restricted the spin and the number of particles to be even. Another related extension is to include single-particle poles. Doing so will affect many of the details, but will not change the global picture. A more interesting generalization is to consider particles with spin, see e.g. [91]. Similarly, it is an open problem to implement the known structure of the UV of the theory, say asymptotic freedom or the CFT data of the UV fixed point, into the S-matrix bootstrap.
Finally, one can wonder if there is anything to be learned from this analysis for the conformal bootstrap. In the latter case, the theory is gapless so naively there is no elastic unitarity. However, CFTs in d > 2 have a twist gap, and multi-twist operators are mapped to the multi-particle states in the AdS dual theory. Therefore, it would be interesting to understand the AdS analog of the various aspects of the present paper more directly. the unit vector n In term of these coordinates, the angular integration in (2.20) reads 2 ) . The above formula is only true in d ≥ 4. In d = 3 we have which can be also obtained as a distributional limit from (A. where in the second step we have used (A.2) with d ≥ 4 and in the third we have used (A.3).

B Useful Identities for Gegenbauer P -and Q-functions
As discussed in the main text, due to the SO(1, d − 1) symmetry, the elastic unitarity kernel and the Mandelstam kernel are diagonal in spin. They take the form (2.36) and (3.6) correspondingly. In this appendix we derive these forms together with the related integrated expression (3.18).
The literature on the properties and identities of the Gegenbauer functions is extensive [93], starting with the Gegenbauer addition formula which dates back to 1893 to Gegenbauer himself [92]. Better suited for us is the integrated form of this identity [93] which in our conventions reads Perhaps the cleanest way to derive the above formula is to use group theoretic techniques [94]. If z 1 and z 2 are the cosines of the polar angles of unit vectors n 1 and n 2 and z is the cosine of the azimuthal angle difference between n 1 and n 2 , then z 1 z 2 + z 1 − z 2 1 1 − z 2 2 = n 1 · n 2 . One can then apply a rotation to make n 2 aligned along the z-axis, as the vector product is invariant under this transformation. The relation between P J (z 2 ) will then involve the SO(d − 1) matrix representation of this rotation, the Wigner D-matrix [95], which for a specific entry is given up to a factor by P (d) J (z 1 ). The integration in z will select this entry and project out the others, so we end up with a closed equation between Gegenbauer polynomials.
Let us change variable to with integration limits or equivalently, The Gegenbauer addition formula then becomes with P d given in (2.21).
Multiplying the above by n with n (d) J given by (2.29). We then arrive at (2.36). When z → 1, we can use (B.6) and the kernel localizes to We can get a similar identity to (2.36) for the Mandelstam kernel, by using the definition of the Gegenbauer function of the second kind, eq. (2.47), to get (3.6).
We can get a hint at the analytic structure of the Mandelstam kernel from (3.6). In terms of η and η the kernel shares the [−1, 1] branch cut of Q (d) J . The situation is more interesting in the z plane. Given that P (d) J (z) is a polynomial in z, analytic everywhere, K d (z, η , η ) can only be singular whenever z is such that the sum in J no longer converges. Indeed, when J → ∞ we have P which signals the singularity of the kernel as deduced in appendix A.
Representation (3.6) makes the symmetries of the kernel manifest. In particular, K d (z, η , η ) is symmetric in its last two arguments, and further obeys where we used Q Finally, let us derive (3.18). We start with integer J and |η 1 |, |η 2 | > 1. We take (B.5) and apply to it . In this way we get where we used the fact that for |η 1 |, |η 2 | > 1 the kernel is analytic in a finite region around [−1, 1]. We now want to blow up this contour to infinity to pick up the branch cut of the Mandelstam kernel (A.9). At infinity we have K d (η → ∞, η 1 , η 2 ) ∼ log η η for d > 3 and ∼ 1 η for d = 3. Given that Q J (η → ∞) ∼ η 3−d−J we get that the integrand goes like ∼ η −J−2 , which gives a null contribution to the arc at infinity for ReJ > −1. In this way we arrive at 1 πi Note that in the main text we included θ(η − η + ) in the discontinuity of the kernel, see (3.10). In the formula above, which is valid for complex (η 1 , η 2 ), the variable η is integrated from η + to ∞ and we can simply use (B.14) Plugging (B.13) and (B.12) back into (B.11) yields (3.18).
The above equation can be continued in spin. Indeed, note that both sides of (3.18) are manifestly analytic in spin J for Re[J] > −1 and coincide for positive integer J. To argue that they coincide for any J we also need to check the growth at infinity. One can check that both sides of (3.18) have a large J asymptotic λ(η + ) −J and therefore the conditions of Carlson's theorem are satisfied.
The large J expansion of the Q-function is given in (6.3) and (6.2). The aim of this appendix is to argue that there are no nonperturbative corrections to this expansion of the form λ(z) −αJ , with α > 1. This fact is used in section 3.5 to derive the Landau curves in the elastic region. We consider separately the cases when the spacetime dimension is even and when it is odd.
Odd dimensions In odd spacetime dimensions the Q-function (2.47) takes a simple form Where P n (J, x) is a polynomial of degree n in x and J. For example, we have This form makes the large J expansion trivial.
between the coefficients of the large J expansion of the Q-functions in even and odd spacetime dimensions. Now supposed Q (d even) had an λ(z) −αJ -type correction. Assuming no cancellations, such a correction would result in an analogous correction in the expansion of Q (d odd) (z 1 ) in the right hand side of (6.22). From the above however it is clear that corrections of this type are absent.

D Gribov's Theorem
It is possible to use elastic unitarity condition continued in spin J (3.15) to constrain the high-energy asymptotic of the discontinuity of the amplitude [96,97].
Consider the following ansatz for the discontinuity of the amplitude where B(s, log t) is a slowly varying function of t that grows slower than a power, e.g. B(s, log t) ∼ (log t) q(s) .
Gribov's Theorem: Let us assume the high energy behavior of the discontinuity (D.1) in the elastic region 4m 2 < s < 16m 2 . Elastic unitarity then implies that Historically, Gribov's theorem excluded the classical picture of diffraction from a black body T The easiest way to prove Gribov's theorem is to note that if the integral (D. 2) diverges f J (s) develops a singularity on the real axis at J = α(s). Taking J = α(s) + where 0 < 1 and real, we get from the Froissart-Gribov formula Elastic unitarity close to the leading Regge singularity J = α(s) then takes the schematic form which can be only consistent if (D.2) holds. Indeed, otherwise we get that the singularity in the RHS of (D.4) does not match the singularity in the LHS of (D.4).
A simple, and physically natural, way out of the contradiction is to assume that α(s) ∈ C. This corresponds to the discontinuity of the amplitude that takes the form T t (s, t) = β(s)t α(s) .
Let us now impose elastic unitarity (3.15) at J = α * (s). We get that the solution is Imα(s) . (D. 6) In d > 3 the consideration above tacitly assumed that Reα(s) > −1. This is related to the fact that in d > 3, as can be seen explicitly from (2.47), Q It would be interesting to understand better properties of the full nonperturbative leading Regge trajectory in the complex s plane, see e.g. discussion in [38] for some common assumptions. The properties of the leading Regge trajectory in the planar theory are relatively well-understood, see e.g. [98,99].
We now take R to just be given by (E.4) as a seed to (E.5) and iterate this equation to find the minimal consistent set of Landau curves in (3.24) that is consistent with this equation.
Note that for s > 4m 2 , R can have additional discontinuities in addition to normal thresholds. Hence, at each step we correct R with the additional thresholds of ρ that are generated through the iteration of (E.5). In the end we check the iteration procedure converges to a closed set of Landau curves.

F Threshold Expansion for Non-Integer J
It is interesting to ask about the continuation of the formula (5.3) for the solution to elastic unitarity to non-integer J. For a related discussion, see e.g. [100]. for |z| > 1 as can be easily seen from (2.47).
Switching in the RHS of (F.1) to the integral ∞ 4m 2 dt we can continue f J (s) to s < 4m 2 . Together with the fact that T t (s, t(z)) is positive and real for 0 < s < 4m 2 we conclude that f J (s) (s−4m 2 ) J is real analytic function of s with a branch cut starting at s = 4m 2 . Moreover, f J (s) (s−4m 2 ) J is positive and real for 0 < s < 4m 2 .
Let us now impose continued in spin elastic unitarity (3.15). It is convenient to rewrite it as follows The general solution to it takes the form Few comments are in order. Let us first discuss how the formula above reduces to (5.3) when J is even integer. In even d it is trivial upon identifyingb J=2k (s) = (s − 4m 2 ) 2k b 2k (s). In odd d the situation is more subtle because in this case e −iπ(J+ d−3 Secondly, note that when |J| → ∞ the ratio e −iπ(J+ d−3

)
sin π(J+ d−3 2 ) is polynomially bounded which is consistent with the expected behavior of J J (s) at infinity. This poses a potential problem in even d for half-integer J and in odd J integer J.
Let us consider even d first. In this case e −iπ(J+ d−3

)
sin π(J+ d−3 2 ) = −i e −iπJ cos πJ develops a pole at odd J which corresponds to zero of f J (s). Therefore unless there is a cancellation between the two terms in the denominator of (5.3) we have that f 1 2 +Z = 0. If this is the case via Carlson theorem we then conclude that f J (s) = 0 identically. Therefore, an infinite number of poles should cancel with the corresponding poles inb J (s). In this way we get G Threshold expansion in J-space: Technical Details In this appendix we collect various results and technical details that are relevant to the inversion of the threshold expansion using the Froissart-Gribov formula performed in section 6.3.
We start with the derivation of (6.22). Consider first d to be even. We would like to evaluate the following integral We would like to do the integral for general J. The strategy is to do the integral first for integer J exactly and then analytically continue it to arbitrary J.
As a first step we note that the integrand can be interepreted as a discontinuity of some simple function For integer n and d this only holds for even d, where the power is half-integer and therefore we have a square-root type discontinuity.
Therefore if we interpret the integrand as a discontinuity T t (s, t) = Here we used that J is integer and dropped the contour at infinity which requires J > n − d−3 2 .
We note that the integral (G.1) satisfies a very simple recursion relation based on the identity ∂ z 1 I(d, n, z 1 ) = −(n − d − 3 2 )I(d, n − 1, z 1 ). (G.4) Therefore we can first compute the integral for n = 0 and then use this differential equation to compute the integral for n > 0.
Let us now find explicitly I(d, 0, z 1 ). To do this let us do the following change of variable 1 t = λ(z 1 ) ≡ z 1 + z 2 1 − 1. (G.5) We can then write where we used the relation between P (G.7) Using the orthogonality property of P The rest we can get trivially using (G.4). Note that the final result holds in any d and for any J.
The basic integral is the following Therefore, it is clear that we get the following result for the integral I(d, n, z 1 ) = π sin πd n k=0 c k,n λ(z 1 ) −2k , (G. 10) where c k,n can be explicitly found.
We have c 0,0 = 1. It is then easy to check that we have the following result I(d, n, z 1 ) = 2 Therefore to apply the formula for n = m we need to deform the contour across α = 0 pole and pick the residue. We then perform the α integration.
For given n this can be quite easily expanded at large J. Important property of this expansion is that higher terms in k have an extra suppression in 1 J . Using this formula one can in principle repeat the analysis of section 6 up to an arbitrary high order in z−z 1 z−1 without spoiling the Regge behavior of the amplitude.

H Keyhole Integrals in Odd d
In this appendix we collect some of the useful integrals in odd d. They appear both in the large J expansion of partial waves, and in the threshold expansion of the double spectral density. A key difference compared to even d is appearance of powers of both logarithm and inverse logarithm of the threshold expansion parameter σ t .

H.1 Partial Wave
In the discussion of the large J expansion we encountered the integral (6.6), and in odd d we introduced a function g n (J) in (6.8) that controls the large J asymptotic behavior of partial waves. Let us compute it explicitly. In our problem n = d−3 2 . The keyhole contour is depicted in figure 9 and it naturally appears when deriving the Froissart-Gribov formula.
For n = 0, 1 (d = 3 and d = 5) it is not necessary to keep the keyhole since the integral (H.1) converges and we can simply write g n=0,1 (log J) = ∞ 0 dz z n 1 log 2 z J + π 2 e −z .

(H.2)
To compute n ≥ 2 it is convenient to slightly modify the integral and use the following recursion relation where we used that e −a∂x g(x) = g(x − a).
By solving the recursion we then then get Together with the result of the previous subsection it allows us to compute the leading threshold contribution to the double spectral density.