The S-matrix Bootstrap IV: Multiple Amplitudes

We explore the space of consistent three-particle couplings in $\mathbb Z_2$-symmetric two-dimensional QFTs using two first-principles approaches. Our first approach relies solely on unitarity, analyticity and crossing symmetry of the two-to-two scattering amplitudes and extends the techniques of [arXiv:1607.06110] to a multi-amplitude setup. Our second approach is based on placing QFTs in AdS to get upper bounds on couplings with the numerical conformal bootstrap, and is a multi-correlator version of [arXiv:1607.06109]. The space of allowed couplings that we carve out is rich in features, some of which we can link to amplitudes in integrable theories with a $\mathbb Z_2$ symmetry, e.g., the three-state Potts and tricritical Ising field theories. Along a specific line our maximal coupling agrees with that of a new exact S-matrix that corresponds to an elliptic deformation of the supersymmetric Sine-Gordon model which preserves unitarity and solves the Yang-Baxter equation.

The bootstrap of the two-to-two S-matrix of the lightest particle in a relativistic unitarity quantum field theory was recently revived in [1][2][3] and extended to particles with flavour in [7][8][9][10]. These works can be seen as gapped counterparts of the conformal bootstrap explorations in [4,5] and [11] (without and with flavour respectively). Here we initiate the bootstrap analysis of S-matrix elements involving different external particles in Z 2 symmetric theories. This multiple amplitude study again mimics a similar recent development in the conformal bootstrap, namely the multiple correlator analysis of the Ising model which famously gave rise to the CFT islands in [12]. 1 We will consider two-dimensional QFTs with exactly two stable particles of masses m 1 and m 2 . We will assume the theory to be parity and time-reversal invariant and both particles to be parity even. For simplicity we will also postulate the existence of a Z 2 symmetry, under which the first particle is odd and the second particle is even. 2 This means that the nonzero three-particle couplings are g 112 and g 222 , which can be defined non-perturbatively in terms of the residues of a pole in a suitable S-matrix element. In the first part of this paper we will analyze all the two-to-two S-matrices of particles 1 and 2 and use crossing symmetry, analyticity and unitarity to explore the space of possible points in the (non-dimensionalized) (g 112 , g 222 ) plane as a function of m 2 /m 1 -see figure 8 on page 19 to get an idea. In order to avoid singularities or Coleman-Thun poles [6], which complicate the analytic structure of the scattering amplitudes, we will restrict ourselves to Note that we allow m 2 < m 1 also.
Under the stated assumptions there are five different physical two-to-two scattering processes as shown in figure 1. These can be grouped either according to the nature of their intermediate states, which can be Z 2 odd or even, or according to whether they are 'diagonal' or not. To wit, for a diagonal process the incoming and outgoing momenta are the same whereas for an off-diagonal process they are different. 3 As is also indicated in figure 1, we call the 12 → 12 diagonal process 'forward' scattering, and the 12 → 12 off-diagonal process 'backward'. 1 While walking hand in hand as illustrated in the previous paragraph, the S-matrix and the conformal bootstrap also have significant differences. In the CFT bootstrap we exclude theories; once excluded, a theory can never be accepted; with better computers we exclude more. In the S-matrix bootstrap of [2,3] we include theories; by constructing explicit solutions to crossing and unitarity some parameters are shown to be allowed; with better computers we include more. (In the S-matrix bootstrap when we impose new physical conditions we will exclude some of the previously found S-matrices though so the process is not as monotonic as in the CFT bootstrap.) The two bootstraps are thus two faces of a same coin, the Yin and the Yang, the darkness and the light, the chaos and the order. For a recent review of the conformal bootstrap state of the art see [13]. 2 In two dimensions theories with fermions and scalars are naturally Z 2 symmetric theories so the setup here applies as well to any theories with scalars and fermions, not necessarily supersymmetric. 3 As explained further below, in two dimensions the scattering angle can take only two values by kinematical restrictions; the outgoing momenta are essentially 'locked' in terms of the incoming momenta. Unlike in higher dimensions, there is therefore no (analytic) function interpolating between forward and backward scattering.  Figure 1: Diagonal processes are those where the incoming and outgoing particles have the same momenta as illustrated in the first row; they are all crossing invariant. The non-diagonal processes in the second row are those for which the final momenta are not the same as the initial momenta. Swapping space and time interchanges the odd and even off-diagonal processes so these off-diagonal processes play a crucial role in connecting these two sectors of different Z 2 charge.
In section 2 we will state in detail the conditions of unitarity, analyticity, and crossing symmetry that these five processes must obey. To guide ideas let us mention two conspicuous facts. First, we note that crossing symmetry flips the s and t axes on the diagram. This relates the two off-diagonal processes and thereby reduces the number of independent amplitudes (i.e. functions of the Mandelstam invariant) to four. Of course, it also imposes a non-trivial constraint on the amplitudes for the diagonal processes. Second, we observe that particle 1 can appear as an intermediate state in all the odd processes and gives rise to a pole in these amplitudes with residue proportional to g 2 112 , whereas particle 2 gives rise to poles in all the even amplitudes with residues proportional to g 2 112 , g 112 g 222 or g 2 222 , depending on the process. These poles can be thought of as our definition of the corresponding couplings. Analytic bound based on the scattering of the lightest odd particle, from [2]. Solid line: Analytic bound arising from the forward (or transmission) scattering of the odd particle against the even particle; it is a much stronger bound. Red dots: The numeric bound obtained from all two-to-two processes as discussed in the main text. The shaded regions represent the allowed regions which nicely shrink as we include more constraints. Any relativistic, unitary, Z 2 invariant theory theory with two stable particles (one odd with mass m 1 and one even with mass m 2 ) must lie inside the darkest blue region.

Quick comparison with single-correlator bounds
As a warm up exercise let us first discuss the three diagonal processes in isolation and explain how the methods discussed in [2,3,15,16] already lead to some constraints on the couplings. Analytic bound based on the scattering of the lightest even particle, from [2]. Red dots: The numeric bound obtained from all two-to-two processes as discussed in the main text. The shaded region represent the allowed region. When the even particle is the lightest, we can solve analytically for the maximal coupling, even considering the full set of amplitudes. When the odd particle is the lightest, the coupling can be bigger, diverging when singularities of the amplitudes corresponding to physical processes collide. This happens at m 2 /m 1 = 2/ √ 3. After this mass ratio the upper bound disappears.
plane is pretty straightforward. For example, the odd process has a two-particle s-channel cut starting at s = (m 1 + m 2 ) 2 and a pole at m 2 1 with residue proportional to g 2 112 , plus the crossed t-channel singularities obtained by swapping s → 2m 2 1 + 2m 2 2 − s. The even processes S 11→11 and S 22→22 have their two-particle s-channel cuts starting at min(4m 2 1 , 4m 2 2 ) and a pole with residue g 2 112 or g 2 222 , again plus the crossed t-channel singularities obtained by swapping s → 4m 2 1 − s or s → 4m 2 2 − s. As for unitarity, notice that the discontinuity across the cut is always positive, but it is bounded from above only for physical s, which means s > 4m 2 1 for S 11→11 and s > 4m 2 2 for S 22→22 . Therefore only for the lightest of the two particles is the discontinuity everywhere bounded from above, whereas for the other particle the discontinuity can be arbitrarily large (but not negative) in the interval between min(4m 2 1 , 4m 2 2 ) and max(4m 2 1 , 4m 2 2 ). We can bound the couplings as follows. First let us bound g 2 112 by using the maximum modulus principle for S forward 12→12 following [2,3,15,16]. We define with h ab (s) ≡ (s − (m a − m b ) 2 )((m a + m b ) 2 − s). The function f 12→12 is free of singularities (since we divided out by functions with poles at the pole location of the amplitudes) and is bounded at the s-and t-channel cuts (since the functions we divided by are phases at those cuts and the amplitude is bounded). Therefore f 12→12 (s) must have absolute value smaller or equal to 1 everywhere, and in particular at m 2 2 and m 2 1 where we can simply read off the maximally allowed couplings in these amplitudes. This leads to a universal upper bound on g 2 112 , which is the solid line in figure 2. The exact same analysis can be used for the elastic amplitude for the lightest of the two particles. If we denote this by , so m = min(m 1 , m 2 ), then the maximum modulus principle for gives a bound on the coupling appearing on S → , which is g 2 2 . This is the dashed line for m 2 > m 1 in figure 2 and the solid line for m 2 < m 1 in figure 3.
Finally we can use the techniques of [2] to also derive a bound on g 222 from the amplitude S 22→22 even when m 2 is not the lightest particle. In this case there is a cut which is not bounded directly by unitarity as depicted in figure 4. As we derive in appendix C, the amplitude with maximal g 2 222 is given by The corresponding bound on g 2 222 is plotted as the solid line in figure 3 for m 2 > m 1 . As the figure shows, the bound actually disappears for m 2 ≥ 2 √ 3 m 1 , which is due the t-channel pole colliding with the s-channel cut in the 22 → 22 process at this mass ratio. This is the simplest instance of a phenomenon we call screening. It is detailed in figure 4 and we will encounter it again below. In the same way we could obtain a bound on g 2 112 from the 11 → 11 process even when m 1 is the heaviest particle. This bound corresponds to the dashed line in figure 2 for µ < 1, and is always less restrictive than the bound from S forward 12→12 . This concludes our discussion of the single-amplitude results. As a preview for the more detailed numerical results presented below, we already marked in figures 2 and 3 in red dots our best numerical values of the coupling obtained from a simultaneous analysis of the full set of two-to-two amplitudes depicted in figure 1. Figure 2 displays a clear improvement over the quick single-amplitude analysis for m 1 / √ 2 < m 2 < √ 2m 1 , with an intriguing kink at m 2 = m 1 . It would be fascinating to find if this kink corresponds to a physical theory. On s 4m 2 2 4m 2 1 two-particle s-channel cut below the physical threshold for this process unitarity/physical region 4m 2 2 m 2 2 = 3m 2 2 (t-channel) single particle pole collides with the (s-channel) two-particle cut (and thus disappears) as m 2 /m 1 & 1.155 Figure 4: Analytic structure of the S 22→22 amplitude (for clarity we do not show the left cut and s-channel pole following from crossing symmetry S 22→22 (s) = S 22→22 (4m 2 2 − s)). If m 2 is not the lightest particle, there is a new feature in the S 22→22 amplitude: a two particle cut starting at s = 4m 2 1 corresponding to the contribution of two particles m 1 . This cut appears before the cut for two particles m 2 at physical energies s ≥ 4m 2 2 where regular unitarity is imposed and the amplitude needs to be bounded. As m 2 grows beyond 2/ √ 3m 1 the t-channel pole corresponding to the exchange of particle m 2 enters the new cut (by crossing symmetry the s-channel pole enters the t-channel cut) so we "lose" this pole. Beyond this point we can no longer bound g 222 since it does not appear in any other diagonal amplitude. This is indeed what we observe in the numerics as illustrated in figure 3. Note that before the bound on g 222 disappears it diverges. This divergence, arising from the collision of the t-channel pole with an s-channel cut is analogous to the divergences in bounds on couplings when s-and t-channel poles collide as already observed in [2]; the dashed line in figure 2 which was taken from [2] diverges at m 2 = √ 2m 1 for exactly this reason.
the other hand, in figure 3 we see no improvement over the single-amplitude results. In fact, in section 3 we will prove that the maximal value of g 222 in the multi-amplitude analysis saturates the single-amplitude analytic bounds just derived. In the same section we will show a more complete picture by considering the entire (g 112 , g 222 ) plane.

QFT in AdS
As shown in [1], there exists a completely orthogonal approach towards the problem of determining the maximal couplings in QFT. Rather than working from the S-matrix, which required analyticity assumptions that in general dimension D are not very well understood, the idea is to consider QFTs on an AdS background. The boundary correlators of such a QFT, which are defined in a similar way as in the AdS/CFT correspondence, behave much like conformal correlation functions in one lower dimension d = D−1. By applying numerical conformal bootstrap methods of [4] one can put a universal upper bound on the three-point couplings of QFTs in AdS. One can then extrapolate this bound to the flat-space limit (by sending all scaling dimensions to infinity), resulting in putative bounds for flat-space QFTs.
In [1,2] this was shown to work extremely well for two-dimensional QFTs: a precise match was found between the single-correlator analysis using the conformal bootstrap, and the single-amplitude analysis that we partially reviewed above.
In this work we continue these explorations. As discussed further in section 4, the Z 2 symmetric setup that we consider is easily translated to a multi-correlator conformal bootstrap problem for QFTs in AdS. In most cases we again find a very good match, and in particular we are able to recover the coupling of the 3-state Potts field theory from the conformal crossing equations. For large-ish mass ratios, however, we will see that the multi-correlator bootstrap appears to be less powerful than even the single-amplitude bootstrap.

