Adding flavour to the S-matrix bootstrap

We explore the S-matrices of gapped, unitary, Lorentz invariant quantum field theories with a global O($N$) symmetry in 1+1 dimensions. We extremize various cubic and quartic couplings in the two-to-two scattering amplitudes of vector particles. Saturating these bounds, we encounter known integrable models with O($N$) symmetry such as the O($N$) Gross-Neveu and non-linear sigma models and the scattering of kinks in the sine-Gordon model. We also considered more general mass spectra for which we move away from the integrable realm. In this regime we find (numerically, through a large N analysis and sometimes even analytically) that the S-matrices saturating the various coupling bounds have an extremely rich structure exhibiting infinite resonances and virtual states in the various kinematical sheets. They are rather exotic in that they admit no particle production yet they do not obey Yang-Baxter equations. We discuss their physical (ir)relevance and speculate, based on some preliminary numerics, that they might be close to more realistic realistic theories with particle production.


Introduction
In this paper we continue the exploration of the space of gapped quantum field theories following [1][2][3] by focusing on two dimensional theories with a global symmetry under which particles transform. We will consider the two-to-two scattering processes of particles transforming in the vector representation of O(N ). As we will see, such processes turn out to be way richer than their analogous counterparts without global symmetry.
As described in much more detail in the next section, kinematically, such two-to-two scattering amplitudes live in the Mandelstam physical sheet or equivalently in the physical strip 0 < Im(θ) < π in terms of the hyperbolic rapidity θ. Direct s-channel processes take place at the lower boundary of this strip, for real rapidity. There, the various amplitudes are bounded by unitarity as |S rep (θ)| 2 ≤ 1 where rep are the various possible representations formed by the two incoming particles. At the upper boundary of the physical sheet we have the crossed t-channel processes. Here is where global symmetry manifests itself rather strikingly as a tension between unitarity and crossing symmetry. Indeed, under crossing transformations, the various individual components are trivially swapped but when we translate that back to the representation basis -for which s-channel unitarity was so conveniently simple -we obtain a very non-trivial mixing of the various representations. In this upper boundary we then have | rep' d rep,rep' S rep' (iπ − θ)| 2 ≤ 1 where d is a constant N -dependent matrix of purely group theoretical origin. Exploring the space of S-matrices with global symmetry thus amounts to studying the space of functions leaving in this strip and bounded at its boundary in such coupled way. It is this fascinating problem which we will consider here. 1 By looking for theories which maximize particular couplings, we will rediscover in this way the two most famous integrable models with O(N ) symmetry -the O(N ) non-linear sigma model and the O(N ) Gross-Neveu model -whose S-matrices were found by Zamolodchikov and Zamolodchikov in their seminal 1979 paper [4]. For N = 2, we rediscovered the S-matrix for the kinks of the Sine-Gordon model. Furthermore, we make contact with a much less known integrable solution [6]. We also found some other analytic solutions whose physical (ir)relevance is discussed below. Some other times, these maximization problems require numerics. Both in the analytic and numerical examples, we will often unveil a very rich structure of infinitely many resonances in the various Riemann sheets for the putative Smatrices with the largest possible physical couplings.
In section 2 we introduce the O(N ) S-matrix setup and we review Zamolodchikovs' Smatrices mentioned above. We also explain the general numerical setup in this section. We then discuss in section 3 how the integrable S-matrices can be rediscovered from these numerics when we impose very specific spectra of bound states (or absence thereof). As we move away from these special points, the S-matrices develop very interesting new features. We present in the same section some numerical results for some non-integrable cases. In section 4 we describe an analytic large N analysis which yields some intuition for the rich mathematical structures found numerically. We then present the analytic solution of various cases along with some general analytic properties for the S-matrices. Finally in section 5 we give some concluding remarks.
Note: While we were finishing the preparation of this article, we became aware of the work [7] which overlaps with section 3 here.

Setup, key examples and numerics 2.1 S-matrices
We consider relativistic particles of mass m transforming in the O(N ) vector representation. Our particles have an O(N ) index i = 1, . . . , N and an hyperbolic rapidity θ parametrizing its energy and momentum as the usual: E ± P = m exp(±θ). Lorentz boosts act as translations in hyperbolic rapidity hence Lorentz invariant quantities depend on difference of rapidities only. An example which will be the focus of this paper is the two-to-two particle S-matrix S kl ij (θ) defined in the standard way 2 |θ 1 , i 1 ; θ 2 , i 2 out = S j 1 j 2 i 1 i 2 (θ = θ 1 − θ 2 )|θ 1 , j 1 ; θ 2 , j 2 in .
(1) Figure 1: Mapping between s and θ variables. The cuts corresponding to the two particle thresholds with branch points at s = 0 and s = 4m 2 get opened in the θ plane. These cuts are square roots, so that going twice around one of the branch points leaves you back where you began. This is exemplified in the red path for the branch point s = 4m 2 (or the regular point θ = 0). The thick lines in black represent possible inelastic thresholds. Unitarity is imposed for physical values of the center of mass energy s ≥ 4m 2 (θ > 0) represented by the green dashed line. The physical sheet(strip) in s(θ) is highlighted in grey. The bound state poles are located in the window s ∈ 0, 4m 2 (θ ∈ [0, iπ]) (blue line). Finally, the points in yellow are related by crossing which acts as a reflection around s = 2m 2 (θ = iπ/2). Notice that in general there are infinitely many sheets.
In terms of Mandelstam invariants, we have s/m 2 = 4 − t/m 2 = 2 + 2 cosh(θ) and u = 0 in two dimensions. The very useful map between s and θ opens up the two particle cut and is depicted in figure 1.
The S-matrix can be coded into three functions corresponding to the three possible O(N ) invariant tensors with four indices K kl ij = δ ij δ kl = , I kl ij = δ l i δ k j = , P kl ij = δ k i δ l j = or, equivalently, to the three irreducible representations arising in the product of two fundamental incoming representations. Explicitly, we have S(θ) = σ 1 (θ)K + σ 2 (θ)I + σ 3 (θ)P = S sing (θ)P sing + S anti (θ)P anti + S sym (θ)P sym .
In the component description crossing symmetry is trivial to impose since under crossing transformations we swap incoming and outgoing particles so that the O(N ) invariant tensors are simply re-shuffled. In the irreducible representation description unitarity is straightforward since by preparing a two-particle state in a definite O(N ) representation we diagonalize (1). Crossing symmetry and unitarity in both descriptions are summarized in tables 1 and 2. Alternatively, one could also work in a basis where crossing is diagonalized as in [8].
We assume here that the particles being scattered are the lightest particles in the theory. There could also be bound states showing up in the S-matrix of two fundamental particles. They would transform in the singlet, anti-symmetric or symmetric traceless representation of O(N ) and show up as a pole in the corresponding channel. For instance, if there is a single bound state transforming in the anti-symmetric representation with mass m BS = 2 cos(λ/2) there is a pole in the corresponding channel as 3 which then, according to crossing, would lead to poles in the t-channel for all components as can be read from table 2, It is because of the interplay between these various channels, also illustrated in the tension between simplifying unitarity and crossing at the same time as summarized in tables 1 and 2, that the O(N ) bootstrap problem is so much richer than the single component case.

