S-matrix bootstrap in 3+1 dimensions: regularization and dual convex problem

The S-matrix bootstrap maps out the space of S-matrices allowed by analyticity, crossing, unitarity, and other constraints. For the $2\rightarrow 2$ scattering matrix $S_{2\rightarrow 2}$ such space is an infinite dimensional convex space whose boundary can be determined by maximizing linear functionals. On the boundary interesting theories can be found, many times at vertices of the space. Here we consider $3+1$ dimensional theories and focus on the equivalent dual convex minimization problem that provides strict upper bounds for the regularized primal problem and has interesting practical and physical advantages over the primal problem. Its variables are dual partial waves $k_\ell(s)$ that are free variables, namely they do not have to obey any crossing, unitarity or other constraints. Nevertheless they are directly related to the partial waves $f_\ell(s)$, for which all crossing, unitarity and symmetry properties result from the minimization. Numerically, it requires only a few dual partial waves, much as one wants to possibly match experimental results. We consider the case of scalar fields which is related to pion physics.


Introduction and Summary
Recently, new insights were found on the old idea [1,2] of solving the S-matrix directly from its analytic structure, symmetries, crossing and unitarity. Although these constraints do not uniquely determine the S-matrix, combining them with a convex maximization program leads to bounds on couplings between particles and their bound states. Furthermore, it was observed that well-known theories such as a subsector of the sine-Gordon model saturate such bounds [3,4]. Similar programs can be carried out for 3+1 dimensional theories [5], and multiple amplitudes [6]. Even in the case where the theory does not have bound states similar ideas can be applied. An interesting example is the 1+1 d O(N ) non-linear sigma model [7] which describes an asymptotically free theory with a dynamically developed mass gap in the infrared. It was exactly solved in [8,9] and more recently revisited with the S-matrix bootstrap approach in [10][11][12]. In particular, in [10] it was argued that one should focus on mapping out the space of allowed S-matrices under the given constraints. From that perspective, the O(N ) non-linear sigma model lies at a special point -a vertex -in the space of allowed theories, and maximizing a linear functional in a convex space generically leads to a vertex. This idea was made more manifest later in [13] where a section of the space was plotted with a clear vertex at the NLSM. Further work on other models [14] showed that sometimes full regions of the boundary correspond to interesting theories if such theories have free parameters. More recently, the study of reflection matrices [15] showed that vertices appear at points where resonances, namely poles in the second sheet cross into the physical region becoming bound states. Since bound states are not allowed, the S-matrix changes form and a discontinuity in the derivative appears.
In this work we concentrate on the dual problem proposed in [13] and further studied in [32,33]. Thinking in terms of the allowed space of S-matrices, the primal approach provides an interior region that becomes larger as the numerics is improved, eventually approaching the full allowed space. On the other hand, the dual problem provides an exterior region that shrinks onto the allowed space of S-matrices as the numerics improves. This is particularly useful in 3 + 1 dimensions where the problem is numerically difficult, since it allows to bracket the boundary between the interior and the exterior. Moreover, we observed that the dual problem requires only a few partial waves much in the same way as the results of experiments where only a few partial waves are available. In summary, the dual problem as described here provides a complementary and in many ways very useful approach to the problem of mapping out the allowed space of S-matrices.
The construction of the dual problem requires several steps. First we need to define the primal problem. For that we use the Mandelstam representation which assumes maximal analyticity and enforces crossing properties for the amplitudes. The unitarity constraints are then imposed on the partial waves with definite angular momentum and isospin. It turns out, however, that the price to pay for a parameterization that has manifest analyticity and crossing is that, it contains fluctuations of the variables that barely modify the partial waves. Therefore, roughly speaking, they are not constrained by unitarity and they do not contribute to the physical partial waves. We solved this problem by regularizing the primal problem, namely putting a bound -a regulator -on the variation of the primal variables. A large enough regulator allows for the amplitudes we need, while suppressing unwanted large fluctuations of irrelevant variables. In addition, the regularization also leads to a well-defined dual problem. Alternatively, one could parameterize the partial waves as analytic functions, but in that case crossing is difficult to impose since it mixes all angular momenta. We have not explored this last possibility further.
Once the primal problem is regularized, it is straightforward to define the dual problem using the tools of conic optimization. This leads to upper bounds on the space of S-matrices allowed by the constraints. As an example, we can see in figure 1 how the maximum (blue curve) is bracketed between the red curves (primal) and the green curves (dual). In this case, the blue curve was obtained by assuming the existence of a pole at threshold whereas the red and green curves did not make such assumption since, in general, we do not know the pole structure of the optimal amplitude. Furthermore, as in 2d, it turns out that one can derive the dual problem just using analytic properties of the amplitudes. For that purpose we introduce what we call the generalized double dispersion relations (4.22) that allows us to relate the values of the amplitude at an unphysical point to the values in the physical region. We do this by defining a dual amplitude (see for example (4.23)) . Both were run without assuming the existence of a pole at threshold. The blue line represents the primal problem once a pole at threshold is included. The black diamonds represent the pioneering results of Lopez and Mennessier [34]. The inset shows more clearly the results near the symmetric point which is the minimum of the blue curve where π 4 F ( 4 3 , 4 3 , 4 3 ) ≃ 2.6613. The factor π 4 is to match the normalization in the literature, e.g. [5]. through a double dispersion relation with support in the physical region. After minimization, the dual amplitude contains all the information necessary to extract the partial waves in the energy range where they saturate unitarity. The construction allows for the possibility of unitarity unsaturation but we did not find that property in the numerical solutions we obtained. It will be interesting to explore this idea further. This paper is organized as follows. In the next section, we review how conic optimization can be used to map out the space of allowed S-matrices. In the following section, we describe first how regularizing the primal problem eliminates irrelevant fluctuations of the primal variables while at the same time provides a well-defined dual problem. In section 4 we rederive the dual problem using arguments from complex analysis, where we introduce the generalized dispersion relation (4.22) and the dual amplitudes. These are analytic functions of two variables that contain the same physical information as the amplitude to the extent that the partial waves can be extracted from it where unitarity constraints are saturated. Its double spectral density has support in the physical region instead of the Mandelstam region. Although the subsequent numerical tests focus on the single flavor case, we have also included the formulation with O(N ) global symmetry in sections 2.2, 3.4 and 4.5, which may be skipped for readers who want to focus on the general logic. In section 5 we describe two different numerical implementations of the primal and dual problems, one based on interpolation points and the other using a basis of functions. We illustrate these methods by doing some simple numerical computations in section 6 and leave a full numerical exploration of the method for upcoming work. We give our conclusions in section 7 and collect a few useful formulas in an appendix.