Outline
The Z 2 symmetric S-matrix bootstrap is fully spelled out in section 2 and analysed numerically in section 3 leading to various bounds on the allowed coupling space for various mass ratios as illustrated in figure 8. In section 3.3 we discuss integrable Z 2 symmetric theories with m 2 = m 1 and how some of them nicely show up at the boundary of the allowed S-matrix space found in the numerical bootstrap. These include a massive deformation of the 3-state Potts model, the super-symmetric Sine-Gordon model and a SUSY breaking integrable elliptic deformation of the super-symmetric Sine-Gordon which seems to be novel as far as we know. Section 4 contains the results from the QFT in AdS analysis. Various appendices complement the main text with further extensions. (For example, the special role of the Tricritical Ising model as a kink in the space of S-matrices is discussed in appendix H.) 2 Multiple amplitudes 2.1 Kinematics of the various Z 2 preserving processes There are six two-to-two processes involving particles m 1 (odd) and m 2 (even) in a two dimensional Z 2 symmetric theory. We also assume time-reversal and parity symmetry. Four of those six are even processes where we scatter either 11 or 22 into either 11 or 22. Of those four, two are trivially related by time-reversal, so we can ignore one of them (say 22 → 11) in what follows. The remaining two processes are Z 2 odd processes where we scatter the odd particle against the even particle obtaining those same two particles in the future. As explained in the introduction this process splits into two possibilities which we call the forward and the backward component, see figures 1 and 5.
In two dimensions, any process depends uniquely on the center of mass energy or equivalently on the Mandelstam invariant If the particles are distinguishable these are two genuinely different processes denoted as the forward or backward process. (They are sometimes also called the transmission and reflection processes.) In higher dimensions, these two scenarios are limiting values of the a single amplitude when the scattering angle tends to θ = 0 or θ = π, but in two dimensions there is no scattering angle and these processes are described by independent functions. As we exchange time and space, i.e. as we analytically continue these processes by swapping t and s we see that the forward process is mapped to itself while the backward process as seen from its crossed channel describes the m a m a → m b m b event. This translates into equations (19) and (20) in the main text.
This in particular means that the other two Mandelstam invariants are completely determined in terms of s. It is important to find the precise relation because crossing symmetry permutes the three Mandelstam invariants and therefore leads to symmetries of the amplitudes M (s) that we need to impose. In a process 5 involving m a m b → m c m d The first equal sign is the two dimensional constraint: in two dimensions p 3 is always a linear combination of the two-vectors p 1 and p 2 and hence the determinant vanishes. In the second equal sign we used the on-shell conditions and momentum conservation. For example c and so on. Evaluated explicitly and combined with the previous linear constraint on the Mandelstam invariants, this can be cast in a nice symmetric form: i . Let us now specialize to the Z 2 preserving cases mentioned above. For the simplest processes corresponding to all equal masses (i.e. for 11 → 11 and 22 → 22) the condition dramatically simplifies into stu = 0 which leads to u = 0 or t = 0 or s = 0. In fact, we can not set s = 0 since by definition we assume s to be constructed from two incoming particles and setting u = 0 or t = 0 is the same up to a simple relabelling of the final particles which we can always do for indistinguishable particles. Hence without loss of generality we can set u = 0 recovering the famous result that elastic scattering of identical particles in two dimensions has zero momentum transfer.
Next we have the processes involving two particles of mass m 1 and two particles of mass m 2 . Here it matters whether the two particles of the same mass are both incoming or if one is incoming and the other is outgoing. Let us start first with the second case so that we can set m a = m d = m 1 and m b = m c = m 2 in agreement with the conventions of figure 5. Then we obtain a nice factorization of the constraint (9) into with two clear solutions: u = 0 corresponding to forward scattering and u = 2m 2 1 + 2m 2 2 − (m 2 1 − m 2 2 ) 2 /s + s corresponding to the more complicated backward scattering. Note that in forward scattering the final momenta are equal to the initial momenta but this is not the case in backward scattering where the momentum transfer is non-zero as highlighted in figure 5.
Lastly we have the even process 11 → 22 where m a = m b = m 1 and m c = m d = m 2 which of course corresponds to a simple relabelling of the previous constraint in which s ↔ u and thus leads, after discarding the s = 0 solution, to whose solutions are u = 1 2 (2m 2 1 + 2m 2 2 ± (4m 2 1 − s) (4m 2 2 − s) − s). In fact, these two solutions are equivalent up to relabelling of the two outgoing particles. Of course, the s ↔ u relation between 11 → 22 and backward 12 → 12 scattering is just crossing symmetry.
All in all we understood that all amplitudes can be thought of as functions of s with the other Mandelstam invariants given by M forward 12→12 (s) : M backward 12→12 (s) : M 11→22 (s) : The above equations allow us to state the crossing symmetry equations which we will impose in the sequel. They are: Note in particular that the last crossing relation plays quite an important role: it connects the even and the odd sectors.
For more on how the above discussion can be related to a similar analysis in higher dimensions see appendix A.

Analyticity, Unitarity and Extended Unitarity
The central hypothesis for the S-matrix bootstrap is that the scattering amplitudes are analytic for arbitrary complex values of s up to so-called Landau singularities [17] corresponding to on-shell intermediate processes. For the amplitudes and mass range discussed in this paper, these singularities in the physical sheet correspond to the possibility of the full a → b scattering process to factorise into two scatterings, first a → c and then c → b. Each on-shell state c of the theory will produce a singularity in the a → b process for s equal to the center of mass energy squared of the state c. This singularity will then proliferate according to its image under crossing transformations, see e.g. (17)(18)(19)(20). The discontinuities around these singularities are governed by the generalized unitarity equations [17], (where the first equality assumes time reversal invariance.) Equation (21) is very powerful and reduces to a number of familiar examples in special cases: • The contribution from one particle intermediate states corresponds to nothing but the usual bound-state poles: there the phase space integral reduces to the energy momentum delta function and the product of amplitudes to the physical three-point couplings, combining to the bound-state pole discontinuity −2πiδ(s − m 2 k )g 12k g 34k . • For real values of s for which there are no on-shell states, (21) reduces to the reality condition ImM 12→34 = 0.
• All of the above are very well known. Indeed, for the lightest two particle states in a given channel, there is nothing more to (21) than bound state poles, real analyticity and unitarity. For heavier external states, however, (21) extends the unitarity relation to the unphysical energy region s < max{(m 1 + m 2 ) 2 , (m 3 + m 4 ) 2 } by keeping the quadratic terms in the unitarity equation that correspond to physical intermediate states of energy √ s. This is what is called extended unitarity.
In our Z 2 symmetric setup and for √ s < min(3m 2 , 2m 1 + m 2 ), 6 we find and for √ s < min(3m 1 , 2m 2 + m 1 ), 2ImM Backward where the denominators come from the phase space factors and θ is the Heaviside step function. For example, if m 2 > m 1 then equation (24) for s > 4m 2 2 is just unitarity for the 22 → 22 process, but for 4m 2 1 < s < 4m 2 2 it is a "new" constraint over the |11 production cut.
Of course, the scattering amplitudes also have cuts and poles corresponding to crossed intermediate processes. The discontinuities around those singularities are governed by the generalised unitarity equations for the crossed scattering, together with the crossing equations (17)(18)(19)(20).
For energies above the three particle threshold, new terms corresponding to three-particle intermediate states should be introduced in the r.h.s. of equations (22)(23)(24)(25)(26) It is useful, however, to keep only the contributions from two-particle intermediate states and replace the full set of equations (22)(23)(24)(25)(26) by a positive semidefinite constraint on the amplitudes. For the Z 2 even sector, by dropping the contributions from intermediate states with three or more particles in (21), we can write in matrix form where takes into account the phase space volume. Note that (27) is saturated for √ s < min(3m 2 , 2m 1 + m 2 ). As discussed in section 3, for the numerical implementation we impose (27) even before multiparticle thresholds, leaving for the computer to achieve saturation where (22)(23)(24) applies. A similar discussion holds for the Z 2 odd sector.
In appendix B, we provide a direct derivation of (27) for √ s > 2 max(m 1 , m 2 ). This derivation elucidates the physical meaning of the matrix M and its relation to transition probabilitues between initial and final states.

Implementation
As discussed in section 2.2, the Z 2 symmetric scattering amplitudes in the mass range (1) are analytic functions in the physical sheet of the the kinematical variable s up to poles corresponding to bound states. This sheet is defined by continuing the amplitudes away from physical kinematics respecting the i prescription and has as its boundaries cuts corresponding to two and higher particle production thresholds which may happen in the s, t and u channels. These can be summarised by expressing the amplitudes through dispersion relations, as illustrated in figure 6. For the case m 1 < m 2 , we obtain with C a→b constant. Equations for the m 1 > m 2 case are obtained by replacing 4m 2 1 → 4m 2 2 in the lower limits of the integrals. Recall that these are the only indepen-  (15), as we move in the s plane, we hit poles and two-particle (as well as multiparticle) thresholds in the s and u channel, but not in the t-channel. This is a consequence of our arbitrary definition of t and u (in the language of appendix A, this comes from choosing for each s a single point in the two-valued hyperbolas of figure 15). To derive the dispersion (31) we start by assuming the amplitude approaches a constant at infinity (but see discussion in the main text) and write the identity where γ is the dotted contour above. We can then neglect the arcs at infinity. The contribution from the arcs around the red singularities correspond to the s-channel pole and s * integral in (31). After changing the integration variable in the remaining terms to u * (s * ) according to equation (15), we find the kernel transformation where we could relabel u * → t * in the second term. Then, after absorbing M ∞ into the constant C 11→22 and using the crossing relation (20) and the discontinuity formula (21) for the pole terms, we obtain the dispersion relation (31). In deriving this relations, see figure 6, we assumed that the scattering amplitudes have no essential singularities at infinity, and in fact approach a constant in this limit, i.e. the S-matrix becomes free. This latter assumption is not crucial nor required: it can be lifted by introducing subtractions as discussed in [2] and the numerical problem of maximising the couplings is not sensitive to this. This is to be expected physically, since the low energy physics of bound state poles should not be much sensitive to the behaviour of the amplitudes at high energies.
To obtain a concrete numerical implementation to the problem, we proceed as follows. First, we define a dispersion grid {x 1 , ..., x M } along the integration domains in (28)(29)(30)(31). We then approximate the discontinuities ImM a→b (x * ) by splines σ a→b (s) 7 linear in between the grid points up to a cutoff point x M , after which we assume the discontinuities decay as ImM a→b (x * ) ∼ 1/x * . 8 With this approximation we can analytically perform the integrals in 7 If m 1 < m 2 , extended unitarity, equations (22)(23)(24), allows for M 22→22 to diverge as 1/ s − 4m 2 1 close to the 4m 2 1 threshold. Due to this, between the first two grid points, we approximate ImM 22→22 ∝ 1/ s − 4m 2 1 . If m 1 > m 2 we should replace 1 ↔ 2 in this discussion. 8 This is similar to the numerical implementation in [2]. We could have parametrised our amplitudes using the ρ variables defined in [3]. These variables provide a cleaner framework for the numerics but, in practice, we find that convergence with the ρ variables is much slower than with the use of discretized dispersion relations. (28)(29)(30)(31) obtaining, in the case m 1 < m 2 and for M 11→11 , as an example, where the functions K i are defined in the appendix A of [2]. Next, we impose (27) along a fine grid over s > min{4m 2 1 , 4m 2 2 } (we impose analogous constraints over analogous ranges for the Z 2 odd channels). Note that we leave for the computer to achieve saturation of (27) before the three-particles thresholds. As shown in appendix I, equation (27)is equivalent to the semidefiniteness constraint and a similar rewriting can be done for the Z 2 odd sector. If we fix α = g 222 g 112 , as well as the masses, then the matrix in the l.h.s. of (33) is linear on the variables {C a→b , g 2 112 , σ a→b (x i )}. The problem of maximising g 2 112 in this space of variables under the positive semidefinite constraint (27) (and equivalent for the Z 2 odd sector) is therefore a semidefinite program and can be solved with, say, SDPB [14]. Details on the numerical implementation, such as parameter settings are available upon request.