Integrable O(N ) S-matrices
In this section we review the integrable S-matrices with O(N ) symmetry obtained by Zamolodchikov and Zamolodchikov in [4]. In this seminal work, starting from the decomposition (2), 3 J = 2 sin(λ) is a trivial Jacobian. Notice that here we define the coupling g 2 rep to be the residue of the S-matrix S rep (s) and not of T rep (s) so we would need one further simple Jacobian if we were to compare the results in this paper with the conventions in [2].
where M is a non-negative integer, α j are free parameters, and where λ GN ≡ 2π N −2 . Here NLSM stands for Non-linear Sigma Model and GN stands for Gross-Neveu for reasons that will become clear shortly. Let us review and highlight a few of the many very interesting features of (7) and (8) 1. The first thing to note is that this solution has a very rich pattern of poles and zeros in the infinitely many copies of the Mandelstam s plane (equivalently, in the infinitely many strips in the θ plane as illustrated in figure 1) as made explicit by the various gamma functions in the functions F a (θ). At the same time note that the product F π−λ GN (θ)F 2π (θ) does not contain any poles (or zeroes) inside the physical strip. The vector in (8) multiplying these F 's also does not have any poles inside the physical strip (only a zero inside the strip at θ = iλ GN and a zero at the boundary of the strip at θ = iπ for some components.) . Integrability implies that unitarity is saturated, i.e. S rep (θ)S rep (−θ) = 1, so that in a given representation a pole at θ has an associated zero at −θ.
Crossing mixes the different representations relating points at θ and iπ − θ (see discussion around (6)).
2. Therefore, potential bound state poles must come from the Castillejo-Dalitz-Dyson (CDD) pre-factor which is the product in (7) for α j ∈ [0, π]. The simplest solution commonly referred to as the 'minimal' solution corresponds to setting M = 0, i.e. to introducing no CDD factors and hence no bound states. 5 Beautifully, there is a theory whose S-matrix is precisely given by this choice: it is the O(N ) non-linear sigma model [4]. The analytic structure of this solution is presented in figure 2.
3. As made explicit in (7), when there are bound states at generic positions θ = iα j (with α j ∈ [0, π]) they are common to all representations (with an important exception to be discussed in the next point). This means that integrable theories produce very degenerate spectra where bound states in different representations come along at once with the same mass. For example, there is no integrable theory where the particles form a single bound state in the singlet channel. (We can view this as a nice feature: the bootstrap of O(N ) symmetric theories with particles in the vector and anti-symmetric representations alone, will necessarily land us outside the integrable world. We will investigate these cases further in section 3.2.)

4.
A simple exception is when we consider a single CDD factor, i.e. M = 1, with α 1 = λ GN . In this case the CDD factor introduces two poles in the physical strip: one at θ = iλ GN and another one at θ = iπ − iλ GN . However, in the symmetric representation the first one is cancelled by the zero explicitly seen in the vector in (8). Therefore we are left with an s-channel pole at θ = iλ GN for the singlet and anti-symmetric representations and a t-channel pole at θ = iπ − iλ GN for all representations.
The former are identified as the bound states of the theory 6 . Since both s-channel poles sit at the same position, the bound states have a common mass m sing = m anti = 2 m cos(λ GN /2). There is also a beautiful physical theory corresponding to this Smatrix: it is the O(N ) Gross-Neveu model [4]. The structure of poles and zeros of this model is depicted in figure 3.

5.
There is a subtlety with the last point. For the Gross-Neveu solution the sign of the residues at the s-channel pole θ = iλ GN is appropriate for N = 7, 8, 9, . . . but opposite of what it should be for N = 5 and the pole disappears altogether from the physical sheet as N = 4 or N = 3. 7 (a) For N = 7, 8, 9, . . . the residues have the proper sign and the bound state interpretation holds perfectly. This case is connected to the N → ∞ limit where things simplify and the theory becomes effectively free. 8 (b) For N ≤ 4 the potential bound state poles leave the physical strip. Actually, for N = 3, 4 the Gross-Neveu and NLSM solutions coincide. 9 This small puzzle is resolved by realizing that the the Gross-Neveu spectrum includes kinks 10 which are the only stable particles below N = 5. So for N = 3, 4 the S-matrix above cannot describe vector particles of the Gross-Neveu model since they do not exist at all for these values of N .
(c) For N = 5 the poles have opposite signs and must be reinterpreted. Actually, the bound states disappear from the spectrum (see equation (6.11) in [4]) so that the poles describe now a pair of kinks instead and correspond to Coleman-Thun poles.
6. Note that at large energy (i.e. s → ∞ or θ → ∞) we have (for any representation) S int rep 1 + ia rep / log(s) so that at very high energy the S-matrices become effectively free but they approach this free limit very slowly, logarithmically so. This logarithm is quite physical, it is a sign of asymptotic freedom of these O(N ) integrable theories.
7. Finally, there is another solution to the factorization (Yang-Baxter) equations which appeared in the appendix of [6] and is much less known. To our knowledge, it is has not been understood if there is a physical theory it corresponds to. We make contact with this solution in section 4. 11 This concludes the review of the N > 2 solution. For N = 2 the solution is even richer, with infinitely many gamma functions and is presented in section 3.1 below. The interpretation of the minimal solution in this case is rather different. It describes the scattering of the kinks and anti-kinks of the sine-Gordon model (which is dual to the massive Thirring model where these kinks correspond to the fundamental fermions).