Physical problem as a primal problem in conic optimization
In this paper we consider 2 → 2 scattering of scalar particles because of its simplicity. Although we do not assume any particular identity for the particles or form of the interaction, it is useful to describe this process in the language of pion scattering because it is a very well-studied example. This type of system has already been studied with the S-matrix bootstrap in [5,16,17,21,22]. In this paper we concentrate on the formulation of the dual problem.

Single pion scattering
The simplest example is single pion scattering, namely a scalar particle of mass m 2 that we set to m 2 = 1. The amplitude for the process is a function of the Mandelstam variables s = ( Its analytic and crossing properties are captured by the Mandelstam representation that can be written as 1 with the kernels The amplitude F (s, t, u) has double jumps on the Mandelstam regions depicted in red in fig.2. The amplitude is crossing symmetric under s ↔ t ↔ u, and only the symmetric part of ρ contributes so we can take ρ(x, y) = ρ(y, x).
The physical partial waves f (s) are defined as which are non-zero only for even . In the integral (2.4), t and u are to be considered as functions of µ and s according to From now on, when integrating over Legendre polynomials of µ we always assume t = t(µ, s), u = u(µ, s) according to (2.5a), (2.5b). This implies that, to compute the partial waves, we only need to know the amplitude in the physical region depicted in green in fig.2 (or any of its crossing conjugates since the amplitude is crossing symmetric). It is convenient to define rescaled partial waves h (s) and the S-matrices S (s): where δ are the (possibly complex) phase-shifts. The unitarity condition is Notice that in the physical region s > 4, 4 − s < t < 0, eq.(2.4) can be inverted as where the sum is over even and s + indicates we evaluate the amplitude on the real axis above the cut [4, ∞). Given the general constraints of crossing, unitarity and analyticity as given by the Mandelstam representation, we want to map out the space of allowed amplitudes. One way to do that is to start from a point, e.g. the free theory S = 1, pick a direction in the space of S-matrices and move along that direction until one of the constraints is violated. This is equivalent to maximizing a linear functional. The space of allowed S-matrices is of course infinite dimensional, but we can concentrate in a lower dimensional subspace. Since unitarity is a convex constraint and we impose linear conditions, the problem is a conic maximization problem [35] that can be solved with standard methods [36][37][38]. A standard functional that we can choose is where (s 0 , t 0 , u 0 ) is a point in the Mandelstam triangle (see fig. 2) where the amplitude is real. So the problem is to find (f 0 , σ(x), ρ(x, y) = ρ(y, x)) that maximizes F P under the unitarity constraint. In convex maximization, an important role is played by the dual problem that we define later and is the main subject of our paper.

O(N ) case
Now we consider a scalar theory with N species of particles and O(N ) symmetry. For the actual pions N = 3. The 2 → 2 amplitude is now F ab,cd = A(s, t, u)δ ab δ cd + A(t, s, u)δ ac δ bd + A(u, t, s)δ ad δ bc (2.12) with A(s, t, u) = A(s, u, t) and a, b, c, d . . . = 1 . . . N . The expression satisfies crossing and isospin symmetry but we still need to impose the unitarity constraint. For a well-defined isospin (I = 0, 1, 2) in the s-channel we get: The Mandelstam representation for the scattering amplitude, including subtractions, is now (α = 1, 2) The amplitude A(s, t, u) has double jumps in the Mandelstam regions depicted in red in fig. 2. The variables now are σ α=1,2 and ρ α=1,2 (x, y) where ρ 2 (x, y) = ρ 2 (y, x) and the rescaled partial waves h I (s) are given by where the even 's are non-vanishing for I = 0, 2 and the odd ones are nonvanishing for I = 1. Keep in mind that t and u are functions of µ, s according to (2.5a), (2.5b). We can also compute a Mandelstam representation for the isospin channels: The kernels follow the pattern in (2.14): α (s, t, u; x) = K α (t, s, u; x) + K α (u, t, s; x) (2.20c) and the same for the double dispersion kernels. The unitarity constraint is: Under these conditions we want to maximize a linear functional F P of the amplitude. A standard example is for a point (s 0 , t 0 , u 0 ) inside the Mandelstam triangle where the amplitude is real and n I are constants. We can also take a general linear functional of the partial waves for some coefficients c I, of which only a finite number we choose to be nonzero. Such functional has a feasible dual problem without the need for regularization, and it can be useful if we believe that the theory we are interested in maximizes certain partial waves.

Primal problem
Both problems can be summarized in the language of convex optimization by using a compact notation as follows. Define a set of variables α n where n denotes a discrete or continuous index. From (2.2) and (2.15), we can read A sum over n is understood as integration when the label is continuous. For example, for the single pion and using (2.2), the primal functional (2.10) can be written as Comparing with (2.2), we can extract the coefficients a n for the maximization of F P = F (s 0 , t 0 , u 0 ): Instead, if we wanted to minimize F P = F (s 0 , t 0 , u 0 ), this is equivalent to maximizing F P = −F (s 0 , t 0 , u 0 ) and therefore we would take the a n in (2.26) with the opposite sign. Other choices of a n would correspond to other functionals, for example those in (2.23). Thus, for any choice of functional, the primal problem can be stated as max αn F P = a n α n , where for the single pion I = 0, even, and for the O(N ) case I = 0, 2, even and I = 1, odd. In (2.27), we have used the compact notation h I n (s) which gives the partial waves as linear combinations of the variables α n . For example, in the case of single pion, we have the following h I n (s) for the variables (2.24) where we used the integrals in the appendix to write the partial waves in terms of the Legendre functions Q as in the Froissart-Gribov form (see e.g.

Regularized primal problem and its dual
In this section we discuss the need and also practicality of regularizing the primal problem by putting an upper bound on the norm of the double spectral density ρ(x, y) in (2.2) or (2.15). This need arises from two different but related issues. One is that highly oscillating functions added to ρ(x, y) make little difference in the physical amplitude, the other is that the dual problem is weakly infeasible 2 . These problems are not related to the Mandelstam representation but to the fact that the unitarity condition is evaluated in the physical region s > 4, 4 − s < t < 0 which is not in the support of the double spectral density, namely the boundary of the analyticity region. For example, we can alternatively do as in [5] and map the analyticity region to unit disks for each Mandelstam variable using with z s ≤ 1 and the same for t and u. Then, in the physical region z s = 1 but z t < 1, z u < 1. If the single pion amplitude is written as a polynomial of some given degree M as then the coefficients ρ nm with both n and m large affect the value of the polynomial very little in the physical region where z t < 1 and z u < 1. Therefore the unitarity condition does not constrain those coefficients because they do not affect the physical amplitudes and therefore they can become large. However since they do not affect the physical amplitudes they are not really relevant. A simple bound such as ρ nm < M reg for some regulator M reg solves that problem by not allowing the coefficients to grow. Related to this, the dual problem is generically not feasible (feasibility depends on the functional chosen) unless we regulate the primal problem in which case the dual is feasible for any primal functional. Let us now consider these issues in more detail using the Mandelstam representation.

Fredholm equations and regularization
A single or double dispersion relation defines the amplitude through a linear operator Given an amplitude F , determining a spectral density ρ is equivalent to solving a Fredholm equation of the first kind with kernel K (see e.g. [39]). For this problem to be well-posed we should require that for any F there exists a unique ρ and that the solution depends continuously on F . It is well-known [39] that this is not true when the kernel K is continuous, the problem of finding ρ is then ill-posed 3 . For example, adding a highly oscillating function to ρ, even with a large coefficient, does not modify F . In the case at hand, it means that the coefficients of such highly oscillating variations of ρ cannot be determined by maximization since they barely modify the functional. On the other hand they also barely modify the partial waves and therefore have no physical interest and should be suppressed. For the case in hand, recall the definition of the amplitude (2.2). The single dispersion part has a singular kernel and therefore the problem is well-posed. In the double dispersion relations, however, if we look for example at the part of the kernel that behaves as it is not singular for physical values of t (i.e. 4 − s < t < 0) because in the region of integration in (2.2) we have y > 4 and then y ≠ t. The same is true for other parts of the kernel in (2.2).
In the case of Fredholm equations of the first kind, the method to make the problem well-posed is to use the so-called Tikhonov regularization. Instead of solving eq. (3.3), we minimize the functional where ⋅ is a norm 4 . The parameter γ provides a regularization. Ideally, the norm should be such that the variations of ρ that affect the partial waves the most are kept and others are suppressed. In practice, the norms that we tried numerically, including the simple square norm, give similar answers. In the following we discuss the dual problem that, if formulated in the usual way, leads to a linear constraint on the dual variables. Again, such constraint can be written as a Fredholm equation of the first kind. In 1+1 dimensions, such equation can be solved, but in 3+1 dimensions, the dual problem turns out to be weakly infeasible. The regularized case has no such problem.
As described in the next section, instead of the regularization in (3.6), it turns out that there is a similar form of regularization that fits perfectly well with convex maximization.

Regularized convex problem and its dual
In this section we consider another motivation to regularize the primal problem. If we dualize a general conic optimization problem, the dual problem might not be feasible [35], i.e. the dual constraints may have no solutions. In this case, regularizing the primal problem solves both issues, makes the dual problem feasible and suppresses unwanted variations of ρ as discussed previously. Let us start with a generic discussion of regularization. Consider real variables v j and the primal problem: where C is a cone. Note that by using extra variables η (s) the unitarity constraint S (s) ≤ 1 can be written as a cone S (s) 2 ≤ η (s) 2 with a linear constraint η (s) = 1. Writing the Lagrangian, we find that where y a are Lagrange multipliers and s ∈ C * where is the cone dual 5 to C. Thus, if we impose on y a and s j the constraints However, in conic maximization, there is no guarantee that the dual problem is feasible, namely that the hyperplane −f − y ⋅ A parameterized by y a actually intersects C * . For some matrices A the hyperplane always intersects the cone, whereas for others it depends on the point f j , namely on the primal functional. For a simple geometric picture one can imagine the future light-cone and lines in various directions and positions that may or may not intersect the cone. Once again, the solution to this issue is to regularize the primal problem. Now we introduce two ways of regularizing that we call the M − and γ−regularizations. In the M −regularization, we define the primal problem as for some norm ⋅ and a regulator M reg . We should introduce also the dual norm ⋅ * defined as With this definition, and using v ≤ M reg from (3.12), we can write 14) The dual constraints can now be solved as namely they just determine µ j . Replacing this result back in (3.14) we get An alternative way of regularization, i.e., the γ-regularization is to con- If we now impose the conditions it follows that subject to the constraints (3.18). The M -and γ-regularizations are actually dual of each other. In the following we choose the M -regularization but the γ-regularization is equally good.

The single pion case and its dual
Now that we have discussed the idea of regularization and the primal problem we are going to construct the dual problem following the formal procedure known from standard convex maximization. In the case of a single pion, the regularized primal problem is the one in (2.27) with the extra condition α ≤ M reg for some norm ⋅ . Thus, we want to solve max αn F P = a n α n where is even. In practice we increase the regulator M reg until a plateau is reached in max αn F P as a function of M reg . The norm ⋅ is largely arbitrary and can be chosen in different ways within reason as discussed below and in section 5.
Let us now write the Lagrangian L = n a n α n + evenˆ∞ 4 where the auxiliary functions k (s) and u (s) satisfy k (s) ≤ u (s). It is then easy to see that L ≥ n a n α n = F P (3.25) since where we have used the unitarity of S (s). Replacing S (s) in terms of the partial waves we get where we defined the scalar product Given the constraint α ≤ M reg , and from the definition of dual norm (3.13), which, together with the condition leads to where we chose the lowest allowed value of u (s), i.e. u (s) = k (s) . Thus the dual problem is ds Im k (s)h ,n (s) .

(3.33)
Although replacing u (s) = k (s) makes the expression more compact, in order to use standard optimization software for the numerical part it is better to keep that part of the objective linear as in (3.27) and require (3.31).
The same can be done with Θ * by introducing a variable ν and requiring Θ * ≤ ν, i.e. the equivalent linear objective used for the numerical part reads with constraints k (s) ≤ u (s), and Θ * ≤ ν. Also, as in the primal problem (3.23), to choose a regulator M reg , we increase its value until we reach a plateau. The Θ n associated with f 0 and σ(x) can be set to zero since that part of the problem is feasible: which can be seen as a choice of norm for those components. As discussed later in (4.54), this cannot be done for the double dispersion relation part For those components a natural dual norm is given by dy Imh (x, y; s)Θ(x, y) , (3.37) that suppresses the Θ(x, y) associated with variations that affect the most the imaginary part of the partial waves. For practical purposes we can choose other norms since the results should be independent of the specific regularization. Since this is a bounded convex problem and, after regularizing, the dual is feasible, the duality gap closes [35] 6 , namely In that case we have −Re(k (s)S (s)) + k (s) = 0. Thus Since the S (s) are the partial waves of the primal problem, they automatically satisfy crossing and analyticity when the minimum is found exactly. On the other hand, as long as the relatively mild constraints (3.35a), (3.35b) are satisfied, the k (s) can be chosen arbitrarily and therefore, specially in the numerical procedure, we can truncate the space of k (s) to a finite number of partial waves and discretized energy s j . In this case the minimum has a gap with the maximum of the primal problem and therefore analyticity and crossing are satisfied only approximately. Notice that taking k (s) = 0 implies that the corresponding S (s) is undefined. Finally, if the duality gap closes we have ⟨α, Θ⟩ = M reg Θ * which allows to compute the primal variables α n .

The O(N ) case and its dual
The O(N ) case is a simple generalization of the previous case. We obtain for the primal problem max αn F P = n I F I (s 0 , t 0 , u 0 ) = a n α n , Again we can set to zero the Θ associated with f 0 and σ α (x): and choose a norm for the double dispersion relation part associated with ρ α (x, y): ds Im k I (s)h I ,α (x, y; s) (3.43) In the previous expressions the sums are over I = 0, 1, 2 and over even for I = 0, 2 and odd for I = 1. The partial waves can be computed as in (3.39): Once again they lead to partial waves that saturate unitarity unless k I (s) vanishes for some values of s in which case the partial waves are undetermined for those s values.

Generalized dispersion relations and the dual amplitudes
In 1+1 d, an important part of the dual problem is that it can be derived using generalized dispersion relations [33]. The main idea is that the functional is the value of the analytic amplitude in the unphysical kinematic region which, by using dispersion relations, can be written in terms of the boundary values, namely, the values in the physical region that are known to obey the unitarity constraints. In 3+1 d, when using double dispersion relations, we are not aware of a similar construction. Therefore, in this section, we derive a generalized double dispersion relation (4.22) which can be used to relate the value of the amplitude at an unphysical point to the values on the physical region, plus an extra term which can be bounded and can be made as small as required by appropriately choosing the function in the dispersion relation.

Dual problem in 1+1 d from dispersion relations
Suppose we have an S-matrix S(s) given by an analytic function with cuts for s > 4 and s < 0. We propose the dispersion relation satisfying crossing: S(s) = S(4 − s). Consider maximizing a functional F P = ReS(s 0 ) with 0 < s 0 < 4 under the unitarity constraint S(s + ) ≤ 1 for all s > 4. We can write where C is a small contour encircling s 0 . By deforming the contour we obtain the dispersion relation (4.1). The observation is that we can replace 1 z−s 0 with any analytic function K(z), as long as it has a simple pole at s 0 with residue one, and no other singularities except possible cuts on the real axis at s < 0 and s > 4. Then we can write For simplicity we assume (anti)-crossing K(s) = −K(4 − s) which requires adding a pole at 4 − s 0 , real analyticity K * (s) = K(s * ), and that K(s) falls sufficiently fast at infinity so that we can deform the contours C to wrap the cuts. After some algebra we find The best bound is obtained from minimizing F D over all possible functions K(z).

Generalized dispersion relations in 1+1 d
The same result can be obtained by using a similar but slightly different formulation of the problem that generalizes directly to 3+1 dimensions. We start by considering analytic functions defined through a dispersion relation 7 Regarding crossing symmetry, notice that g(4 − x) = ±g(x) implies G(4 − z) = ∓G(z). For any such function we can use contour integration assuming that we can drop the contribution from the arc at infinity to derivê Take now two such functions: H(z) given by and G(z) as written in (4.5). We havê where the last identity follows from the fact that G(z)H(z) is analytic in the upper half-plane and also in the lower half-plane as was the case in (4.7). Using and (4.6), we then find Now consider a crossing symmetric analytic function S(z) whose discontinuity is given by where θ is the Heaviside step function, and for g(z) we take such that g(4 − x) = g * (x) or G * (4 − z) = −G(z * ). In this case (4.11) reduces to F P =ReS(s 0 ) In 1+1 d we have two options. We can regularize the problem by requiring ρ < M reg and then obtain the bound where Θ(x) = Re G(x − ). In this case the dual functional contains k(x) which is not the boundary value of an analytic function, but a spectral density as in (4.13). Another possibility is to impose Re G(x − ) = 0 for x > 4 instead. In this case, using (4.6) and (4.13), one can compute the imaginary part of k(x) in terms of its real part. The result is Since κ(ξ) ∈ R and κ(4 − ξ) = κ(ξ), the functioñ defines what we call the dual amplitude. It is anti-crossing-symmetric, real analytic, has poles with residue one at s 0 and 4 − s 0 and cuts on (−∞, 0) ∪ (4, ∞). Its boundary value above the cut (4, ∞) is k(x) from where we can compute the S-matrix. Therefore we go back to the previous case. There is however an important difference between the two cases: in the first case, k(x) is not the boundary value of an analytic function but its jump across the cut. Therefore, if it vanishes on a segment, it does not necessarily vanish everywhere. This is important since vanishing of k(x) allows unitarity of S(x) in the same region not to be saturated. Therefore imposing a restriction on ρ(x) (such as ρ ≤ M reg ) allows the duality gap to close without implying unitarity saturation at all energies.

Generalized dispersion relations in 3+1 d and the dual amplitude
In the case of 3+1 dimensions we have amplitudes that depend analytically on two independent Mandelstam variables (s, t). We assume they obey double dispersion relations with support on the real axis of s and t. We are going to ignore subtractions, since the single dispersion relations can be treated as the 1+1 d case. Let us consider analytic functions of the form and define the double jump for x, y ∈ R: For any function such as G(w 1 , w 2 ), analytic in the upper and lower halfplanes (excluding the real axis) and assuming we can drop the arc at infinity, we getˆ+ where we used (4.21a), (4.21b) to drop all terms where the integral over either x or y vanished. We call (4.22) a generalized dispersion relation since, as we will see below, it can be used to relate the values of the analytic amplitude at an unphysical point to the values in the physical region. Using the identity (4.22), we can now extract from the analytic amplitude its value on the physical region s > 4, 4 − s < t < 0, by introducing an analytic function, i.e., the 3+1 d dual amplitude It is important to note that we do not takek(x, y) to be real. Now, take an amplitude H st (s, t) to be given by so that ∆ 12 H st = −4ρ st , and assume that ρ st (x, y) vanishes 9 in the physical region x > 4, 4 − x < y < 0. Then (4.22) gives (4.26) In the physical region s > 4, 4 − s < t < 0, the values H st (s + , t) can be written in terms of its partial waves as in (2.9): (4.28) where we have defined the dual partial waves k (s) as and as usual, t in the integrand is a function of µ, s as in (2.5a). We can invert this last relation as When the minimum of the dual problem is achieved, the dual partial waves are directly related to the physical partial waves as seen in (3.39). Taking the real part in (4.28) we find For comparison with (3.36) it is convenient to define namely the coefficient of the partial waves: With these coefficients we can write Replacing in (4.23) and after some algebra we obtain from (4.32) which is the same as formula (3.36). Suppose now that, we have an amplitude defined by and the dual amplitude wherek(s, t) is the same function as in (4.23) and u in the integrand should be understood as u = 4 − s − t. The same procedure leads to where the partial waves are defined in the usual manner (2.4) and (2.6). Thus, for an amplitude of the form Now we are ready to consider an amplitude of the generic Mandelstam form eq. (4.43) can be written in the form of (4.41) with We then have where K tu is given by At this point we might want to set Θ a (x, y) = 0 for a = 1, 2, 3 and x > 4, y > 4 so that we remove the ρ a 's from eq. (4.51). This is however not possible. By interchanging the order of integration in (4.23) one can write x − s + (4.53) The left hand side defines an analytic function of y with a cut on the (−∞, 0] line. If such function equals the right hand side for y > 4, then, by analytic continuation they are equal everywhere. However the left hand side is finite when y = t 0 and therefore has no pole at y = t 0 . Thus the equality is impossible and the dual problem where Θ 1 = 0 is infeasible. To elaborate this further, let us map the plane with a cut on the (−∞, 0] line to the unit disk using The region y ≥ 4 maps to the segment [−1, 0] and t 0 to a point w 0 on the positive real axis inside the disk. The function 1 y−t 0 on the right-hand side of (4.54) maps to the continuous function 1 y(w)−t 0 for w ∈ [−1, 0] and therefore can be approximated arbitrarily close by a polynomial of high enough degree. Since a polynomial is an analytic function, when mapped back to y, we can represent it as a Cauchy integral and thus write it as in the left-hand side of (4.54). That means that we can make both sides of (4.54) arbitrarily close. However, once again, if we want to make them exactly equal we need a polynomial of infinite degree, namely a series expansion. We know that such series expansion has to diverge at w 0 and does not define an analytic function inside the disk. This means that the dual problem is weakly infeasible. It is also true that for some particular primal functionals, namely some particular right-hand side of equation (4.54), the problem can be solved and in that case the dual is feasible. One simple example is a functional written in terms of a few partial waves as in (2.23). In that case, if we choosek I (s) = −ic I (s) in (3.43) then Θ α (x, y) = 0. Furthermore, in that case, the functional F D is finite since only a finite number of the c I (s) are non-vanishing. This is an important distinction, since the functional (2.10) can be written as in (2.23) but with an infinite number of coefficients as explained in appendix A. On the other hand, given that regularizing the primal problem is also important from the perspective of section 3, we are going to always consider regularized primal problems.

The dual amplitude for single pion
For the single pion all ρ a (x, y) are the same and symmetric. Therefore we get and Θ 1,2,3 (x, y) are those in (4.47a)-(4.47c). We also assumed ρ ≤ M reg and S (s) ≤ 1. The dual amplitude is where we used the kernel (2.3b). Replacing Θ a (x, y) in (4.57) one can see that Θ(x, y) is the same as the one in (3.36) and therefore only contains k (s) with even. The same is true for the first term since in (4.46), for this case, h (s) = 0 if is odd. Thus, we rederive the dual problem using the generalized dispersion relations. Although the result is the same as with the standard conic optimization method in section 3.3, the main result of this section is the existence of the dual amplitude (4.58), an analytic function of two variables that also captures all information on the physical partial waves. Namely, from its double jump k(s, t) one can compute the partial waves f (s) through eqs.(4.30) and (3.39). Thus, it contains the same physical information as F (s, t, u). A caveat is that the partial waves extracted this way from (3.39) always saturate unitarity, unless the dual partial waves k (s) vanish for some range of s. In such case the partial waves cannot be obtained from the dual amplitude and the corresponding range of energy allows the possibility of unitarity unsaturation. We leave more explorations in this direction for future work.

The dual amplitudes for the O(N ) model
For the O(N ) model the primal functional is Each F I can be treated independently, therefore we are led to the dual amplitudes (4.60) and the same for K su,I (x, y) and K tu,I (x, y). Then In this problem there are only two independent densities: ρ 1 (x, y), and ρ 2 (x, y) = ρ 2 (y, x). This allows us to write with the dual amplitudes (α = 1, 2) with an implicit sum over I and the kernels K I α are those in (2.20a)-(2.20c). Notice that the dual amplitudes depend on the point (s 0 , t 0 , u 0 ) at which we evaluate the usual amplitude. Each dual amplitude is associated with one primal variable, namely with one spectral density and has a double jump on the physical regions depicted in green in fig.4. After minimization, the double jump determines the physical partial waves from (4.29) and (3.44).

Numerical implementation
For the primal problem we have to compute the primal functional and the partial waves to impose unitarity. For that purpose, we map the region 4 ≤ s ≤ ∞ to ξ ∈ [0, π] using s = 4 cos 2 ξ 2 . (5.1) Through a redefinition of the constant f 0 and spectral density σ(x), we can rewrite the amplitude as To obtain this result we used and renamed σ 4 cos −2 ξ 2 → σ(ξ). The first term is a constant that can be absorbed into f 0 . Doing the same for ρ(x, y) we get extra terms that can be absorbed into σ(x) and f 0 . To compute the partial waves it is useful to define where, as usual, t and u are functions of µ, s as in (2.5a) and (2.5b). Now that all functions are defined in a finite interval [0, π] we can proceed to discretize the problem by evaluating the functions on a one-dimensional grid or by choosing a basis of functions. Choosing a basis leads to the method employed in [5] and we will discuss it later. First we consider the case of a discrete set of interpolation points since it gives a new alternative method that has some practical advantages.

Interpolation points
The most straight-forward approach is to discretize ξ by choosing equally spaced points: The primal functional becomes and the partial waves are where K(ξ j 1 , s + j ) should be understood as The matrixK j 1 j 2 implements the principal part integral To implement this integral we extend g(ξ 1 ) to the range −π ≤ ξ 1 ≤ π by assuming that it is anti-symmetric g(−ξ 1 ) = −g(ξ 1 ) so that the integrand is symmetric. For the spectral densities σ(ξ 1 ), ρ(ξ 1 , ξ 2 ) this agrees with them being the imaginary part of the amplitude that changes sign across the cut. Now we compute namely the discretized kernel that relates the real and imaginary part of an analytic function at the boundary of the unit disk. The dual problem can be constructed from this by discretizing the dual partial waves, namely using variables and computing the dual functional as For each primal variable there is a Θ component given by (3.33): where the coefficients are those in (5.8a)-(5.10c). Now we require which can be thought as taking Θ and Θ j 1 to have infinite norm unless they vanish. For the double dispersion we define In this way we have the dual of the numerical primal and therefore the minimum is equal to the maximum of the primal. It is a very useful way to evaluate the dual since it requires fewer partial waves. To obtain the actual dual that is above the maximum we have to restrict the space of dual partial waves. The way to accomplish this is to define a set of coefficients and expand the dual partial waves as k (ξ) = n k n e inξ (5.20) The coefficients k ln are taken real. Thus k has a symmetric real part and antisymmetric imaginary part under ξ → −ξ as corresponds to the dual of a real analytic function as in (3.39). In our approach we evaluate the integrals numerically being careful to keep enough interpolation points to correctly describe the functions. As a practical rule we have that, if we have N max coefficients then we should have at the very least M = 4N max interpolation points. So, the final form of the dual problem is and we further impose the constraints 0 = a + n Im[k n θ n ] (5.23a) 0 = a j 1 + n Im[k n θ n;j 1 ] (5.23b) The coefficients θ (a) are given by which can be efficiently evaluated using the Fast Fourier Transform.

Discrete basis of functions
It is useful both conceptually and also for numerical purposes to consider the same problem in a discrete basis of functions by taking The simplest functions to use as a basis are (n ≥ 1) To this set it is convenient to add an extra function φ 0 (x) = 1 √ x−4 that allows the inclusion of a divergence at x = 4 (or y = 4). Here we only do that as a check since the purpose is to obtain the existence of a pole from the optimization. We then have

Let us defineΦ
where as always, in the integrand, t and u are functions of µ, s as in (2.5a), (2.5b). Now we can write the partial waves ( even) as The functional is F P = a 0 f 0 + n a n σ n + nm a nm ρ nm (5.33) For example if we take F P = F (s 0 , t 0 , u 0 ) then a 0 = 1 , (5.34) a n = Φ n (s 0 ) + Φ n (t 0 ) + Φ n (u 0 ) , (5.35) The dual problem is which come from setting to zero the coefficients Θ associated with f 0 and σ n . For the double dispersion relation we have where we used the usual square norm but others lead to the same results. The constraints (5.38b) are a small set of contraints so the space of functions k (s) is largely unconstrained. The full set of dual constraints should include Θ nm = 0 that, as mentioned cannot be imposed since it has no solutions. However, if M reg ≫ 1, minimizing the functional generically results in functions k (s) such that Θ nm ≃ 0. In the O(N ) isospin case, we proceed in the same way. The (t ↔ u) symmetric amplitude A(s, t, u) in (2.12) is written as using the same basis functions as in (5.29). The primal variables are the σ (1,2) n , ρ (1,2) nm with the only condition that ρ where (we omit the argument s in all functions): and

44b)
where for I = 0, 2 we only have even and for I = 1 only odd. Finally we choose a functional The dual problem is now min k I (s) notice that, in constructing the dual problem we can directly use the formulas (3.41) since they do not depend on how we choose the variables α n to define the primal problem. Notice that once again in the dual problem the functions k I (s) are independent. Crossing only appears in the dual problem through the coefficients A ,nm (s), B ,nm (s). Finally, the dual problem thus defined is the dual of the numerical primal problem and then its minimum should agree with the maximum of the primal. To obtain the actual convex dual we have to restrict the space of variables k (s) by introducing once again a set of real coefficients to parameterize the dual variables: The number of coefficients k p should be much smaller than the number of coefficients in representing the amplitude (5.26). The reason is that now, those coefficients are taken as a way to evaluate the integrals in the dual problem and cannot do that accurately if the functions k (s) fluctuate more rapidly than the functions in the basis (5.26).