Results for any m 2 /m 1
For each mass ratio m 2 /m 1 and for each coupling ratio α ≡ g 222 /g 112 we can now look for the maximum value of g 112 . By varying all parameters we obtain a nice 3D plot which is presented in appendix D; by contrast, in this section we will restrict ourselves to showing only 2D plots that each correspond to a fixed value of m 2 /m 1 . For example, at equal masses m 2 /m 1 = 1 we have figure 7 which shows the upper bound as a function of α. Although holding α fixed is convenient for the numerics (as explained above), it is sometimes more useful to visualize the allowed space of couplings (g 112 , g 222 ) instead. To do this we simply multiply the α axis in the numerics by g 112 , and in this way we can represent the same m 2 /m 1 = 1 data as in figure 9. Applying the same mapping to other mass ratios in the range m 1 / √ 2 < m 2 < √ 2m 1 we furthermore obtain the panels shown in figure 8. (As explained below, the results for m 2 < m 1 / √ 2 are somehow trivial due to screening.) For the most part, the numerical bounds in these figures significantly improve the bounds single amplitude bounds derived in the introduction which set the box sizes.
The most remarkable feature of figures 7 and 8 is the existence of a pronounced maximum of g 112 , which is attained for a non-trivial value of the ratio α = g 222 /g 112 . In particular, for equal masses this maximum (point B in figure 7) is a clear kink in the boundary of the allowed region. It would be remarkable if there is a physical theory sitting close to this kink. As shown in figure 10, such a theory should not be integrable.
Sometimes the numerical red dots in figures 8 approach the solid black lines. When this happens the full numerical bounds saturate the analytically derived diagonal bounds. We : Maximum coupling g 112 as a function of α = g 222 /g 112 for a Z 2 symmetric theory with an odd and an even particle both with the same mass. Solid black: bounds from single amplitude analytics. Red: bounds from multiple amplitudes numerics. The interesting points A, B, C, D are discussed in more detail in the next section. Multiplying the α axis by g 112 we convert this plot into a plot of the allowed coupling space (g 112 , g 222 ), see figure 9. see that this happens for very small g 112 9 and when we approach the boundaries of the mass range m 1 / √ 2 < m 2 < √ 2 (for some small values of α). This is not surprising: when g 112 → 0 we decouple the odd and even particles. Since there would be no poles in any amplitude but in M 22→22 , the bound would reduce to the single amplitude bound coming from the 22 → 22 process and yielding This explains the analytic bound in figure 3. In the second case, when we approach the 9 The fact that the numerical points do not exactly touch the vertical lines in panels (a)-(d) in figure 8 when g 112 0 is due to numerical convergence. In that region it would be more sensible to ask for the computer to maximise g 222 instead of g 112 . This would lead to numerical saturation of the vertical lines.  boundary of the mass range, we expect screening to be very important since the extended unitarity region becomes quite large. The poles in the M 11→22 component can now be almost perfectly screened, see also appendix E.2, allowing for the diagonal amplitude bounds on g 112 to be saturated. We omitted panels for m 2 < m 1 / √ 2 since in this range we can have perfect screening for any value of g 112 /g 222 , so that the multiple amplitudes bounds in the (g 112 , g 222 ) plane coincide with the rectangular frame derived from diagonal processes.
Note also that there are no vertical walls in the last row of panels in figure 8 since for m 2 > 2 √ 3 m 1 there are no longer analytic bounds on g 222 from the 22 → 22 amplitude. As the extended unitarity region in 22 → 22 becomes bigger, it becomes increasingly more effective at screening the pole, until at m 2 /m 1 = 2/ √ 3 the s (t) channel 22 → 11 production threshold collide with the t (s) channel pole as discussed in figure 4. After this mass ratio, the discontinuity across the cut can completely screen the bound state pole implying that its residue can be arbitrary. This is indeed nicely backed up by our numerics as seen in the last two panels in figure 8 where we see that g 222 becomes unbounded at this mass range.
Finally, we can also look for the maximum value of either coupling (g 112 or g 222 ) leaving the other coupling arbitrary. In other words, how tall (g 112 ) and wide (g 222 ) are the darker allowed regions in (8) where the allowed coupling live. Once plotted for various mass ratios, this gives figures 2 and 3 in the introduction.
Each optimal S-matrix at the boundary of the allowed coupling space is numerically seen to saturate the extended unitarity equations (22)(23)(24)(25)(26). This means that the scattering of two particles of type 1 or 2 can never lead to multi-particle production. Processes such as 11 → 222 are forbidden. When dealing with 2D S-matrices, in particular extremal examples saturating unitarity such as the ones stemming from this numerical computation, we are commanded to look for integrable field theories. For m 2 = m 1 , these are only possible if the inelastic amplitudes M 11→22 = M Backward 12→12 = 0 but no S-matrices we found satisfy this condition 10 so the boundary S-matrices we find for m 2 = m 1 can at most be close to those describing good physical theories. This still leaves the possibility of finding interesting physical theories along the equal mass line m 2 = m 1 . 11

(Surprises at) the m 2 = m 1 line
Indeed, nice surprises are to be found in the m 2 = m 1 line depicted in the two equivalent figure 9 (depicting the space of allowed couplings) and figure 7 (for the maximum coupling g 112 as a function of the coupling ratio α ≡ g 222 /g 112 ). Although equivalent these two figures highlight different aspects of this interesting line so it is worth having both in mind.
As concluded in the last section, the line m 2 = m 1 is where our hope lies if we are to match the S-matrices we obtained numerically with physical integrable theories. This necessary condition is not sufficient. For an extremal S-matrix to correspond to an integrable theory it should also obey the factorization conditions encoded in the so-called Yang-Baxter equations [22,23]. In figure 10 we see how our extremal S-matrices fail to satisfy these conditions as we move along the allowed coupling region (by sweeping α).
Before unveiling which analytic S-matrices we successfully identified along the m 1 = m 2 line let us go over these numerics in some detail: We observe that for large negative α Yang-Baxter is violated until we reach α = −1 (i.e. when the couplings are equal up to a sign, g 222 = −g 112 ) at which point Yang-Baxter is beautifully satisfied. This point is isolated; immediately to the right of α = −1 Yang-Baxter fails again. It is curious to note that this special point -our first candidate for a physical integrable theory -marked with an A in figures 9 and 7 looks absolutely innocent there, without any apparent kink features. As we increase α further into positive values we reach point B for α 0.76 where the coupling g 112 is maximal. As seen in figures 9 and 7 this point is a nice kink. Since it does not obey Yang-Baxter, however, this can hardly correspond to a physical theory. As we increase α further we reach α = +1 marked with a C in figures 9 and 7 where again something interesting happens. At that point something goes unstable as far as testing of Yang-Baxter goes indicated by the shower of points in figure 10 from α = 1 until somewhere around α 1.2. Once this mess settles we observe a nice line where Yang-Baxter is satisfied throughout! A particularly nice point along that line is point D located at α = √ 3. That point actually is the one furthest from the origin in figure 9, in other words, it maximizes g 2 112 + g 2 222 and as described below it corresponds to a nice known physical theory.
Now we unveil what we found about these points. In short (setting m 1 = 1 here): • Point A is a massive deformation of the three state Potts Model.
• Point B is yet to be identified. We do not know the analytic form of the corresponding S-matrix; Since it does not obey YB it can at most be close to a physical theory.
• Point D is (an analytic continuation of the lightest breather S-matrix of) the supersymmetric Sine-Gordon model.
• There is a line going from point C at α = 1 all the way to α = ∞ where the optimal S-matrix is given by an elliptic deformation of the super-symmetric Sine-Gordon. We are unaware of a physical theory with this S-matrix. Point C is the tip of the elliptic deformation where it becomes hyperbolic. Here For comparison recall that the analytic diagonal bounds were |g max 112 | = |g max 222 | = 4.56. We will now slowly build up towards those conclusions. The first observation, reviewed in appendix F, is that the full numerical optimization problem can actually be diagonalized and solved exactly for α = ±1, that is when the two physical couplings are the same up to a sign. For α = −1 the result is the S-matrix of the massive deformation of the three-state Potts model [19] 12 From this solution we read off g 112 = −g 222 = 3 √ 3m 2 1 which matches perfectly with point B in figures 9 and 7. In appendix G, we briefly review the 3 state Potts field theory.
As also explained in appendix F, the point α = +1 is the other point where we can find a clever change of basis to diagonalize our problem and compute the maximal couplings analytically to find g 112 = +g 222 = 6 √ 3m 2 1 which again matches perfectly with point C in figures 9 and 7. What we also observe in the process of deriving that analytic solution is that the S-matrix saturating this bound is not unique; there are zero modes. This is probably the explanation of the shower in figure 10. These zero modes are probably only present for α = +1 but in the vicinity of this point there is probably still some small numerical remnant thereof. We thus expect the shower in figure 10 to be nothing but a zero-mode related numerical artifact; the true solution to the optimization problem probably obeys 12 Here we rotated the one particle basis from [19] √ 2 , so that the charge conjugation operator is diagonalized. This operator is to be interpreted as the Z 2 symmetry generator. In the |A , |A † basis the S-matrix is diagonal and that is why we can solve this point exactly, see Appendix F.
Yang-Baxter for any α > 1. Yet, since this seems to be a zero mode issue, we expect the coupling as predicted by the numerics to still be correct. We will soon provide very strong evidence for these claims.
Point D for α = √ 3 is a potentially interesting point if we interpret the Z 2 symmetry as fermion number and think of particles 1 and 2 a Majorana fermion and a boson respectively. Then the condition g 222 /g 112 = √ 3 would follow for theories where these two particles are part of a N = 1 supersymmetry multiplet, see also [26]. Inspired by this -and by [26] -we tried to compare the optimal S-matrices at g 222 /g 112 = √ 3 to those of the lightest breathers of the super-symmetric sine-Gordon theory. 13 Beautifully, although we only impose the SUSY condition at the level of the couplings, we see that SUSY emerges at the level of the full S-matrix elements and indeed the optimal S-matrix saturating our bounds at point D is an analytical continuation of the lightest breather supermultiplet of the super-symmetric Sine-Gordon! Unfortunately, while we are able to check this to very convincing numerical accuracy we have no analytic derivation of this statement. For completeness, here are some super SG formulae [21].
The lightest breather supermultiplet SSG S-matrix S and where γ is fixed so that the bound state mass is equal to 1. That is γ = 2π/3. Note that even though the overall scalar factor in the S-matrix is invariant under γ = 2π/3 ↔ γ = π/3, the matrix part is not. This is the sense in which our S-matrix is an analytic continuation of SSG. (compare with SG, in which picking m b = 1 instead of m b = √ 3 only leads to an 13 Strictly speaking we are comparing with an analytic continuation of that S-matrix since our bound-states have mass equal to the external particles while the next-to-lighests breathers of super-symmetric sine-Gordon have mass bigger than √ 2 times that of the external particles. In [2] the usual bosonic sine-Gordon S-matrix was identified as the theory with the largest coupling in the S-matrix of the lightest particle with a single bound-state of mass m b . When m b > √ 2 this is kosher but as m b < √ 2 (and in particular for m b = 1) we also need to extend the definition of the SG S-matrix beyond its original mass range. In that case it amounted to multiplying the S-matrix by −1. Here the situation is morally the same but the modification ends up a bit more complicated. This means that here -as there -we do not know a physical theory and we can only write an exact S-matrix that saturates the bound. 14 In the basis |11 , |12 , |21 , |22 so that the second (third) element on the second row is the forward (backward) 12 → 12 component, for instance. overall minus sign). More generally m b /m 1 = 2 cos γ/2 and the physical mass range for SSG is 2 > m b > √ 2.
Two-particle Z 2 symmetric solutions of the Yang-Baxter equations are classified [18]. It is natural that if an extension of SSG exists with an extra parameter, that it is given by an elliptic solution of the Yang-Baxter equations. In fact, examining the classified solutions we see that the only good candidate for being promoted to an S-matrix with all the symmetry properties we have and that reduces to SSG in the trigonometric limit is solution 8VII of [18], which is equivalent to where we normalised by the 12 → 12 forward component so that comparison with SSG is easier. = ±1 from YB. The s-channel poles of 12 → 12 forward and backward correspond to the flow of a particle of type one, and therefore the residues of this two amplitudes at the θ = 2πi/3 pole must coincide. This fixes sn(γω|κ) = sn( 2πi 3 ω|κ). Crossing symmetry together with the fact that in the trigonometric limit κ → 0 we have to recover SSG fix ω = − i π K(κ) where K is the complete elliptic integral of the first kind (more precisely, crossing gives ω = − i π (2n + 1)K(κ) with n an integer. The κ → 0 limit fixes n). This completely fix a crossing symmetric matrix structure up to one free parameter, κ, which hopefully is unconstrained. There is a miracle going on. For our amplitudes, we have that and, moreover, res If (38) is a candidate of matrix structure of the S-matrices saturating the numerical bounds, we must have the same relation between the respective components. It turns out that this holds automatically for any κ after all the conditions above are imposed. Otherwise this would fix κ = 0 and we would conclude that there are no Yang-Baxter deformations respecting the symmetries and spectrum of our problem. So all we need to do now is to unitarize and introduce the poles. Note that and g(θ) ≥ 1 for θ ∈ R. Therefore, as follows from [2], to unitarize ED we just need to multiply it by while to introduce the poles, we multiply by CDD pole with direct channel pole at 2πi/3, At the end of the day, a candidate for a unitary, crossing symmetric, integrable deformation of the supersymmetric sine-Gordon reads Points D and C with α = √ 3 and α = 1 would now correspond to κ → 0 and κ → 1 respectively. As κ → −∞, α → ∞. We can now compute the couplings associated to the elliptic deformation (44), cross our fingers and compare those couplings with the numerics of figures 9 and 7. The elliptic deformation analytic results are the solid chartreuse 15 lines in those figures. The agreement could not be better. Note that the agreement goes all the way to point C and since by construction the elliptic solution obeys Yang-Baxter, the shower in figure 10 should indeed be a simple zero-mode related numerical artifact. To make precise the elliptic notation used here, we present in appendix J a representation of this S-matrix in Mathematica friendly notation, ready to be copy pasted so the reader can more easily explore this exotic solution.
An obvious question is whether this elliptic deformation corresponds to a nice physical theory. 16 Since the supersymmetric sine-Gordon we encountered here is not a totally kosher theory but an analytic continuation thereof, it is natural to first extend this analysis to the mass range were the super-sine Gordon breather lives and to study its elliptic deformation for those more physical set of parameters. This is beyond the scope of this paper and is currently being investigated in [26].
The m 2 = m 1 line was indeed full of surprises.