Numerical setup
Following the analytic structure of the S-matrix, one can write discretized dispersion relations like the ones used in [3] to get numerical bounds on the couplings. In the present case with global symmetry O(N ), there are two main differences: there are three functions S rep (s) coupled by crossing and the high energy behaviour requires the use of subtractions.
Crossing symmetry tell us that the vector S = (S sing , S anti , S sym ) satisfies (setting m = 1 In other words, although each element is very small, the phase space is very large leading to a final net result for the singlet channel. This is a general expected feature of large N theories: because of the large phase space, the singlet component should dominate in this limit. We will re-encounter this clearly in a large N analysis in section 4.1 below. 9 It is simple to see that the CDD factor with α 1 = λ GN = 2π/(N − 2) equals one for those values of N. 10 These transform in the spinor representation of O(N ) and were studied in [9,10]. The Gross-Neveu spectrum includes as well a tower of anti-symmetric tensors whose S-matrices were obtained in [11]. 11 We became aware of the work [6] after we published version 1 of the present paper. We thank A.B. Zamolodchikov and S. Komatsu for various inspiring discussions of closely related S-matrices which eventually led us to this reference.
as can be read off from table 2 in section 2.1. Therefore the t-channel poles and the t-channel discontinuities of the vector S are related to their s-channel counterparts through this same matrix (while in [2] they were simply identical.) As for the slow decay at infinity this is solved by starting our dispersion relation derivation with the identity 12 where the integral goes over a small circle around a point s inside the physical sheet. Because we divided by s − 2 the integrand decays very fast at infinity so we can blow up the contour safely 13 and in this way end up with a contour around the S-matrix poles and the s− and t− channel multi-particle cuts. Using the matrix d in (9) to relate t-channel processes to s-channel ones we finally end up with the dispersion relation where the various couplings and bound state masses for the n BS bound states are contained in the vectors where rep n is the representation under which the n-th bound state transforms. For instance, for the non-linear sigma model we have n BS = 0 so we would need no pole vectors. For the Gross-Neveu model we have n BS = 2 and these two bound states have the same mass and transform in the singlet and anti-symmetric representation respectively so that in this case we would have two pole vectors pole GN In the dispersion relations (11) we have trivialized crossing taking into account the O(N ) symmetry of the problem and possible large energy behaviour. The remaining ingredient 12 Note that the subtraction constants S(2) satisfy the crossing condition S(2) = d · S(2) which fixes one of the constants in terms of the other two (e.g. we can eliminate S sing (2) by writing S sing (2) = 1 2 [(N + 2)S sym (2) − N S anti (2)]). 13 If the S-matrix approaches a constant at infinity this single subtraction suffices. Note that with this subtraction the integrand vanishes at 1/(x 2 log(x)) at infinity for the O(N ) symmetric integrable S-matrices of the previous section. Without the subtraction we would have S(x)/(x − s) ∼ 1/(x log(x)) which would not decay fast enough to allow us to blow up the contour.
is to impose unitarity. In order to do this numerically, we follow [2] and choose a grid of n grid points s j > 4 (j = 1, . . . , n grid ) in which we are to impose unitarity and discretize the densities ρ → ρ j = ρ(s j ) with a linear interpolation between these points so that we can explicitly perform the integrals in (11) and obtain discretized dispersion relations. 14 In this way we impose n grid unitarity constraints |S d (s j )| ≤ 1, where the superscript d denotes the approximated S-matrix using the discretized dispersion relations.
With this ansatz at hand, we can now choose a bound state spectrum and maximize one of the variables (or combinations of the variables) in the dispersion relation. We could maximize bound state couplings g 2 n as originally proposed in [2] or the value of the S-matrix at a symmetric point s = t = 2 as in the pion toy models in [3] or we could impose a zero at a given value for some component (i.e. add a resonance) and maximize its slope as in [12]. In the next section we will start with a few such maximization questions (some of which are mixed versions of the ones just mentioned) which lead to the integrable S-matrices discussed previously. Then we will move away from these integrable models and find new S-matrices whose physical (ir)relevance we shall speculate about.

Reproducing integrable models
In this section we show how we can reproduce known integrable theories with O(N ) global symmetry. In general, we expect the numerical maximization to reproduce the physical S-matrices for the appropriate mass spectrum. First, we consider a simple problem which reproduces free theory. Then we focus on the case with bound states and N > 2 which gives rise to the Gross-Neveu S-matrix. After that we discuss the maximization procedure leading to the non-linear sigma model. Finally we present the case with N = 2 and bound states which reproduces the S-matrix for the kinks of sine-Gordon (or massive Thirring model).

Free theory
Consider a theory with no bound states. Then the various S-matrix components have no poles inside the strip and are therefore bounded by their values at the two boundaries of the strip (i.e. the lower boundary at θ ∈ R and the upper one at θ ∈ iπ + R). At the lower boundary of the physical strip we have unitarity which states that in any of the three representations we have |S rep (θ)| ≤ 1 while in the upper boundary these functions can be as large as allowed by the crossing relations in table 2 which relate that upper boundary with the value of the same functions in the lower boundary. For the symmetric component these equations tell us that we can have modulus as large as 1 in the upper boundary which is the same bound as in the lower boundary. 15 Therefore the largest this component can be, anywhere inside the strip is 1 and this is attained for S symmetric = 1 which in turn leads to S rep = 1 for all representations. We see that the theory with no bound states and the largest absolute value for the symmetric component anywhere in the physical strip is the free theory. This is something we can check numerically; it obviously works and is indeed a nice albeit trivial check of our code.

Gross-Neveu
As described in section 2.2, there are two bound states appearing in the 2 → 2 S-matrix of the fermions of the Gross-Neveu model. These bound states appear in the singlet and antisymmetric representations and have a common mass m GN = 2 cos π N −2 . In principle we could maximize any of the two corresponding cubic couplings g sing , g anti , but it turns out that it is the maximization of g anti which reproduces the Gross-Neveu S-matrix.
Considering the spectrum above in the dispersion relations (11) and implementing the numeric maximization of g anti as explained in section 2.3 we find perfect agreement with the integrable solution where F a (θ) is defined in (8). We performed the numerics for various values of N . In figure 4 we show some numerical results against the solution (13) for N=7,11,15. The analytic structure of one of the representations is depicted in figure 5. There is a zero at s = 0 (in green), a bound state pole at s = m 2 GN (in dark blue) and a t-channel pole in lighter blue. The green, blue and red dashed lines are respectively the regions below the left cut, between s ∈ [0, 4m 2 ] and above the right hand cut where we impose unitarity. (b) Numerical results for the corresponding regions in the θ = α + iβ plane wit N = 7 (m GN = 2 cos(π/5)) and n grid = 75. The top panel in green depicts the zero at θ = iπ (s = 0). The middle panel in blue shows the bound state and t-channel poles. In the red bottom panel we see that unitarity is saturated (solid curve representing the absolute value of the function).

Non-linear sigma model
Another famous integrable model we would like to reproduce is the non-linear sigma model. As discussed in section 2.2, this model has no bound states and its S-matrix differs from the Gross-Neveu one by an overall CDD factor. Instead of bound states, the parameter λ GN in (8) labels the position of a zero in the symmetric representation. In the numerics, we impose this zero S sym (θ = iλ GN ) = 0 as well as the absence of bound states by setting all cubic couplings to zero. Inspired by the Gross-Neveu case, we maximize the effective quartic anti-symmetric coupling given by S anti (iπ/2). With these conditions we are able to reproduce the solution (8) as exemplified in figure 6 for N = 6.
As we decrease N , the zero in the symmetric representation moves towards the upper boundary of the physical strip, reaching it at N = 4 and leaving the physical strip for N < 4. It would be interesting to understand what is the condition one should impose in order to reproduce the non-linear sigma model for N < 5.
In the beautiful work [8] it was pointed out that the space of unitarity, crossing symmetric relativistic S-matrices with O(N ) symmetry has a very rich geometric structure with special cusps. Provided we point more or less towards one such cusp, any maximization functional will push us there. They could identify the non-linear sigma model as a cusp in this theory space, without imposing the zero at S sym (θ = iλ GN ). It would be fascinating to explore more throughly the geometry of this S-matrix space, in particular for the various other maximization problems considered in the following sections.

Sine-Gordon kinks
Another important result derived in [4] is given by the case N = 2 corresponding to the integrable S-matrix for the sine-Gordon kinks and antikinks. 16 The kink A and antikink A can be packed into a doublet A = A 1 + iA 2 (Ā = A 1 − iA 2 ). The solution has a free parameter γ 17 and reads where the prefactor is given by The S-matrix (14) has a very rich structure. As we vary the parameter γ, the bound state spectrum changes. The bound state masses are given by m n = 2m sin(nγ/16) with n < 8π/γ and they alternate between the anti-symmetric (n odd) and singlet representations (n even). 17 As explained in [4], this parameter is related to the coupling in the sine-Gordon Lagrangian L sG = 1 2 (∂φ) 2 + m 2 0 /β 2 cos(βφ) through γ = β 2 /(1 − β 2 8π ).
The bound state with n = 1 is the sine-Gordon breather whose S-matrix was rediscovered in [2]. For γ ≥ 8π all bound states disappear from the spectrum (i.e. the corresponding poles leave the physical strip).
We shall focus on two regimes which comprehend 1 < m 1 < 2: • For 4π ≤ γ < 8π there is a single bound state in the anti-symmetric representation whose mass takes values √ 2 < m anti < 2 and is given by m anti = 2 sin(γ/16). However, imposing this bound state and maximizing g anti is not enough to reproduce the integrable model. Taking inspiration from the NLSM discussed in the previous section, we shall impose one of the zeros in the solution (14). One zero which is always inside the physical strip for this range of γ is S sing [i(−π + γ/8)] = 0 which can easily be seen in the first component of (14).
• For 8π/3 < γ < 4π there is another bound state in the singlet representation with mass m sing = 2 sin(γ/8), as well as the previous bound state which for this range takes values 1 < m anti < √ 2. For the numerics, we keep maximizing the anti-symmetric coupling and impose one zero inside the physical strip, this time at S sing [i(−π + γ/2)].
In figure 7 we show the maximum coupling g anti obtained numerically for both regimes which nicely matches the analytic solution (14). The S-matrices agree perfectly as well.
In [2] it was shown that the sine-Gordon S-matrix for the lightest breathers is quite special. It has a single bound state flowing (the second lightest breather) and of all theories with a single bound state it is the one with the largest possible coupling. Here we see that the kink S-matrix of sine-Gordon also comes about from a maximization procedure as the theory with O(2) symmetry with the largest possible coupling and appropriate resonances (zeroes). These two are related. As mentioned earlier, breathers are themselves bound states of kinks so we can obtain the breather S-matrix by fusing the kink S-matrix. In that sense, the kink S-matrix is the most fundamental one and it is nice to see it come up naturally here. There are also intermediate objets such as the scattering of kinks and breathers. It would be very interesting to look for those as well.

Away from integrable points
In this section we present some of the numerical results obtained for non integrable spectra. The simplest possibility is to consider theories with a single bound state transforming in one of the three representations and maximize its associated cubic coupling g rep . Figure 8 depicts the maximum couplings as we vary the bound state mass m rep and N . In these numerics, we can also read of the corresponding S-matrices. We observe -as in [2] -that they all saturate unitarity (and thus admit no multi-particle production). That is, for s > 4m 2 they all satisfy |S rep | = 1 within our numerical precision. Also as in [2] we observe that these couplings vanish as the bound state becomes weakly coupled as expected. It should be possible (and interesting) to analyse this perturbative corner analytically. 18 In the next section we will show how these maximal coupling curves (and others) can be reproduced analytically at large N .
What is not shown in figure 8 is how interesting the corresponding S-matrices are. By plotting them in the complex θ plane we found remarkably rich structures of infinite poles and zeros in the various θ-strips. 19 An example is depicted in figure 9.
Interesting as they might be, we will not be able to compare the theories saturating the bounds in figure 8 with any known theory. The reason is that since they have no particle production, they should probably be integrable. But since they have a single bound state they cannot satisfy Yang-Baxter. So, if they exist they must be something more exotic. We will speculate further on the physics of such theories in the discussions.
To make contact with known integrable examples we can consider one bound state in the anti-symmetric representation and another one in the anti-symmetric representation with equal mass and maximize the cubic coupling to the anti-symmetric one. When the mass passes by the value m GN = 2 cos( π N −2 ) we do obtain the Gross-Neveu model as seen in the previous section. Away from this point, the bound is deformed smoothly but the S-matrices are quite different, with a very exotic analytic structure akin to that illustrated in figure 9. The maximum couplings we obtained for various N and bound state masses are shown in figure 10. The divergence at m = √ 2 comes from the collision of s and t-channel poles as in [2]. In this case they screen each other by simply setting g 2 sing = N 2 g 2 anti . Below, we will revisit some of the cases found numerically here. To do so, first we will get some further insight from a large N analysis in section 4.1 before attacking the full finite N case in section 4.2.

Large N
In this section we will consider large N . We will be able in this way to proceed analytically to a much greater extent and get some insights about the general analytic structures to be expected in the finite N optimization numerics.
The key place where N enters and where large N presents some simplifications is in 19 Note that we explore the S-matrices across all the Mandelstam cuts while our numerical ansatz is originally designed to parametrize the physical sheet alone. We can do this because unitarity is always saturated from 4m 2 until the first inelastic threshold which leads to an exact functional relation S rep (θ)S rep (−θ) = 1 which we can combine with crossing to visit any strip in θ, i.e. any sheet in s. See e.g. [14] for a similar analysis for a theory with a single particle.
At real θ all S-matrices saturate unitarity; they are phases so they are all order O(N 0 ) functions, even at large N . At the upper boundary of the physical strip, at θ ∈ iπ + R, we can use these crossing equations to see what values the various S-matrices can take. If we can bound the functions in the two boundaries of the physical strip then we can use the maximum modulus principle to bound the functions everywhere inside the strip. That is the main idea which we will explore.
The singlet component can be very large in the upper boundary of the strip. It can be of order O(N ) since it is a linear combination of phases and two of them are weighted by huge coefficients, highlighted in blue in (16). The anti-symmetric and symmetric components will be at most of order O(N 0 ) at the upper boundary since they are linear combinations of phases but the prefactors in (16)(17) are at most of order O(N 0 ). We can imagine the values inside the strip to smoothly interpolate between their two boundaries. Then we would conclude that anti-symmetric and symmetric representations would remain small while singlet would be the much larger dominant contribution. This is of course what we expect in a large N theory since there is a lot of phase space to scatter into when we form a singlet. Technically, in large N vector theories, dimmer-like colour loops dominate. Indeed, in all the examples we will now explore, it is by analyzing the singlet component that we can learn a great deal very easily.

Maximum singlet effective quartic coupling
Suppose we are after the S-matrix with the largest singlet component in the middle of the physical strip, that is at θ = iπ/2, in a theory with no bound states (i.e. no poles inside the physical strip). The value at the middle of the strip can only be as large as its value at the boundary of the strip by the maximum modulus principle. At the upper boundary, the singlet S-matrix can have at most absolute value N/2 + N/2 = N if the phases S sym and S anti are opposite to each other (at least for real θ) so that the terms in blue in (16) add up. So all in all, we have a function without poles which satisfies |S singlet (θ)| ≤ 1 at real θ and |S singlet (iπ + θ)| ≤ N at the other boundary of the physical strip. One function which clearly saturates these bounds is S optimal singlet (θ) = e θ iπ log(N ) which is indeed a phase for θ ∈ R, equal to N times a phase for θ ∈ iπ + R and real analytic for purely imaginary θ. This is the optimal solution to our problem since the ratio S singlet θ)/S optimal singlet (θ) has absolute value smaller or equal to 1 at both boundaries and thus is at most 1 in the middle of the strip. Hence Now we can test this against the numerics described before. Performing them for various values of N leads to results in figure 11 which beautifully matches with our prediction (20) as we extrapolate to infinite N .
We can also learn about all the S-matrices themselves and their analytic properties. As we saw before, for the singlet component to be as large as possible the symmetric and antisymmetric components ought to be opposite of each other for real θ so that the two blue terms in (16) would add up constructively. In this case, where we impose no poles that force these components to be different, it is natural to assume that S sym = −S anti everywhere. (This is confirmed by the numerics at any finite N .) Then, (16) immediately leads to Note that the last term can be dropped inside the physical strip and its boundary except at the upper boundary where both terms are small and of order O(N −1 ). So, summarizing, we S anti (θ) = −S sym (θ) (24) valid anywhere inside and on the boundary of the physical strip. (But not elsewhere!) At the upper boundary we get S sym (iπ − θ) ∝ 1 N sin( θ π log(N )) so we see that this solution exhibits an infinite amount of zeros (resonances) right at the upper boundary of the physical strip in the symmetric and anti-symmetric channels. We can again compare this expectation with the finite N numerics to see that these zeros are there indeed, see also figure 15 below. These zeroes are a generic feature of the solutions in fact and arose already in another example, see figure 9 above. The appearance of infinite resonances and periodicity along the real θ direction were encountered previously in other models like the ones studied in [15] and [16]. In the context of the S-matrix bootstrap resonances were recently discussed in [12].
Finally, note that while (22)(23)(24) are valid only in the physical strip, we can extend the solution to any value of θ since we can combine crossing symmetry with unitarity saturation which reads S rep (−θ) = 1/S rep (θ). For example, using unitarity we see that these infinitely many zeros at the upper boundary of the physical sheet become poles on the lower boundary of the so called mirror sheet Im(θ) ∈ [−π, 0]. Then we can use crossing to learn about the analytic structure of this solution in the strip above the physical strip where Im(θ) ∈ [0, 2π] and then unitarity again etc. In this way, we can deduce the full analytic structure of this function to be as given as discussed in section 4.2.1 below. We see that at large N the problem simplifies so much that we cannot only solve the maximization problem but also get key insights about the analytical structure of the various S-matrix elements. Indeed, assuming such structure persists at finite N we will be able to guess the full analytic solution for this maximization problem below.