A few numerical results
The main idea for the numerical approach is that, in the primal we parameterize the amplitude using the interpolation points as described in 5.1, or as linear combinations of a set of functions as in 5.2. As we increase the size of this set, the numerical maximum increases approaching the true maximum.
In the dual, we write the dual partial waves as linear combinations of a set of functions and compute the minimum. As we increase the size of the set, the minimum decreases approaching the maximum of the primal problem from above. To check the method, we choose a known but challenging problem: maximizing the value of the amplitude for single pions at a point in the Mandelstam triangle, for example at the symmetric point s 0 = t 0 = u 0 = 4 3. This is numerically challenging because the maximum is attained by a function that has a pole at threshold (s = 4) as studied in [5]. There, it was noticed that the convergence to the maximum in terms of the size of the basis is slow unless one explicitly includes a pole at threshold. For that reason the dual problem we study here is important since it provides an upper bound that is above the maximum. Note that this is all done without assuming the existence of a pole at threshold. In fact one can use the same no-pole coefficients for the primal and dual. Of course, although the coefficients are the same, they enter differently in the primal and the dual leading to lower and upper bounds. For the actual optimization we use standard software [36,37]. We considered first the amplitude at the symmetric point F = F ( 4 3 , 4 3 , 4 3 ). In fig.3 we see how, increasing the regulator M reg the maximum increases until it reaches a plateau that cover several orders of magnitude of the regulator (notice the logarithmic scale in the horizontal axis). That plateau is taken as the maximum. For larger values of the regulator the maximum can increase further or fluctuate randomly due to numerical errors, that part should be ignored. The dual problem reaches the plateau always above the actual maximum and the primal always below. As shown in fig. 3, in the primal problem we used M = 20, 30, 40, 60, 100, 160 interpolation points for the single and double dispersion relations where unitarity constraints are imposed. For the dual we used 512 interpolation points and 16 ≤ N coeff ≤ 128 and only dual partial waves up to max = 4. It is interesting to extrapolate the result to a very large number of interpolation points or coefficients. This is done in fig.4 where one can see that for an infinite number of coefficients (vertical axis) the maximum converges to a number very close to the value obtained assuming the existence of a pole. The dual seems to be doing better without needing extrapolation.
At the plateau of figure 1, we can compute the physical partial waves h (s) from the dual partial waves k (s) using eq. (3.39). As mentioned before, the dual formulation does not assume the threshold pole as the maximum is attained, and can therefore detect such a pole as a result of the optimization. The corresponding phase shift δ (s) can then be computed using the definition (2.7). In figs. 5 and 6, we show the S-wave h 0 (s) and its phase shift δ 0 (s) obtained from the dual partial waves and compare it with the primal result when a pole is included in the setup, and find perfect agreement.
To show the power of the method further, we consider F = max F (s 0 , s 0 , 4− 2s 0 ), with 0 < s 0 < 4 whose result have been displayed in fig.1 in section 1. The primal data points were obtained by using basis functions as described in section 5.2. The lower curve is done with 8 × 8 coefficients for the double dispersion, 100 coefficients for the single dispersion and checking unitarity up to max = 10 at 100 values of s. The upper red curve is 10 × 10 coefficients for double dispersion and 300 for single, checking unitarity at 300 points for max = 20. To check we also run 100 interpolation points (for single and double dispersion relation) which agreed with the lower curve. This is be- as the regulator M reg is increased for the primal (red) and dual (green) problems. The horizontal axis gives the regulator in a logarithmic scale and up to a horizontal shift to plot the curves together. As explained in the main text, the primal improves by increasing the number of interpolation points, as indicated with the red upward arrow, and the dual improves by increasing the number of coefficients, as indicated with the green downward arrow. We do not asume the existence of a pole at threshold which would give the horizontal blue line. The values are multiplied by a factor π 4 to match the normalization in the literature, e.g. [5]. changes for the primal (red) and the dual (green). The vertical axis correspond to infinite number of interpolation points for the primal or coefficients for the dual. We do not assume the existence of a pole at threshold which would give the horizontal line. The values are multiplied by a factor π 4 to match the normalization in the literature, e.g. see [5].  cause the pole at threshold is in the s-wave and can be reproduced from the single dispersion relation with enough coefficients. In general running 100 interpolation points should be better. The dual (green points) was run using 512 interpolation points and using partial waves up to max = 12. Finally, the blue line is the result of including a pole 11 . We can see that the primal and the dual bracket that line. Even if we did not have the blue line, we are still assured that the program is converging from above and below. The black diamonds are the pioneering results of Lopez and Mennessier in [34] where upper bounds were found by a somewhat different method. When s 0 = t 0 → 0 and s 0 = t 0 → 4, numerically we observe a divergence F max ∼ 1 √ s 0 and F max ∼ 1 √ 4−s 0 respectively and we also find a minimum at the symmetric point s 0 = t 0 = u 0 = 4 3 where π 4 F ( 4 3 , 4 3 , 4 3 ) ≃ 2.6613. This value at the symmetric point agrees with the result in [5], we just observe that it is the minimum of the blue curve.