QFT in AdS
In the previous section we have numerically explored the space of scattering amplitudes that allow for a Mandelstam representation and we found examples of amplitudes that appear to maximize couplings subject to the unitarity constraints. As explained in footnote 1, these extremal coupling constants are not true 'upper bounds': although our numerical results appear to have converged, a numerically more refined ansatz will find slightly larger values.
An orthogonal approach to the extremization of three-point couplings in field theories was developed in [1]. The idea is to consider a field theory in an AdS background and investigate the 'boundary' correlation functions that are so familiar from the AdS/CFT correspondence. In our setup gravity is non-dynamical and this translates into the absence of a stress tensor among the set of boundary operators. Nevertheless it is natural to claim [1] that these correlation functions obey all the other axioms of a unitary CFT, including crossing symmetry, making them amenable to a numerical bootstrap analysis as in [4]. In this way 15 Chartreuse, of selcouth beauty, is a colour half-way between yellow and green. 16 We thank Davide Gaiotto for illuminating discussions on related topics. any general constraints on CFT data directly imply corresponding constraints for QFTs in AdS, and by extrapolating these results to the flat-space limit we can get constraints on flat-space QFTs as well. (For a gapped QFT in AdS 2 scaling dimensions and masses are related as m 2 R 2 = ∆(∆ − 1) and therefore the flat-space limit is typefied by sending all scaling dimensions ∆ → ∞ whilst keeping ratios fixed, The QFT in AdS approach uses CFT axioms to provide rigorous upper bounds, at least modulo our extrapolation procedures. It does not assume analyticity or any particular behavior at large complex energies, and the unitarity constraints are phrased in terms of reflection positivity rather than in terms of probabilities. And yet it was shown in [1] that it provides upper bounds on flat-space couplings that are numerically equal (up to three significant digits in some cases) to the extremal couplings obtained with the S-matrix bootstrap methods. In this section we demonstrate that this striking equivalence was not just a fluke by employing once more the QFT in AdS approach to reproduce some of the previous S-matrix bootstrap results from conformal crossing equations.

Setup
We will consider four-point functions of operators on the real line, which we think of as the boundary of an AdS 2 space with curvature radius R. There are two distinguished operators φ 1 and φ 2 of dimensions ∆ 1 and ∆ 2 , which correspond to the two single-particle states of the setup described above -in particular it is understood that Besides the assumed Z 2 symmetry, under which φ 1 is odd and φ 2 is even, we will also assume that the QFT is parity invariant and that φ 1 and φ 2 are both parity even. 17 The OPEs are, schematically, Here the (non-)appearance of φ 1 and φ 2 on the right-hand sides is dictated by Z 2 symmetry. The other operators are meant to correspond to multi-particle states for the QFT in AdS and their minimal scaling dimension mimicks the beginning of the two-particule cuts in the corresponding scattering amplitudes. The parity properties are dictated by the parity of the operators on the left-hand side. We should add that the OPE coefficients λ ijk are related to bulk couplings with the unsightly relative normalization coefficient [27] where ∆ ijk = ∆ i + ∆ j − ∆ k . This relation was explained in [1].
In one dimension conformal transformations preserve operator ordering modulo cyclic permutations. This leads to the following non-equivalent four-point functions and we will numerically analyze the lot of them. 18 Our recipe follows that of [1] with minor variations. Suppose that we wish to obtain a bound on g 2 112 (in units of m 1 ) for a given coupling ratio α = g 222 /g 112 and mass ratio µ = m 2 /m 1 . We then proceed as follows: 1. Choose a ∆ 1 . Then set ∆ 2 = µ∆ 1 and also fix the ratio 2. A single conformal bootstrap analysis of the four correlators listed above now yields a numerical upper bound on λ 2 112 . Our multi-correlator bootstrap analysis is very similar to the one introduced in [12] where it was successfully applied it to the threedimensional Ising model. The systematics of our analysis (normalizations, conformal block decompositions, functionals) can be found in appendix K. The bound so obtained also depends on the number of derivatives of the crossing equations that we analyze and this introduces a new parameter Λ, so we write where we use (46) to pass from (λ 2 112 ) max to (g 2 112 ) max .
3. Upon repeating step 2 for various Λ one finds that (g 2 112 ) max depends significantly on Λ. To obtain an estimate of the bound that we would obtain if we could analyze all the crossing equations, i.e., if we possessed infinite computational resources, we extrapolate the results for various Λ to estimate In practice we do this by fitting a polynomial through data points ranging from Λ = 32 up to Λ = 140. 19 Examples of this extrapolation are shown in figure 26 on page 68. This limit provides our estimate for the best possible upper bound for a QFT in AdS with two particles with masses determined by ∆ 1 and µ and bulk coupling constant ratio given by α. 18 An interesting observation is that the φ 1 φ 2 φ 1 φ 2 correlator does not feature the identity operator. A conformal bootstrap analysis of this correlator in itself therefore does not give any bounds whatsoever because it lacks an overall normalization. This is completely different from the forward 12 → 12 amplitude which we have seen can give a meaningful bound on g 112 . However we will shortly see that the ensemble of correlators does give numerical results that mostly agree with the ensemble of amplitudes. 19 In [1] we were able to obtain results up to Λ = 200 or Λ = 300 for the different scenarios. The multi-correlator analysis of this paper is numerically more demanding, even more so because the rho series expansion [28] for conformal blocks with large unequal dimensions converges much more slowly. 4. We view ∆ 1 as a proxy for the AdS curvature radius. We therefore repeat steps 1 to 3 for a number of different values of ∆ 1 and once more extrapolate to infinite ∆ 1 to obtain a result on the flat-space coupling: This is the coupling we can compare with the flat-space S-matrix bootstrap analysis.
Appendix K.5 contains technical details of the extrapolation procedure.
For finite Λ and ∆ our results provide rigorous upper bounds on the three-point couplings for any QFT in AdS that obeys the stated assumptions. Once we begin the extrapolations we introduce errors that are hard to quantify and this is an unavoidable drawback of our method. Nevertheless we will soon see, as was the case in [1], that the extrapolated bounds appear to accurately reproduce the S-matrix bootstrap results in most cases.

Results
The numerical algorithm outlined above is computationally demanding. For a single µ and α we need about 10 different values of ∆ 1 and for each of these we need about 15 different values of Λ, implying about 150 multi-correlator bootstrap runs. We have therefore chosen a few representative values of µ and α to demonstrate both the feasibility of the multicorrelator conformal bootstrap approach to scattering processes and the match with the flat-space S-matrix bootstrap results.

Results for equal masses
Our first plot is for µ = 1 so we have two particles of equal masses. In figure 11 we overlay the QFT in AdS results (isolated data points) with the S-matrix bootstrap region shown before in figure 9. The black frame again indicates the single-amplitude bounds, which are in fact equal to the single-correlator bounds found in [1]. We have performed a multi correlator QFT in AdS analysis for ratios α = g 222 /g 112 equal to +1, 0, −1 and −8/3 and in all cases we find reasonably good agreement with the multiple amplitude S-matrix bootstrap result. For α = −8/3 our bound comes out somewhat higher than the value reached by the S-matrix bootstrap. This might be due to our extrapolation procedure, which also makes it difficult to put error bars on the QFT in AdS points, but it might also be a consequence of the finite truncation level in the S-matrix bootstrap. It is of course reassuring that the S-matrix bootstrap (Yin) always gives lower values than the conformal bootstrap (Yang).
For the data points in figure 11 the extrapolation is standard, i.e., as outlined above and elaborated on in appendix K, but the data point requires a comment. In that case we found that the maximal squared coupling (g 2 112 ) max [µ, α, ∆ 1 , Λ] from the multi-correlator analysis is always numerically equal to one half of the corresponding maximal squared coupling obtained from the single-correlator analysis -even for finite ∆ 1 and Λ so before any extrapolations. Therefore, rather than doing a detailed multi-correlator analysis, we just plotted one half the single-correlator result. This factor of one half is understandable: using a change of operator basis similar to the one described in appendix F, one finds that the multi-correlator problem effectively becomes that of two decoupled single-correlator problems which each feature a squared coupling that is rescaled by a factor 2. 20

Results for α = −1
Our next result is shown in figure 12, where we have assumed g 222 /g 112 = −1. We will discuss in turn the black curve, the red shaded region, and the green (new) data points.
The black curve corresponds to the best single-correlator bound for the given mass ratio. It is actually made up of two parts: for m 2 > m 1 it is the bound obtained from the 1111 four-point function, whereas for m 2 < m 1 it is the bound obtained from the 2222 fourpoint function. These single-correlator bounds were already obtained in [1] and were shown to agree with the single-amplitude analysis of [2].
In red we show the multi-amplitude results obtained with the methods discussed in section 3.1. It is again made up of different parts: for 1/ √ 2 < m 2 /m 1 < √ 2 we can use the numerical 20 It is essential here that α = 1 so λ 222 = λ 112 .  Figure 12: The green data points show the maximal values of log(g 2 112 ) provided by the multicorrelator QFT in AdS analysis, now for α = g 222 /g 112 = −1 and as a function of the mass ratio m 2 /m 1 . For comparison we have also added the single-correlator QFT in AdS bounds from [1], as solid lines, as well as the S-matrix bootstrap data, in red. The plot naturally splits into three regions. First of all, for m 2 /m 1 < 1/ √ 2 there is screening and the multi-amplitude bound reduces to the single-amplitude bound. As shown, the correlator bounds nicely follow this behavior. Moving rightward, for 1/ √ 2 < m 2 /m 2 < √ 2 we find a respectable match between the multi-correlator and the multi-amplitude data, in particular we again recover the three-state Potts field theory at m 2 = m 1 . For √ 2 < m 2 /m 1 there are Landau singularities and the multi-amplitude analysis becomes complicated. However we know that the multi-amplitude bound must lie at or below the single-amplitude bound from S forward 12→12 , meaning that it must end up somewhere in the striped region. The multi-correlator analysis, on the other hand, appears unable to improve on the weaker 1111 single-correlator bound.
analysis and we take the α = −1 slice from figure 16. For m 2 /m 1 < 1/ √ 2 we have screening and the multiple amplitude analysis does not give stronger results than the analysis of the single amplitude S 22→22 which, as we stated before, agrees with the 2222 single-correlator bound. For √ 2 < m 2 /m 1 there are Landau singularities and a more sophisticated analysis is necessary to obtain multiple amplitude results, but we do know that the maximal coupling from the multiple amplitude analysis can only lie below the single amplitude bounds. In particular, it must lie below the bound obtained from S forward 12→12 , which was given as the solid line in figure 2 in the introduction and here yields the striped region in figure 12. 21 Finally, the new data points obtained from the multi-correlator conformal bootstrap are indicated in green. The data points are obtained from a standard extrapolation, as before, whereas for the data points the numerical multi-correlator analysis gave identical results to the numerical single-correlator analysis for all ∆ and Λ. The extrapolation will therefore trivially equal the single-correlator result as is indicated in the plot. 21 We explain in appendix L that Landau singularities do not appear in M forward 12→12 for any mass ratio in the range 0 < m 2 /m 1 < 2, so the corresponding single-amplitude bound should be perfectly valid.
For all points with m 2 /m 1 < √ 2 our extrapolated QFT in AdS results lie at or just above the S-matrix results. As for the previous plot, the small finite difference might be either due to our extrapolation procedure but also due to the S-matrix bootstrap results not yet having converged. The points with m 2 /m 1 > √ 2 are more puzzling. The striped domain, as we mentioned, arises from an analysis of S forward 12→12 alone and so a multi-amplitude analysis (with Landau singularities and all) will only be able to land somewhere in that domain. Unfortunately this single-amplitude bound does not seem to be picked up by the multi-correlator analysis at all. 22 It would be nice to know why this is the case: are we missing constraints to be imposed for the QFT in AdS bootstrap? 23 Alternatively, is there maybe a 'phase transition' where by pushing to a very high number Λ of derivatives the single-and multi-correlator analysis begin to differ? Such a phenomenon would obviously invalidate our large Λ extrapolations and might therefore resolve the puzzle. It would be somewhat analogous to the observations discussed in appendix K.4, see in particular figure 25, where we explain that taking different Λ's for different crossing equations leads to nonsmooth behavior.
Of course, as we mentioned at the beginning of this section, at a technical level the conformal bootstrap analysis looks completely different from the S-matrix bootstrap. We are confident that both analyses yield valid constraints on three-point couplings, but besides physical intuition there was no a priori guarantee that these constraints had to be exactly the same. From this perspective the aspect most in need of an explanation in figure 12 is the quantitative match between the results for m 2 /m 1 < √ 2 (and similarly for all points in figure  11 and the results of [1]) rather than the discrepancy in the other points. Either way, the precise connection between conformal correlators and scattering amplitudes warrants further investigation.

Discussion
We have demonstrated the feasibility of the multiple-amplitude bootstrap for the lightest two particles, and shown that it gives stronger bounds compared to the simpler S-matrix bootstrap of only the lightest particle. Clearly one expects to get increasingly stronger bounds by considering more and more scattering processes. It is interesting to consider how such results could converge to an 'optimal bound' that we would obtain by considering the entire S-matrix, as follows.
In all our numerical experiments it turns out that the unitarity condition at all energies is (numerically) saturated in the subspace we work on. To illustrate this, consider for example the various possible outcomes from scattering the lightest particle in our setup -let us say 22 See also footnote 18. 23 At least the multi-correlator result, which is a hard upper bound, is above or equal to the S-matrix bootstrap results in all cases, so the results are not in direct conflict with each other.
that it is the Z 2 odd particle -against itself. Probabilities must add up to 1 so where the red and blue inequalities follow trivially from probability positivity as indicated in dark green. In [2] we effectively considered only the weakest red inequality. The theories which lie on the boundary of this space turn out to saturate this inequality for all energies. This would mean that any other process has zero probability, and a theory saturating the bounds of [2] must therefore have zero particle creation or transmutation since the only allowed process is elastic scattering. This is possible for special cases like integrable theories, but generically it cannot be the case. Therefore, by including also the constraints of the other processes we should get better bounds and this is indeed what we have observed in this paper: our improved bounds are due to the stronger constraint given by the lower blue inequality.
However, we once more found that the optimal solutions now saturate this new condition for all values of the energy, so the remaining processes for the extremal S-matrix are again all zero. 24 In other words, theories lying on the boundary of the new space are theories where particles 11 can continue into 11 or transmute into 22 but we still get zero probability for all other processes that have a different final-state particle content. We expect this pattern to continue by including more and more processes in the game, i.e. we will continue to observe unitarity saturation within the subspace we consider, no matter how large. As we increase the size of our truncation we will hopefully asymptotically approach an optimal bound, but we are unlikely to hit a non-integrable theory at our boundary if we consider only a finite number of processes. 25 This gives our S-matrix bootstrap approach a more asymptotic flavor than the numerical conformal bootstrap.
Indeed, for our setup with two particles with mass m 1 and m 2 we find no integrable theories if the masses are different, whereas if they are equal then there are exciting physical theories at our boundary: the three-states Potts model and (an analytic continuation of) the super-symmetric sine-Gordon model. We also find a full segment around the supersymmetric sine-Gordon theory which seems to obey all the necessary factorisation requirements to be an integrable theory; it would be very interesting to see if that is the case. (More on this in [26].) We also discussed how the same bootstrap results can be obtained from AdS, using the setup first discussed in [1]. Putting a gapped Z 2 symmetric theory into an AdS box induces a one dimensional Z 2 symmetric conformal theory in its boundary which we can analyze by numerical conformal bootstrap methods. To make contact with flat space, we take this box to be large which corresponds to large scaling dimensions on the boundary. As it happens the numerical conformal bootstrap results become rather weak at large scaling dimensions and this makes it computationally quite challenging to obtain reliable results. This is the main drawback of the AdS approach. Fortunately there are interesting and potentially very helpful developments on this front: according to [29], convergence can be much improved by a smarter choice of functionals (see also [30,31]). It would be very interesting to explore this further.
An important advantage of the AdS approach, on the other hand, is that it requires no subtle assumptions about the various analytic properties of scattering amplitudes. The AdS box is thus a literal black box from which we can get beautiful S-matrix bootstrap results even when the analytic properties of such amplitudes might be less obvious.
A good example of such a subtle assumption is extended unitarity, which we have seen is crucial for our multiple-amplitude bounds. Recall that this is a generalisation of usual unitarity which controls the analytic behaviour of scattering amplitudes for unphysical energies, below physical thresholds. For example, when we scatter the next-to-lightest particle against itself we have a two particle cut associated to the lightest particle starting at s = (2m lightest ) 2 , before the physical two particle cut at s = (2m next-to-lightest ) 2 , and extended unitarity governs the discontinuity of scattering amplitudes in the segment between those two values. Extended unitarity is built into perturbation theory [17] but it is not straightforwardly justified non-perturbatively. However, the fact that our QFT in AdS approach exactly reproduces the flat space results provides strong evidence for the validity of the extended unitarity assumption.
Relatedly, it is puzzling that the AdS bounds for m 2 > √ 2m 1 are so much weaker than the 12 → 12 forward flat space extremal coupling, see figure 12. Either the AdS numerics did not converge yet or perhaps there is something deeper to be learned there. It would be very interesting to extend our flat space analysis beyond the mass range (1) into the m 2 > √ 2m 1 domain. Here we would need to include so-called Coleman-Thun singularities in our setup. An example is shown in figure 13.
Another very important new ingredient which is not unrelated to the anomalous cuts arising in the extended unitarity region and which appeared in this work is the phenomenon of screening. Amplitudes involving heavier particles can sometimes produce discontinuities in some of these anomalous cuts which cancel, i.e. screen, other singularities such as physical poles corresponding to bound-state stable particles. That is, these discontinuities can be tuned so that the amplitudes can often be quite small in the physical energy region where experiments are done. This mechanism is not only possible but it is actually realised by some amplitudes which lie at the boundary of the truncated 2 → 2 multiple correlator S-matrix space. It is natural to expect that similar phenomena would also be realised at the boundary of the full S-matrix space. If this boundary is physically significant then screening would be an interesting way for nature to hide strong couplings from the observer.

Let us conclude with some interesting open problems and future directions.
One open problem for which we now have all tools to explore concerns the tricritical Ising field theory. This theory is obtained by deforming a conformal minimal model with two relevant deformations, see appendix H, and is integrable if one of the deformations is set to zero. It would be very interesting to bootstrap this theory when both parameters are non-zero. This is a particularly nice case study because the next to lightest particle here is well below √ 2 times the mass of the lightest particle -see table 1 for its value at the integrable point -which means we can readily apply all the methods developed here. The only modification would be to include further poles in the ansatz corresponding to the other stable particles this theory has -see again table 1. So here is a homework exercise: consider a line in the mass ratio parameter space which passes by the masses of the integrable theory. Something remarkable should happen: In the anomalous cut of the 22 → 22 amplitudes we should see a peak developing as we approach the integrable theory. This peak is going to become a new stable particle in the integrable theory with a mass we know. Seeing this peak show up in detail would be great, as it would constitute a "discovery" of a new particle through the S-matrix bootstrap. Of course, more interesting still would then be to move away from the masses of the integrable theory and explore the full tricritical Ising field theory, nonintegrable and all. The previously discovered sharp peak -typical of an integrable theorywould now be smoothened out and correspond to a nearly stable resonance. Because there are so many masses and couplings this would be a challenge numerically, albeit a worthwhile one.
Another open problem which we could now easily address is the problem of multiple amplitudes without Z 2 symmetry. We would now also include amplitudes such as 11 → 12 which have amusing 2D kinematics by themselves. The Ising field theory perturbed by both magnetic field and temperature would be a perfect case study for this case.
A much more challenging but very interesting open problem would be to extend the multiple amplitude analysis to higher dimensions (as in [3] and [10]). The Z 2 spontaneously broken phase of the φ 4 model in 3D, for instance, seems to have a single stable bound state of mass m 2 1.8m 1 [32]; it would be fascinating to try to bootstrap this S-matrix.
Finally, another frontier in 2D would be to delve into the multiple particle S-matrix bootstrap. Can we tame scattering of 2 particles into 3, 4, . . . final particles? There are two obvious obstacles. One is the analytic structure of these amplitudes. They depend now on more kinematical variables and have a huge plethora of Landau singularities; it is unclear if we can characterise them fully. The other challenge is even more basic: can we close up the system of equations? Suppose we consider a basis of initial and final states with both two and three particles. Then we need to deal with the 3 → 3 amplitudes. But those, by crossing, are related to 2 → 4. By unitarity we would then need to include four particles in the final and initial states as well. But then we are forced to consider 4 → 4 processes which are now related by crossing to 3 → 5 and 2 → 6 scattering and so on. It seems we are suddenly obliged to consider any number of final particles at once which of course would be computationally completely infeasible. Hopefully we can find a suitable truncation scheme. Along these lines, perhaps we could first try to get some inspiration from the AdS side. Some of the necessary higher point conformal blocks are well known in 1D [33], so can we use this to devise a 1D CFT bootstrap numerical problem dual to the very intimidating flat space multiple particle bootstrap? Even if very challenging numerically, this would prove of extreme conceptual value.

A Two dimensions and higher dimensional kinematics
The Mandelstam plane provides us with a very useful depiction of the (real sections of the) interrelations between the three Mandelstam variables u, s, t. The first important object in this plane is the Mandelstam triangle.
Consider a two-to-two process involving particles with momenta p a , p b , p c , p d associated to particles of mass m a , m b , m c , m d . We define the three Mandelstam invariants s = (p a + p b ) 2 , t = (p a + p c ) 2 and u = (p a + p d ) 2 . If particles p a , p b are the two incoming particles then √ s is the centre of mass energy of the scattering process. The same is true if the incoming particles are particles p c , p d . In either of these cases the process can only be physical if we have enough energy to produce both the initial and final state, so for s ≥ max((m a + m b ) 2 , (m c + m d ) 2 ). Of course, the same scattering amplitudes can describe other channels. 26 If p a , p c (or p b , p d ) are the two incoming particles then √ t is the centre of mass energy of the scattering process and similarly for u so the physical conditions in those cases would read t ≥ max((m a + s t u  Figure 14: The Mandelstam triangle is the region where all Mandelstam variables are below their corresponding two-particle threshold: To be in a physical region we thus need to be in the pink region. This is necessary but not sufficient. We need to have enough energy but we also need to scatter at a real angle. For instance for incoming particles p a , p b we can easily compute the scattering angle to find . For any left hand side between −1 and +1 corresponding to a real angle, and for any s in the physical range this equation determines a physical t. (Of course, u = m 2 a +m 2 b +m 2 c +m 2 d −s−t is automatically fixed.) The set of physical s and t determined in this way determine the physical region in the s-channel. The other channels are treated similarly with cos(θ ac ) = RHS of (54) m b ↔mc, s↔t , cos(θ ad ) = RHS of (54) m b ↔m d , s↔u If all masses are equal (to m) then (54) reduces to the famous relation so that the physical region in the s-channel is simply s > 4m 2 and 4m 2 −s ≤ t ≤ 0 represented by the top blue region in figure 15a. In two dimensions the angle ought to be 0 or π so that we have either t = 0 or t = 4m 2 − s, i.e. u = 0. These two conditions (t = 0 and u = 0) are equivalent if the external particles are indistinguishable so in that case we can pick either; in the main text we took u = 0. Note that these two conditions are nothing but the boundary of the darker blue region. Similarly we could study all other channels which we can simply obtain by relabelling the Mandelstam variables in (56). The two extra physical regions are the other two blue regions in the same figure 15a. Note that the boundary of all these blue regions can be written concisely as stu = 0 which is nothing but the constraint obtained in the main text from the two dimensional constraint (9) and represented by the dashed lines in figure 15a.
Finally, we come to the more interesting case where the external masses are only pairwise equal: m a = m b = m 1 and m c = m d = m 2 . This same configuration can describe the 11 → 22 process (with √ s being the centre of mass energy), and the 12 → 12 processes (with √ t being the centre of mass energy). For the first case we use (54) to get .
The physical region corresponding to −1 ≤ cos θ 11→22 ≤ +1 is now a more interesting curved region, represented by the darker blue region in figure 15b. Again, in two dimensions we can only have backward or forward scattering so t must saturate one of these inequalities. If the particles of the same mass are indistinguishable then both solutions are equivalent as before.
We can thus take θ 11→22 = 0 without loss of generality leading to and reproduce in this way the results below (11) in the main text. In the t-channel particles m a = m 1 , m c = m 2 are incoming so we are scattering a odd particle against an even particle. Then we use the first relation in (55) which now reads Note that in our convention it is √ t (and not √ s) who is the center of mass energy in this channel. 27 The physical region corresponding to | cos(θ 12→12 )| ≤ 1 is now represented by the lighter blue region in figure 15b. The two dimensional conditions that the angle is 0 or π are now quite different. The former corresponds to forward scattering and is obtained by s = 0 (obtained by equating the RHS of (59) to +1) while the later corresponds to backward scattering and yields the more involved relation (obtained by equating the RHS of (59) to −1.) Note that this relation is nothing but (58) if we solve for s. In other words, these two configurations are simply related by crossing symmetry s ↔ t.
To summarize: the boundary of the physical regions are now given by the black solid line and by the black dashed curve in figure 15b. Crossing u ↔ t at s = 0 relates the left to the right of the straight line -leading to condition (19). Crossing symmetry also relates the top to the bottom branch of the hyperbolic looking curve -reflected in equation (20). In two dimensions these two curves (the hyperbola and the straight line) are independent while in higher dimensions they are smoothly connected (by moving in angle space).

B Unitarity and final state probabilities
The S-matrix is defined by the expansion of in-states in terms of out-states The in-states and the out-states are both a complete basis of the Hilbert space. Let us start by discussing the physical meaning of the diagonal unitarity equations (22) and (24). These follow from the statement that the state |A in above is normalized. However, due to the continuum of states this is a bit subtle. The trick is to contract the state with itself but with different momenta Here, the state |A in represents a state with the same particle content but different momenta. We use the standard normalization for the inner products: where the product runs over each particle in the state |B in . Unitarity then reads where J A,B is the jacobian defined by The natural physical interpretation is that is the probability of the in-state |A in end up in the out-state |B out .
For two particle states |A in = |12 in and |B out = |34 out , equation (65) reduces to with E i = m 2 i + p 2 i and This gives We conclude that, for two particle states, the transition probabilities are given by For example, for the initial state |11 in , we can write This equation is equivalent to (22) in the energy range where only 2 particle states are available. To show this we contract (61) with a generic out-state out C|, to find and use the standard definition of the amplitude M : where P A denotes the total momentum of the state |A in . This relates S with M . In the particular case of two particle states one finds where the first term is present (and equal to 1) if and only if the initial and final states are the same. This allows us to rewrite (72) in terms of M , which should be compared with (22). Notice that the physical derivation given here is not valid for 2 min(m 1 , m 2 ) < √ s < 2 max(m 1 , m 2 ) because the two particle state of the heavier particle is not available. This is the regime of extended unitarity where we must use (22).
There is also a more intuitive derivation of the full matrix form of the unitarity constraints (27). Consider the Z 2 even sector for simplicity. The matrix of inner products of the states {|11 in , |22 in , |11 out , |22 out }: must be positive semi-definite. In fact, for the range of energies (73) where there are only two particle states, the rank of this matrix must be 2 because both the in and the out states are complete basis. This is a very intuitive way to derive the unitary constraints. However, one must be careful with Jacobian factors that relate different delta-functions. Factoring out in 11|11 in and using (61) and (65), we can define a positive semi-definite matrix (without delta-functions) It is convenient to rescale the states |22 → ρ 22 ρ 11 |22 so that all diagonal entries become 1. This leads to the following positive semi-definite matrix: Notice that using (76), the 4 × 4 matrix V can be written as where M and ρ are the 2 × 2 matrices defined in (27). We will now show that the condition V 0 is equivalent to (27) for √ s > 2 max(m 1 , m 2 ). First notice that the eigenvalues of the hermitian matrix I − V take the form (−λ 2 , −λ 1 , λ 1 , λ 2 ). 28 Then, the condition V 0 implies that λ 2 i < 1. On the other hand, if we compute explicitly the square of I − V, we find Therefore, the eigenvalues of S † S must be less than 1. Equivalently, we can say that and (27) follows. 29 In the extended unitarity region 2 min(m 1 , m 2 ) < √ s < 2 max(m 1 , m 2 ) this derivation does not apply but we can still use (27).

B.1 Bounding ImM 11→22
The discontinuity of the amplitude M 11→22 does not have well defined sign. Furthermore, in the region 2 min(m 1 , m 2 ) < √ s < 2 max(m 1 , m 2 ) below the physical regime, one could worry that the discontinuity could be very large and lead to screening. However, the generalized unitarity equations (22)- (24) forbid this phenomena in the range of masses we consider.
where the label = 1 or = 2 stands for the lightest particle. Taking the modulus square of equation (85) and using the other two equations, we find Therefore the size of ImM 11→22 is related to the positive discontinuities ImM 11→11 and ImM 22→22 , which for this reason are bounded. In fact, in our numerical procedure we impose unitarity as the set of inequalities (27), which in particular implies This leads to which bounds ImM 11→22 in our setup.

B.2 Phase shifts
In the Z 2 even sector, it is trivial to define diagonal phase shifts The Z 2 odd sector is slightly more interesting. In this sector, it is convenient to use states of definite parity, Then, we define the phase shifts by In our numerical algorithm, we impose unitarity in the odd sector by the following positive semi-definite condition It is instructive to see what this implies for the phase shifts δ ± 12 . Using (92), we conclude that e 2iδ ± 12 (s) where Thus |b| ≤ a and we recover the usual unitarity inequality for the phase shifts The goal of this appendix is to prove that (4) is the amplitude with maximal coupling g 2 222 compatible with crossing symmetry and unitarity. To this end, it is convenient to define and .
Furthermore, q(s) is analytic in the s-plane minus the s-channel cut (4m 2 1 , +∞) and the t-channel cut (−∞, 4m 2 2 − 4m 2 1 ). In the extended unitarity region 4m 2 1 > s > 4m 2 2 , we have 30 Im where we assumed m 2 1 < m 2 2 < 4 3 m 2 1 . We conclude that maximizing g 2 222 is equivalent to maximizing q(m 2 2 ) subject to q(s) = q(4m 2 2 − s), Im q(s) ≤ 0 for 4m 2 1 > s > 4m 2 2 and |q(s)| 2 ≤ 1 for s > 4m 2 2 . To prove that the optimal solution is given by q(s) = 1 it is useful to change to the coordinate such that the unit disk |z| < 1 covers (the half Re s > 2m 2 2 of) the physical sheet (see figure 1 of [3] for more details). In these coordinates, the optimization problem translates to maximizing q(z = 0) subject to |q(z)| ≤ 1 for |z| = 1 and Im q(z) ≤ 0 for 0 < z(4m 2 1 ) ≡ z 0 < z < 1, with q(z) analytic on the unit disk minus the cut from z 0 to 1. Then Cauchy's theorem leads to Clearly, the optimal solution is given by q(z) = 1 (and Im q(z) = 0 inside the unit disk). Figure 16 presents the numerical results for the maximum value of g 112 for each mass ratio m 2 /m 1 and for each coupling ratio α = g 222 /g 112 . By changing axes and looking at different (d) Figure 16: Bounds on g 2 112 following from the multiple amplitudes analysis, as a function of α and m 2 /m 1 . These bounds should hold for any Z 2 symmetric quantum field theory with two particles in the spectrum, 1 being odd and 2 being even. They improve the bounds derived from individual amplitudes, corresponding to the red and black surfaces. As can be seen from the various angles (a)-(d), the bound surface has many interesting features that are described in the main text.

D 3D Plots
sections we obtain figures 8 and 9 from the main text. Many of its features were discussed in section 3.2. The black and red surface correspond to the bounds coming from diagonal processes and are, respectively, the translated versions of the horizontal and vertical solid lines in figure 8. Figure 16 has a ridge where the coupling g 112 is maximal for each mass ratio. That maximal value set an upper bound for the question: how big can the coupling g 112 be in a Z 2 symmetric theory with only two stable particles? In figure 2 we depict this maximum g max 112 (m 2 /m 1 ) (a similar analysis can be done for g 222 , leading to figure 3). We see that this maximal coupling approaches the analytic bound derived from the diagonal 12 → 12 component as the mass approaches the boundary of the mass range (1) and is otherwise significantly stronger, specially when the particles are mass degenerate where we observe a nice kink feature in figure 2.
This kink has a cute geometrical interpretation in the full 3d plot in figure 16: The top of the ridge meets a valley at m 1 = m 2 . The valley is a kink for any α. For equal masses there is no extended unitarity region and that renders the numerics way more manageable. This is why we can afford to have so many points along the valley as clearly seen in the figure.
There is one more motivation for resolving this valley region very finely: It is the natural place to look for interesting physical theories. Indeed, each optimal S-matrix in the surface of figure 16 saturates the extended unitarity equations (22)(23)(24)(25)(26). This means that the scattering of two particles of type 1 or 2 can never lead to multiparticle production. Processes such as 11 → 222 are forbidden. When dealing with 2D S-matrices, in particular extremal examples saturating unitarity such as the ones stemming from this numerical computation, we are commanded to look for integrable field theories.

E.1 Invisibility Cloak Toy Model
In this appendix we highlight the importance of not leaving any densities unconstrained as they can lead to very efficient screening thus invalidating any possible bounds. To this purpose consider as a toy model the function where 0 < a < b. We think of the first term, the pole at the origin, as a target which we would like to screen. The region [a, b] where the density term is defined is denoted as the screening region and that second term is denoted as the screening term. We can think of it as an invisibility cloak whose role is to make the full function small in pre-defined regions.
To make it concrete suppose we want the target to be screened in a region to the right of the screening region (as we would usually associate to an invisibility cloak a la Harry Potter ) but also -thus making it much more challenging -in a region between the target and the screening region and even in another region to the left of the target! The point we want to make here is that this screening is trivial to achieve if we put no bounds on the density ρ. For that purpose it suffices to consider a discretized version of the problem f (x) = 1 x + grid  figure 17. Of course this only works because the density is unconstrained here otherwise the amount of possible screening is limited. In the screening region between a and b the function we get is pretty huge and wild. In other words, the invisibility cloak is working hard so that spectators in the blue screened regions see nothing.

E.2 Screening in our setup
The phenomenon of screening happens in our setup for m 2 < m 1 / √ 2. In this case, the optimal bounds on the couplings g 112 and g 222 are just the same as the ones obtained from the single amplitude analysis. In fact, the off-diagonal amplitudes M 11→22 that saturate these bounds vanish in the physical region. To understand how this is possible it is convenient to The discontinuities in this region can not be huge. Hence the discontinuities in this region can not be huge either first understand why the bounds improve in the region m 1 / √ 2 < m 2 < m 1 .
For m 1 / √ 2 < m 2 < m 1 , the amplitude M 11→22 can not vanish generically so the diagonal amplitude M 22→22 can not saturate unitarity since some production of 11 is inevitable. To see this more precisely note that M 11→22 has poles and an extended unitarity region where the discontinuity is bounded by the diagonal M 11→11 and M 22→22 components as in (89). On the one hand, M 22→22 is bounded by unitarity because this amplitude has no extended unitarity region. On the other hand, the discontinuity of the 11 → 11 amplitude is positive in the extended unitarity region and therefore is bounded by unitarity in the physical region. This is depicted in figure 18. To summarize: we see that some screening is possible but it can not lead to a vanishing 11 → 22 amplitude in the physical region and thus to unitarity saturation for 11 → 11. This is why the multiple amplitude analysis had to improve the bound obtained from a purely diagonal analysis.
This same analysis also explains why for m 2 < m 1 √ 2 the diagonal bound is optimal in our setup. This is because in this range we have a collision between the s-channel and t-channel cuts corresponding to intermediate production of two of the lightest particles in the scattering of the heaviest particle 11 → 11, see figure 19a. The s-and t-channel discontinuities can now be huge as long as they cancel each other and do not lead to a violation of unitarity for the 11 → 11 component in the physical region. If they are unbounded, then there is a region in the 11 → 22 non-diagonal component where this amplitude can also be unbounded. (Note that if the cuts overlap then the imaginary part in the right hand side of (22) should be understood as the s-channel discontinuity.) If the amplitude can be huge with both signs in a finite segment then it can screen as illustrated in the previous section and can thus  Figure 19: The s-and t-channel discontinuities come with opposite signs. As such, when there is an overlapping region as in the 11 → 11 amplitude represented at the top, they can both be very large as long as their sum remains bounded and does not lead to violation of unitarity in the physical region represented by the black solid lines at the top. Furthermore the discontinuity of M 11→22 (at the bottom) in that same kinematical region is bounded by the ImM 11→11 and has no definite sign. Therefore, in the region where ImM 11→11 is unbounded, the amplitude M 11→22 can take advantage of screenning. kill 11 → 22 in the physical region (the backward 12 → 12 amplitude, related by analytic continuation to the 11 → 22 process can also be killed of course). In other words, for m 2 < m 1 √ 2 we can set to zero all non-diagonal processes without violating any of our physical constraints! Therefore the full numerical plots are expected to coincide with the analytical diagonal bounds. This is indeed what we observed in our numerics.
Would be great to find a way to improve our bounds for m 2 < m 1 For m 2 > √ 2m 1 there is a similar screening phenomena that occurs. Furthermore, in this mass range there are other Landau singularities known in two dimensions as higher pole Coleman-Thun singularities [6]. Indeed, if the mass m 2 > √ 2m 1 an on-shell diagram as in figure 13 will produce a double pole. Its residue, as seen in the figure, is governed by the coupling (to the fourth power) and by the 2 → 2 on-shell S-matrix of the lightest particle. These are all objects which we are already manipulating and it should thus be  Figure 20: We can realize the screening mechanism with a large number of resonances with masses m 2 A ∈ [4m 2 2 , 4m 2 1 − 4m 2 2 ] in crossing related pairs so that their contribution is moderate in the diagonal channels due to s-and t-channel pair-wise cancelations while, at the same time, having the potential to screen the non-diagonal processes where no t-channel poles show up (since the resonances are taken to be Z 2 even) and where the s-channel poles can have arbitrary sign as this is not a reflection symmetric process.
possible to tame these singularities if we properly understand how to deal with the inherent non-linearities. We look forward to reporting on this interesting problem in the future.

E.3 Multiple resonance toy model
One might wonder if the screening mechanism is a numerical artifact or could actually be realized in a reasonable QFT. Here we provide an example pointing towards the later provided we accept some fine tuning.
Consider a theory where m 1 > √ 2m 2 . The dangerous screening region is the region where the s and t channel cuts overlap in the 11 → 11 amplitude, i.e. for s ∈ [4m 2 2 , 4m 2 1 − 4m 2 2 ] as described in the previous section. Suppose we have many extra Z 2 even particles m 3 , m 4 , . . . , m 2N with m 2 A in that screening region range 32 and suppose further that for each particle with mass squaredm 2 A in that region there is a particle with mass squared where is a small parameter. Assume further that the couplings scale with this small parameter as while g 12A = 0 since the extra particles are even. Finally, assume that the couplings g 11A for two particles related as in (105) are the same up to small corrections. In this scenario, several interesting things might happen, including screening: • The particles would not appear in the 22 → 22 channel since they would come as poles with residues of order 2 → 0 or in the 12 → 12 channel since they are Z 2 even.
• The particles would appear in the 11 → 11 amplitude. Each particle contributes a huge amount since each coupling is of order 1/ 2 . However, because of the condition (105), for each s-channel pole there is a corresponding nearby t-channel pole which nearly cancels it, leading to a finite O( 0 ) result, see figure 20a. If the mass degeneracy if very tiny the couplings could be huge and still lead to a very good cancelation, compatible with unitarity for this 11 → 11 amplitude.
• The particles would also appear as s-channel poles in the 11 → 22 amplitude since g 11A g 22A = O( 0 ), see figure 20b. Note that this product could be very large and can take any sign. Hence it could lead to screening if N is large as explained in the toy example in section E.1. This screening could then lead to 11 → 22 being very small in the physical region.
Of course we could ask: "who ordered these extra particles?" No one did but since it is a priori consistent to add them, the numerics will take advantage of them and add them whenever useful for the optimization goal.
F Solvable Points at m 1 = m 2 For m 1 = m 2 it is sometimes possible to rotate the one-particle basis as so that the off-diagonal amplitudes have no poles. If that is the case, we can consistently set these amplitudes to zero and allow for the diagonal processes to saturate unitarity. For the spectrum considered in this paper, the poles terms in the Z 2 basis are (108) A straightforward brute force analysis in Mathematica shows that a change of basis as in equation (107) can diagonalise (108) only if in which case the rotated poles terms become of the form Diag where now the poles terms reduce to Diag It is then straightforward to apply the maximum modulus principle as in section 1 to obtain bounds the optimal couplings. For the (109) scenario, we conclude that a and e are fixed while b = c = d = 0, since the diagonal processes must saturate unitarity. This solution corresponds to the S-matrix of the 3-states Potts model at T = T c . Note that in the rotated basis we no longer have s ↔ t symmetry for diagonal processes in this case, since the external particles no longer diagonalise the charge conjugation operator, see appendix G. On the other hand, for (110) we conclude only that a is fixed and b = c = 0 (note that after changing basis b is no longer related to d by crossing), while d and e correspond to a zero modes that do not affect the optimal coupling. A particular choice of d and e lead to the hyperbolic limit of the elliptic deformation of the supersymmetric Sine-Gordon model, discussed in section 3.3.

G 3-state Potts field theory
The 3-state Potts model in two dimensions has a continuous phase transition described by a non-diagonal minimal model with central charge c = 4/5 and a global permutation symmetry S 3 . This conformal field theory (CFT) contains 3 relevant scalar operators invariant under Z 2 ⊂ S 3 (see [25] for a nice introduction to this topic). This allows us to define a family of Z 2 symmetric QFTs with action where τ , h and h are relevant couplings and we used the notation of [24]. The scaling dimensions of the relevant operators are ∆ = 4 5 , ∆ σ = 2 15 and ∆ Ω = 4 3 . The purely thermal deformation (τ = 0 and h = h = 0) preserves the S 3 symmetry and leads to an integrable QFT. This theory has only two stable particles with the same mass m transforming as a doublet of S 3 . These particles are usually described in a basis |A , |A † where the Z 2 acts as charge conjugation C|A = |A † with C 2 = 1 [19,20]. In this basis, the S-matrix is diagonal, i.e. S AA→A † A † = S Backward AA † →AA † = 0 and where we used the rapidity θ to parametrize the Mandelstam invariant s = 4m 2 cosh 2 θ 2 . Notice that, in this basis, the Yang-Baxter equations are trivially satisfied.
In this paper, we work in the eigenbasis of the Z 2 global symmetry generated by charge conjugation, In this basis, we find Using equation (76) we can obtain the expressions (36) for the connected scattering amplitudes that maximize g 2 112 for m 1 = m 2 and g 222 = −g 112 (point A in figure 9). The magnetic deformations h and h in (112) preserve Z 2 and therefore must be compatible with our bounds, at least for small magnetic deformations hτ − 14 1 that do not give rise to more stable particles. In fact, the mass spectrum of these theories (with h = 0 and h = 0) has been studied in [24] using the Truncated Conformal Space Approach. The authors observed that the degeneracy between the two particles is lifted for h = 0. It would be interesting to study the cubic couplings and the S-matrices of this 2-parameter family of Z 2 symmetric QFTs and compare them to our bounds.

H Tricritical Ising (cusp)
The tricritical Ising model in two dimensions has a continuous phase transition described by a diagonal minimal model with central charge c = 7/10 and a Z 2 (spin flip) symmetry. This conformal field theory (CFT) contains 2 relevant scalar operators invariant under Z 2 (see [25] for a nice introduction to this topic). This allows us to define a family of Z 2 symmetric QFTs with action where τ , τ are relevant couplings. The scaling dimensions of the relevant operators are ∆ = 1 5 and ∆ = 6 5 . The purely thermal deformation (τ = 0 and τ = 0) leads to an integrable QFT. This theory has seven particles but only four of them have masses below the continuum of multi-particle states. The masses of these particles are given in table 1.
In this paper, we constrained the space of 2D S-matrices with Z 2 symmetry by requiring unitarity and analyticity for the two-to-two S-matrix elements involving only the two lightest particles {m 1 , m 2 } as external states. We restrained ourselves to study the subset of theories for which only m 1 and m 2 themselves appeared as bound states in these matrix elements. This excludes the tricritical Ising field theory from our analysis in the many body of this paper. However, one can easily relax this restriction. We shall not do a full numerical study of the multiple amplitude bootstrap in this more general setup. We will just derive the analytic bounds that follow from the amplitudes S 11→11 and S Forward 12→12 . For concreteness, we consider a theory with Z 2 symmetry and a mass spectrum as in the table 1. 33

Particle
Mass 2m cos 5π/18 ≈ 1.29 m even m 3 2m cos π/9 ≈ 1.88 m odd m 4 2m cos π/18 ≈ 1.97 m even + possibly extra particles with masses bigger than m 1 + m 2 .  . One question that could be asked is: what values for the pair (g 2 112 , g 2 114 ) are allowed by the unitarity constraint |S 11→11 | ≤ 1? To answer this question, it is useful to introduce the variables 33 Any masses within the range 0 < 4m 2 1 − m 2 4 < m 2 2 < 2m 2 1 and m 2 1 < 2m 2 1 + 2m 2 2 − m 2 3 < m 2 1 + m 2 2 would lead to qualitatively the same conclusions regarding single amplitude bounds as below but the precise locations of cusps and edges on the bounds does depend on the masses. If we deviate away from these mass constraints then we we would have to redo the analysis. This is the same discontinuous nature of the bounds already observed in [2], see e.g. figures 10 and 11 therein. where we now think of f as a function of w 4 by inverting equation 121. Since w(f, α) is an increasing function of f , to maximise f (w 4 (z 2 )) = g 2

112
(g max 112 ) 2 is equivalent to maximise g(w 4 (z 2 )). Moreover, since g is an automorphism of the disk, unitarity implies |g(w 4 )| ≤ 1 for w 4 on the unit disc. Finally, g(0) = 0. Now recall Schwartz Lemma: Lemma 1 Schwartz Lemma: Let D be the unit disk and g : D → D be a holomorphic map such that g(0) = 0 and |g(w)| ≤ 1 on D. Then |g(w)| ≤ |w|. Moreover, if the inequality is saturated for any non-zero point in D, then g(w) = aw with |a| = 1.
We conclude that under the extra constraint g 2 114 = α(g max 114 ) 2 , the maximal value for g 2 112 is given by the solution of the following algebraic equation on S 11→11 : After performing the equivalent exercise for the maximisation of g 2 114 under g 2 112 = β(g max 112 ) 2 , where β ∈ [0, 1], we obtain the allowed space of (g 2 112 , g 2 114 ) as depicted in the figure 23.
One can play the same game for the S Forward 12→12 matrix element, which contains bound state poles corresponding to particles m 1 and m 3 (and therefore constrain the space (g 2 112 , g 2 123 ). One can then combine both results into a 3D plot of allowed triplets (g 2 112 , g 2 114 , g 2 123 ) compatible with |S 11→11 | ≤ 1 and |S 12→12 | ≤ 1. This is figure 24.
It would be interesting to perform a multiple amplitudes analysis for this setup and explore the space of masses (m 1 , m 2 , m 3 , m 4 ) in the vicinity of the values of table 1. Notice that in such an analysis, the values in table 1 would be single out by the condition that the multiple amplitude bounds saturate the single amplitude bounds of figure 24 at the tip of the spear. That is because only for those particular masses can the multiple poles in the off-diagonal channels coincide and cancel (à la integrable bootstrap), allowing for the off-diagonal amplitudes to vanish and for the diagonal processes to saturate unitarity, as in the boundary of the yellow surface in figure 24. It would also be interesting to see if the sub-leading (non-integrable) deformation τ of the tricritical Ising model leads to any feature of the bounds.

I Numerical optimization as a SDP
As discussed in section 3.1, once we fix α = g 222 g 112 , our discretised ansatz for the amplitudes depends only on the variables η = {g 2 112 , C a→b , σ a→b (x i )}. To reduce the maximisation problem to an SDP, all we need to do is to write the extended unitarity constraint (27) as a semidefinite constraint linear on η, as in (33). The purpose of this appendix is to prove the equivalence of (27) and (33) or, explicitly, Proof of ⇒: This becomes the r.h.s of (127) if we pick v = −ρM w.
Proof of ⇐: where we used the r.h.s. of (127) in the second inequality.

J Elliptic Deformation
Bellow one can find the definition of the supersymmetric sine-Gordon elliptic deformation S-matrix in Mathematica friendly notation. Note that to compute the S-matrix for physical θ ∈ R, one must be careful and take the appropriate principal value around the singularity at θ = x in the integrand of IntR.

K Conformal computations K.1 Crossing symmetry in one dimension
We consider correlation functions of primary operators φ i (x i ) in a one-dimensional boundary field theory. The weight of the operator φ i will be denoted as ∆ i . We will work in the Euclidean theory on R ∪ {∞}. The conformal group is a two-fold cover of P SL(2, R); its elements act by the usual fractional linear transformations on the positions, with ad−bc = ±1. The elements with negative determinant involve the parity transformation x → −x. The Jacobian of this transformation is (ad − bc)(cx + d) −2 and the fields transform as with P φ ∈ {0, 1} dictated by the parity of φ. 34 In a correlation function we should remember that the parity operation also reverses the operator ordering.
In a suitable basis the two-point functions take the familiar form The two-point function of a parity odd operator is negative, so the associated norm φ|φ is positive. 35 The operator product expansion reads where we assume that x < 0. If we act with the parity operator on both sides then we find and therefore the reflected OPE coefficients between primaries are with σ 12k := P 1 + P 2 + P k mod 2 .
The previous relation in particular implies that parity odd operators can never appear in the OPE of two identical operators. We will work in a basis where two-point functions are diagonal. The structure of three-point functions then dictates that and so the OPE coefficients transform either in the trivial representation (if σ 12k = 0) or the sign representation (if σ 12k = 1) of the permutation group S 3 . Notice that, even if parity is broken, the cyclic symmetry is always preserved. This for example means that C 12k = C k12 , whereas C 12k and C 21k are not always the same.
where we added an absolute value sign so the expression is unambiguous also for negative values of its argument.
Let us briefly discuss operator ordering. Correlation functions with operator orderings that are cyclic permutations of each other are directly related, as follows from covariance under the (orientation-preserving) inversion x i → −1/x i . Furthermore, parity symmetry dictates that To see the complications that arise if we just swap two adjacent operators, let us swap operators 1 and 2. A simple relabeling leads to the block decomposition and assuming a parity invariant theory this is equal to where we used a standard hypergeometric transformation formula, valid for x < 1. The factor (−1) P k and the absolute value sign in the definition (141) imply that G 2134 and G 1234 are not, in general, related in an obvious manner. The symmetries and non-symmetries altogether leave us with three independent four-point functions from which the others follow. We can take these to be 2134 , 1234 , and 1324 , which respectively correspond to x < 0, 0 < x < 1, and 1 < x. We will here be interested in just the second of these correlators.
With the ordering fixed there are only two OPE channels. For the 1234 ordering we gave the first one above; the crossed channel OPE reads with y = x 23 x 41 /x 24 x 31 = 1 − x. Crossing symmetry then takes the form 36 Let us consider all correlation functions of two parity-even operators φ 1 and φ 2 . We then 36 One foolproof way to obtain this expression is to relabel the operators in the original expression (138) and then use a conformal transformation to relate 2341 to 1234 . To verify that directly fusing operators 2 and 3 together in (138) gives the same OPE limit requires that C 1p4 = C p41 which we proved previously. find, in a diagonal operator basis, the following set of non-trivial crossing equations: where the operators labeled k are parity even, whereas the operators labeled p can be either parity even or parity odd.

K.2 Transition to convex optimization
In matrix form, the I'th crossing symmetry equation can be written as: and with We can act with a functional on each equation, and then add all of them. This yields where F I p (without an argument) is shorthand for the function of ∆ p obtained by acting with the corresponding component of the functional.
If we single out the identity (with C 110 = C 220 = 1 and C 120 = 0) and the external operators from the sums then we obtain For a feasibility study we can normalize the functionals on the unit operators, giving and furthermore demand positive semidefiniteness of the three square matrices listed above. We can also get bounds on products of OPE coefficients. For example, if we set and then maximize/minimize F 1 0 + 2F 6 0 + F 2 0 , we obtain a lower/upper bound on C 2 111 . More precisely, if we extremize and the result is positive then crossing cannot be solved. If the result is negative then the absolute value of the result is our upper (for minimization) or lower (for maximization) bound.

K.3 Setup with Z 2 symmetry
In the previous section we discussed the general one-dimensional conformal bootstrap analysis for two operators. Let us now specialize to the case discussed in the main text, so we consider a QFT in AdS 2 with a Z 2 symmetry and only two stable parity even particles; a Z 2 odd one created by an operator φ 1 and a Z 2 even one created by an operator φ 2 .
With this additional symmetry F 3 = F 4 = 0 automatically which leaves four non-trivial correlation functions of φ 1 and φ 2 . In these correlators we should also label the internal operators with an even/odd quantum number depending on the OPE channel in which they appear. Our notation above is already adapted to this situation: in equation (147) the operators labeled k are necessarily parity and Z 2 even, whereas the operators labeled p are parity even or odd but always Z 2 odd. (Operators that are Z 2 even but parity odd do not feature in this set of correlation functions.) With the exception of φ 1 and φ 2 themselves we take their dimensions to lie above the following 'gap' values: gap index in (156) even even min(2∆ 1 , 2∆ 2 ) k odd even ∆ 1 + ∆ 2 p, P p = 1 odd odd ∆ 1 + ∆ 2 p, P p = −1 These gaps are precisely the two-particle values for a QFT in AdS, reflecting our assumption that there are no further stable single-particle states.
Going then through the same logic as before we write the crossing equations in matrix form as p,Pp=1 with (note some redefinitions with respect to the previous formulae) and with After acting with the functional and singling out the important operators again we find If we want to find bounds in the (C 112 , C 222 ) plane we can set C 222 = βC 112 and set a normalization With this normalization our problem takes the form of a semidefinite program and we can proceed using the well-known numerical bootstrap methods of [4,12] and the specialized solver of [14]. Further technical details are available upon request from the authors.

K.4 Functionals and derivative combinations
As in previous works, our choice of functionals are linear combinations of derivatives of the crossing equations at z = 1/2, so our functionals α can be written as  Figure 25: Bound on the maximal coupling (without any extrapolations) as a function of Λ 5 = Λ 6 for fixed Λ 1 = Λ 2 = 80. We see that the multi-correlator bound does not improve over the singlecorrelator bound (in black) for a range of values of Λ 5 = Λ 6 , and this may lead one to believe that no improvement is possible whatsoever. On the other hand, increasing Λ 5 = Λ 6 above the natural value given by Λ 1 = Λ 2 does not seem to lead to further improvements for the range of values that we tested. The case shown here has ∆ 1 = ∆ 2 = 6.58895 and g 222 /g 112 = 1, but other cases look similar: we start with a plateau, then a kink (not necessarily at Λ 1 /2) marks the start of a downward trajectory (which is not necessarily this linear), and then another kink at Λ 5 = Λ 1 leads to a second plateau where the bound is constant.
The coefficients α n should be thought of as the main degrees of freedom to be fixed by the optimization procedure. Furthermore, in a multi-correlator study we can act with a different functional on each of the crossing equations in (147), so the degrees of freedom are α I n , I ∈ {1, 2, 5, 6}, n ∈ {0, 1, 2, . . . Λ} .
where, as before, I denotes the different crossing equations and we used that the third and fourth equation are identically zero by our Z 2 symmetry assumptions.
Of course the functions that are odd around z = 1/2, like the first and second crossing equation in (147), contribute only about Λ/2 non-trivial degrees of freedom since the α n for odd n are meaningless (and in fact should be set to zero to prevent numerical instabilities). This reduces the scaling of the number of components to 3Λ rather than 4Λ.
In our earliest attempts we thought it unfair to first two crossing equations that they could only contribute half as many functional components as the fifth and sixth equation. Therefore we decided to cut off the sum in (161) at Λ/2 for the fifth and sixth equation but at Λ for the first and second equation. Such an egalitarian approach, however, would have led to completely different and incorrect results. To illustrate this we plot in figure 25 the upper bound on the coupling as a function of Λ 5 = Λ 6 , i.e. the cutoffs for the fifth and sixth equation, whilst holding fixed every other parameter, including the cutoffs Λ 1 = Λ 2 for the first and second equation. Clearly for Λ 5 = Λ 1 /2 the multi-correlator bound offers exactly no improvement over the single-correlator bound. So, if we had continued to work with the egalitarian cutoffs then we would erroneously conclude that no improvement would have been possible over the single-correlator bound! Only when increasing Λ 5 beyond Λ 1 /2 do we begin to see a gradual improvement, which stops as abruptly as it started at Λ 5 = Λ 1 .
Notice that the plot is drawn for the data point in figure 11, which has the special property that the multi-correlator bound turns out to be exactly half that of the single-correlator bound. We however observed very similar behavior, including the kinks and stabilization, also for other data points where the final multi-correlator bound is completely non-trivial.

K.5 Extrapolations
An example of the extrapolation procedure outlined in section 4.1 is shown in figure 26. One important subtlety not mentioned in the main text is that we extrapolate the log-ratio of the multi-correlator and the single-correlator result. That is, for every multi-correlator optimization run we also ran a single-correlator optimization run (for either 1111 or 2222 ) and so we get raw data that we can denote (g 2 112 ) max,multi [µ, α, ∆ 1 , Λ] and (g 2 112 ) max,single [µ, α, ∆ 1 , Λ]. We found that direct extrapolation of the multi-correlator bounds led to a relatively large dependence on our fitting procedure, whereas extrapolation of the log-ratio log(g 2 112 ) max,multi [µ, α, ∆ 1 , Λ] − log(g 2 112 ) max,single [µ, α, ∆ 1 , Λ] could, as shown, be done with relatively low-degree fits. Since we know that the singlecorrelator bounds match the analytic single-amplitude bounds with large accuracy [1], we can add (the logarithm of) this known answer to our extrapolated log-ratio to obtain a better estimate of the flat-space multi-correlator bound.
The first 8 plots show the raw data and subsequent extrapolations to infinite Λ. The three curves correspond to fits with a polynomial in Λ −1 of degree 3 (in blue), degree 4 (in orange, mostly coinciding with blue), and degree 2 (in green). In the fits we did not include the (three) data points with Λ < 60. We observe a rather small difference between the different extrapolations, and we have checked that these fits give good predictions for the high Λ raw data points if we exclude one or more of them by hand.
The final plot in figure 26 collects all the extrapolated points to infinite Λ. The data points line up nicely, providing evidence for a small non-systematic error in our first extrapolation. We have fitted a linear function in ∆ −1 to the degree 3 (in blue) points, leading to the extrapolated value of −1.446 for the log-ratio of this data point. From this plot it is clear that there is no meaningful difference if we had extrapolated the degree 4 (in orange) points instead.
We employed exactly the same extrapolation procedure, including the choice of the degree of the fitting functions, for all the other data points shown in the main text.  figure 12. Vertically we plot the logarithm of the maximal coupling, but with the logarithm of the single-correlator bound subtracted, so δ(log(g 2 max )) = log(g 2 max,multi ) − log(g 2 max,single ). In the first 8 plots we show our raw data in black, and the curves correspond to three different extrapolations to infinite Λ. In the final plot we collected the Λ → ∞ extrapolations and the single red line represents our ∆ 1 → ∞ extrapolation. Our final answer gives δ(log(g 2 max )) ≈ −1.446 in the flat-space limit, and adding the single-correlator bound 3.673 gives the 2.226 plotted in figure 12. Figure 27: An on-shell scattering process must lie inside the convex hull defined by the vertices attached to external lines. This is due to momentum conservation: an on-shell internal particle that wandered outside the convex hull would never be able to move back inside, since there would be no external momentum available to kick it back in. There are 3 possible convex hulls for 2 → 2 amplitudes: a quadrilateral, a triangle, or a line.

L Landau Singularities in 12 → 12 Scattering
Landau singularities are associated to diagrams representing particle interactions with all lines on-shell and momentum conservation at each vertex [6,17]. We claimed in the text that for the forward 12 → 12 scattering there ought to be no new such singularities (to be added to the bound state poles already there) whereas for the 12 → 12 backward component we claimed that when m 2 > √ 2m 1 we do find such new on-shell processes. It is easy to convince oneself of that with a few pictures. For that purpose, we adapt here some beautiful discussion in chapter 18 of the book [35] by Bjorken and Drell, translating it to our two dimensional case of interest.
Start with a diagram representing a two-to-two on-shell process. Each of the four external particles eventually encounter a vertex (several can meet at the same vertex). Consider those vertices containing an external line. They define a convex hull which can be a quadrilateral, a triangle or a line as shown in figure 27. Next we draw all other internal lines to complete the Landau diagram. The first important claim is that all those lines must lie inside the convex hulls just defined. That is simply because of momentum conservation: if they got out of the convex hull there would be no external particle to kick them back in (and the diagram would never close)! With this in mind we can go on to analyse the three possible convex hulls in turn.
Consider first the quadrilateral case and assume first that at each of the four external vertices we have a cubic vertex. Then, momentum conservation and on-shellness constrains the angles at those cubic vertices 37 . For example, for the case at hand with two particles 2⇡/3 x x = 2 arccos m 2 2m 1 < ⇡ 2 if m 2 > p 2m 1 Figure 28: On-shellness and energy-momentum conservation fix angles in cubic vertices. Importantly, when m 2 > m 1 √ 2 there are acute angles in cubic couplings and the analysis of Landau singularities becomes more intricate.
of different masses we have the angles in figure 28. For m 2 < √ 2m 1 we realize right away that the angles are obtuse so it is generically impossible to form a quadrilateral! What about using other vertices when particles first interact? Well, that would be even worse as the total opening angles would be even greater in that case. If m 2 > √ 2m 1 then one of the angles is smaller than π/2 so we have a better chance of finding extra singularities and, indeed, we can sometimes form such diagrams, as the one in figure 13, but not if two external particles are of type 2 and two particles are of type 1 as illustrated in figure 29(a). Hence, there are no "quadrilateral convex hull" Landau diagrams for 12 → 12 scattering.
What about "triangular convex hull" singularities? When m 2 > √ 2m 1 , as explained in figure 29(b), those are indeed present in the backward component, but not in the forward one. So the backward amplitude should have extra Landau poles but the forward amplitude should not. As such, the bound on g 2 112 derived from the forward component should hold even for m 2 > √ 2m 1 . As discussed in the main text, this bound is not captured by the QFT in AdS bootstrap, see figure 12.
We did not discuss the case where the convex hull is the line: those are just the usual singularities such as the bound state poles and all the production cuts which open for multiparticles at rest, and thus moving parallel to each other, along the convex hull line. from their absence in the euclidean region after a causality argument, as detailed in [35].  Figure 29: (a) There are no possible "quadrilateral convex hull" Landau singularities in the 12 → 12 amplitude even when m 2 > √ 2m 1 . To understand this, consider a vertex with an external leg attached. The total opening between internal edges on such vertex must be greater than x (y) for an external even (odd) leg. Hence, for the 12 → 12 process, the total internal angles at the four external vertices must be greater than 2(x + y) ≥ 2π, so that it is inconsistent to have a closed singular diagram inside the convex hull. (b) The same is true for "triangular convex hull" diagrams unless two odd particles meet at the same external vertex, as on the bottom right diagram -as in other cases the total opening at external vertices would be greater than π. For the forward component there would be no momentum transfer through such vertex, so that by on-shellness and energy-momentum conservation the internal opening at this vertex would be at least π, once again making it impossible to have a singular Landau diagram. For the backward component, on the other hand, the momentum transfer through the vertex, and hence the opening angle, can be arbitrary. In turn, this amplitude will have extra Landau singularities for m 2 > √ 2m 1 .