One bound state
Next let us consider the problem of large N theories with a single bound state whose coupling we maximize.
First assume that the bound state is in the singlet representation. From the crossing relation (16) we learn two things. The first one is that the t-channel pole residue is smaller than the s-channel pole by a factor of 1/N so that at large N we can effectively drop it and keep a single pole at θ = iλ (the s-channel pole for a bound state of mass m BS = 2m cos(λ/2)). The second thing was discussed in the last section: to make the singlet component as large as possible it is optimal to anti-align the symmetric and anti-symmetric components; in  that case the two terms in blue in (16) add up so that the singlet component has maximal magnitude N at the upper boundary. The punch-line is that we simply should multiply the solution (19) of the previous section by a function which introduces a single pole at θ = iλ, is real for purely imaginary θ and is a phase in both boundaries of the physical strip. Such function is We thus conclude that the optimal solution is given by and therefore max g 2 sing = (2 sin λ) × (2N λ/π , sin λ) where the first factor is a simple Jacobian since we define the coupling as the residue in s. This indeed matches beautifully with the large N extrapolation of our numerics as indicated in figure 12. We could also, as in the previous section, further determine the other channel S-matrices. The logic is exactly as before and leads to At the boundary of the strip this solution exhibits again a rich pattern of resonances which one can indeed observe in the numerics at finite N .
As a second example, consider a single pole in the anti-symmetric representation or in the symmetric representation. From (25) we see that an s-channel pole in either of these representations induces a t-channel pole in the singlet component whose magnitude is N/2 larger. In the right hand side of (25) the anti-symmetric and symmetric representations show up with opposite prefactors at large N however the residue in the anti-symmetric case comes with a further minus sign -recall (12). Therefore, we see that both problems lead to the very same maximization problem as far as the singlet component is concerned: that of finding a t-channel pole in the singlet channel to be as large and negative as possible. Following the same logic as in the previous case we thus conclude that the optimal solution for the singlet representation would read Computing the residue and multiplying by −2/N and the usual Jacobian we then obtain the optimal large N solution which as explained above coincides for these two problems: max g 2 anti single pole anti = max g 2 sym single pole sym = 8 sin 2 (λ)N −λ/π .
In figures 13(a) and 13(b) we see that the numerics for these two cases, although different at finite N , beautifully converge towards this universal prediction as we extrapolate them to infinite N .
Note that in all three cases discussed in this section, even when we wanted to maximize couplings in other channels, the key was to study the singlet component at large N . This is in line with the intuition advocated at the beginning of this section.