Conclusions
The S-matrix bootstrap provides a method to investigate the space of Smatrices allowed by the general constraints of analyticity, unitarity, crossing, and global symmetries. When restricted to a subsector such as 2 → 2 scattering the space can be efficiently mapped by using conic maximization numerical tools. In particular the space is convex. At the boundary of the space sometimes distinguished points such as vertices can be found that define interesting S-matrices which correspond to physically important theories. Since the numerics requires discretizing the problem, an approximation is made. One can control the approximation by running a primal problem that reduces the number of variables on which the S-matrix depends. In that case we map an interior region of the actual allowed space that grows as we increase the number of variables. Alternatively, one can do an approximation in the dual convex problem, in which one restricts the space of dual variables, namely of Lagrange multipliers that enforce the constraints. In that case one obtains an exterior region that contains the allowed region of S-matrices. By run-ning both problems we bracket the region between an interior and an exterior one [13]. This is particularly useful in the case of a challenging numerical problem such as the scattering of scalar particles in 3 + 1 dimensions. The problem is that the maximum is attained by an amplitude that has a pole at threshold [5]. One would like to reproduce those results without assuming the existence of a pole. In this paper we constructed explicitly the dual problem for scattering of scalars in 3 + 1 dimensions, both for a single scalar and the more general O(N ) model. Numerically, when applied to the single scalar case we obtain a good bracketing of the maximum without assuming the existence of a pole. Another interesting property is that the dual problem allows for the possibility of unitarity unsaturation, however, in practice, it is not clear when that would happen. In summary, this work provides a new valuable tool that can be efficiently used to map out the space of S-matrices in higher dimensional theories.

Acknowledgements
We are very grateful to Lucia Cordova, Harish Murali, Joao Penedones and Pedro Vieira for discussions, as well as comments and suggestions on the draft. We would also like to thank the referee from JHEP for many interesting questions and suggestions. In addition, M.K. is very grateful to the DOE that supported in part this work through grants DE-SC0007884, DE-SC0019202 and the QuantiSED Fermilab consortium, as well as to the Keck Foundation that also provided partial support for this work.

A Useful formulas
In (2.28c) and (2.28c), using results from [41], the integrals over µ were done analytically in terms of the Legendre functions Q usinĝ where, t and u are functions of µ and s given by (2.5a) and (2.5b). Another useful formula is [41] ∞ =0 (2 + 1)P (x 1 )Q (x 2 ) = 1 It is interesting to notice that, using this last identity one can see that solves the constraint Θ(x, y) = 0 in (3.36). However the argument of the Legendre polynomial is larger than one and therefore the functions k (s) are real and negative and grow in absolute value with . Therefore the functional F D in (3.33) is divergent. Also the sum in (4.30) is divergent and (A.4) does not lead to a well defined dual amplitude in agreement with the observation after (4.54).