Two bound states with same mass
Let us consider a final large N example where the theory has two bound states with equal mass, one in the singlet representation and another one in the anti-symmetric representation and we are maximizing the anti-symmetric coupling. This case is relevant since it contains the Gross-Neveu model for the particular value of the masses m sing = m anti = 2 cos(π/(N − 2)). As before, we want to find the large N optimal solution for the singlet representation and use (25)   bound state in the singlet representation, we should have an s-and a t-channel pole in that component. However, in contrast to (27), we cannot omit the t-channel pole because the anti-symmetric representation has also an s-channel pole, i.e. when evaluating (25) at θ = iλ we need to keep the first two terms on the right. We thus conclude that the optimal solution is given by 20 S optimal From (25) we get the following large N equation between the residues which after the replacement of the optimal solution (32) leads to the maximum coupling max(g 2 anti ) = (2 sin λ) × N −λ/π − N −2+λ/π 4|tan λ| .
The term N −2+λ/π in the previous equation is subleading for small λ but becomes relevant as we reach λ = π (e.g. as the mass of the bound states approaches zero). The divergence at m 2 BS = 2 given by the screening of singlet and anti-symmetric poles is nicely reproduced by the factor of |tan λ| in (34). In figure 14 we confirm that this result nicely matches with the large N extrapolation of the numerics discussed in section 2.3.

Finite N
In this section we discuss the analytic solutions at finite N found in two cases: the maximization of S sing (iπ/2) without bound states, and a single bound state in the singlet representation. As we will see, these solutions preserve some of the features found at large N in the previous section and nicely match the numerical results at finite N .

Maximization/minimization singlet channel without bound states
Consider a theory without bound states where we maximize S sing (iπ/2). As discussed in section 4.1.1, at large N the symmetric and anti-symmetric representations antialign so that the singlet channel is as large as possible. Somewhat surprisingly, our numerics show that in the present case the equality S anti (θ) = −S sym (θ) holds at finite N . We shall impose this condition in our derivation below. Another important feature observed at large N is the presence of a sequence of zeros in S sym (θ) at the upper boundary of the physical strip. For each of these zeros, unitarity+crossing implies there is an infinite sequence of zeros and poles in higher sheets. Since this is a general property of O(N ) S-matrices, let us carefully explain the logic leading to this structure.
Let us examine the consequences of a zero in the physical strip situated at some value θ * for a given representation. From unitarity we see that there should be a pole in the θ ∈ [−iπ, 0] strip at −θ * for the same representation. Through the crossing equations, this pole implies the presence of more poles at iπ + θ * for all representations in the θ ∈ [iπ, 2iπ] strip. We can use unitarity again to observe that all representations should have a zero at −iπ − θ * . Furthermore, crossing tells us that there are zeros as well at 2iπ + θ * . 21 Continuing this logic we are naturally lead to the functions F a (θ) defined in (8) and present in the integrable solutions of section 2.2. These functions are explicitly unitary and encode some of the properties the S-matrices should satisfy from crossing.
Going back to the case at hand, we assume (based on numerical observations) that the zeros at θ ∈ iπ + R are periodically spaced by π 2 /ν, where ν is a parameter to determine. Following these zeros through crossing and unitarity as in the argument above we arrive at the following ansatz Replacing this ansatz into the crossing equations (16)(17)(18) fixes the remaining parameter ν to be N = 2 cosh(ν). The analytic structure of this solution is described in figure 15. 22 Another way to write the S-matrices in (35) is sinh ν π (iθ + (2l − 1)π) sinh ν π (iθ − 2lπ) sinh ν π (iθ − (2l − 1)π) sinh ν π (iθ + 2lπ) 21 It is easy to see in the crossing equations (16)(17)(18)) that a single pole (zero) at θ * in one of the representations does (not) imply extra poles (zeros) at iπ − θ * . If we have more than one pole at θ * , the functions may conspire in such a way that their residues cancel so that there is no pole at iπ − θ * in a given representation. In contrast, zeros at θ * in all representations does always imply zeros at iπ − θ * for all representations. 22 For N = 2 the solution (35) reduces to S(θ) = − π−iθ π+iθ , 1, −1 F π (−θ). Figure 15: Analytic structure of the solution to the maximization of the singlet representation S sing (iπ/2) without bound states presented in (35). The anti-symmetric and symmetric representations are the same up to a minus sign. The horizontal spacing between zeros and poles is π 2 /ν, where ν is given by N = 2 cosh(ν). The central poles and zeros in dark green form the structure that remains in the N = 2 limit. The minimization problem differs by an overall minus sign.
which allows one to make direct contact with the large N results of section 4.1.1 rather trivially. Indeed, at large N we have ν log(N ) so that all arguments in these expressions are very large. As such, inside each θ strip we can simplify all sinh's to one of its exponentials. In this way we see that inside any θ strip any component is given by a simple power of N times e ∓iθν/π . For instance, in the physical sheet we recover (22) and (23) (recall that inside the strip we can drop the last term in (23)). The above representation is also more suitable for numerical evaluation, since when truncating the product to some l max we keep the poles and zeros closer to the physical strip. In fact, it was using this representation with l max = 50 that we generated the yellow curve in figure 11.
If instead we want to minimize S sing (iπ/2) we would be led to the same solution with an overall minus sign.
Given that the solution (35) saturates unitarity, it is natural to ask if it also obeys Yang-Baxter equations. As one can check it does satisfy factorization and in fact was written before in the appendix of [6]! 23 To our knowledge, a physical model corresponding to this integrable S-matrix is yet to be identified. 24

One bound state in singlet representation
Let us now derive the finite N solution for the case of a theory with a single bound state in the singlet representation.
The first thing we notice from the numerics is that, as in the previous example, the large N equality S anti (θ) = −S sym (θ) holds at finite N . Secondly, there are also zeros at the upper boundary of the physical strip. This time however there is no zero at θ = iπ and what was a single sequence of zeros in the previous example splits now into two strings starting at θ = iπ ± ζ, where ζ is a parameter to determine. In the following we will assume that S anti (θ) = −S sym (θ) and that the spacing between zeros in each string is again given by π 2 /ν which is consistent with numerical observations (see green structure in figure 16). 25 Let us now consider the implications of a single pole at θ = iλ in the singlet channel corresponding to the bound state. For this we use the crossing relations simplified by the condition S anti (θ) = −S sym (θ). One of the equations reads If we evaluate this equation at θ = iλ, we conclude that the t-channel poles should cancel each other so that S sym (iλ) remains finite (recall that there is only one s-channel pole for the singlet representation). In other words, their residues must be the same. One can check using unitarity that this implies the existence of a double zero (pole) for S sym (2iπ − iλ) (S sym (−2iπ + iλ)). 26 By the crossing+unitarity logic described before, the double pole will propagate to higher sheets creating the pattern depicted in blue in figure 16. Let us stress the general property of the O(N ) S-matrices we just derived: whenever a pole appears only in the singlet representation there will be a tower of double poles and zeros in higher sheets for all representations.
The next thing to note is that the value S(0) fixes S(iπn) for any n ∈ Z through crossing. In particular, if we read from the numerics the value S(0) we can compute the sign of the Smatrices at S(iπn). In the present case we have S(0) = (1, 1, −1) which leads to alternating 23 There, equation A.5 should have an infinite product instead of a sum. 24 This integrable S-matrix bears strong resemblance with the one in [17] where also S anti (θ) = −S sym (θ) (or equivalently σ 2 (θ) = 0 in the decomposition (2)). In that case there is also a periodicity related to N , which takes values |N | < 2 and for N = 0 this solutions beautifully relates to self-avoiding polymers, see also [18,19]. The relative coefficients of the various S-matrix channels in [6] and [17] differ by a simple analytic continuation ν → iµ while the overall factor is related in a more involved way. It would be interesting to explore further the connections between these two solutions. 25 The fact that the there are two copies of strings of zeros at θ = iπ±ζ ±nπ 2 /ν follows from real analyticity S * (s) = S(s * ) or S * (θ) = S(−θ * ). 26 We have that S sing,sym (θ) ∼ r θ−(iπ−iλ) and S sym (iπ + θ) = 1  The double crosses and circles represent double zeros and poles, respectively. The parameters µ i specify the position of the new zeros and poles appearing in the singlet representation. The first zeroes at the upper boundary of the physical strip depicted in green are located at θ = iπ ± ζ and the rest are spaced by multiples of π 2 /ν. signs every iπ for all representations. For S sing (θ), this means that as we move in the imaginary line segment θ ∈ [0, iπ], we should start with a positive function which then has the s-and t-channel poles and end with a negative value. Since each pole changes the sign of the function, we conclude that there should also be a zero in this line (another pole would imply there are other bound states) at some position θ = iµ 1 − iπ. For S anti (θ) = −S sym (θ) the single t-channel pole is consistent with the change of sign. Following the new zero through crossing and unitarity -keeping in mind the double poles and zeros produced in higher sheets-we arrive at the orange structure in figure 16.
If we repeat the exercise in the second strip θ ∈ [iπ, 2iπ] we easily see the need for a new zero in S sing (θ) (yellow cross in figure 16). This is actually the case in every strip, so that we end up with infinite parameters µ n which label the position of the first new pole in the  Table 3: Values for the parameters ζ and µ i≤4 with N = 5 and m 2 sing = 3.1 (λ = 0.98), obtained by the method explained in appendix A.
n−th strip.
Collecting the analytic structure discussed above and depicted in figure 16, we write the ansatz for the S-matrix where µ i and ζ are parameters to determine in function of λ and N and we have defined The remaining parameters µ i and ζ can be fixed by using crossing symmetry. As explained in details in appendix A, the system of equations we need to solve has the nice form which can be solved numerically very efficiently and to arbitrary accuracy. An example of the values obtained by solving the above equations is given in table 3. We suspect (39) should have a nice and simple physical interpretation; it would be very interesting to find it.
The cases discussed above satisfied S anti (θ) = −S sym (θ) which greatly simplified our task of finding a solution to unitarity and crossing equations. For all other cases we studied this equality no longer holds, so that finding the full analytic solution is more complicated. Although we do not have the analogue of (39) for other cases, we were able to obtain (at least some of) the analytic structure of other solutions. For example, considering bound states in the singlet and anti-symmetric representations with the same mass, one finds -for the mass range m sing = m anti < m GN -the poles and zeros depicted in figure 17. This case exhibits an ever richer structure with pairs of zeros inside the physical strip and patterns that repeat each other in a fractal way as we move outside the physical strip. Other maximization problems lead to S-matrices of roughly the same complexity. It would be nice to obtain the full analytic solution in these cases.

Discussion
We analysed the space of O(N ) symmetric two-dimensional relativistic quantum field theories whose lightest particles transform in the fundamental representation. We carved out part of this space by looking for 2 → 2 candidate S-matrices which maximize various couplings (be them cubic couplings given by the physical residues of the S-matrix or effective four point couplings given by the S-matrix amplitudes evaluated in some symmetric point in the physical sheet).
In the boundary of this theory space we identified known integrable S-matrices describing the excitation scattering in the O(N ) non-linear sigma model (NLSM), the O(N ) Gross-Neveu model, the kinks of the Sine-Gordon theory or the yet unknown model for solution (35) first discovered by Hortacsu, Schroer and Thun in [6]. The NLSM S-matrix was beautifully identified to lie at a cusp of the space of unitarity and crossing invariant S-matrices in [8]. It would be very interesting to connect our various bounds to their proposed picture for the S-matrix space and its cusps. We studied, for instance, the other cusp found in that work and found that its analytic structure is reminiscent of the general pattern found here, with infinitely many resonances at the boundary and inside the physical sheet.
We can already anticipate that the less known and quite exotic integrable solution (35) is also a cusp. This is yet another reason to look for the corresponding physical theory. We expect it should describe some planar processes as the self-avoiding polymers of [17] since the only allowed color structures are S kl ij (s) = σ 1 (s) + σ 3 (s) .
It would be very interesting to explore the possible connections with the models described in [18] and [19].
So far, we encountered already at least three cusps in the space of theories where integrability holds: The NLSM found in [8], the free theory of section 3.1 and the more exotic solution of [6] in section 4.2.1. It would be fascinating to understand mathematically the emergence of Yang-Baxter equations at the cusps of this space. While natural from a physical point of view, the emergence of these cubic equations remains a completely mathematical mystery in this context. Can we find other solutions to YB (perhaps some never discovered before) for other symmetry groups?
Also at the boundary of this physical space we identified even more exotic S-matrices. These saturate unitarity so that where the sum is over the three possible representations rep = (singlet, anti, sym) . These theories have an extremely rich pattern of resonances with infinitely many of them in the various sheets of the Mandelstam plane as explored in the main text (see for example figures 16 and 17). 27 Because they saturate unitarity, these optimal solutions admit zero particle production so that S exotic 2→3 (s) = S exotic 2→4 (s) = · · · = 0 just like the integrable models. However, contrary to those these exotic S-matrices do not obey the Yang-Baxter relations.
Can we have proper physical theories with no particle production and no Yang-Baxter? Probably not. Since Yang-Baxter is not satisfied, the 3 → 3 S-matrix cannot consistently factorize into a product of two-to-two S-matrices (like for an integrable theory). If this S-matrix is not factorized (and hence supported on an additional conservation δ function consistent with this factorization), it should have a regular part. This regular part can be analytically continued by crossing into S 2→4 which thus should not be strictly zero [22,23]. These exotic S-matrices, with strictly zero particle production are thus likely unphysical.
Having said that, perhaps we need not completely discard these S-matrices altogether. After all, while we can probably not accept a strictly zero particle production in a theory whose two body S-matrix does not obey the factorization conditions, can this particle production be very small? In other words, could these S-matrices, with their very rich analytic properties and their infinitely many resonances, be close enough to real physical S-matrices which simply have little particle production?
To properly analyze this question we could try to bootstrap not only the two body Smatrix but also multi-particle components, such as S 2→4 and S 3→3 , altogether and figure out the fate of these S-matrices in this more complete setup. This is a fascinating but very involved problem which we are currently attacking together with Alexandre Holmrich. Further intuition about the exotic S-matrices might be found in the study of the "special form factors" introduced in [24] properly generalized to a setup with O(N ) symmetry. What will these tell us about deformations around the NLSM or the Gross-Neveau model for instance? Can this shed light over the exotic S-matrices we found for masses near those of the Gross-Neveu model?
In the meantime, we can preform some phenomenological games and add some particle production by hand. More precisely, we can redo the numerics modifying the unitarity constraints to read |S rep (s)| 2 ≤ f rep (s) where f rep (s) = 1 − σ (abs) rep is equal to 1 until we reach inelastic thresholds and decays to zero for larger values of s. Of course we do not know what these functions are so the best we can do we can do in this simple exercise is to set them to some arbitrary simple functions with these properties and see what one gets. If all functions are taken to be the same (for example, we could set f rep (s) = exp(−α(s−16m 2 ) 2 θ(s−16m 2 )) starting at the 4 particle threshold) then the new optimization problems are actually trivially related to the ones in the main text as S with particle production (θ) = F (θ)S main text (θ) where F (θ) is given by formula (11) in [2]. In particular, because the solutions are simply dressed by this overall function, all the rich pattern of zeros and poles which we found are still present with particle production. More realistically, different channels will produce particles at different rates so things will be less trivial. At large N we can repeat most of the analysis in section 4.1 for fixed imposed particle production. We can also easily run numerics with modified unitarity. The conclusion of this simple analysis -using some randomly chosen particle production functions -is that the various resonances are still there but their positions move around slightly; some acquire some imaginary part, some stay on the t-channel cut. For low energy the S-matrices are comparable to the unphysical ones we got without any particle production provided particle production is not large or changes too abruptly. This is perhaps not so surprising since particle production is a high energy effect of sorts while our maximization probes mostly the low energy region of the S-matrices.
If the infinitely many resonances present in the integrable solution (35) and more exotic S-matrices are indeed there for some physical theory (albeit in some slightly shifted positions in the second case), it would be fascinating to unveil their physical origin. grant 2016/01343-7 for partial financial support. Now we only need to consider i max +1 different values of n to solve the system of equations. The values obtained in table 3 were found by evaluating the effective equations (45) for n = 1, . . . , 5 and i max = 4. We performed this exercise for various m sing and N and found always perfect agreement with the position of zeros and poles in our numerics.