# Scattering amplitudes from dispersive iterations of unitarity

Piotr Tourkine, Alexander Zhiboedov

► To cite this version:

Piotr Tourkine, Alexander Zhiboedov. Scattering amplitudes from dispersive iterations of unitarity. Journal of High Energy Physics, 2023, 11, pp.5. 10.1007/JHEP11(2023)005. hal-04053799

HAL Id: hal-04053799

https://hal.science/hal-04053799v1

Submitted on 15 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Scattering amplitudes from dispersive iterations of unitarity

Piotr Tourkine$^a$ and Alexander Zhiboedov$^b$

$^a$Laboratoire d'Annecy-le-Vieux de Physique Théorique (LAPTh),
CNRS & Université Savoie Mont-Blanc,
9 chemin de Bellevue, 74940, Annecy, France

$^b$CERN, Theoretical Physics Department,
1 Esplanade des Particules, CH-1211 Geneva 23, Switzerland

E-mail: piotr.tourkine@lapth.cnrs.fr, alexander.zhiboedov@cern.ch

**ABSTRACT:** We present and numerically implement a computational method to construct relativistic scattering amplitudes that obey analyticity, crossing, elastic and inelastic unitarity in three and four spacetime dimensions. The algorithm is based on the Mandelstam representation of the amplitude and iterations of unitarity. The input for the iterative procedure is given by the multi-particle double spectral density, the S-wave inelasticity, and the value of the amplitude at the crossing-symmetric point. The output, obtained at the fixed point of the iteration of unitarity, is a nonperturbative scattering amplitude. The amplitudes we obtain exhibit interesting features, such as non-zero particle production, intricate high-energy and near the two-particle threshold behavior. Scattering amplitudes obtained by initializing the iteration process with zero (or small) multi-particle input end up close to saturating the S-matrix bounds derived by other methods. There is a version of the iterative algorithm that is directly related to Feynman diagrams: it effectively re-sums infinitely many two-particle reducible planar Feynman graphs in the $\phi^4$ theory, which remarkably produces a unitary nonperturbative scattering amplitude function. Finally, we discuss how the algorithm can be further refined by including multi-particle unitarity.

**KEYWORDS:** Nonperturbative Effects, Scattering Amplitudes

ARXIV ePRINT: 2303.08839
# Contents

## 1 Introduction and summary of results 2

1.1 Basic idea of the iteration algorithm 6

1.2 Summary of results and plan of the paper 8

## 2 Algorithm 9

2.1 Basic concepts 10

2.2 Iteration algorithm 13

2.3 Numerical implementation 15

2.3.1 Discretization 15

2.3.2 Numerical integrals 18

2.3.3 Convergence 19

2.3.4 Computation times 20

2.4 Fixing $\eta_{MP}$: analyticity in spin up to $J = 0$ and Feynman graphs 20

2.5 Graph interpretation and summability 21

## 3 Toy-model amplitudes: zero double discontinuity 24

3.1 $d = 3$ 25

3.1.1 The algorithm 25

3.1.2 Numerical implementation: hotstart near the two-particle threshold 27

3.2 $d = 4$ 34

3.2.1 The algorithm 34

3.2.2 Numerical implementation: “hotstart” in the Regge limit 36

3.3 Proof of existence of the amplitudes in $3 < d < 4$ 40

3.4 $d > 4$: nonrenormalizability and unitarity 42

## 4 Quasi-elastic amplitudes in $d = 3$ 43

4.1 2QE amplitude for $\lambda = -3\pi$ 43

4.1.1 Single and double spectral functions 43

4.1.2 Partial waves, impact parameter 45

4.1.3 Amplitude 45

4.2 2PR amplitude 47

4.3 Coupling dependence 47

## 5 Quasi-elastic amplitudes in $d = 4$ 49

5.1 2QE amplitude for $\lambda = 5\pi$ 49

5.1.1 Single and double spectral functions 49

5.1.2 Partial waves, impact parameter 49

5.1.3 Amplitude 50

5.2 2PR amplitudes 53

5.3 Coupling-dependence 53
<table>
  <tbody>
    <tr>
      <td>6</td>
      <td>Multi-particle double discontinuity: a case of acnode in d = 4</td>
      <td>54</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Computing the s<sub>67</sub> integral</td>
      <td>56</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Analytic continuation in t</td>
      <td>57</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Nonperturbative lesson</td>
      <td>59</td>
    </tr>
    <tr>
      <td>6.4</td>
      <td>Atkinson scattering-from-acnode</td>
      <td>60</td>
    </tr>
    <tr>
      <td>7</td>
      <td>Low energy observables</td>
      <td>61</td>
    </tr>
    <tr>
      <td>7.1</td>
      <td>(&tau;<sub>0,0</sub>, &tau;<sub>0,1</sub>) plots in d = 3, d = 4</td>
      <td>63</td>
    </tr>
    <tr>
      <td>7.2</td>
      <td>c<sub>2</sub> plots</td>
      <td>65</td>
    </tr>
    <tr>
      <td>7.3</td>
      <td>Grid dependence</td>
      <td>67</td>
    </tr>
    <tr>
      <td>8</td>
      <td>Mandelstam representation and quantum gravity</td>
      <td>68</td>
    </tr>
    <tr>
      <td>9</td>
      <td>Conclusions and open directions</td>
      <td>71</td>
    </tr>
    <tr>
      <td>9.1</td>
      <td>Summary</td>
      <td>71</td>
    </tr>
    <tr>
      <td>9.2</td>
      <td>Open directions</td>
      <td>73</td>
    </tr>
    <tr>
      <td>A</td>
      <td>Definitions</td>
      <td>77</td>
    </tr>
    <tr>
      <td>B</td>
      <td>Grid generation</td>
      <td>78</td>
    </tr>
    <tr>
      <td>C</td>
      <td>Elastic amplitudes in d = 2: product of CDD zeros</td>
      <td>81</td>
    </tr>
    <tr>
      <td>D</td>
      <td>Analytic Regge limit</td>
      <td>82</td>
    </tr>
    <tr>
      <td>E</td>
      <td>Proof of convergence à la Atkinson for toy-model in 3 &lt; d &lt; 4</td>
      <td>83</td>
    </tr>
    <tr>
      <td>E.1</td>
      <td>Bounding the norm ||&rho;||</td>
      <td>84</td>
    </tr>
    <tr>
      <td>E.1.1</td>
      <td>Bounding the principal value integral</td>
      <td>85</td>
    </tr>
    <tr>
      <td>E.2</td>
      <td>Proving the contracting property of the map</td>
      <td>87</td>
    </tr>
    <tr>
      <td>E.3</td>
      <td>Generalization to d = 3 and d = 4</td>
      <td>89</td>
    </tr>
    <tr>
      <td>F</td>
      <td>Perturbative &phi;<sup>4</sup></td>
      <td>90</td>
    </tr>
    <tr>
      <td>F.1</td>
      <td>d = 4</td>
      <td>90</td>
    </tr>
    <tr>
      <td>F.2</td>
      <td>d = 3</td>
      <td>91</td>
    </tr>
    <tr>
      <td>G</td>
      <td>Aks graph inelasticity</td>
      <td>93</td>
    </tr>
    <tr>
      <td>G.1</td>
      <td>d = 3</td>
      <td>93</td>
    </tr>
    <tr>
      <td>G.2</td>
      <td>d = 4</td>
      <td>96</td>
    </tr>
    <tr>
      <td>H</td>
      <td>Full coupling dependence of d = 4 2QE and 2PR partial waves</td>
      <td>99</td>
    </tr>
    <tr>
      <td>I</td>
      <td>Impact parameter space and the Froissart-Gribov formula</td>
      <td>100</td>
    </tr>
    <tr>
      <td>J</td>
      <td>Acnode evaluation</td>
      <td>102</td>
    </tr>
  </tbody>
</table>
# 1 Introduction and summary of results

The S-matrix bootstrap arises from the delicate tension between the relativistic concept of causality and the quantum-mechanical notion of unitarity [1–3]. Causality manifests itself through analyticity of scattering amplitudes, as well as crossing symmetry which relates different parts of the complex scattering energy plane to different scattering processes. Unitarity relates different scattering processes to each other in a non-linear fashion.<fnref n="1" /> At its infancy, the S-matrix bootstrap almost exclusively attempted to describe scattering of hadrons which turned out to be a hard problem, still unsolved today. Its multifarious modern incarnation pursues a broader goal of exploring the space of relativistic, causal, quantum-mechanical theories (see e.g. the recent reviews [4, 5]). The full characterization of this space includes the tasks of finding its boundaries, and locating in it theories of physical interest, such as QCD. These tasks come with the general belief that if a theory resides on the boundary, it is conceivable that it can be “solved”. This phenomenon is familiar in the domain of the conformal bootstrap [6], and 2d S-matrices [7]. To this day, it remains an open question as to whether there exist solvable, physical theories whose S-matrices can be efficiently calculated nonperturbatively in spacetime dimension $d > 2$. More generally, the S-matrix bootstrap techniques can be combined with inputs from experiment, perturbative calculations, lattice simulations, or other methods [8], to pinpoint a given physical theory.

Various approaches to the S-matrix bootstrap enforce various amounts of the basic constraints of analyticity, crossing and unitarity. The more analyticity, crossing and unitarity is put in, the more stringent the constraints on the space of amplitudes are. The novelty of the techniques that we develop and implement in this paper is that

* we enforce the constraints of two-particle (elastic and inelastic) unitarity fully;

* we put in certain constraints dictated by multi-particle unitarity, e.g. the support of the multi-particle double spectral density;

* we control the amplitudes up to arbitrarily high energies and for any scattering angles (or impact parameters).

We implement these constraints by solving unitarity iteratively, while keeping analyticity and crossing symmetry manifest at every step. The basic idea of the iteration algorithm goes back to the early days of the S-matrix bootstrap and the works of Mandelstam [9–11], Chew and Frautschi [12], Ter-Martirosyan [13] and others.<fnref n="2" /> The idea was significantly sharpened and further developed by Atkinson in a series of papers [14–17], and various follow-ups [18–23], where the *existence* of analytic, crossing-symmetric, unitary amplitudes, and the convergence of the iteration algorithm based on the Mandelstam representation, were rigorously proven using functional analysis methods. A pedagogical introduction to these ideas can be found in the lecture notes by Atkinson [24].

<footnote n="1">¹The simplest example being the optical theorem which equates the imaginary part of the forward amplitude to the probability to scatter in an arbitrary final state.</footnote>

<footnote n="2">²In the lecture notes [11], Mandelstam mentions in this context a preprint by Kenneth Wilson, but we did not manage to find it.</footnote>
To the best of our knowledge, this approach to the S-matrix bootstrap has never been implemented in $d > 2$.<fnref n="3" /> In other words, while the equations have been written long time ago, no solutions have been produced. This is the problem we address in the present paper using numerical methods.

The precise goal of this paper is to apply these techniques to the simplest setup: $2 \to 2$ scattering amplitude of identical scalar particles in three and four spacetime dimensions. The amplitudes functions<fnref n="4" /> that we construct here satisfy maximal analyticity, crossing, elastic and inelastic unitarity. They exhibit many interesting features, such as clearly identifiable particle production (consistent with the Aks theorem [27]),<fnref n="5" /> the detailed expected analytic structure coming from the Landau curves, and a nontrivial fixed angle and Regge behavior [29]. They are specified by a certain multi-particle production data which is the subject of our choice, and which serves as an input for the algorithm. To describe a given physical theory, such as QCD, one then needs to develop a model for the multi-particle production data. We describe all these notions in detail below.

Before delving into the details of these techniques, let us mention that an equivalent, maybe more accessible way to state what our algorithm achieves is a crossing-symmetric unitarization of scattering amplitudes via a re-summation of a class of two-particle reducible planar graphs, whose vertices are made of the quartic coupling, and an additional effective vertex which accounts for multi-particle physics, see figure 1.<fnref n="6" /> In general, these graphs are not Feynman diagrams, but signify the discontinuity structure of the amplitude in the spirit of generalized unitarity [31–33], see also [34]. For a specific, restricted choice of the input parameters to be described below, our algorithm literally re-sums a class of Feynman diagrams, allowing to make direct connection with perturbation theory at low energies.

We now explain our set-up. Firstly, we work with a gapped theory with a single stable particle of mass $m$. We parameterize the two-to-two scattering amplitude $T(s, t)$ using the Mandelstam representation [10]

$$ 
\begin{aligned}
T(s,t) &= \lambda + B(s,t) + B(s,u) + B(t,u), \\
B(s,t) &= \frac{1}{2} \int_{4m^2}^{\infty} \frac{ds' \rho(s')}{\pi (s' - s_0)} \left( \frac{s - s_0}{s' - s} + \frac{t - s_0}{s' - t} \right) \\
&\quad + \int_{4m^2}^{\infty} \frac{ds' dt' (s - s_0)(t - t_0)\rho(s', t')}{\pi^2 (s' - s)(t' - t)(s' - s_0)(t' - t_0)},
\end{aligned} 
\quad (1.1) 
$$

where $\rho(s)$ is called the single spectral density, $\rho(s,t) = \rho(t,s)$ is called the double spectral density,<fnref n="7" /> and we set $s_0, t_0, u_0$ to the crossing-symmetric value, $s_0 = t_0 = u_0 = \frac{4m^2}{3}$. The

<footnote n="3">³In $d > 2$, we could only find one attempt to numerically implement a simplified unitarity-inspired set of equations in $d = 4$ by Boguta [25]. In $d = 2$ the algorithm was successfully implemented in our previous work [26], as a warm-up for the present paper.</footnote>

<footnote n="4">⁴We call them amplitude functions to emphasize that we do not discuss an algorithm to generate the S-matrix beyond the $2 \to 2$ amplitude. In the rest of the paper, to avoid cluttering we will simply call them amplitudes.</footnote>

<footnote n="5">⁵The amplitudes that we obtain have much larger inelasticity compared to the very weak asymptotic bound obtained in [28].</footnote>

<footnote n="6">⁶See e.g. [30] for a review of standard approaches to unitarization.</footnote>

<footnote n="7">⁷We sometimes interchangeably use the term single/double spectral function.</footnote>
subtraction constant $\lambda$ has dimensionality $m^{4-d}$ and we will always measure it in the units of external particle mass $m$. The representation (1.1) makes analyticity and crossing symmetry

$$T(s, t) = T(t, s) = T(u, t), \quad s + t + u = 4m^2, \qquad (1.2)$$

manifest.

Let us briefly comment on the assumptions that go into the Mandelstam representation. Firstly, the analyticity domain of the amplitude given by the Mandelstam representation, the so-called maximal analyticity in which the only singularities of the amplitude are those required by unitarity, i.e. the amplitude is analytic in $(s, t, u) \in \mathbb{C} \setminus \{\text{unitarity cuts}\}$. Maximal analyticity goes beyond the axiomatic analyticity domain proven from the axioms of QFT [2, 35]. It is however a natural domain of analyticity for the scattering of lightest particles, where it is believed that no “anomalous thresholds” arise, see [34, 36, 37] for recent discussions. Maximal analyticity plays a central role in our analysis since we use it to analytically continue unitarity relations, as reviewed for example in [34, 38]. Secondly, the Mandelstam representation with $N$ subtractions assumes that the amplitude is polynomially bounded, $|T(s, t)| < |s|^N$ for any $t$ on the physical sheet. For our purposes, the physical sheet can be defined by analytic continuation from the kinematical region $0 < s, t, u < 4m^2$ to the complex plane without crossing any unitarity cuts.<fnref n="8" />

The nontrivial task therefore is to impose unitarity on (1.1). To discuss unitarity, it is convenient to introduce the partial wave expansion of the amplitude. In $d=4$ we have<fnref n="9" />

$$T(s,t) = 16\pi \sum_{J=0, J \text{-even}}^{\infty} (2J+1)f_J(s)P_J\left(1 + \frac{2t}{s-4m^2}\right), \qquad (1.3)$$

$$S_J(s) = 1 + i \frac{\sqrt{s-4m^2}}{\sqrt{s}} f_J(s), \qquad (1.4)$$

where $P_J(z)$ are the usual Legendre polynomials.

It is convenient to distinguish the following unitarity conditions:

**Elastic unitarity.** For $s$ between the elastic and inelastic thresholds, we have:

$$|S_J(s)| = 1, \quad 4m^2 \le s < s_{\text{MP}}, \qquad (1.5)$$

where $s_{\text{MP}}$ is the minimal energy required for multi-particle production. This condition exists only in the gapped theories, where production of two and three or more particles are separated by energy. Even though it is most simply formulated for positive integer $J$, elastic unitarity also holds for complex $J$ with $\text{Re}J > J_0$, where $\lim_{|t| \to \infty} |T(s, t)| \lesssim |t|^{J_0}$ for $s_{\text{MP}} > s \ge 4m^2$, see e.g. [40] and [38].

<footnote n="8">⁸One might also wonder, why instead of the unproven Mandelstam representation not to iterate crossing-symmetric dispersion relations, see e.g. [39]? The reason is that in the latter case the required analytic properties of the amplitude are not manifest.</footnote>

<footnote n="9">⁹Here we present the formula in $d=4$ not to clutter the notations. In the main text and appendix A we present $d$-dimensional formulas.</footnote>
The *Mandelstam equation*, first derived in [9], is an expression of elastic unitarity at complex *J* and it takes the form<fnref n="10" />

$$ \rho(s,t) = \frac{(s - 4m^2)^{\frac{1}{2}}}{(4\pi)^2 \sqrt{s}} \int_{z_1}^{\infty} d\eta' \int_{z_1}^{\infty} d\eta'' \theta(z - \eta_+) \frac{T_t(s + i\epsilon, t(\eta'))T_t(s - i\epsilon, t(\eta''))}{\sqrt{(z - \eta_ -)(z - \eta_ +)}}, \quad 4m^2 \le s < s_{\text{MP}} \qquad (1.6) $$

where $\eta', \eta''$ are cosines of the complexified scattering angles that parameterize the two-particle unitarity cut, $\eta_{\pm} = \eta'\eta'' \pm \sqrt{\eta'^2 - 1}\sqrt{\eta''^2 - 1}$, $t(\eta) = \frac{s-4m^2}{2}(\eta - 1)$, $z_1 = 1 + \frac{8m^2}{s-4m^2}$, and $T_t \equiv \frac{T(s,t+i\epsilon)-T(s,t-i\epsilon)}{2i}$ is the t-channel discontinuity of the amplitude. Via the Froissart-Gribov formula, see e.g. [38, 40] and appendix A, (1.6) guarantees elastic unitarity for spins Re[$J$] > $J_0$, where $J_0$, called the *Regge intercept*, is mentioned above. For the amplitudes discussed in the present paper $J_0 = 0$.

## Inelastic unitarity.

$$ |S_J(s)| \le 1, \quad s \ge s_{\text{MP}}, \qquad (1.7) $$

expresses the fact that the total probability of the initial state of two-particle with angular momentum *J* to go to the two-particle final state is less or equal than one. This constraint can be restated as a positive semi-definiteness condition and was successfully implemented first in [41], and then in many follow-ups.

The condition (1.7) implies at large *J* unitarity for scattering at fixed impact parameter *b*

$$ |S_{\frac{\sqrt{s-4m^2b}}{2}}(s)| \le 1, \quad s \gg m^2, \qquad (1.8) $$

which we also check.<fnref n="11" />

### Multi-particle unitarity.
The conditions above are the simplest consequences of unitarity of the S-matrix $\hat{S}^\dagger\hat{S} = \hat{1}$. Naively they exhaust constraints that can be imposed purely at the level of the $2 \to 2$ scattering amplitude. While it is true if we use unitarity for physical kinematics only, it is not true for the analytic continuation of unitarity away from the physical region. Such an analytic continuation is particularly natural for amplitudes that satisfy maximal analyticity, such as the ones considered in the present paper, or which are common in the literature starting from [41]. In this case multi-particle unitarity manifests itself at the level of the two-to-two scattering amplitude in various ways.

First of all, in this paper, we assume the asymptotic Hilbert space to be the Fock space of multi-particle states made of a single stable particle. This puts natural constraints on the analytic structure of the amplitude. For example, this implies the existence of a set of branch points when $s, t, u = n^2m^2$ (with $n \in \mathbb{N}, n \ge 2$), known as normal thresholds. Another, maybe slightly less familiar fact, is that when combined with maximal (or extended) analyticity this structure also implies that the two-to-two amplitude exhibits a

<footnote n="10">¹⁰Loosely speaking, one obtains it by taking another discontinuity of the usual Im $T \sim \int |T|^2$ elastic unitarity equation.</footnote>

<footnote n="11">¹¹Here, we recall that $\frac{\sqrt{s-4m^2}}{2}$ is the spatial momentum, so $J = \sqrt{s-4m^2b}/2$ is just the standard definition of the angular momentum.</footnote>
complicated but computable structure of the Landau curves, which are the thresholds for the double spectral density $\rho(s, t)$, see e.g. [34, 42, 43]. The double discontinuity of the amplitude that develops along these Landau curves can be sometimes further expressed in terms of the two-to-two scattering amplitude. The simplest example of this type is the Mandelstam equation described above that expresses the double spectral density above the leading Landau curve in terms the square of the discontinuity of the amplitude (1.6), and follows from analytic continuation of elastic unitarity. An example of a condition which follows from the analytic continuation of multi-particle unitarity is the extended regime of validity of the Mandelstam equation below the leading multi-particle Landau curve identified in [34]. More generally, the double spectral density $\rho(s, t)$ that is developed across the leading multi-particle Landau curve can be again expressed in terms of the $2 \to 2$ amplitude leading to extra unitarity constraints that go beyond (1.5) and (1.7). We do not implement these type of constraints in the present paper and leave this task for future work. These extra constraints are nonperturbative analogues of the familiar generalized unitarity relations that we have already mentioned above.

The three unitarity conditions above are nonperturbative. In contrast, in perturbation theory they are only satisfied order by order in the expansion parameter. However, perturbation theory computations guarantee that multi-particle unitarity is satisfied order by order in the coupling as well. This goes well beyond the current nonperturbative S-matrix analysis which focuses almost exclusively on the two-to-two scattering amplitude.

Let us now explain how we solve unitarity by means of an iterative procedure.

## 1.1 Basic idea of the iteration algorithm

The amplitude (1.1) is fully specified by the triad $(\lambda, \rho(s), \rho(s, t))$. We keep $\lambda$ fixed at all times, whereas $\rho(s)$ and $\rho(s, t)$ are subject to iterations using unitarity.

Let us start with the double spectral density $\rho(s, t)$. Given an amplitude that satisfies elastic unitarity, we explain later that it can always be written as

$$ \rho(s,t) = \rho_{\text{el}}(s,t) + \rho_{\text{el}}(t,s) + \rho_{\text{MP}}(s,t), \quad \rho_{\text{MP}}(s,t) = \rho_{\text{MP}}(t,s), \quad (1.9) $$

where $\rho_{\text{el}}(s, t)$ satisfies the Mandelstam equation (1.6), and $\rho_{\text{MP}}(s, t)$ has nontrivial support only in the multi-particle region $s, t \ge s_{\text{MP}}$. In fact, multi-particle unitarity implies that it has support in a smaller region $s, t \ge s_{\text{MP,LC}}(t)$ given by the position of the leading multi-particle Landau curve, see e.g. [34, 42]. In the iteration algorithm considered in this paper we take $\rho_{\text{MP}}(s, t)$ to be given and fixed, and $\rho_{\text{el}}(s, t)$ is computed iteratively using the Mandelstam equation (1.6) that guarantees elastic unitarity of partial waves with $\text{Re}[J] > 0$. Intuitively, $\rho_{\text{MP}}(s, t)$ controls physics which is “multi-particle in all channels.” In the language of Feynman diagrams, it is captured by the two-particle irreducible graphs, see [34] for a detailed explanation.

Because the Mandelstam equation involves only the double-discontinuity, we also need to enforce a separate unitarity condition to make sure that elastic unitarity is satisfied for the $J = 0$ partial wave. Without loss of generality, we can write it as follows

$$ 1 - |S_0(s)|^2 = \eta_{\text{MP}}(s), \quad (1.10) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Diagrams that generate $2 \to 2$ scattering amplitude in $\phi^4$ theory up to three loops. Diagrams inside a dashed frame are those for which the double spectral density is non-zero. In our iteration scheme the tree-level contact diagram (magenta) has the meaning of the fixed value of the amplitude at the crossing-symmetric point $\lambda \equiv T(\frac{4m^2}{3}, \frac{4m^2}{3})$. The black diagrams are then generated by iterations of two-particle unitarity and crossing-symmetric Mandelstam representation. The “frog” diagram (green) describes multi-particle production in the S-wave, which serves as an input for the algorithm $\eta_{\text{MP}}(s)$. The “open envelope” diagram (blue) is the first one that contributes to the multi-particle double spectral density $\rho_{\text{MP}}(s, t)$ which we also take as an input. In this way the iterative algorithm is specified by the triad $(\lambda, \eta_{\text{MP}}(s), \rho_{\text{MP}}(s, t))$. Given multi-particle data, the nonperturbative amplitude is effectively generated by iterations of two-particle unitarity gluing and dispersive crossing-symmetrization via the Mandelstam representation. Importantly, imposing multi-particle unitarity would lead to extra unitarity relations capturing both the frog and the open envelope diagrams.</figcaption></figure>

where $0 \le \eta_{\text{MP}}(s) \le 1$ characterizes particle production in the $J=0$ partial wave. We use (1.10) to compute $\rho(s)$. For $\eta_{\text{MP}}(s)$ we consider two schemes: a) it is *given* and *fixed*; b) it is expressed in terms of $(\lambda, \rho_{\text{MP}}(s, t))$. We describe both approaches in the text.

In the most general case, the iteration process is therefore *initialized* by the triad $(\lambda, \eta_{\text{MP}}(s), \rho_{\text{MP}}(s, t))$ which are kept fixed, and $\rho(s)$ and $\rho(s, t)$ are then computed iteratively using (1.9) and (1.10).

Loosely speaking, we can say that we are reconstructing *scattering* from *production*, specified by $\eta_{\text{MP}}(s), \rho_{\text{MP}}(s, t)$, and *subtraction*, specified by $\lambda$.<fnref n="12" /> At the fixed point of the iteration process the scattering amplitude $T(s, t)$ automatically satisfies:

* maximal analyticity

* crossing

* elastic unitarity

<footnote n="12">¹²As we explain below, one role of subtraction in our dispersion relations is to fix the value of the amplitude at the crossing-symmetric point to $\lambda$.</footnote>
We then separately check that the amplitude also satisfies

• inelastic unitarity

This last condition is not guaranteed by the algorithm, but for the class of amplitudes considered in this paper we will find that inelastic unitarity is always satisfied.

## 1.2 Summary of results and plan of the paper

In this paper, we describe a numerical implementation of the iteration algorithm presented above and report on a class of nonperturbative amplitudes in $d = 3$ and $d = 4$ which we produced through this procedure. Specifically, we have constructed three different classes of amplitudes depending on the input $(\lambda, \eta_{\text{MP}}(s), \rho_{\text{MP}}(s, t))$:

• *toy-model amplitudes*, for which the double discontinuity is zero and scattering in S-wave is purely elastic. These serve a pedagogical purpose and also speed up the full algorithm thanks to a “hot-start” procedure explained below.

• *two-particle quasi-elastic amplitudes* (2QE) correspond to the input $(\lambda, 0, 0)$. This algorithm is characterized by the fact that the $J = 0$ partial wave is purely elastic at all energies. From the point of view of our methods, these are truly minimal amplitudes driven solely by their value at the crossing-symmetric point.

• *two-particle reducible amplitudes* (2PR) correspond to the input $(\lambda, (2.24), 0)$, where we made an extra assumption that the imaginary part of partial waves is analytic in spin all the way to $J = 0$. This fixes $\eta_{\text{MP}}(s)$ in terms of the double spectral density, see (2.24), which in turn is fixed by $\lambda$. We expect that this scheme is equivalent to re-summation of fully two-particle reducible planar Feynman diagrams computed in the renormalization scheme $\lambda = T(\frac{4m^2}{3}, \frac{4m^2}{3})$.

• amplitudes with $\rho_{\text{MP}}(s, t) \neq 0$. All known physical theories are of this type, and we only initiate the exploration of this space in the present paper with a specific example using the “acnode graph” [44].

After we have constructed the amplitudes we can explore their physical properties. These include particle production, near-threshold behavior, the Regge limit, fixed-angle high energy scattering, elastic and inelastic cross sections. We find that they exhibit non-zero but small particle production, therefore we can call them *quasi-elastic amplitudes*.<fnref n="13" /> We can also place various amplitudes in the space of allowed couplings explored recently in [45, 46]. All amplitudes constructed in the present paper live very close to the boundary of the allowed region. The most “extremal” amplitudes are the 2QE amplitudes.

Both 2QE and 2PR amplitudes have $\rho_{\text{MP}}(s, t) = 0$. In the old S-matrix bootstrap literature, the condition $\rho_{\text{MP}}(s, t) = 0$ was called “strip approximation” and it was applied to pion scattering, see e.g. [47]. Based on our results, we expect this approximation might

<footnote n="13">¹³Let us emphasize that, while the third class of amplitudes we construct have an inelastic input, we still refer to it as “quasi-elastic”, because the overall inelasticity, as measured in the partial waves, is small and vanishes at high energies.</footnote>
work reasonably well for scalar, $\phi^4$-like theories, but not for confining gauge theories for which $\rho_{\text{MP}}(s, t) \neq 0$ should be very important.

The plan of the paper is as follows:

* In section 2 we introduce the iteration algorithm in detail. We also discuss its numerical implementation. The key ingredient in our work is the Mandelstam representation of the amplitude that allows us to combine unitarity and crossing in a elegant way.

* In section 3, we introduce a simplified version of the iteration algorithm where we construct quasi-elastic amplitudes with zero double discontinuity, $\rho(s, t) = 0$ and elastic S-wave scattering $\eta_{\text{MP}}(s) = 1 - |S_0(s)|^2 = 0$, in spacetime dimensions $3 \le d \le 4$. This section is mostly self-contained and can provide an accessible, less technical way to delve into the formalism.

* In sections 4 and 5 we numerically construct examples of what we call quasi-elastic amplitudes. These are amplitudes which exhibit non-zero, but small particle production. We construct both 2QE and 2PR amplitudes in $d=3$ and $d=4$ in a finite range of $\lambda$ around 0. We explicitly check that the amplitudes satisfy both elastic and inelastic unitarity. We analyze their behavior at high energies both at fixed angles and in the Regge limit.

* In section 6 we derive the double discontinuity of the acnode graph in four dimensions. Following the work of Gribov and Dyatlov [48], we analytically continue the three-particle phase space integral from negative to positive $t$. We then proceed to use the acnode double spectral density as a proxy for $\rho_{\text{MP}}(s, t) \neq 0$, and we construct the corresponding amplitude functions.

* In section 7 we look at some low energy observables, to make comparison with the parameter space analyzed recently in [45, 46]. We find that the quasi-elastic amplitudes constructed in this paper all lie close to the boundary of the allowed region. We find that the 2QE amplitude is the most “extremal.”

* In section 8 we consider the question whether the Mandelstam representation holds in gravity. In this context, we revisit the old argument by Cerulus and Martin [49], which connects the Regge limit, the fixed angle scattering and the polynomial boundedness required for the validity of the Mandelstam representation.

* In section 9 we conclude and collect thoughts on open problems and future directions.

## 2 Algorithm

In this section we first review the key concepts and equations that serve as our starting point. We then present details of the numerical implementation of the iteration procedure.
## 2.1 Basic concepts

**Subtracted Mandelstam representation.** We parameterize our amplitude using the double-dispersive representation called Mandelstam representation [10], introduced in (1.1). This rewriting of the Mandelstam representation involves one subtraction. Subtractions are used to write down dispersion relations when the arc of the Cauchy contour at infinity does not vanish. For instance, let $f(z)$ be a function analytic in the complex plane apart from a branch cut on the real axis starting from $z = z_1$ and extending to infinity. If $f(z)$ goes to a constant $a$ at infinity, the usual dispersion integral that would give $f(z) = \frac{1}{2i\pi} \int_{z_1}^{\infty} \frac{\mathrm{Disc}f(w)}{z-w} dw$ cannot be written because the arc of the Cauchy contour around infinity does not vanish. A dispersion relation can only be written for the subtracted function $f(z) - a$. If the function has a polynomial growth $f(z) \sim z^n$, $n + 1$ subtractions need to be used. The same result can be achieved by dividing instead of subtracting, and obtaining dispersion relations for $\frac{f(z)}{z^{n+1}}$, in terms of $n + 1$ undetermined constants.

The Mandelstam representation is a double dispersive representation: we consider the dispersion relation for the discontinuity of the function and then we plug it into the usual dispersion relation. Similarly, double dispersion relations can be written.

The subtractions constants are unknowns, which are not determined by the dispersion relation. For us, eq. (1.1) has one unknown, $\lambda$, the value of the amplitude at the crossing symmetric point,

$$ \lambda = T\left(\frac{4m^2}{3}, \frac{4m^2}{3}\right). \qquad (2.1) $$

which is usually taken to represent the quartic coupling of the amplitude. We can turn this to our advantage, as it allows us to input the nonperturbative value of the coupling for any amplitude, thus we can explore in a controlled manner the space of couplings.

Furthermore, non-zero $\lambda$ acts as a source for the iterations, and creates amplitudes which are necessarily non-trivial.

Note finally that (1.1) allows to describe a class of functions that go at worst to a constant at infinity, but we do not put in the value of this constant. The behavior at infinity, $T(s,t) \to 0$ or $T(s,t) \to \text{const} \neq 0$ is dynamically generated by the algorithm.

**Mandelstam equation.** The Mandelstam equation is the analytic continuation in $t$ of unitarity in the elastic strip in $s$. It relates the double discontinuity of the amplitude to the square of the single discontinuity of the amplitude. Equivalently, the Mandelstam equation expresses elastic unitarity for complex spins $J$. Schematically, it is obtained by taking a second discontinuity in the standard $\mathrm{Im}T \propto |T|^2$ formulation of elastic unitarity. It is convenient to introduce the following notation

$$ \mathcal{T}(s, z) \equiv T(s, t(z)), \quad t(z) = -\frac{s - 4m^2}{2}(1 - z). \qquad (2.2) $$

For a full modern account of the Mandelstam equation and elastic unitarity, see [38]. For $s$ in the elastic strip, and $4m^2 < t < \infty$, we have

$$ \rho(s,t) = \frac{(s - 4m^2)^{\frac{d-3}{2}}}{4\pi^2(4\pi)^{d-2}\sqrt{s}} \int_{z_1}^{\infty} d\eta' \int_{z_1}^{\infty} d\eta'' \mathcal{T}_t^{(+)}(s, \eta') \mathcal{T}_t^{(-)}(s, \eta'') \tilde{K}_d(z, \eta', \eta'') \qquad (2.3) $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Domain of support of the double spectral density $\rho(s, t)$ in a theory with $Z_2$ symmetry ($\phi^4$-type interaction) with $s_{MP} = 16m^2$. The Landau curves which separate the region $\rho(s, t) = 0$ from $\rho(s, t) \neq 0$ are given in (2.5).</figcaption></figure>

where $\eta', \eta''$ are cosines of the complexified scattering angles, $\mathcal{T}_t$ is the $t$ channel discontinuity of $\mathcal{T}$ and $z_1 = 1 + \frac{4m^2}{s-4m^2}$. The upper index $\pm$ specifies the direction of real $s$ from the complex plane, namely $\pm i\epsilon$. The double discontinuity of the amplitude is only non-zero above the so-called Landau curves and vanishes below them:

$$s \le \sigma(t) : \quad \rho(s, t) = 0 \qquad (2.4)$$

where $\sigma(t)$ is a function given by the union of the leading Landau curves in two channels, see solid-orange curve in figure 2,

$$\sigma(t) = \min \left( \frac{4t}{t - 16m^2}, \frac{16t}{t - 4m^2} \right) \qquad (2.5)$$

The kernel $\tilde{K}_d$ is of purely kinematic origin and is given by, in $d > 3$ dimensions:

$$\tilde{K}_d(z, \eta', \eta'') = \frac{4\pi^{\frac{d+1}{2}}}{\Gamma(\frac{d-3}{2})} \Theta(z - \eta_+) \frac{(z^2 - 1)^{\frac{4-d}{2}}}{(z - \eta_-)^{\frac{5-d}{2}} (z - \eta_+)^{\frac{5-d}{2}}}. \qquad (2.6)$$

where $\eta_{\pm} = \eta'\eta'' \pm \sqrt{(\eta'^2 - 1)(\eta''^2 - 1)}$. In $d > 3$, the kernel reduces the double integral $\int_{z_1}^{\infty} dz' \int_{z_1}^{\infty} d\eta''$ to a domain contained below a hyperbola-shaped curve define by $\eta', \eta'' > z_1$ and $\eta_+ < z$, see figure 3. In $d = 3$ we get

$$\tilde{K}_3(z, \eta', \eta'') = 4\pi^2 \delta(z - \eta_+) \frac{\sqrt{z^2 - 1}}{\eta_+ - \eta_-}, \qquad (2.7)$$

where the delta-function localizes the 2d integral to a one-dimensional segment.
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Typical integration domain of the Mandelstam equation (2.3). The integration range for the complexified angles is $\infty > \eta', \eta'' \ge 1 + \frac{4m^2}{s-4m^2}$ and $\eta_+ \le z$, where recall that $\eta_+ = \eta'\eta'' + \sqrt{(\eta'^2-1)(\eta''^2-1)}$.</figcaption></figure>

In terms of $\rho(s)$ and $\rho(s, t)$, the t-channel discontinuity $T_t(s, t)$ is given by

$$T_t(s, t) = \rho(t) + \int_{4m^2}^{\infty} \frac{ds' (s-s_0)\rho(s', t)}{\pi (s'-s)(s'-s_0)} + \int_{4m^2}^{\infty} \frac{du' (u-u_0)\rho(u', t)}{\pi (u'-u)(u'-u_0)}. \quad (2.8)$$

Equation (2.3) is defined for $4m^2 \le s < s_{MP}$. Let us now define a function $\rho_{el}(s, t)$ by the r.h.s. of eq. (2.3) for all values of $s \ge 4m^2$. Considering a purely multi-particle function $\rho_{MP}(s, t)$ supported in the quadrant $s, t \ge s_{MP}$, we have now defined with precision the model of the amplitude written in (1.9), which is the model that we use. For us, as we said, $\rho_{MP}(s, t)$ will be treated as an input, and we will typically define it above some realistic multi-particle Landau curves, following the analysis of [34].

It is important to note that no assumption was made in writing (1.9): any double spectral density of an amplitude that satisfies elastic unitarity can be written in the form (1.9). The key aspect of this formula is the fact that adding the double spectral density needed to satisfy elastic unitarity in the t-channel, namely $\rho_{el}(t, s)$, restores crossing symmetry and at the same time does not spoil elastic unitarity in the s-channel. This happens due to the fact that $\rho_{el}(t, s)$ has only support for $s \ge s_{MP}$. In this way we can have both unitarity and crossing!

**J = 0 unitarity.** The Mandelstam equation is a dispersive version of elastic unitarity. Indeed, by integrating both sides of (2.3) against the Legendre Q-function and using the Froissart-Gribov formula for partial waves, see appendix A and [38], one recovers $|S_J| = 1$ for complex $J$ as long as $\text{Re}[J] > J_0$, where $J_0$ is the Regge intercept of the amplitude. The last condition is required for validity of the Froissart-Gribov representation of partial waves. In this paper, we consider amplitudes with $J_0 = 0$. Therefore, the Mandelstam equation does not guarantee elastic unitarity of $J = 0$ partial wave and we need an extra unitarity constraint to impose it. It was described in the introduction and is much simpler
to state than the Mandelstam equation. We simply enforce the $J = 0$ partial wave unitarity equation (1.10) as follows:

$$2\mathrm{Im}f_0(s) = \frac{(s - 4m^2)^{\frac{d-3}{2}}}{\sqrt{s}} |f_0(s)|^2 + \frac{\sqrt{s}}{(s - 4m^2)^{\frac{d-3}{2}}} \eta_{\mathrm{MP}}(s), \quad (2.9)$$

where $\eta_{\mathrm{MP}}(s) \neq 0$ only in the multi-particle region $s \ge s_{\mathrm{MP}}$. Now, to be practically implemented in terms of the data of our amplitude, and in particular $\rho(s)$ and $\lambda$ (but also $\rho(s, t)$), we need to perform the partial wave projection to $J = 0$ of the Mandelstam representation eq. (1.1). We display the lengthy expression of the result of the projection in (A.7) in general dimension $d$.

To unpack (2.9) a bit more, and prepare it for the iterations explained below, we keep $|f_0(s)|^2$ in the right-hand side implicit, but spell out the left-hand side using (A.7) in terms of $\rho(s)$, $\rho(s, t)$. We thus obtain the following intermediate step equation, in spacetime dimension $d$:

$$\rho(s) + 2 \int_{4m^2}^{\infty} \frac{dt'}{\pi} \rho(s, t') \left( \frac{1}{t_0 - t'} + \frac{1}{t'} {}_2F_1 \left( 1, \frac{d-2}{2}, d-2, \frac{s-4m^2}{t'} \right) \right) \\ = 4(16\pi)^{\frac{d-3}{2}} \Gamma\left(\frac{d-1}{2}\right) \left( \frac{(s-4m^2)^{1/2}}{\sqrt{s}} |f_0(s)|^2 + \frac{\sqrt{s}}{(s-4m^2)^{\frac{d-3}{2}}} \eta_{\mathrm{MP}}(s) \right). \quad (2.10)$$

Here, just like in (A.7), in integer dimensions $d = 3, 4, 5, \dots$, the hypergeometric function simplifies to functions with logarithmic or square roots singularities, depending on whether $d$ is even or odd, respectively.<fnref n="14" /> Note finally that although $\lambda$ is not visible in this is equation, it is explicitly entering $f_0(s)$ in the r.h.s., see (A.7).

## 2.2 Iteration algorithm

We are now well equipped to explain the iteration algorithm. As we already mentioned above, the basic idea behind it goes back to the S-matrix program pursued in the 60's [9-11, 13-17].

The idea is to input fixed values of the triad $(\lambda, \eta_{\mathrm{MP}}(s), \rho_{\mathrm{MP}}(s, t))$, and construct the amplitude $T(s, t)$ that satisfies analyticity, crossing and elastic unitarity, that is to say obtain the functions $\rho(s)$ and $\rho(s, t)$ which satisfy the equations (2.3) and (2.10). The algorithm does not guarantee inelastic unitarity (1.7) as an output. It has to be checked separately on the fixed point of the iteration algorithm. In this sense, inelastic unitarity constrains a *posteriori* the space of inputs $(\lambda, \eta_{\mathrm{MP}}(s), \rho_{\mathrm{MP}}(s, t))$.

The unitarity equations can be formally rewritten in terms of a functional $\Phi$, whose elements are already defined above, in the following form:

$$(\rho(s), \rho(s, t)) = \Phi_{\mathrm{MP},\lambda}[\rho(s), \rho(s, t)] \quad (2.11)$$

where we specified the input triad elements as an index of $\Phi$ to insist on the fact that these are external, fixed parameters.

<footnote n="14">¹⁴The $_2F_1$ function simply comes from performing the partial wave projection of the t- and u-channel cuts in the Mandelstam representation.</footnote>
We solve these equations by iterating the map $\Phi$: initialising on a starting point $\rho^{(0)}(s), \rho^{(0)}(s,t)$, we define a sequence of functions $\rho^{(n)}(s), \rho^{(n)}(s,t)$, $n = 1, 2, \dots$ by

$$ (\rho^{(n+1)}(s,t), \rho^{(n+1)}(s)) = \Phi_{\text{MP},\lambda} [\rho^{(n)}(s), \rho^{(n)}(s,t)] \quad (2.12) $$

If the map converges, the fixed point, $(\rho^{(\infty)}(s,t), \rho^{(\infty)}(s)) = \Phi_{\text{MP},\lambda}[\rho^{(\infty)}(s), \rho^{(\infty)}(s,t)]$ gives by definition a scattering amplitude which satisfies all of the equations above, therefore satisfies analyticity, crossing and elastic unitarity (with inelastic unitarity checked separately).

**Iteration steps.** The iteration algorithm proceeds in the following steps, in order:

1. Single discontinuity iteration, from (2.10),

$$ \begin{aligned} \rho^{(n+1)}(s) = & 4(16\pi)^{\frac{d-3}{2}} \Gamma\left(\frac{d-1}{2}\right) \left( \frac{(s-4m^2)^{1/2}}{\sqrt{s}} |f_0^{(n)}(s)|^2 + \frac{\sqrt{s}}{(s-4m^2)^{\frac{d-3}{2}}} \eta_{\text{MP}}(s) \right) \\ & - 2 \int_{4m^2}^{\infty} \frac{dt'}{\pi} \rho^{(n)}(s,t') \left( \frac{1}{t_0-t'} + \frac{1}{t'} \right) {}_2F_1 \left( 1, \frac{d-2}{2}, d-2, \frac{s-4m^2}{t'} \right), \end{aligned} \quad (2.13) $$

where $f_0^{(n)}(s)$ is obtained from formula (A.7), adding $(n)$ iteration counters exponents to the spectral $\rho$ functions.

2. Double discontinuity iteration

$$ \rho_{\text{el}}^{(n+1)}(s,t) = M_d[T_t^{(n)}, T_t^{(n)}](s,t), \quad (2.14) $$

$$ \rho^{(n+1)}(s,t) = \rho_{\text{el}}^{(n+1)}(s,t) + \rho_{\text{el}}^{(n+1)}(t,s) + \rho_{\text{MP}}(s,t). \quad (2.15) $$

where $M_d[\cdot, \cdot]$ is the quadratic functional of the Mandelstam equation, which is exactly the r.h.s. of eq. (2.3), with $T_t^{\pm}$ replace by $T_t^{\pm,(n)}$, so we do not reproduce the expression here.

We simply rewrite the $t$-channel single disc, to clarify that we do not iterate any function with that step, but simply reconstruct the single disc. from the double disc.:

$$ T_t^{(n)}(s,t) = \rho^{(n)}(t) + \int_{4m^2}^{\infty} \frac{ds' (s-s_0) \rho^{(n)}(s',t)}{(s'-s)(s'-s_0)} + \int_{4m^2}^{\infty} \frac{du' (u-u_0) \rho^{(n)}(u',t)}{(u'-u)(u'-u_0)}. \quad (2.16) $$

In practice, for the initialization step, we set $\rho^{(0)}(s) = \rho^{(0)}(s,t) = T_t^{(0)}(s,t) = 0$. The condition $T^{(0)}(s,t) = \lambda$ is automatically generated by the subtracted representation we use, and it enters $|f_0|^2$ in the r.h.s. of (2.13) at step 1.

Let us now imagine a situation where this iterative process converges and we reach the fixed point. The physical content of (2.13) at the fixed point is then nothing but unitarity for the $J = 0$ partial wave (2.9). Elastic unitarity for partial waves with $J > 0$ is satisfied at the fixed point due to (2.14). Crossing symmetry is manifest due to (2.15), thanks to adding $\rho_{\text{el}}^{(n+1)}(t,s)$ which restores crossing, and at the same time does not spoil elastic unitarity since it has only support for $s \ge s_{\text{MP}}$.

Lastly, inelastic unitarity for $J > 0$ partial waves (1.7) is not guaranteed to hold at the fixed point of the iteration. It presents a nontrivial constraint on the triad $(\lambda, \eta_{\text{MP}}(s), \rho_{\text{MP}}(s,t))$ which we check below.
## 2.3 Numerical implementation

Let us now describe how we numerically implemented this procedure.

### 2.3.1 Discretization

The numerical implementation is rather simple in essence. We discretize the functions $\rho(s)$, $\rho(s, t)$, $T_t(s, t)$ on grids and define linear interpolants as the functions we iterate. To perform an iteration, we evaluate the updated functions of step $n + 1$ on the grid points, and obtain the iterated interpolants in this way.

For the discretization, we need three grids:

* a one-dimensional grid for $\rho(s)$, that spans $s \ge 4m^2$,

* a two-dimensional grid for $T_t(s, t)$ that spans the whole quadrant $s, t \ge 4m^2$ (the single discontinuity of the amplitude is non zero after the first threshold, this is standard),

* a two-dimensional grid for $\rho(s, t)$ that spans the portion of the quadrant $s, t \ge 4m^2$ which is above the Landau curves.

We start by mapping these domains to $[0, 1]$ and the unit square $[0, 1]^2$ via the following change of variables:

$$ (s, t) \rightarrow (x, y) = \left( \frac{4m^2}{s}, \frac{4m^2}{t} \right) \qquad (2.17) $$

By doing so we effectively introduce a UV and IR cutoffs related to the sampling of the functions close to 0 and 1, which we discuss below. We provide more details on the explicit construction of these grids in appendix B.

**Grid for $\rho(s, t)$**. As was said above, $\rho(s, t)$ is non-zero only above the union of the two Landau curves ($\frac{16s}{s-4}$, $\frac{4s}{s-16}$). In the $x, y$ variables, these curves become lines, defined by

$$ \left( \frac{1-x}{4}, 1-\frac{x}{4} \right), \qquad (2.18) $$

and, as a function of $x, y$, the double spectral function $\rho(x, y)$ has support below these lines, as is depicted in figure 4.

The various high energy limits of the amplitude are now described as follows: the Regge limits correspond to the segments $x = 0, y = 0$. The fixed (imaginary)$^{15}$ angle region corresponds to $x, y \rightarrow 0$ with $x/y$ fixed. The threshold region, which matters in $d = 3$ but is trivial in $d = 4$, is mapped to $x \rightarrow 1$ and $y \rightarrow 1$. The grids have to sample accurately these various regions and for this reason, we choose grids with exponentially many points up to some cut-off.

After experimenting a lot with various types of grids, we settled for a default “fine” $(s, t)$-grid with 11445 points, with cutoff at $x = 10^{-12}$, which provided a good interplay between the need to describe inelastic effects to very high energies while keeping manageable computation times. We give more detail on this grid in appendix.

<footnote>$^{15}$Recall that we are in the unphysical kinematics quadrant $s, t \ge 0$.</footnote>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** The region of the nontrivial support of the double spectral density in the $x, y$ variables. It is given in terms of the location of the leading Landau curves (2.4).</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 5.** Linear interpolation of a slowly decaying function, $-1/\log(x)$ in this example. The first grid point is at $10^{-9}$, and then the interpolating function sharply connects to zero. The dispersive integrals used in the iteration process have $1/s$ kernels which suppresses the integration at the cutoff scale and renders this effect innocuous, if the cutoff is high enough. Note also that the garland-shape of the piecewise-linear interpolating function is due to the log-log scale.</figcaption></figure>

**A note about the cutoff.** For the class of amplitudes considered in this paper, we always have $\rho(x), \rho(x, y) \to 0$ at high energies. The cutoff mentioned above is to be understood in the following sense. Above the cut-off, even though we use linear interpolation, the functions are dynamically generated by the iteration process. Between the cut-off and zero, we continue to use linear interpolation in $x, y$, which means, in the $s, t$ variables, that we force a $1/s$ or $1/t$ decay of the functions. An illustration of this phenomenon is given in figure 5.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Fine mesh for the double spectral density (11445 points). Red dots are points $(x_i, y_i)$ on which the functions is evaluated, while the interpolation is linear on each facet.</figcaption></figure>

One last important detail about the 2d grid is that when we compute $\rho_{\text{el}}^{(n)}(x_i, y_j)$ in step 2 of the iteration, it is only non-vanishing below the Landau curve $y = (1-x)/4$, thus we only need to discretize the wedge between $y=0$ and $y=(1-x)/4$, for $0 < x < 1$. The full double-discontinuity is obtained after crossing-symmetrizing, $\rho(x, y) = \rho_{\text{el}}(x, y) + \rho_{\text{MP}}(x, y)$ and is correctly defined below the union of both Landau curves, as it should.

The fine grid with which we generated most of the data is represented in figure 6. Later in section 6 we comment on the influence of density of the grids on the precision of our result. This grid has essentially has a subleading influence on the significant digits of the low-energy observables described in section 7. However, this grid is essential to produce high energy quantities, such that inelasticity in the partial waves, accurately.

In general, all of the errors we face in this paper are errors coming from the “trapezoidal rule”, familiar from approximating one-dimensional integrals with piecewise linear functions.

**Grid for $T_t(s, t)$**. The single discontinuity $T_t$ is defined on the whole quadrant $[0, 1]^2$. Therefore, we need another grid to sample this domain. For the latter we simply chose a uniformly-spaced grid of points at $(x_i, y_j) = (i/N, j/N)$ with $i, j = 0, 1, \dots, N$. In practice, this appeared sufficient. We used a regular square grid with $72 \times 72 = 5184$ points. This grid has essentially no influence low-energy observables nor high energies. This comes from the fact that $T_t(s, t)$ is dominated by $\rho(s)$, which has its own grid.

**Grid for $\rho(s)$**. Finally, the $t$-independent simple spectral function $\rho(s)$ is discretized on a 1-dimensional grid which also samples logarithmically the extremities $x \to 0, 1$ so as to have more information on threshold and Regge behavior. We also took a cut-off $10^{-12}$ in order to have consistency between our grids and 91 points in total. The explicit definition of the grid is given by the following sequence of points:

$$x_0 = 0, \quad x_1 = 10^{-12}, \quad \dots, \quad x_N = 1. \qquad (2.19)$$

In section 7, we produce some estimation of the convergence of low-energy observables as we densify this grid, and find it to be the major source of error. Again, this error appears to be of the simple trapezoidal-rule type, as it scales with the square of the inverse of the number of points, $1/N_{\text{grid-points}}^2$.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** An example of the integration domain for the Mandelstam equation (2.3) in the Regge limit for $s = 5m^2, t = 10^6 m^2$.</figcaption></figure>

### 2.3.2 Numerical integrals

Once we have defined the interpolating functions, we can run the algorithm. At step (n+1), we compute numerically all the right-hand sides in the equations defined in section 2.2. These numerical integrals need to be handled with care. One important point relates to all the dispersion integrals. The real part of these is computed using the standard principal value prescription (the imaginary part is just a delta function and poses no problem). These principal values can be improved in a standard way, using the following identity

$$ \text{P.V.} \int_0^1 \frac{f(x')}{x' - x} dx' = \int_0^1 \frac{f(x') - f(x)}{x' - x} dx' + f(x) \log \frac{1-x}{x}. \qquad (2.20) $$

and variations thereof when the singular point lies at the extremity and for two-dimensional integrals.

Another source of numerical intricacies come from the deformation of the integration domain in the r.h.s. of the Mandelstam equation in the Regge limit. In a very deep Regge kinematics, the two-dimensional hyperbola-shaped domain of integration of the Mandelstam integral becomes very skewed and requires care in the numerical evaluation, see figure 7. For these integrals, we have found that we obtained reliable results using the standard **GlobalAdaptive** method of Mathematica, supplied with various **PrecisionGoal**, which we would in general keep to 5 to keep good runtimes and control on the numerics. We did not give much importance to the estimated error of the numerical integrations, because we have a precise and independent way to check our results, which is to measure violation of the imposed unitarity condition, and we have found that it is satisfied to a very good precision, e.g. $1 - |S_0| \sim 10^{-10}$ for the 2QE amplitudes.

We plan to release the code in the near future,<fnref n="16" /> once we have optimized it fully. In the meantime, we make available with the arXiv submission a notebook with a few of the most representative examples of our amplitudes.

<footnote n="16">¹⁶Definitely before the earth reaches a warming of +1°, which is unfortunately bound to happen soon.</footnote>
<table><caption>**Table 1.** Approximate range of convergence of the algorithm. “Back-reacted” refers to the amplitudes with non-zero double spectral density constructed in section 4 and section 5.</caption><thead><tr><th></th><th>Model</th><th>λ<sub>min</sub></th><th>λ<sub>max</sub></th></tr></thead><tbody><tr><td rowspan="2">d = 3</td><td>toy-model</td><td>-5π</td><td>0</td></tr><tr><td>back-reacted</td><td>-3π</td><td>0</td></tr><tr><td rowspan="2">d = 4</td><td>toy-model</td><td>0</td><td>25π</td></tr><tr><td>back-reacted</td><td>0</td><td>5π</td></tr></tbody></table>

### 2.3.3 Convergence

In this program we can discuss two types of convergence. The first one refers to the convergence of the fixed-point iteration with respect to the number of iterations, whereas the second one is the convergence of the discretized solution to the continuum limit solution as we take the grid size and cut-offs to zero.

**Convergence of the discretized algorithm.** Fixed-points can be attractive or repulsive. Our iterations converge in some range of the input parameter space, which depends on the size of the inputs $(\lambda, \eta_{MP}(s), \rho_{MP}(s, t))$. To first approximation, we find that the iteration algorithm converges when the input is small enough. Physically, it corresponds to interactions which are not too strong. In $d = 3$, we found convergence for $-3\pi \le \lambda \le 0$, and in $d = 4$ for $0 \le \lambda \le 5\pi$. Convergence ranges are recapitulated in table 1, for the toy-model of section 3 and the full algorithm of sections 4, 5.

Let us next talk about the speed of convergence. A sequence $(f)_{n,n\ge0}$ is said to have rate of convergence $q$ if there exists a non-zero constant $\mu$ such that

$$ \lim_{n \to \infty} \frac{\|f_{n+1} - f_*\|}{\|f_n - f_*\|^q} = \mu. \qquad (2.21) $$

where we stay generic about the norm and the space in which $f_n$ lives. If $q = 1$, the convergence is called *linear*, this is the case for fixed-point iterations in general, and what we observe here, see figure 10. Gradient methods such as Newton-Raphson have *quadratic* convergence, which we observed explicitly in $d = 2$ [26], where we implemented both methods. As can be seen from (2.21), *linear convergence* implies that the approach to the fixed point is *exponentially fast*.

**Convergence in the continuum limit.** Let us next discuss convergence of the algorithm to the full, smooth, solution to the equations, as the grid size goes to zero.

First of all, let us notice that no continuum limit proofs for the existence of solutions to the equations described in the present paper exist up-to-date. In [16], Atkinson considered a $d = 4$ amplitude given by the Mandelstam representation with one subtraction, but when imposing unitarity he used a modified set of equations which effectively introduces a UV cutoff that trivializes the Regge limit. We describe this phenomenon in a simplified setting in appendix E. When trying to remove this cutoff, the parameter space in which convergence of the algorithm can be proven shrinks to zero size, unless one restricts to $3 < d < 4$. Note
finally that in the original paper, Atkinson works with no subtractions and in this case no such cut-off is needed, [14, 15], but this setup does not allow to describe $\phi^4$-like interactions, which do not decay in the Regge limit in $d=4$ and in all kinematics in $d=3$.

As a consequence, the only way for us to describe non-discretized solutions in $d=3$ and $d=4$ is empirical: we generate solutions on different grids of increasingly smaller minimum spacing and observe convergence to a smooth solution.

Next, when taking the continuum limit, we observed two qualitatively different behaviors: either the solutions were simply becoming more and more smooth and valid in a larger range of energies, or we were losing the solutions. We observed that the continuum limit is smooth as we are removing all cutoffs for $\lambda < 0$ in $d=3$, and $\lambda > 0$ in $d=4$. For the opposite signs of couplings, positive in $d=3$ and negative in $d=4$, we observed that the algorithm would converge on a finite grid with small number of points, and start to diverge as we took the cutoff to infinity and grid spacing to zero. In $d=3$, the continuum limit solution in the class of functions that we considered requires $\lambda < 0$ via the simple dispersive argument, see (3.17). In $d=4$, we believe that the divergence of the algorithm is related to the existence of the Landau pole, which causes the theory to become strongly coupled as we send the UV cutoff to infinity. For us, it has manifested itself through the fact that the amplitude becoming bigger in the UV eventually caused a divergence of the algorithm.

In conclusion, we find that we can define a continuous limit while maintaining the contracting property of the algorithm.

### 2.3.4 Computation times

To close this section, let us just comment briefly on the computation times. We chose our typical grids for $\rho(s, t)$ 11,455 points and $\rho(s)$ with 100 points, respectively, such that computation times were manageable. As a the data plotted in this paper could be ran on a modern laptop in less a week. Furthermore, the codes are easily parallelized, since the iteration proceeds by evaluating integrals of the same functions $\rho^{(n)}(s)$, $\rho^{(n)}(s, t)$ at all the grid points $(s_i, t_j)$ and are completely independent of each other. To produce a typical dataset, with 10 different couplings, and 10 iterations (with hotstart to achieve good precision), can be done in $O(10h)$ with 50 CPUs, using mathematica's native lightweight grid-computing environment. Given how easy it is to parallelize this code, which essentially performs thousands of independent integrals at each step, we expect a huge speed-up if it were adapted to run on a GPU.

## 2.4 Fixing $\eta_{MP}$: analyticity in spin up to $J=0$ and Feynman graphs

In the discussion above we considered the S-wave inelasticity $\eta_{MP}(s)$ and $\rho_{MP}(s, t)$ to be independent inputs. However, there is a simple situation, where the two are not independent: consider an amplitude for which discontinuity $T_s(s, t)$ vanishes in the Regge limit $t \to \infty$ for any $s$. In this case we can use the Froissart-Gribov formula for the imaginary part of the partial wave

$$ \mathrm{Im} f_J(s) = \frac{2N_d}{\pi} \int_{z_1}^{\infty} dz (z^2 - 1)^{\frac{d-4}{2}} Q_J^{(d)}(z) \rho(s, t(z)). \qquad (2.22) $$
all the way to $J = 0$. Recall that we can always decompose the double spectral density as follows

$$ \rho(s, t) = \rho_{\text{el}}(s, t) + \rho_{\text{el}}(t, s) + \rho_{\text{MP}}(s, t). \qquad (2.23) $$

Plugging this decomposition into the Froissart-Gribov formula (2.22), we get the desired relationship between $\eta_{\text{MP}}(s)$ and $\rho_{\text{MP}}(s, t)$:

Analyticity in spin up to $J = 0$:

$$ \eta_{\text{MP}}(s) = \frac{(s - 4m^2)^{\frac{d-3}{2}}}{\sqrt{s}} \frac{4N_d}{\pi} \int_{z_1}^{\infty} dz (z^2 - 1)^{\frac{d-4}{2}} Q_0^{(d)}(z) (\rho_{\text{el}}(t(z), s) + \rho_{\text{MP}}(s, t(z))), \qquad (2.24) $$

where we also used the fact that $\rho_{\text{el}}(s, t)$ satisfies the Mandelstam equation and therefore is responsible for elastic scattering. In the numerical iterations, we use the single spectral function $\rho(s)$, so to implement (2.24) it is convenient to rewrite it as follows

$$ \rho_{\text{MP}}(s) = \frac{4\Gamma(\frac{d-1}{2})}{\pi^{3/2}\Gamma(\frac{d-2}{2})} \int_{z_1}^{\infty} dz (z^2 - 1)^{\frac{d-4}{2}} Q_0^{(d)}(z) (\rho_{\text{el}}(t(z), s) + \rho_{\text{MP}}(s, t(z))), \qquad (2.25) $$

where $\rho_{\text{MP}}(s)$ is the multi-particle single spectral density.

We call the algorithm in which we fix the S-wave particle production using (2.25) two-particle reducible (2PR). The reason is that it has simple relationship to Feynman diagrams in $\phi^4$ theory, which we discuss now. Let us start by setting $\rho_{\text{MP}}(s, t) = 0$: this corresponds to throwing away all the graphs which are not two-particle irreducible, e.g. the open envelope graph in figure 1. We can then generate the amplitude by iterations which is graphically summarized in figure 8. We also set $\eta_{\text{MP}}(s)$ to be given by (2.24). We expect that this procedure corresponds to re-summing all the two-particle reducible Feynman graphs generated by two-particle unitarity and crossing. In figure 1 these are graphs depicted in black. This expectation is based on the assumption that analyticity in spin of inelasticity up to $J = 0$ is a true property of this subset of Feynman graphs, something that we verified explicitly up to three loops.

In the full theory we do not expect (2.24) to hold. The reason is that we get diagrams similar to the frog diagram, see figure 1. On one hand, this diagrams contributes to $\eta_{\text{MP}}(s)$. On the other hand, it does not have any double spectral density. Interestingly, the frog diagram is just a first in a series of planar melonic diagrams of this type.

Such diagrams have been recently studied in the context of the SYK model, see e.g. [50, 51], and it would be very interesting to explore them in the context of multi-particle scattering. Similarly, understanding unitarity structure of these graphs nonperturbatively is an interesting problem.

## 2.5 Graph interpretation and summability

After having exposed this whole procedure, the reader might legitimately ask: why would this iteration procedure converge? After all, iterating unitarity is very similar in essence to the standard perturbative computations, and perturbation theory is known to produce an asymptotic series due to a factorial growth of the total number of Feynman graphs.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Graphical representation of the iteration process for the double spectral density. The first three graphs in the r.h.s. are generated by the Mandelstam double integral and the Mandelstam equation. Graphs in the second row are generated through the addition of the crossed term $\rho_{el}(t, s)$ in the iteration process for $\rho(s, t)$. Finally, the last term specifies the addition of production multi-particle term $\rho_{MP}(s, t)$ at every iteration.</figcaption></figure>

Two simple arguments indicate that our iteration scheme should be convergent. The iteration can be schematically presented as a simple graphical recursion relation of figure 8. The first element in favor of convergence is that this recursion only generates planar graphs, and the planar sector is known to be convergent in some cases, see e.g. [52]). Secondly, a more precise counting can be made. The number of terms $N_{n+1}$ at the $n+1$-st step can be seen to be roughly given in terms of the number of terms $N_n$ at the $n$-th step by

$$N_{n+1} = 6(N_n)^2 + 1. \qquad (2.26)$$

Compared to Feynman graphs at step $n$, we generate graphs with $N_{n-1}^2$ loops, as opposed to $N_{n-1} + 1$ loops which we would have in perturbation theory by adding a loop. Therefore, the relation between $n$ and the number of loops $L$ is not linear, but rather exponential:

$$L \approx 2^n. \qquad (2.27)$$

From (2.26), we can estimate that $N_n \simeq e^{\log c_0 L} \simeq c_0^L$ where $c_0$ depends on the initialization. This does not exactly tell us the number of $L$-loop graphs yet it motivates that the growth should be a power-law, rather than a factorial which implies a finite radius of convergence of the algorithm. A slightly more refined expectation is that we have a series of the type $\sum_L (\lambda_p c_0)^L$, where $\lambda_p$ is the strength of scattering at energy scale $p$. For the series to converge therefore we expect that the theory should stay weakly coupled at all scales. All the amplitudes that we construct in this paper are indeed of this type.<fnref n="17" />

A closely related problem was solved in [54], and later in [55] using matrix model techniques. The generating function for the number of two-particle fully two-particle reducible planar graphs was found to be given by the following generating function

$$\Gamma(g) = \frac{1}{2} \left( 1 - g - \sqrt{1 - 6g + g^2} \right) = g + 2g^2 + 6g^3 + 22g^4 + \dots \qquad (2.28)$$

<footnote n="17">¹⁷As a side remark, let us notice that our split of Feynman diagrams into the multi-particle graphs which serve as a building block for loops of two-particle reducible graphs is reminiscent of the “loops of loops” and “hard” vs. “easy” graphs in the discussion in [53].</footnote>
<table><caption>**Table 2.** The number of non-automorphic graphs contributing to $2 \to 2$ scattering in the $\phi^4$ theory as a function of number of vertices. We distinguish two-particle fully reducible graphs which are automatically generated by the iteration algorithm considered in the present paper and all graphs. Given that the number of 2PR graphs grows like a power, we see that most of the graphs (which add up to the factorial growth) contribute to the multi-particle input. Similarly, if we take a given multi-particle graph and add it to the iteration algorithm, it will generate the number of extra graphs that grow like a power. We thank Miguel Correia for producing the counting presented in this table.</caption><thead><tr><th># of vertices</th><th>1</th><th>2</th><th>3</th><th>4</th><th>5</th><th>6</th><th>7</th><th>8</th><th>9</th><th>10</th></tr></thead><tbody><tr><td>All graphs</td><td>1</td><td>1</td><td>2</td><td>8</td><td>26</td><td>124</td><td>627</td><td>3794</td><td>25306</td><td>188746</td></tr><tr><td>2PR graphs</td><td>1</td><td>1</td><td>2</td><td>6</td><td>19</td><td>74</td><td>310</td><td>1451</td><td>7130</td><td>35398</td></tr><tr><td>Planar graphs</td><td>1</td><td>1</td><td>2</td><td>7</td><td>23</td><td>95</td><td>411</td><td>2005</td><td>10214</td><td>54873</td></tr></tbody></table>

Here the coefficients describe the number of different graphs at order $g^n$. For example, $2g^2$ stands for horizontal plus vertical one-loop bubble; $6g^3$ comes from horizontal and vertical two-loop bubbles as well as four different orientations of the bubble plus triangle diagram, etc. From this expression we immediately see that the number of graphs indeed grows like a power, and the radius of convergence of the small $g$ expansion is finite $g < g_*$, and is controlled by the smallest root of $\sqrt{1 - 6g + g^2}$

$$g_* = 3 - 2\sqrt{2} \approx 0.171573. \qquad (2.29)$$

Given that the actual perturbation theory expansion parameter involves extra numerical factor, e.g. $\frac{\lambda}{16\pi^2}$ in $d=4$, we expect that the radius of convergence of our algorithm to be bigger than given by (2.29). This is indeed what we observed in practice.

It is also interesting to list the number of non-automorphic graphs both two-particle fully reducible and total number [56], see table 2. As expected the number of full graphs grows much faster very quickly only a small fraction of diagrams comes from two-particle iterations, whereas most of the graphs contribute to $\eta_{MP}(s)$ and $\rho_{MP}(s, t)$. When one accounts for automorphisms $(1, 1, 2, 6, \dots)$ from table 2 becomes $(1, 2, 6, 22, \dots)$ from (2.28). Generalizing the estimate above for the number of graphs, from the two-particle reducible to $n$-particle reducible, we still get a power-law behavior (though with a larger power). The factorial growth thus originates from the fact that $n$ is unbounded from above.

A few comments are in order. Although our summation involves only planar graphs and fewer terms than expected, we still obtain a fully nonperturbative amplitude in the end, which is generated by the inputs $(\lambda, \eta_{MP}(s), \rho_{MP}(s, t))$. These inputs parameterize the space of amplitude functions: navigating the space of inputs formally maps to navigating the space of outputs, i.e. the space of amplitudes. In particular, the actual physical amplitudes within this space. Practically, this construction allows us to formally ignore the problem of re-summation of all Feynman graphs (which should result in $\eta_{MP}(s)$ and $\rho_{MP}(s, t)$), and yet be able to describe nonperturbative physical amplitudes.
### 3 Toy-model amplitudes: zero double discontinuity

Having described the full set of equations and the algorithm in detail, we now turn to a simplified version of the problem. The amplitudes we consider are obtained by discarding the double spectral function of the Mandelstam representation, and assume the form:

$$T(s, t) = \lambda + \int_{4m^2}^{\infty} \frac{ds'}{\pi} \frac{\rho(s')}{s' - \frac{4m^2}{3}} \left( \frac{s - \frac{4}{3}m^2}{s' - s} + \frac{t - \frac{4m^2}{3}}{s' - t} + \frac{u - \frac{4m^2}{3}}{s' - u} \right). \quad (3.1)$$

Thanks to the subtraction, we still have that

$$T\left(\frac{4m^2}{3}, \frac{4m^2}{3}\right) = \lambda. \qquad (3.2)$$

This representation is manifestly crossing symmetric and has maximal analyticity. However, given zero double spectral density, unitarity can only be imposed in the $J = 0$ sector.<fnref n="18" /> It takes the form

$$|S_0(s)|^2 = 1 - \eta_{\text{MP}}(s), \qquad (3.3)$$

where $1 \ge \eta_{\text{MP}}(s) \ge 0$ for $s \ge s_{\text{MP}}$.

However, in this section, we describe in detail the case $\eta_{\text{MP}}(s) = 0$, so we impose that the scattering in the S-wave is elastic at all energies:<fnref n="19" />

$$|S_0(s)| = 1, \quad s > 4m^2. \qquad (3.4)$$

We further recall that in this paper we work with the assumption that we have no bound states below the two-particle threshold. Below, we refer to the amplitudes that satisfy (3.4) and (3.1) as toy-model amplitudes.

This toy-model serves three main purposes. Firstly, it is a simpler pedagogical setup to learn the general techniques of this paper. Secondly, we actually use this algorithm in the first step of generating our solutions, to reach a first fixed-point which is close to that of the solution to the full set of equations. We refer to this procedure as a hotstart. Lastly, these amplitudes are interesting in their own right, and they provide a sort of analogue of the CDD two-dimensional integrable S-matrices, see e.g. [57] and appendix C.

Note finally that for the problem discussed in this section, all integrals are one-dimensional. This implies that all the methods, including the gradient-based methods developed in our former work [26] can be applied to implement a more efficient iteration algorithm. The use of gradient-guiding would be expected to increase the space of parameters in which the iteration converges. The reduction of numerical integrations to tensorial operations also speeds up significantly the process. This would consequently render hotstart procedure described here very efficient. We did not attempt this here and leave this improvement for future work.

<footnote n="18">¹⁸Later we study the full solution, with nonzero double spectral density, which restores unitarity for all spins. This will be the subject of the following sections.</footnote>

<footnote n="19">¹⁹Notice that existence of such S-matrices in $d > 2$ is consistent with the Aks theorem [27], which only requires particle production for partial waves with $J \ge 2$.</footnote>
Let us now begin the detailed study of this toy-model. We split our discussion according to the number of spacetime dimensions. The situation in $d = 2$ is reviewed in appendix C. We analyze in detail the new cases $d = 3$ and $d = 4$, which turn out to be quite subtle, because of threshold and Regge behavior, respectively, and the iteration algorithm requires extra care to converge in practice.

Note finally that, as we commented on before, the situation is the simpler in $3 < d < 4$, and for this case we can even present an argument about the existence of the toy-model amplitudes directly in the continuum limit, along the lines of the original work by Atkinson [24], see appendix E. We do not attempt to solve the problem in $d > 4$.

## 3.1 $d = 3$

In three dimensions, the near two-particle threshold region requires extra care. We start the discussion by explaining the equations to be solved iteratively, derive the expected two-particle threshold behavior, and then present the results of our numerical implementation. The section is meant to be, to a large extent, self-contained.

### 3.1.1 The algorithm

Our starting point are amplitude functions that admit once-subtracted Mandelstam representation with zero double discontinuity (3.1). We start with discussing S-wave unitarity (3.4), which is the sole unitarity condition to be satisfied. In section 2, we referred to appendix A for the derivation of eq. (2.10). Here, the absence of double-discontinuity makes this derivation much less cumbersome, and to keep the discussion self-contained, let us derive the analogue of (2.10) in a fully explicit manner (with explicit $f_0$). To start, we expand the amplitude in three-dimensional partial waves$^{20}$

$$T(s,t) = 16 \sum_{J=0}^{\infty} \frac{1}{1 + \delta_{J,0}} f_J(s) \cos(J\theta), \quad t = -\frac{s-4m^2}{2}(1-\cos\theta), \qquad (3.5)$$

$$f_J(s) = \frac{1}{16\pi} \int_0^{2\pi} d\theta \cos(J\theta) T\left(s, -\frac{s-4m^2}{2}(1-\cos\theta)\right). \qquad (3.6)$$

The unitarity equation that we are solving, $|S_0(s)|^2 = 1$, takes the form

$$2\mathrm{Im}f_0(s) = \frac{1}{\sqrt{s}}|f_0(s)|^2, \quad s \ge 4m^2, \qquad (3.7)$$

where we used that, in $d = 3$,

$$S_J(s) = 1 + \frac{i}{\sqrt{s}} f_J(s). \qquad (3.8)$$

We are interested in finding an amplitude, or rather its discontinuity $\rho(s)$, that solves (3.7) iteratively. Using (3.6), we can compute the spin zero partial wave in terms of

<footnote>$^{20}$We follow the conventions of [38] which are reviewed in appendix A.</footnote>
the spectral density $\rho(s)$. We plug (3.1) into (3.6) to get

$$f_0(s) = \frac{1}{8} \left( \lambda + \int_{4m^2}^{\infty} \frac{ds'}{\pi} \frac{\rho(s')}{s' - \frac{4m^2}{3}} \left( \frac{s - \frac{4}{3}m^2}{s' - s} + K_0^{(d=3)}(s', s) \right) \right), \\ K_0^{(d=3)}(s', s) = 2 \left( \frac{s' - \frac{4m^2}{3}}{\sqrt{s'\sqrt{s} + s' - 4m^2}} - 1 \right). \qquad (3.9)$$

From this equation, it follows readily that

$$\mathrm{Im} f_0(s) = \frac{1}{8} \rho(s), \qquad (3.10)$$

$$\mathrm{Re} f_0 = \frac{1}{8} \left( \lambda + \mathrm{P.V.} \int_{4m^2}^{\infty} \frac{ds'}{\pi} \frac{\rho(s')}{s' - \frac{4m^2}{3}} \left( \frac{s - \frac{4}{3}m^2}{s' - s} + K_0^{(d=3)}(s', s) \right) \right). \qquad (3.11)$$

In this way we obtain a closed equation for $\rho(s)$ given by

$$\rho(s) = \frac{1}{16\sqrt{s}} \left[ \left(\rho(s)\right)^2 + \left(\lambda + \mathrm{P.V.} \int_{4m^2}^{\infty} \frac{ds'}{\pi} \frac{\rho(s')}{s' - \frac{4m^2}{3}} \left( \frac{s - \frac{4}{3}m^2}{s' - s} + K_0^{(d=3)}(s', s) \right) \right)^2 \right] \qquad (3.12)$$

We start the iteration process by setting $\rho^{(0)}(s) = 0$ and at step $n + 1$, we define $\rho_{n+1}(s)$ as a function of $\rho_n(s)$ using the r.h.s. of (3.12), as was explained in section 2.

One can check by inspection that the first two iterations reproduce exactly the two-loop computation of the amplitude in the $\frac{\lambda}{4!}\phi^4$ theory in a renormalization scheme defined by (3.2). At step 0, we have

$$T^{(0)}(s, t) = \lambda. \qquad (3.13)$$

Plugging this into the unitarity relation (3.12), we get

$$\rho^{(1)}(s) = \frac{\lambda^2}{16\sqrt{s}}. \qquad (3.14)$$

Plugging this into the Mandelstam representation, we get $T^{(1)}(s, t)$, from which the two-loop $\rho^{(2)}(s)$ can be computed using unitarity, see appendix F.2 for details. Starting from three loops, or equivalently at order $\lambda^4$, we would have to introduce the double spectral density in order to produce further actual $\phi^4$ theory graphs, as well as $\eta_{\mathrm{MP}}(s)$. This will be the subject of next sections.

Let us note the following features of the iteration process. First, we see that after the first iteration the spectral density $\rho^{(1)}(s)$ goes to a constant at the two-particle threshold $s = 4m^2$. Upon performing the dispersive integral, this leads to the appearance of logarithms $\int_{4m^2}^{\infty} \frac{ds'}{s' - s} \sim \log(s - 4m^2)$, and therefore we expect perturbation theory to break down close to the two-particle threshold due to large logarithms $(\lambda \log(s - 4m^2))^k$. Second, in the Regge limit, we find that $\rho^{(1)}(s) \sim s^{-1/2}$, and by analyzing the iteration equations it is easy to convince oneself that this behavior is preserved, and should hence be the behavior of the solution.21

<footnote>21We will see later that in the back-reacted solution, this Regge behavior, a fixed Regge pole, is not compatible with unitarity and indeed gets dressed by logarithmic corrections.</footnote>
The nonperturbative behavior of the amplitude close to the two-particle threshold is actually known and it takes the following form [58, 59]

$$ \rho(s) = \frac{4}{\sqrt{s}} \left[ \left( b_0(\lambda) + O\left(\frac{s}{4m^2} - 1\right) + \frac{1}{2\pi\sqrt{s}} \log\left(\frac{s}{4m^2} - 1\right) \right)^2 + \frac{1}{4s} \right]^{-1} . \quad (3.15) $$

The solutions we found matched very well this behavior, and we could observe empirically that $b_0(\lambda) \sim \frac{1}{\lambda}$ for small $\lambda$.<fnref n="22" /> The form of the partial wave is dictated by unitarity near the threshold, see formula (5.6) in [38]. In this way the nonperturbative spectral density at the fixed point actually vanishes, albeit at a slow rate, at the two-particle threshold:

$$ \text{Threshold:} \quad \lim_{s \to 4m^2} \rho(s) \approx \frac{32\pi^2 m}{\log^2\left(\frac{s}{4m^2} - 1\right)} \to 0. \quad (3.16) $$

This is in contrast with the result of the first iteration (3.14). Numerically, we observe that this tension creates a divergence in the iterations if we try to solve the discretized system on a grid which goes too far in the IR (close to the two-particle threshold).

To resolve this tension, we used a set of successive grids, that extend further and further near the threshold. In this way, we find that the iteration algorithm converges and smoothly matches to the expected asymptotic behavior (3.16). We describe this in detail below.

Let us also notice, that given (3.16), unitarity (3.7) implies that $Re f_0(4m^2) = 0$. This condition leads to the following sum rule

$$ \lambda + \int_{4m^2}^{\infty} \frac{ds' \rho(s')}{\pi s' (s' - 4m^2)(s' - \frac{4m^2}{3})} = 0. \quad (3.17) $$

This sum rule is an interesting example of a situation where the “coupling constant” $\lambda$ is dispersive, even though the amplitude $T(s, t)$ goes to a constant (different from $\lambda$) at infinity. On the other hand, for the spectral density $\rho(s)$ in the Regge limit we get

$$ \text{Regge:} \quad \lim_{s \to \infty} \rho(s) \sim s^{-1/2} \to 0. \quad (3.18) $$

Let us also comment briefly on the sign of $\lambda$. The standard $\phi^4$-theory with a repulsive potential bounded from below corresponds in our conventions to $\lambda \le 0$. This is the case for which we find our iterations to converge. Positive $\lambda > 0$, on the other hand, are excluded by the sum rule (3.17) and the fact that unitarity implies that $\rho(s) \ge 0$ in the toy-model. Consistently, we find that our algorithm diverges in this case as we try to remove the IR cutoff.

### 3.1.2 Numerical implementation: hotstart near the two-particle threshold

We now turn to the description of the numerical implementation of this algorithm, and present some explicit results. We use the setup described in section 2.3, and interpolate the single-spectral function $\rho(x)$ on grids $x_0 = 0, \dots, x_N = 1$.

Then we perform iterations of unitarity on this interpolant. We evaluate the value of the spectral density at the grid point $\rho_{n+1}(x_i)$ at $(n+1)$-th iteration using the piece-wise linear $\rho_n(x)$ plugged into the unitarity equation (3.12) after performing the change of

<footnote n="22">²²The same behavior was also argued to hold in the full $\phi^4$ theory in [58, 59].</footnote>
variable $s' \to 4m^2/x'$. We also impose the following boundary conditions for the spectral density $\rho(s)$ which are conserved by the iteration process

$$ \rho(0) = \rho(1) = 0. \qquad (3.19) $$

In the limit of grid spacing becoming small, we observe that the function describes a smooth curve, and we naturally expect the function to approach the true continuous solution.

In practice, the simple algorithm described above requires care because of the threshold behavior, as was explained in the previous subsection. The unitarity equation (3.12) fixes the behavior close to $x = 0$ and $x = 1$ to be of particular form. For the example, close to the two-particle threshold $x = 1$, we expect the universal behavior (3.16) and not a simple linear behavior. Moreover, as discussed above, initializing the iteration process with $\rho_0(x) = 0$ leads to large logarithms close to the threshold, which destabilize the iteration process.

Luckily, there is an efficient way to circumvent the difficulty above and find the desired solution in an arbitrary large range of energies using a version of the hotstart described above.

In the present case, we use the hotstart in the following way. The iteration process is not defined on one fixed grid, but on a family of grids, with increasing near two-particle threshold density of points and higher-energy cut-off. We initialize the iteration process of a finer grid using the fixed point solution of the previous, rougher, grid. In this way, and proceeding with sufficient caution, we are able to build functions which can be extended arbitrarily close to the threshold, and probe arbitrarily far the slow logarithmic decay when $x \to 1$.

This procedure provides a guidance in the fixed-point iteration procedure. The fixed-point on a fine grid is not reachable by a standard iteration because of the large logarithms that destabilize the iterations. However, the guidance provided by the hotstart at each improvement of the grid induces a little “kick” in the right direction which allows to push to arbitrary high energies the problem and can converge to a solution up to arbitrarily high cutoffs. We then find that the solution converges to the expected universal behavior (3.16).

Let us illustrate the procedure described above with an example. Let us say we take $\lambda = -5\pi$, we find a minimum grid first point $x_{N-1}$ ($x_N = 1$). For concreteness we can take $x_{N-1} = 1 - 10^{-4}$. We run the iteration algorithm starting from $\rho_0(x) = 0$ and we find that it converges after $n \sim 10$ iterations. Then, we push the cut-off higher and set $x_{N'-1} = 10^{-8}$. On this new grid, with possibly more elements than the first ($N' \ge N$) the iteration would not have converged starting from $\rho_0 = 0$, but it does so if we set $\rho_0(x)$ to the converged value of the previous round of iterations. In this way we can empirically push the cut-off to values arbitrarily close to 1.<fnref n="23" /> An example of this procedure is provided in figure 9. Note that, as was said above, if we start directly on a cut-off which is too close to the threshold, for instance $10^{-16}$, the algorithm does not converge. Note also that since the Regge limit is trivial, no problem is found near $x \to 0$ and we do not need any caution there.

<footnote n="23">²³For instance, we could easily reach cut-offs of order $1 - 10^{-100}$. A small technical details is that we encountered a bug with mathematica's interpolation routines which would stop working at precision of order `$MachinePrecision` and so we had to interpolate $\rho(1 - x)$ near 0 instead of $\rho(x)$ near 1.</footnote>
<figure><img src="image_10.png" /><figcaption>**Figure 9.** A typical example of the hotstart iterations. Each solid line is the interpolating function $ρ_n(x)$ at the final stage of its iteration ($n = 8$ in this example). Without these intermediate steps, the algorithm would not converge. The list of $x_1$'s (closest point to 1) is listed on the right. The gray line is the threshold behavior given in (3.15) with $b_0(λ) = 8/λ$ determined experimentally.</figcaption></figure>

We present various results regarding the amplitudes we obtain in the list of figures that follows:

• Figure 9 illustrates the hotstart procedure described above,

• in figure 10 we show the linear convergence of the algorithm on a typical process,

• in figure 11 we plot the resulting spectral densities $ρ(x)$ for various values of the coupling $λ$. In figure 12 we compare our results with the perturbative calculation.

• in figure 13 we plot the spin zero partial wave $S_0$ as a function of energy $x$ for $λ = -3π$.

• in figure 14 we plot the maximum value of $\text{Im}S_0$ as a function of $λ$. A linear interpolation suggests that our algorithm would violate unitarity at $λ \sim -6π$.

• in figure 15 we plot the absolute of the spin $J = 2$ partial wave and observe that it violates inelastic unitarity (1.7). This is not surprising given that we set the double spectral density to zero.

In summary, we were able to construct the desired amplitudes in $d = 3$ for $-5π \lesssim λ \le 0$. These amplitudes exhibit the universal near two-particle threshold behavior (3.16) and the Regge limit (3.18). They also show interesting behavior for fixed angle scattering, see figure 16.
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Convergence of the algorithm with hotstart, represented as the convergence of $|S_0(x)|$ to 1. Every 8 steps (the total number of iteration per given grid size), a small glitch is generated by the hotstart but is soon erased (the last step has 20 iterations). Overall, the process induces nice linear convergence, as expected in a fixed-point iteration (the definition of linear convergence was given around eq. (2.21)).</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 11.** Endpoint of the iteration for various couplings. Solid COLORS = numerics, DASHED = fit with universal solution with $b_0 = 8/\lambda$. This demonstrates agreement with universal threshold behavior in the near-threshold limit.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** $\rho(x)$ vs perturbation theory for toy-model in $d = 3$. We observe that at small $\lambda$ the matching between $\rho(s)$ and the perturbative results derived in appendix F is good for all $x$. At higher $\lambda$'s the difference becomes significant.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 13.** We plot the spin zero partial wave $S_0$ as a function of energy. The scattering becomes trivial, $S_0 = 1$, both at the two-particle threshold, $x = 1$, and at high energies, $x = 0$.</figcaption></figure>
<figure><img src="image_15.png" /><figcaption>**Figure 14.** Left panel: Im$S_0$ as a function of energy for various $\lambda$. Right panel: maximum of Im$S_0$ of the fixed-point solution as a function of $\lambda$. For $\lambda \simeq 6\pi$, the solution would appear to exceed the unitarity bound, which is consistent with the fact that in our method $\lambda \lesssim 6\pi$ is the convergence bound. Note further that the bound is reached inside the elastic band $\frac{1}{4} \le x \le 1$, so inelasticity could only help indirectly. It would be interesting to understand the mechanism by which amplitudes with larger couplings $\lambda > 6\pi$ can be generated, as in [45].</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Figure 15.** $S_2$ of the toy-model amplitudes in $d=3$. Unitarity violation is manifest ($|S_2| > 1$), but is not imposed in the model so it is not surprising. This is cured in the full algorithm later when the double spectral density is added, see e.g. figure 27.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>**Figure 16.** Here we plot the absolute value of the toy model amplitude, $|T(s, t)|$, at several different fixed-angles as a function of energy $\frac{s}{m^2}$ in $d=3$. This plot exhibits that the amplitude has two asymptotic constant regimes: Regge, and fixed-angles, which are constants. The plot also shows that, as a finite angle is slowly turned on from the strict Regge limit, the amplitude leaves Regge to join fixed angles at energies given by $s = (\theta/\pi)^{-2}$ (or impact parameters $b \sim \frac{1}{m}$). These two asymptotic constants persist in the back-reacted model, and the amplitudes are indistinguishable by eye, so the reader is invited to think of this plot as the representation of the fully back-reacted quasi-amplitudes themselves.</figcaption></figure>
## 3.2 $d=4$

Next we repeat the same exercise in $d=4$. The region that requires care is now the Regge limit, $x \to 0$.

### 3.2.1 The algorithm

The dispersive representation (3.1) stays intact, however the unitarity equation changes in an important manner. We start again with the partial wave expansion, which now reads

$$T(s,t) = 16\pi \sum_{J=0,J \text{-even}}^{\infty} (2J+1)f_J(s)P_J\left(1+\frac{2t}{s-4m^2}\right), \quad (3.20)$$

$$f_J(s) = \frac{1}{32\pi} \int_{-1}^{1} dz P_J(z) T(s, t(z)), \quad (3.21)$$

where $P_J(z)$ are the ordinary Legendre polynomials. The unitarity equation $|S_0(s)|^2 = 1$ can be written as

$$2\operatorname{Im}f_0(s) = \frac{\sqrt{s-4m^2}}{\sqrt{s}} |f_0(s)|^2, \quad s \ge 4m^2, \quad (3.22)$$

where we used that $S_0(s) = 1 + i \frac{\sqrt{s-4m^2}}{\sqrt{s}} f_0(s)$ in $d=4$.

Our goal is again to solve (3.22) iteratively. For this purpose it is convenient to compute the spin zero partial wave in terms of the spectral density $\rho(s)$. We therefore plug (3.1) into (3.21) to get

$$f_0(s) = \frac{1}{16\pi} \left( \lambda + \int_{4m^2}^{\infty} \frac{ds'}{\pi} \frac{\rho(s')}{s' - \frac{4m^2}{3}} \left( \frac{s - \frac{4}{3}m^2}{s' - s} + K_0^{(d=4)}(s', s) \right) \right), \quad (3.23)$$

$$K_0^{(d=4)}(s', s) = 2 \left( \frac{s' - \frac{4m^2}{3}}{s - 4m^2} \log \frac{s' + s - 4m^2}{s'} - 1 \right). \quad (3.24)$$

We start the iteration by setting $\rho^{(0)}(s) = 0$. The first two iterations are then nothing but a two-loop computation in the $\frac{\lambda}{4!} \phi^4$ theory. Indeed, we have

$$T^{(0)}(s, t) = \lambda. \quad (3.25)$$

Plugging this into the unitarity relation we get

$$\rho^{(1)}(s) = \frac{\lambda^2}{32\pi} \sqrt{1 - \frac{4m^2}{s}}, \quad (3.26)$$

where we used that $\rho(s) = 16\pi \operatorname{Im}f_0(s)$. This is exactly the single discontinuity of the one-loop bubble. In appendix F, we explain the next order of the calculation.

The location of the “trivial” and “nonperturbative” regimes in $d=4$ is exchanged compared to $d=3$. More precisely, we see that the near-threshold behavior of $\rho^{(1)}(s) \sim \sqrt{s-4m^2}$ is preserved under the iterations of unitarity. On the other hand, the Regge limit receives large logarithmic corrections under iterations because $\int^{\infty} \frac{ds'}{s' s' - s} \sim \log s$.
Thus, the nonperturbative spectral density at the fixed point has a square-root singularity at the threshold

$$ \text{Threshold: } \lim_{s \to 4m^2} \rho(s) \sim \sqrt{s - 4m^2} \to 0. \qquad (3.27) $$

In the Regge limit the amplitude exhibits the following behavior which is non-obvious from the perturbation theory point of view

$$ \text{Regge: } \rho(s) \underset{s \to \infty}{\sim} \frac{32\pi^3}{9 \log^2 s} \to 0. \qquad (3.28) $$

Let us now present the nonperturbative solution of the Regge limit for this toy model. It turns out that the expansion is naturally organized in terms of $y = \log s$ and it takes the following form (see appendix D for the derivation)

$$ \rho(e^y) = \frac{32\pi^3}{9} \frac{1}{y^2} \left( 1 - \frac{8}{3} \frac{\log y}{y} + \frac{8 + \frac{\pi^2}{3} - \frac{32}{9} \log y + \frac{16}{3} (\log y)^2}{y^2} + \dots \right) \\ + \frac{b_\infty(\lambda)}{y^3} \left( 1 + \frac{4}{3y} - \frac{4 \log y}{y^4} + \dots \right) \\ + \frac{27}{128\pi^3} \frac{b_\infty^2(\lambda)}{y^4} + \dots, \qquad (3.29) $$

where asymptotically, around $s = \infty$, $b_\infty(\lambda)$ is a free parameter. Empirically, we observed that $b_\infty(\lambda) \sim \frac{1}{\lambda}$ at small coupling. We comment more on this below, when we discuss the numerics.

It is also very interesting to understand the Regge limit of the full amplitude in this case. Given that $\rho(s) \to 0$ at large $s$, we can rewrite (3.1) as

$$ T(s,t) = \left( \lambda - 3 \int_{4m^2}^{\infty} \frac{ds'}{\pi} \frac{\rho(s')}{s' - \frac{4m^2}{3}} \right) + \int_{4m^2}^{\infty} \frac{ds'}{\pi} \rho(s') \left( \frac{1}{s' - s} + \frac{1}{s' - t} + \frac{1}{s' - u} \right). \qquad (3.30) $$

Moreover, consistency with unitarity, or $\text{Re}f_0(\infty) = 0$, in the Regge limit leads to the following sum rule

$$ \lambda - 3 \int_{4m^2}^{\infty} \frac{ds'}{\pi} \frac{\rho(s')}{s' - \frac{4m^2}{3}} = 0, \qquad (3.31) $$

which, together with (3.28), guarantees that $\text{Re}f_0(s) \sim 1/\log s$ at large energies. Note that (3.31) immediately implies that $\lambda \ge 0$ thanks to unitarity $\rho(s) \ge 0$. Therefore, interactions in this model are attractive. It is interesting to notice that $\lambda$ is dispersive (expressible through the discontinuity of the amplitude), despite the fact that in the Regge limit $T(s,t) \sim \text{const}$.

Let us next discuss the sign of $\lambda$ and the behavior of the amplitude at fixed angles $T(s, \theta) \equiv T(s, -\frac{s-4m^2}{2}(1 - \cos \theta))$. In $d=4$, $\lambda < 0$ corresponds to theory which has the Landau pole in the UV. From the point of view of the unitarity iterations, it corresponds to the fact that $\rho(s)$ grows in the UV, which eventually destabilizes the convergence of iterations. On the other hand, for $\lambda > 0$, which from the Lagrangian point of view corresponds to the potential which is unbounded from below, the theory is asymptotically free.
This is the sign for which our algorithm converges and the Regge limit is the one observed above. Let us recall that the running of the coupling in the $\phi^4$ theory is given by

$$ \lambda_p = \frac{\lambda}{1 + \frac{3\lambda}{16\pi^2} \log \frac{p}{m}}, \qquad (3.32) $$

where we adopted our sign convention for $\lambda$ which differs from the one in the standard textbooks (see e.g. Peskin and Schroeder [60], eq. (12.82)). Below we will see that (3.32) correctly captures the qualitative behavior of the fixed angle scattering amplitude as a function of energy $T(s, \theta)$. In particular $T(s, \theta) \sim \frac{32\pi^2}{3} \frac{1}{\log \frac{s}{m^2}}$ at high energies. The amplitude we obtain also exhibits the nontrivial Regge limit as $\theta \to 0$, and transitions between Regge and fixed-angle at $\frac{s}{m^2} \sim 1/\theta^2$. Finally, let us discuss the behavior at $\theta = 0$. In this case we find that $T(s, 0) \to \text{const} \sim \lambda$ when $s \to \infty$.

### 3.2.2 Numerical implementation: “hotstart” in the Regge limit

The numerical implementation of the iterations in $d = 4$ goes through the same basic steps as in three dimensions except that the two-particle threshold and Regge regions get interchanged. In terms of difficulty the threshold region is now trivial, with a square root fall-off, and the Regge region is non-trivial, with $1/(\log s)^2$ behavior leading order behavior, with fixed coefficient, given in (3.29). We show the results of our fixed-point iteration as a function of the coupling $\lambda$ in figure 17.

In $d = 3$, it was easy to work out the $\lambda$-dependence of the asymptotic behavior. In $d = 4$, the question is more subtle, see the discussion around (3.29). In particular, a curious phenomenon arises, when solving asymptotically unitarity with a series expansion in $1/(\log s)^n$. The term $1/(\log s)^3$ is in the kernel of the r.h.s. of unitarity equation, therefore is not canceled and would appear to produce unitarity violation. We have been able to find one correction that allows to generate a $1/(\log s)^3$ term which can then be fixed to remove the $1/(\log s)^3$ term in the l.h.s., which adds a $\log(\log(s))$ term in the numerator, $1/(\log s)^3 \to (a + b \log(\log(s)))/(\log s)^3$. This means that the coupling-dependence of our solution is governed by the $\frac{\log(\log(s))}{(\log s)^3}$ term, and therefore we cannot test it unless we go to super-exponential grids, which is not feasible for our numerics. Nevertheless, assuming that $b_\infty \sim 1/\lambda$, as in 3d, a consistent picture emerges. At large lambda, all curves are on top of each other. At smaller lambda, curves start to differ, but the contribution of the $\frac{\log(\log(s))}{(\log s)^3}$ term is simply not detectable before very high scales.

We present our results according to the following list:

* In figure 17, we show various $\rho(s)$ functions for different couplings, up to $\lambda = 25\pi$.

* In figure 18, we compare these $\rho(s)$ to perturbation theory up to one and two loops.

* In figure 19, we display the S-wave of the amplitude at $\lambda = 5\pi$. We choose this intermediate coupling because it coincides with the maximum coupling we can reach with the full algorithm later. In addition, the back-reaction created by $\rho(s, t)$ is very small in $S_0$, thus the reader is invited to think of this partial wave as that of the amplitude coming from the full amplitude.
<figure><img src="image_18.png" /><figcaption>**Figure 17.** $\rho(s)$ for various couplings in $d = 4$ (color on the right). The Regge limit ansatz is given in (3.29) (dashed line). The only place where our one-parameter enters the ansatz is in front of the $\log(\log(s))$ term, which makes coupling dependence impossible to detect at these scales. For $\lambda \gtrsim 2\pi$, all the dotted curves are essentially on top of each other.</figcaption></figure>

• In figure 20, we display the growth of Im $S_0$ with $\lambda$, to illustrate that our algorithm reaches the end of its convergence range when Im $S_0 \sim 1$. Since amplitudes with larger couplings are generated by other methods, [45, 46], this means that our algorithm needs to be adapted in order to capture these amplitudes.

• In figure 21, we display the unitarity violation observed in $S_2$. Again, in the toy-model we do not impose $|S_J| \le 1$ for $J > 0$ so this is expected.

• In figure 22, we plot the toy-model resulting scattering amplitudes. At low energies, we simply have $T(s,t) \sim \lambda$. At high energies, at fixed angles $(s, t \to \infty, s/t \text{ fixed})$, our amplitude asymptotes to the effective running coupling $\lambda_p$, given in (3.32) in the UV. In the Regge limit, $(s \to \infty, t \text{ fixed})$, the amplitude goes to a constant, which satisfies the sum rule given in (3.31). Increasing the scattering angle $\theta$ from 0 (Regge) to $O(1)$, we observe numerically that the amplitude remains of the Regge-type up to scattering energies $s$ of order $1/\theta^2$, and then transitions to the fixed-angle regime. In the impact parameter space this transition is associated to $b \sim m$.
<figure><img src="image_19.png" /><figcaption>**Figure 18.** Comparison of the nonperturbative numerical solution for the toy-model in $d = 4$ vs perturbation theory: good at low energies, bad at high energies. x-axis is the variable $x = 4m^2/s$, solid line is numerics and dashed is perturbation theory. This plot looks very similar to the one in $d = 3$, except with Regge and two-particle threshold regions interchanged.</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>**Figure 19.** Spin zero partial wave $S_0$ in $d = 4$ for the toy-model amplitude at $\lambda = 5\pi$.</figcaption></figure>
<figure><img src="image_21.png" /><figcaption>**Figure 20.** Left panel: Im$S_0$ for various $\lambda$ in $d = 4$ the. Right panel: maximal value of Im$S_0$ for various $\lambda$. Right: extrapolation of the maximum value of Im $S_0$, which would seem to hit 1 around $40\pi$. Consistently, our algorithm starts to struggle around $25 - 30\pi$.</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 21.** Violation of unitarity in the spin two partial wave $S_2$ in the $d = 4$ toy model. Recall that unitarity implies that $|S_2| \le 1$.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>**Figure 22.** $|\mathcal{T}(s, \theta)|$ as a function of the scattering angle $\theta$ in $d = 4$. We see that the amplitude decays as $1/\log s$ and transitions between two asymptotes: the Regge asymptote (red dashed curve) and the fixed-angle asymptote (black dashed curve), given by the beta function, and defined in (3.32). The transition occurs at energies $s/m^2 \sim \theta^{-2}$ (or impact parameters $b \sim 1/m$, see footnote 32). We have not tried to compute $a_0$ and $a_1$ and obtained them by fitting.</figcaption></figure>
### 3.3 Proof of existence of the amplitudes in $3 < d < 4$

In the discussion above, we have constructed the desired amplitudes numerically. A natural question arises: what happens in the continuum limit? Our explorations of the grid support the idea that nothing dramatic happens and the amplitudes we have constructed exist in the space of continuous functions. This expectation can be made rigorous in $4 > d > 3$ following the argument of Atkinson [24]. The idea is to argue that iterations of unitarity converges to a unique fixed point in the space of Hölder continuous functions. By doing so, we will see that going to $d = 3$ and $d = 4$ requires an extension of the proof, which we currently do not have.

Let us summarize the main steps of the proof leaving the technical details to appendix E. We introduce $s = 4m^2/x$ so that $x \in [0, 1]$ and consider a class of $\rho(x)$ such that

$$ \rho(0) = \rho(1) = 0. \qquad (3.33) $$

We also consider $\rho(x)$ to be Hölder-continuous

$$ |\rho(x_1) - \rho(x_2)| \le \xi|x_1 - x_2|^\mu, \quad 0 < \mu < \min\left(\frac{d-3}{2}, \frac{4-d}{2}\right), \qquad (3.34) $$

where $\mu < \frac{d-3}{2}$ will arise from imposing unitarity close to the two-particle threshold, whereas $\mu < \frac{4-d}{2}$ comes from unitarity in the Regge limit. We introduce the following norm in this functional space

$$ ||\rho|| = \sup_{0 \le x_1, x_2 \le 1} \frac{|\rho(x_1) - \rho(x_2)|}{|x_1 - x_2|^\mu}. \qquad (3.35) $$

One can check that such functions form a complete, normed, linear space, i.e. the Banach space, see [24].

We use S-wave unitarity to write down an iteration procedure

$$ \rho' = \Phi[\rho], \qquad (3.36) $$

see (E.7) for the explicit form.

Given a bound on the norm $||\rho|| \le B$ we would like to show that $||\rho'|| \le B'$ and impose $B' \le B$ to make sure that the iterated spectral density stays in the same space of functions. We then would like to show that the map is actually contracting meaning that

$$ ||\rho'_2 - \rho'_1|| \le k||\rho_2 - \rho_1||, \quad k < 1. \qquad (3.37) $$

Given (3.37), the fixed point of the unitarity equation (3.36) $\rho_* = \Phi[\rho_*]$ is unique and iterations of unitarity converge exponentially fast

$$ ||\rho_* - \rho_N|| \le 2B \frac{k^N}{1-k}. \qquad (3.38) $$

In appendix E we derive $B'$ and $k$ explicitly in terms of $(\lambda, \mu, d, B)$. Here we simply plot the maximal value of $|\lambda|$ as a function of $d$ for which we can prove that the map to be contracting, see figure 23.
<figure><img src="image_24.png" /><figcaption>Figure 23. The maximal absolute value of the coupling $|\lambda|$ for which we can proof convergence of the iterations as a function of the number of spacetime dimensions $d$.</figcaption></figure>

The characteristic feature of this plot is that the range of $\lambda$ for which the convergence is proven shrinks to zero as $d \to 3$ or $d \to 4$. How is it consistent with our explicit results described above? Intuitively, the reason is clear. The amplitudes there involved large amount of cancelations between various terms in the amplitude, see (3.17) and (3.31). Such cancelations are not taken into account in the proof in appendix E, or in the original Atkinson proofs, where we bound each term separately.

Similarly, we observed both in $d = 2$, see [26], and throughout the whole parameter space explored in this work that the iterations perform better than what is expected from the proofs.

There is another way to understand how we managed to find the fixed point in $d = 3$ and $d = 4$. Effectively, what our algorithm does is it modifies the unitarity equation (due to the cutoff at some maximal or minimal energy) as follows. Let us for concreteness discuss the case of $d = 4$.<fnref n="24" /> The basic idea is to modify the unitarity equation

$$ \rho'(x) = \frac{h(x)}{32\pi} (1-x)^{\frac{1}{2}} \left( (16\pi \text{Re}f_0[\rho](x))^2 + \rho(x)^2 \right), \quad (3.39) $$

where we can for example choose (the precise form is not essential, as long as $h(0) = 0$)

$$ h(x) = \theta(x - x_0) + \frac{x}{x_0} \theta(x_0 - x). \quad (3.40) $$

The role of $h(x)$ is to effectively soften the behavior of $\rho'(x)$ in the Regge limit $x \to 0$. With this modification we can consider a class of Holder-continuous functions with $\mu = 1/2$ and run the argument of the previous section. All the bounds will now depend on $x_0$.

Physically, having such a cut-off is not desirable. Indeed, in a physical theory in the presence of S-wave particle production we have

$$ \rho'(x) = \frac{1}{32\pi} (1-x)^{\frac{1}{2}} \left( (16\pi \text{Re}f_0[\rho](x))^2 + \rho(x)^2 \right) + \rho_{\text{MP}}(x), \quad \rho_{\text{MP}}(x) \ge 0. \quad (3.41) $$

<footnote n="24">²⁴The case of $d = 3$ is essentially identical, where the Regge limit is replaced by the near two-particle threshold region.</footnote>
It is therefore curious to see how we effectively solved this problem in the numerical implementation of the algorithm.

Effectively, when adopting the hotstart for a family of grids, we considered iteration-dependent modification of the unitarity equation

$$h_n(x) = \theta(x - x_n) + \frac{x}{x_n} \theta(x_n - x), \qquad (3.42)$$

where we chose $x_n$ such that $\lim_{n \to \infty} x_n = 0$. We observed convergence and smooth extrapolation of the solutions to the UV as we remove the cutoff $x_n$, where we could smoothly match to the analytic UV complete solution of the model.

### 3.4 $d > 4$: nonrenormalizability and unitarity

Above we discussed the cases $d = 3, 4$ and $3 < d < 4$. Here we briefly comment what happens if we try to extend the same ideas to $d > 4$, and we also make contact with the notion of renormalizability. The main difference comes from the form of the unitarity equation that takes the form

$$2\mathrm{Im}f_0(s) = \frac{(s - 4m^2)^{\frac{d-3}{2}}}{\sqrt{s}} |f_0(s)|^2, \quad s \ge 4m^2. \qquad (3.43)$$

Compared to $d \le 4$, the phase space pre-factor $(s-4m^2)^{\frac{d-3}{2}}/\sqrt{s}$ now grows in the Regge limit as $s^{\frac{d-4}{2}}$. This leads to the fact iterations of unitarity with the starting point $T^{(0)}(s, t) = \lambda$ lead to polynomially growing amplitudes with the power growing with the number of iterations.<fnref n="25" /> In particular, the dispersion relations with a finite number of subtractions will eventually break down. Of course, this is nothing but the fact that $\frac{\lambda}{4!}\phi^4$ is nonrenormalizable in $d > 4$, and here we just see it as a simple consequence of unitarity. Therefore, simple iterations of unitarity starting from $T^{(0)}(s, t) = \lambda$ will not work in this case.

However, as hopefully became clear from the discussion of the $d = 4$ case, this is not a problem per se. Indeed, even in $d = 4$ to get a convergent amplitude we effectively had to “UV complete” the amplitude by supplementing it with the correct nonperturbative Regge limit (namely that $\rho(s) \to 0$), which then made the iterations convergent upon increasing the UV cutoff. One can imagine using the same strategy in $d > 4$. For example, we can try to “eikonalize” the leading order amplitude, which should also be relevant for applying the methods discussed here to gravity. We do not explore the case of $\phi^4$ in $d > 4$, or other nonrenormalizable theories such as gravity here, and leave this interesting problem for the future.

Finally, note that a related problem occurs in 2d [26], where to describe a certain class of CDD amplitudes (with unequal number of zeros and poles), we had to specify boundary conditions, and use the Newton-Raphson method to force the iteration to preserve them.

<footnote n="25">²⁵At the first step, we get $f_0^{(1)} \sim s^{(d-4)/2}$, which then gives $f_0^{(2)} \sim s^{(d-4)/2}s^{d-4}$, and so on.</footnote>
<figure><img src="image_25.png" /><figcaption>**Figure 24.** Maximal deviation of the absolute value of $J = 0$ partial wave from 1 as a function of number of iterations $n$ for $\lambda = -3\pi$ in $d = 3$. We observe linear convergence, see (2.21).</figcaption></figure>

# 4 Quasi-elastic amplitudes in $d = 3$

In this section, we present the results of the numerical implementation of the iteration algorithm presented in section 2 in $d = 3$ spacetime dimensions. The key difference compared to the previous section is that here $\rho(s, t) \neq 0$, and that we can now implement unitarity for partial waves with $J > 0$.

We start with a presentation of the amplitudes with a reasonably large coupling, $\lambda = -3\pi$, which is close to the boundary of the range of convergence of our algorithm. Then, we move to describe the coupling dependence of certain elements of the amplitude, partial waves, etc. We also distinguish two iteration schemes for $(\lambda, \eta_{\text{MP}}(s), \rho_{\text{MP}}(s,t))$ introduced earlier: 2QE given by $(\lambda, 0, 0)$; 2PR given by $(\lambda, (2.24), 0)$. We start the section with the 2QE case.

## 4.1 2QE amplitude for $\lambda = -3\pi$

Below we present results for the 2QE amplitude in $d = 3$ and $\lambda = -3\pi$. Recall that this algorithm is characterized by the fact that the $J = 0$ partial wave is purely elastic at all energies.

To illustrate the convergence of the algorithm, we first plot the maximal deviation of the absolute value of the $J = 0$ partial wave from 1 in figure 24. We observe linear convergence and reach the precision $\sim 10^{-9}$ after $\sim 25$ iterations. For smaller couplings, the convergence is much faster and usually only about 5 iterations are needed.

### 4.1.1 Single and double spectral functions

Next we present various features of the resulting amplitude. The single spectral function $\rho(s)$ is depicted in figure 25. It exhibits logarithmic behavior close to the two-particle threshold $s \to 4m^2$, and decays in the Regge limit $s \to \infty$ in a way which is distinctively different from $\sqrt{x}$, as observed from results of section 3.1.

The double spectral density $\rho_{\text{el}}(s, t)$ is depicted in figure 26. It peaks close to $x = 1$, and decays in the Regge limit $x \to 0$. Recall that $\rho(x, y) = \rho_{\text{el}}(x, y) + \rho_{\text{el}}(y, x)$, where
<figure><img src="image_26.png" /><figcaption>**Figure 25.** Left panel: single discontinuity $\rho(x)$ as a function of energy $x$. The results for $\rho(s, t) \neq 0$ obtained in this section are depicted in red, whereas the toy-model results of the previous section are shown in blue (dashed). Middle panel: the two-particle threshold region $x \sim 1$, we do not observed any visible difference compared to the toy-model result. Right panel: the Regge region $x \sim 0$. Here we see that $\rho(s)$ is significantly corrected compared to the toy model.</figcaption></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 26.** Double discontinuity $\rho_{\text{el}}(x, y)$ as a function of energies in $d = 3, \lambda = -3\pi$. Left panel: the blue line represents the leading Landau curve to the left of which $\rho_{\text{el}}(x, y) = 0$. Right panel: we take the fixed energy slice $\rho_{\text{el}}(0.5, y)$ to exhibit a nontrivial emergent Regge limit as $y \to 0$. By the dashed line we exhibit the double discontinuity coming from the three-loop Aks graph, see figure 61.</figcaption></figure>

$\rho_{\text{el}}(x, y)$ satisfies the Mandelstam equation (1.6) and has nonzero support only above the leading Landau curve, namely $y \le \frac{1-x}{4}$.

In contrast to the toy-model consideration, we do not have an analytic solution to the Regge limit of the $\rho(x)$ and $\rho(x, y)$ for the 2QE amplitudes. Due to the Gribov's theorem [29], it cannot be a simple power, e.g. $x^{1/2}$. Indeed, this is what we observe in the right panel of figure 25. It would be very interesting to try to solve the Regge limit of the 2QE model analytically.
<figure><img src="image_28.png" /><figcaption>**Figure 27.** The absolute value of $|S_J|$ as a function of energy $x$ in $d = 3$, $\lambda = -3\pi$. We see that the amount of scattering in the higher spin partial waves quickly decreases as a function of spin $J$. Note that there is a dynamically emergent scale in the problem where each partial wave peaks at $x \sim 0.01$, which corresponds to $\frac{s}{4m^2} \sim 10^2$. The dashed line is at $x = \frac{1}{4}$, or, equivalently, $s = 16m^2$, which separates the elastic region $x \ge \frac{1}{4}$ from the inelastic one $x < \frac{1}{4}$.</figcaption></figure>

### 4.1.2 Partial waves, impact parameter

Next we analyze the partial waves of the amplitude. We find that the spin-zero partial wave, the S-wave, is essentially identical to what we got in the toy model, see figure 13. In particular, the S-wave scattering decays both in the Regge limit $x \to 0$ and near the two-particle threshold $x \to 1$.

In contrast the structure of the higher spin partial wave is quite different. We plot the absolute value of the higher spin partial waves $J > 0$ in figure 27. We see that this time they obey unitarity. Moreover, the departure from 1 is the signal of particle production. The appearance of particle production is expected from the Aks theorem [27], but here we see the quantitative amount of particle production necessary for elastic unitarity and crossing. As expected, the amplitude is *quasi-elastic*, in other words, the amount of particle production is quite small. In the physical $\phi^4$ theory we expect particle production to be bigger, since in that case $\eta_{\text{MP}}(s)$ and $\rho_{\text{MP}}(s, t)$ are non-zero, essentially, due to multi-particle unitarity. We shall see below, when we analyze the coupling dependence, the magnitude of this production is of order $\lambda^4$ at small lambda, and deviates from it at larger couplings.

We can also consider scattering at fixed impact parameters $b \equiv \frac{2J}{\sqrt{s-4m^2}}$ where we recall that $p = \sqrt{s-4m^2}$ is the spatial momentum. As expected, the amplitude decreases exponentially fast as a function of impact parameters, see figure 28.

Finally, we can analyze unitarity in the partial waves at higher spins. In agreement with general arguments [38, 61] it becomes mostly inelastic, see figure 29.

### 4.1.3 Amplitude

Finally, we can look at the amplitude itself. We plug our expressions for the single and double spectral function into the Mandelstam representation, and compute the resulting amplitude numerically.
<figure><img src="image_29.png" /><figcaption>**Figure 28.** The absolute value of $1 - |S(s, b)|$ as a function of impact parameter $b$ in $d = 3$, $\lambda = -3\pi$. As expected interactions quickly decay as a function of impact parameter, in a Yukawa-like, exponential fashion, see appendix I. Here the dashed line stands for the fit $\sim e^{-2bm}/bm$.</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 29.** Higher-$J$ inelasticity in $d = 3$ as a function of spin $J$. Elastic scattering corresponds to 0 on this plot, meaning that at large spin scattering is weak but mostly inelastic as it goes to 1. This property of partial waves is a direct consequence of elastic unitarity and crossing, see [38] for the derivation of this fact.</figcaption></figure>

Since the double spectral function is numerically small compared to $\rho(s)$, by eye, the amplitude looks indistinguishable from that obtained in the toy-model, so we do not represent them again here and refer the reader to figure 16. In particular, it also asymptotes the running coupling at high energies for non-zero angles.

Next we plot the forward scattering by setting $t = 0$. Via the optical theorem, this is related to the cross-section. For the real part we find figure 30, and for the imaginary part figure 31. Interestingly, we see that at high energies scattering becomes *mostly inelastic* albeit weak.
<figure><img src="image_31.png" /><figcaption>**Figure 30.** We plot the real part of the forward amplitude $T(s, 0)$ for $s \ge 4m^2$ in $d = 3$, $\lambda = -3\pi$. We see that it monotonically decreasing and it approaches a constant at infinite $s$.</figcaption></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 31.** We plot the imaginary part of the forward amplitude $T(s, 0)$ in $d = 3$, $\lambda = -3\pi$. Via the optical theorem it is trivially related to the scattering cross section. We plot three different cross sections: total, elastic, inelastic. Inelastic part of scattering in this model comes fully from the crossed terms $\rho_{el}(y, x)$ in the double spectral density. We see that at high energies scattering becomes more and more inelastic.</figcaption></figure>

## 4.2 2PR amplitude

Next consider the 2PR iteration scheme. In practice, it means that instead of setting $\eta_{MP}(s)$ to zero, we iterate it using (2.24). Essentially, the only visible difference for this scheme is that spin zero scattering is not purely elastic anymore, see figure 32. We also plot the corresponding “multi-particle” spectral density.

## 4.3 Coupling dependence

Here we present some elements of the coupling-dependence of the quasi-elastic amplitudes in $d = 3$, in figure 33. In $d = 3$, we control all couplings from small to the largest we can reach, contrary to $d = 4$ where we cannot control the small couplings, as explained below. The partial waves look smooth, and depend at leading order on $\lambda^4$, which corresponds to 3-loop Aks graph, and receive corrections from higher orders.
<figure><img src="image_33.png" /><figcaption>**Figure 32.** The spin zero partial wave in the 2PR scheme, $d = 3$, $\lambda = -3\pi$. Left panel: the absolute value which is strictly below 1 signifies particle production. Right panel: the multi-particle spectral density $\rho_{\text{MP}}(s)$ as given by (2.25) at the fixed point of the iteration algorithm.</figcaption></figure>

<figure><img src="image_34.png" /><figcaption>**Figure 33.** Inelasticities in $S_J$ for $J = 0$ (2PR, top left) and $J = 2$ (QE, top right). The partial waves are renormalized by a factor proportional to $\lambda^4$, which is the leading order effect for particle production. In these plots, $a_0 = 0.021$ and $a_2 = 0.00021$. The bottom row shows the maximum value of $1 - |S_J|$ renormalized by $\lambda^4$: at small lambda, $a_0$ and $a_2$ are fitted to go to 1, at larger $\lambda$ the inelasticity deviates from pure $\lambda^4$ effects.</figcaption></figure>
# 5 Quasi-elastic amplitudes in $d = 4$

In this section, we present results for the quasi-elastic amplitude in four dimensions for $\lambda = 5\pi$. The result was obtained by hot-starting on the corresponding solution to the toy-model (that had zero double discontinuity). Why this procedure helps can be easily understood: it is because $\rho_{el}(s, t)$ is numerically small, and thus its back reaction on $\rho(s)$ is not big and thus it is worth coming close to the final answer first for $\rho(s)$ with the toy-model.

## 5.1 2QE amplitude for $\lambda = 5\pi$

To evaluate the convergence of the algorithm, we can measure how fast the absolute value of the S-wave converges to one. Recall that in the quasi-elastic case, we force the S-wave to be purely elastic at all energies. We have found that within 20 iterations we have reached the precision accuracy of our computation which for $|S_0(x)|$ is $\sim 10^{-7}$, see figure 34.

We now present various aspects of the amplitude after 20 iterations.

### 5.1.1 Single and double spectral functions

The single discontinuity is shown in figure 35. As in 3d, we observe that the spectral density stays identical to the toy-model one at low energies, $x \to 1$, and deviates from it at high energies, $x \to 0$.

Next we plot the double discontinuity in figure 36. Recall that $\rho(x, y) = \rho_{el}(x, y) + \rho_{el}(y, x)$, where $\rho_{el}(x, y)$ satisfies the Mandelstam equation (1.6) and has nonzero support only above the leading Landau curve, namely $y \le \frac{1-x}{4}$.

### 5.1.2 Partial waves, impact parameter

Consider next the structure of the partial waves as a function of energy $S_J(x)$, where $x = \frac{4m^2}{\hbar^2}$. The scattering is dominated by the spin $J = 0$ partial wave which is by indistinguishable from the corresponding toy-model result shown in figure 19.

<figure><img src="image_35.png" /></figure>

**Figure 34.** Maximal deviation of the absolute value of S-wave from 1 in $d = 4$ as a function of number of iterations $n$. We observe linear convergence and reach the precision $\sim 10^{-7}$ after $\sim 20$ iterations.
<figure><img src="image_36.png" /><figcaption>**Figure 35.** Single discontinuity $\rho(x)$ as a function of energy $x$ in $d = 4$, $\lambda = 5\pi$. Here we focus on the Regge limit $x \to 0$ which reveals some emergent nontrivial behavior that differs from the toy-model result.</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>**Figure 36.** Double discontinuity $\rho_{el}(x, y)$ as a function of energies in $d = 4$, $\lambda = 5\pi$. Left panel: the blue line represents the leading Landau curve to the left of which $\rho_{el}(x, y) = 0$. Right panel: we take the fixed energy slice $\rho_{el}(0.5, y)$ to exhibit the nontrivial emergent Regge limit as $y \to 0$. The dashed line represents the result from the three-loop Aks diagram, see figure 1 the third diagram in the third row. We see good agreement at low energies, $x \to 1$, and different behavior at high energies, $x \to 0$.</figcaption></figure>

Scattering in the higher partial waves $S_{J>0}$ is much weaker, nevertheless it is nonzero in agreement with the expectation based on elastic unitarity and the Aks theorem [27]. We present a first few partial waves in figure 37.

We can also consider scattering at fixed impact parameters $b \equiv \frac{2J}{\sqrt{s-4m^2}}$. As expected the amplitude decreases quickly as a function of impact parameters, see figure 38.

### 5.1.3 Amplitude

Next let us consider scattering at fixed angle scattering. We find that the fixed angle amplitude looks essentially the same as the corresponding plot for the toy model, see figure 22. We find that the amplitude acquires its maximum value at $s = 4m^2$ and closely follow the one-loop running of the coupling $\lambda(s) = \frac{\lambda}{1 + \frac{3\lambda}{16\pi^2} \log \frac{\sqrt{s}}{\sqrt{4/3m}}}$ in $\phi^4$ theory. Therefore, the scattering amplitude that we constructed describes interactions that weaken in the UV, namely asymptotic freedom.
<figure><img src="image_38.png" /><figcaption>**Figure 37.** The absolute value of $|S_J|$ as a function of energy $x$ in $d = 4$, $\lambda = 5\pi$. We see that the amount of scattering in the higher spin partial waves quickly decreases as a function of spin $J$. The dashed line is at $x = \frac{1}{4}$, or, equivalently, $s = 16m^2$, which separates the elastic region $x \ge \frac{1}{4}$ from the inelastic one $x < \frac{1}{4}$. Note that there is a dynamically emergent scale in the problem where each partial wave peaks.</figcaption></figure>

<figure><img src="image_39.png" /><figcaption>**Figure 38.** The absolute value of $1 - |S(s, b)|$ as a function of impact parameter $b$ in $d = 4$, $\lambda = 5\pi$. As expected interactions quickly decay as a function of impact parameter, in a Yukawa-like, exponential fashion, see appendix I. Here the dashed line stands for the fit $\sim e^{-2bm}/(bm)^2$.</figcaption></figure>

Finally, let us consider the forward scattering (or zero angle scattering). It is convenient to consider the real, figure 39, and imaginary parts, figure 40, separately.

Let us summarize the salient features of the nonperturbative scattering amplitude that we have constructed:

* it obeys elastic and inelastic unitarity (in addition to crossing and analyticity). Let us clarify what we mean by this exactly. We impose elastic unitarity numerically on a finite set of grid points. For example, for $|S_0|$ we observed that violations of elastic unitarity are $\sim 10^{-7}$. For inelastic unitarity we have checked both a few low spin partial waves, as well as scattering at various impact parameters and energies.
<figure><img src="image_40.png" /><figcaption>**Figure 39.** We plot the real part of the forward QE amplitude $T(s, 0)$ in $d = 4$. At low energies it start at the values around $\lambda$ and it goes to a constant at high energies. The constant is $\sim \frac{\lambda}{3}$ which is related to the sum rule (3.31) in the toy model. Here we simply observe that the same behavior continues in the full back-reacted model.</figcaption></figure>

<figure><img src="image_41.png" /><figcaption>**Figure 40.** We plot the imaginary part of the forward amplitude, $\text{Im}\,T(s, 0)$. Via the optical theorem it is trivially related to the scattering cross section. We show three different cross sections: total, elastic, inelastic. The inelastic part of scattering in this model comes fully from the crossed terms $\rho_{\text{el}}(y, x)$ in the double spectral density. We see that at high energies scattering becomes more and more inelastic.</figcaption></figure>
<figure><img src="image_42.png" /><figcaption>**Figure 41.** Spin zero partial wave $S_0(x)$, and $\rho_{\text{MP}}(x)$ in the 2PR scheme, $d = 4$, $\lambda = 5\pi$, as a function of energy. We see that the spin zero partial wave takes the same characteristic shape, as the higher spin partial waves depicted in figure 37.</figcaption></figure>

* it is elastic at any energies for the S-wave projection. This is what we called quasi-elastic amplitude in this paper.

* it has non-zero particle production for partial waves with $J > 0$. This is related to the fact that the double spectral density of the amplitude has a correct support above the leading Landau curves. In fact, we see that as energy grows scattering becomes more and more inelastic.

* it exhibits asymptotic freedom. By this we mean that the scattering amplitude at fixed angles decays logarithmically as a function of energy.

* it has a nontrivial Regge limit. For example, for $t = 0$ we have found that the real part of the amplitude goes to a constant, whereas the imaginary part decays logarithmically in energy. The nontrivial Regge limit of the amplitude is forced upon us by elastic unitarity [29].

Finally, we believe that the plots presented in this section will not change as one takes the continuous limit of the amplitude. Our belief is based on experimenting with various grids and observing stability of the results presented here, but of course an actual proof using the fixed point methods in the space of continuous functions is highly desirable.

## 5.2 2PR amplitudes

Next we consider the same amplitudes but in the 2PR scheme which means that scattering in the $J = 0$ partial wave is not purely elastic. Structurally and for the couplings $\lambda$ analyzed in this paper, the 2PR amplitudes look very similar to the 2QE amplitudes. The main difference is that now $S_0$ has nonzero inelasticity figure 41:

## 5.3 Coupling-dependence

We find again that inelasticity in $S_0$ and $S_2$ is driven by $\lambda^4$ at low energies and deviates from it at larger energies, see figure 42.
<figure><img src="image_43.png" /><figcaption>**Figure 42.** All the inelasticities, renormalized by leading order $(\lambda/(16\pi^2))^4$ Aks graph effect, in $d = 4$ for $J = 0$ (top left, 2PR models) and $J = 2$ (top right, 2QE models) for couplings $\lambda \in [\pi; 5\pi]$. The higher waves $J > 0$ are visually indistinguishable in 2PR and 2QE models so we display a mixture of the two models. The rescaling by $\lambda^4$ of $1 - |S_J|$ allows to simultaneously show various couplings on the same plot, for otherwise only the $5\pi$ curve would be visible and the other would appear completely flat, due to the approximate $\lambda^4$ dependence. Bottom: maximum value of $(1 - |S_J|)$ renormalized by $\lambda^4$, and a an extra coefficient $b_J$ such that at $\lambda = \pi$, the curves go to one. This choice is arbitrary and simply allow to show that, in $d = 4$, contrary to $d = 3$, the coupling dependence is approximately identical in $J = 0$ and $J = 2$. We used empirical fitting to find $b_0 = 9.5$ and $b_2 = 1.33$. Finally, importantly, we did not display smaller couplings $\lambda$ because the numerics appeared not trustable, possibly due to a grid/cut-off effect. We show the full set of data in appendix H.</figcaption></figure>

## 6 Multi-particle double discontinuity: a case of acnode in $d = 4$

In the previous sections we worked with the approximation where the multi-particle double spectral density was set to zero, $\rho_{\text{MP}}(s, t) = 0$. While this approximation might be a good leading order calculation, e.g. for $\phi^4$ theory in $d = 3$, it is not exactly true in physical theories. Moreover, in confining gauge theories, for example in QCD, we expect $\rho_{\text{MP}}(s, t)$ to be sizable and important [62]. Therefore it is important to understand better the role of $\rho_{\text{MP}}(s, t)$ and its effect on the physical amplitude.

Here we consider the simplest graph which has three-particle cut both in the $s$ and $t$ channel, namely the acnode graph, see figure 43.<fnref n="26" /> If we are to run the iteration process in the presence of bound states it would be the simplest graph that develops $\rho_{\text{MP}}(s, t)$.

<footnote n="26">²⁶This graph is also known as the four-point kite graph, see e.g. [63, 64] for the recent analysis of this graph using geometric methods.</footnote>
<figure><img src="image_44.png" /><figcaption>**Figure 43.** The acnode graph produces the leading contribution to $\rho_{MP}$ in a theory with both the $\phi^3$ and $\phi^4$-type vertices.</figcaption></figure>

For the $\phi^4$-type theory (or $\mathbb{Z}_2$ symmetric scattering) the first analogous graph is the open envelope, see figure 1.

We would like to compute $\rho_{MP}(s, t)$ for the acnode. This problem was analyzed by Gribov and Dyatlov in a short, insightful, but not fully explicit paper [48]. Moreover, their results have not been checked to the best of our knowledge. Here we fill these gaps, by writing down an explicit formula for the double discontinuity and checking that it correctly reproduces the single discontinuity upon doing a dispersive integral.

Let us describe the derivation that goes through several steps. The starting point is the three-body unitarity integral for the discontinuity of the acnode that takes the form

$$2T_s = \frac{1}{3!} \int \prod_{i=5}^{7} \frac{d^3 \vec{q}_i}{(2\pi)^3 (2E_{\vec{q}_i})} (2\pi)^4 \delta(p_1 + p_2 - \sum_{i=5}^{7} q_i) T(p_1, p_2; q_5, q_6, q_7) T^*(q_5, q_6, q_7; p_3, p_4), \qquad (6.1)$$

where $s = (p_1 + p_2)^2$ and $t = (p_1 - p_3)^2$. The explicit formulas for the amplitudes that enter into the unitarity integral take the form

$$T(p_1, p_2; q_5, q_6, q_7) = -\frac{1}{t_{15} - m^2}, \qquad (6.2)$$

$$T^*(q_5, q_6, q_7; p_3, p_4) = -\frac{1}{t_{46} - m^2}, \qquad (6.3)$$

where $t_{15}^2 = (q_5 - p_1)^2$ and similarly for $t_{46}$.

The unitarity integral can be simplified to the following form, see appendix J for details and the definitions of various parameters,

$$2T_s = \frac{1}{32s} \frac{1}{3!} \frac{1}{(2\pi)^3} \int^{(\sqrt{s}-m)^2} \frac{ds_57 ds_{67} \theta(1-(z_{56}^0)^2)}{4|\vec{p}_1||\vec{p}_4||\vec{q}_5||\vec{q}_6|} \int_{-1}^{1} dz_{26} \frac{1}{\sqrt{K(z_{15}^0, z_{56}^0, -z_{26}) K(z, z_{26}, z_{46}^0)}} \qquad (6.4)$$

This agrees with the expression in [48] (up to 1/3! due to the identical particles that is not imposed in that paper). The expression above is amenable for numerical evaluation. It also admits the following dispersive representation

$$T_s(s, t) = \int_{t_{LC}(s)}^{\infty} \frac{dt' \rho_{ac}(s, t')}{\pi (t' - t)}, \qquad (6.5)$$
where $t_{LC}(s)$ is the position of the acnode Landau curve. A convenient representation of the acnode Landau curve is given by the following formula [44]

$$s(\phi) = 5 + 4 \cos \phi + 2 \left( \frac{3}{2} + \cos \theta + \cos \phi \right) \frac{\sin \phi}{\sin \theta}, \qquad (6.6)$$

$$t(\phi) = 5 + 4 \cos \phi + 2 \left( \frac{3}{2} + \cos \theta + \cos \phi \right) \frac{\sin \theta}{\sin \phi}, \qquad (6.7)$$

where $\theta + \phi = \frac{\pi}{3}$ and $0 \le \phi \le \frac{\pi}{3}$.

Our task below will be to derive a formula for $\rho_{ac}(s, t')$ in (6.5) starting from the unitarity integral (6.4). For simplicity below we set $m = 1$.

## 6.1 Computing the $s_{67}$ integral

A convenient strategy is to start with the $s_{67}$ integral that can be easily done explicitly, see e.g. [65]. Indeed, one can check that the integral takes the following form

$$\int_{s_{67}^{(1)}}^{s_{67}^{(2)}} \frac{ds_{67}}{\sqrt{(s_{67} - s_{67}^{+})(s_{67} - s_{67}^{-})}}. \qquad (6.8)$$

From now on we switch $z_{26} \to t_{26}$ via the linear map

$$t_{26} = \frac{1}{2} \left( z_{26} \sqrt{\frac{s-4}{s}} \sqrt{(s-1)^2 + s_{57}^2 - 2(s+1)s_{57}} - s + s_{57} + 3 \right). \qquad (6.9)$$

To present the result for the single discontinuity $T_s$ it is convenient to introduce the following variables

$$\begin{align*}
\tilde{t}_{26}^{\pm} &= (\sqrt{s_{57}} \pm 1)^2, \\
t_{26}^{(1,2)} &= 1 + \frac{s_{57}}{2} \mp \frac{\sqrt{3}}{2} \sqrt{s_{57}(4-s_{57})}, \\
t_{26}^{\mp} &= \frac{\pm \frac{2\sqrt{s^2-(s_{57}+2)s+(s_{57}-1)^2}\sqrt{t\sqrt{s+t-4}}}{\sqrt{s}} + (-s+s_{57}+1)t}{4-s} + 1, \\
t' &= \frac{1}{2} \left( -s+s_{57} - \sqrt{\frac{s-4}{s}} \sqrt{(s-1)^2 + s_{57}^2 - 2(s+1)s_{57} + 3} \right), \\
t'' &= \frac{1}{2} \left( -s+s_{57} - \sqrt{\frac{s-4}{s}} \sqrt{(s-1)^2 + s_{57}^2 - 2(s+1)s_{57} + 3} \right). \tag{6.10}
\end{align*}$$

Note that they are functions of $s, t, s_{57}$ only.

In terms of these variables the discontinuity of the amplitude takes the following compact form

$$2T_s(s,t) = \frac{1}{768\pi^3\sqrt{s(s-4)}} \int_4^{(\sqrt{s}-1)^2} ds_{57} \int_{t'}^{t''} dt_{26} \frac{1}{\sqrt{(t_{26}-t_{26}^-)(t_{26}-t_{26}^+)} \times \frac{\log \frac{t_{26}-s_{57}+1+\sqrt{(1-4/s_{57})(t_{26}-t_{26}^+)(t_{26}-t_{26}^-)}}{t_{26}-s_{57}+1-\sqrt{(1-4/s_{57})(t_{26}-t_{26}^+)(t_{26}-t_{26}^-)}}}{\sqrt{(t_{26}-t_{26}^+)(t_{26}-t_{26}^-)}}. \qquad (6.11)$$
<figure><img src="image_45.png" /><figcaption>**Figure 44.** Deformation of the $t_{26}$ contour relevant for the double spectral density computation. The original contour goes along the real axis $t'' \ge t_{26} \ge t'$. As we increase $t$ the singularity of the integrand $t_{26}^{-}$ enters the integration contour and drags it away from the original contour. When $t_{26}^{-}$ collides with other singularities after that the integral develops a singularity, due to the pinching of the contour.</figcaption></figure>

## 6.2 Analytic continuation in $t$

Next we would like to analytically continue the expression above from negative to positive $t$. We note that the last two integrals in (6.11) look too hard to be computed explicitly, but it is possible to understand qualitatively how the integration contour deforms during the continuation.

Let us introduce the function

$$F(s, t, s_{57}) = \int_{t'}^{t''} dt_{26} \frac{1}{\sqrt{(t_{26} - t_{26}^{-})(t_{26} - t_{26}^{+})}} \frac{\log \frac{t_{26} - s_{57} + 1 + \sqrt{(1-4/s_{57})(t_{26} - \tilde{t}_{26}^{+})(t_{26} - \tilde{t}_{26}^{-})}}{t_{26} - s_{57} + 1 - \sqrt{(1-4/s_{57})(t_{26} - \tilde{t}_{26}^{+})(t_{26} - \tilde{t}_{26}^{-})}}}{\sqrt{(t_{26} - \tilde{t}_{26}^{+})(t_{26} - \tilde{t}_{26}^{-})}} \quad (6.12)$$

In terms of $F$, we have that

$$2T_s(s, t) = \frac{1}{768\pi^3\sqrt{s(s-4)}} \int_4^{(\sqrt{s-1})^2} ds_{57} F(s, t, s_{57}). \quad (6.13)$$

We now would like to understand the singularities of $F(s, t, s_{57})$ as we increase $t$. One can trace the motion of the singularities of the integrand and observe that as we increase $t$, $t_{26}^{(-)}$ enters the integration contour and drags it, see figure 44.

The singularities of $F(s, t, s_{57})$ arise when the deformed contour gets pinched, which happens for $s_{57}^{\pm}(s, t)$ such that

$$s_{57}^{\pm}(s, t) : \quad t_{26}^{-}(s, t, s_{57}) = t_{26}^{(2)}(s, t, s_{57}). \quad (6.14)$$

Let us denote these singularities as $s_{57}^{\pm}(s, t)$ (they can be written explicitly). Another relevant for us singularity is given by

$$s_{\max}(s, t) : \quad t_{26}^{-}(s, t, s_{57}) = \tilde{t}_{26}^{+}(s_{57}). \quad (6.15)$$

Next we analyze the last integral (6.13) as we increase $t$. We can follow the motion of singularities $s_{57}^{\pm}(s, t)$ which is depicted in figure 45. The pinch occurs on the real axis
<figure><img src="image_46.png" /><figcaption>**Figure 45.** Deformation of the $s_{57}$ contour relevant for the double spectral density computation. The original integration integral goes along the real axis with $s_{57} \ge 4m^2$ which is the condition that the invariant mass of a pair of on-shell particles should be greater than $4m^2$. As we increase $t$ to the positive values, the singularity of the integrand $s_{57}^{+}(s, t)$ enters the integration contour and drags it in the region $s_{57} < 4m^2$. The Landau curve develops when the integration contour is pinched between the two singularities of the integrand $s_{57}^{+}(s, t)$ and $s_{57}^{-}(s, t)$. As we further increase $t$ these singularities move into the complex plane, leading eventually to the contour depicted in figure 46.</figcaption></figure>

when

$$ \text{Landau curve: } \quad s_{57}^{+}(s,t) = s_{57}^{-}(s,t). \qquad (6.16) $$

One can check that this is indeed the correct location of the Landau curve for the acnode, namely it coincides with (6.6). Note that in terms of the original unitarity integral it involves a pinch in $s_{57}$ away from the physical values of this invariant $s_{57} \ge 4m^2$. In fact from the point of view of $T_{2\to 4} \sim T_{2\to 2} \times T_{1\to 2}$ it involves analytic continuation to the second sheet!

To compute $\rho(s, t)$ we need to understand what happens as we continue further. The singularities $s_{57}^{\pm}$ then recede in the complex plane. In the $t_{26}$-plane we take the discontinuity across $t_{26}^{-}(s, t, s_{57}) = t_{26}^{(2)}(s, t, s_{57})$ which is just logarithmic. As a result the relevant integral takes the form

$$ \begin{aligned} \mathcal{I}(s, t, s_{57}) &= \int_{t_{26}^{(2)}}^{t_{26}^{-}} dt_{26} \frac{1}{\sqrt{(t_{26} - \tilde{t}_{26}^{+})(t_{26} - \tilde{t}_{26}^{-})}} \frac{1}{\sqrt{(t_{26} - t_{26}^{+})(t_{26} - t_{26}^{-})}} \\ &= -2 \frac{(\mathcal{F}(\arcsin \chi_1^{1/2} |\chi_2| - \mathbf{K}(\chi_2)))}{\sqrt{(\tilde{t}_{26}^{+} - t_{26}^{-})(t_{26}^{+} - \tilde{t}_{26}^{-})}}, \end{aligned} \qquad (6.17) $$

$$ \chi_1 = \frac{(t_{26}^{(2)} - \tilde{t}_{26}^{-})(\tilde{t}_{26}^{-} - \tilde{t}_{26}^{+})}{(t_{26}^{(2)} - \tilde{t}_{26}^{+})(\tilde{t}_{26}^{-} - \tilde{t}_{26}^{-})}, \quad \chi_2 = \frac{(\tilde{t}_{26}^{-} - \tilde{t}_{26}^{-})(t_{26}^{+} - \tilde{t}_{26}^{+})}{(\tilde{t}_{26}^{-} - \tilde{t}_{26}^{+})(t_{26}^{+} - \tilde{t}_{26}^{-})}. \qquad (6.18) $$

In the expression above $\mathcal{F}(\phi|\chi)$ is the elliptic integral of the first kind, given by EllipticF[$\phi, \chi$] in Mathematica, and $\mathbf{K}(\chi)$ is the complete elliptic integral of the first kind, given by EllipticK[$\chi$] in Mathematica.

We then need to integrate $\mathcal{I}(s, t, s_{57})$ across the contour in the $s_{57}$ plane which is given in figure 46, and we obtain finally

$$ \rho_{\text{acnode}}(s,t) = \frac{1}{3(8\pi)^3} \frac{2\pi}{\sqrt{s(s-4)}} \left( \int_{s_{57}^{(-)}}^{s_{57}^{(+)}} \frac{ds_{57}}{i} \mathcal{I}(s,t,s_{57}) - 2 \int_{4m^2}^{s_{\text{max}}(s,t)} ds_{57} \text{Disc}_{s_{57}} \mathcal{I}(s,t,s_{57}) \right), \qquad (6.19) $$
<figure><img src="image_47.png" /><figcaption>Figure 46. The final integration contour for $\mathcal{I}(s, t, s_{57})$ in the $s_{57}$ complex plane which computes the double spectral density $\rho(s, t)$. The horizontal part is given by $\text{Disc}_{s_{57}}\mathcal{I}(s, t, s_{57})$.</figcaption></figure>

where $\text{Disc}_{s_{57}}\mathcal{I}(s, t, s_{57}) = 2\frac{K(\frac{1}{1-\chi_2})}{\sqrt{1-\chi_2}\sqrt{(\tilde{t}_{26}^+-t_{26})(t_{26}^+-\tilde{t}_{26})}}$. We attach a Mathematica notebook together with the arXiv submission that numerically computes $\rho_{\text{acnode}}(s, t)$ using the formula above.

There are several nontrivial tests for the result (6.19). First of all, it should be crossing-symmetric

$$ \rho_{\text{acnode}}(s, t) = \rho_{\text{acnode}}(t, s). \qquad (6.20) $$

This is highly non-obvious from (6.19), but we have checked numerically that it is indeed the case. Second, it should reproduce the single discontinuity (6.4) upon doing the dispersive integral (6.5), which we again have tested numerically and found perfect agreement. We plot the double spectral density in figure 47.

## 6.3 Nonperturbative lesson

What does the computation of the double discontinuity of the acnode graph teach us about the nonperturbative scattering amplitudes? The starting point in this case is the nonperturbative relation, which includes the three-particle amplitude for $s > 9m^2$ as in (6.1). We can then analytically continue the unitarity relation in $t$ to find the double discontinuity in terms of the analytically continued $T_{2\to3}$ amplitudes, as reviewed for example in [34]. In the case of the acnode the relevant singularity was due to one the one-particle poles $\frac{1}{t_{15}-m^2}$ and $\frac{1}{t_{16}-m^2}$. In the theory with the cubic coupling, these will still be present nonperturbatively. The residue of the pole however, instead of being just a constant will involve a nonperturbative $T_{2\to2}$ amplitude. Therefore, some of the formulas that rely on the detailed form $T_{2\to2}$ will change, however the most striking feature of the acnode double discontinuity calculation is that it probes $T_{2\to2}$ on the second sheet, see figure 45! This fact has implications for closing the system of equations in the iteration scheme.
<figure><img src="image_48.png" /><figcaption>Figure 47. Acnode double spectral density given in $768\pi^3 \times (6.19)$ evaluated for $(x, y) = (\frac{4m^2}{s}, \frac{4m^2}{t})$, such that $0 \le x, y \le 1$. The nontrivial support of the double spectral density is controlled by the Landau curve (6.6) and is localized in the region $s, t > 9m^2$, or, equivalently, $0 \le x, y \le \frac{4}{9}$.</figcaption></figure>

The nonperturbative version of the equation (6.19) schematically will presumably take the form

$$ \rho_{\text{MP}}(s,t) \stackrel{?}{=} \int ds_{57} ds_{67} K(s,t, s_{57}, s_{67}) T_{2 \to 2}^{(2)}(s_{57}) (T_{2 \to 2}^{(2)}(s_{67}))^*, \quad (6.21) $$

where $T_{2 \to 2}^{(2)}$ stands for the $2 \to 2$ scattering amplitude analytically continued to the second sheet through the elastic unitarity cut. This analytic continuation and the structure of the amplitude on the second sheet was recently discussed in [66]. Note that finding the amplitude on the second sheet continued through elastic unitarity cut might also be possible via an iteration process, see [66], therefore equations of the type (6.21) do provide some interesting multi-particle extension of the current $2 \to 2$ S-matrix bootstrap program.

## 6.4 Atkinson scattering-from-acnode

To explore the effect of $\rho_{\text{MP}}(s,t)$ on the amplitudes constructed in the previous sections, we next construct amplitude functions for which, we set

$$ \rho_{\text{MP}}(s,t) = c_a 768\pi^3 \rho_{\text{acnode}} \left( \frac{16s}{9}, \frac{16t}{9} \right), \quad (6.22) $$

where $c_a$ is a constant, $\rho_{\text{acnode}}(s,t)$ is given by (6.19), and we rescaled the arguments to make the acnode Landau curve asymptote to $16m^2$.

We therefore initialize our iteration algorithm with the following input $(\lambda, 0, (6.22))$ which generates an amplitude parameterized by $(\lambda, c_a)$. Let us present the results for such amplitudes in $d=4$. In essence, the effect is to introduce more inelasticity to our amplitudes, see figure 48.

We use the results for these amplitudes below in section 7 for the low energy observables and navigate the space of allowed theories.
<figure><img src="image_49.png" /><figcaption>**Figure 48.** $S_2$ absolute value for the 2QE+acnode $\rho_{\text{MP}}(s, t)$ scheme in $d = 4$, $\lambda = 3\pi$. Left panel: we see the expected growth of inelasticity as we increase $c_a$, given in (6.22). Right panel: the same plot, zoomed in to compare the 2QE iteration scheme with $\rho_{\text{MP}}(s, t) = 0$, versus $c_a = 1$.</figcaption></figure>

<figure><img src="image_50.png" /><figcaption>**Figure 49.** We plot the imaginary part of the amplitude (red), proportional to the total cross section, versus the inelastic cross section (green). The dashed curves correspond to $\lambda = 3\pi$, $d = 4$, 2QE amplitude. Other four curves are given by introducing the acnode $\rho_{\text{MP}}(s, t)$, as in (6.22), with $c_a = 1, 5, 10, 15$.</figcaption></figure>

# 7 Low energy observables

In this section, we look at the comparison of our approach to the current main approach, initiated in [41, 57, 67], and we compare with results in with the recent works [45, 46]. In these works, the authors investigated in particular some bootstrap bounds on some low energy observables, akin to nonperturbative QFT couplings.

These approaches are extremely powerful to derive bounds, but do not provide actual S-matrices whose various properties like the presence of LC or fine details of unitarity like elastic unitarity or Aks production can be tested.

Below, we consider the same low energy couplings and are able to make very precise statements about theories that live extremely close to the boundary. Following [45]
and [46], we thus look at two low-energy pairs of observables, $(\lambda, c_2)$ and $(\tau_{0,0}, \tau_{1,0})$. They are defined as:

$$ 
\begin{aligned}
c_2 &\equiv \frac{1}{4} \partial_s^2 T(s,t)|_{s=t=\frac{4m^2}{3}}, \\
\tau_{0,0} &= 2^{d-1} m^{d-4} f_0(2m^2), \\
\tau_{1,0} &= 2^{d-1} m^{d-2} \partial_s f_0(2m^2).
\end{aligned} 
\qquad (7.1) $$

Both $(\lambda, c_2)$ and $(\tau_{0,0}, \tau_{1,0})$ probe the amplitude at low energies and very similar almond-shape, two-sided bounds on the values of these parameters were derived in the aforementioned works, see [45, figure 3] and [46, figure 2].

While our methods allow us to resolve the fine structure near the boundary, they do not allow us to explore the whole parameter space of allowed theories discussed in these papers. Our amplitudes are located near the origin $(0,0)$ on both plots (on the scale set by the size of the almond/butterfly/leaf shape of allowed couplings). More precisely, we find that our amplitudes are located very close to the lower-bound on $c_2$ for given $\lambda$, or $\tau_{1,0}$ for given $\tau_{0,0}$.

We are thus interested in the following question: which amplitudes among the ones that we have constructed are more extremal? Here, by extremal we imagine amplitudes that minimize $c_2$ for given $\lambda$, or $\tau_{1,0}$ for given $\tau_{0,0}$. We address this question in two steps. First, we consider 2QE and 2PR amplitudes for which $\rho_{MP}(s,t) = 0$. We find that QE amplitudes are more extremal. Second, we turn on $\rho_{MP}(s,t)$ by taking as a proxy the double spectral density of the acnode graph. We find, that adding $\rho_{acnode}(s,t)$ moves the amplitude inside the allowed region, or away from extremality.

Looking at the plots and observing that the iterated amplitudes end up close to the boundary of the allowed region, it is tempting to ask: how close to the boundary are they? This question requires higher precision calculations, and thus a detailed and systematic analysis of our algorithm as we take the continuum limit. Leaving this for future work, we have chosen a particular 2QE amplitude for which the primal bootstrap data is available from [46], namely $\lambda = 2.012434211\pi$, $d = 4$

$$ \frac{c_{2,\min}^{\text{primal}}}{32\pi} \approx 0.00000414818. \qquad (7.2) $$

It would be also interesting to explore how close to this value the dual bounds on this observable [46, 68, 69] can be brought.

Constructing the corresponding 2QE amplitude for the same value of $\lambda$ we have gotten that

$$ \frac{c_2^{2\text{QE}}}{32\pi} \approx 0.0000041484 - 0.00000002 = 0.0000041482, \qquad (7.3) $$

where we indicated the contribution to $c_2$ coming from the single and double discontinuity of the amplitude. In particular, the contribution of particle production to $c_2$ is $\sim 1 \times 10^{-9}$.

To get this value we have explored how our results are affected by changing various grids and cutoffs. Therefore with the current precision we cannot definitively answer the question whether the 2QE amplitude is truly extremal.

Also let us compare the results above with the two-loop result in $\phi^4$, see appendix F,

$$ \frac{c_2^{\text{two-loop}}}{32\pi} \approx 0.000041475. \qquad (7.4) $$
It would be interesting to do the three-loop computation of $c_2$ since it is expected to be even closer to the numbers quoted above.

More conceptually, to understand if $T_{2QE}(s, t)$ is extremal, we would need to study $\frac{\delta c_2}{\delta \eta_{MP}(s)}$ and $\frac{\delta c_2}{\delta \rho_{MP}(s,t)}$ for infinitesimal deformations in the functional space of multi-particle data. More generally, we can write down two extra equations for fixed $\lambda$

$$ \text{Extremal } c_2: \quad \frac{\delta c_2}{\delta \eta_{MP}(s)} = 0, \quad \frac{\delta c_2}{\delta \rho_{MP}(s,t)} = 0, \qquad (7.5) $$

and we try to use them to solve for $(\eta_{MP}(s), \rho_{MP}(s, t))$. In general we observed that turning on inelasticity tends to move the amplitude inside the allowed region. This motivates the expectation that setting $\eta_{MP} = 0$ is extremal. We do not have a similar argument for $\rho_{MP}(s, t)$. Recall that particle production in the spin $J > 0$ sector $\eta_J(s) \equiv 1 - |S_J|^2$ is given by the Froissart-Gribov formula

$$ \eta_J(s) = \frac{(s - 4m^2)^{\frac{d-3}{2}}}{\sqrt{s}} \frac{4N_d}{\pi} \int_{z_1}^{\infty} dz (z^2 - 1)^{\frac{d-4}{2}} Q_J^{(d)}(z) (\rho_{el}(t(z), s) + \rho_{MP}(s, t(z))) . \quad (7.6) $$

Minimizing production in $J > 0$ sector is not the same as setting $\rho_{MP}(s, t) = 0$. In particular, to make contact with amplitudes constructed in [45, 46], it would be interesting to try iterating $\rho_{MP}(s, t)$ in a way that cancels as much particle production as possible in some finite energy (and spin) range. We do not explore this possibility in the present paper. A fixed point of such a procedure would be an extremal amplitude. We leave exploration of this possibility for future work.

## 7.1 ($\tau_{0,0}$, $\tau_{0,1}$) plots in $d=3, d=4$

Here we provide the plots for the $(\tau_{0,0}, \tau_{0,1})$ observables in $d=4$, see figure 50 and figure 51, and in $d=3$, see figure 52 and figure 53.

<figure><img src="image_51.png" /><figcaption>Figure 50. We plot the $(\tau_{0,0}, \tau_{1,0})$ plane in $d=4$. It was observed in [45] that the space of allowed values is very well approximated by the one-loop perturbative result. We see that on this plot all our amplitudes: toy-model amplitudes with zero double discontinuity, 2QE and 2PR all lie on the boundary of the allowed space and are not resolvable.</figcaption></figure>
<figure><img src="image_52.png" /><figcaption>**Figure 51.** The same plot as above still in $d = 4$, only now we zoomed on the point close to $\tau_{0,0} \simeq 1.5$. On this plot we still cannot resolve the 2QE and 2PR amplitudes, but we see that the toy model amplitude, as well as the amplitudes with $\rho_{\text{acnode}}(s, t)$ turned on have higher values of $\tau_{1,0}$, and are thus less extremal.</figcaption></figure>

<figure><img src="image_53.png" /><figcaption>**Figure 52.** We plot the $(\tau_{0,0}, \tau_{1,0})$ plane in $d = 3$. It was observed in [45] that the space of allowed values is well-approximated by the one-loop result. We see that on this plot all our amplitudes: toy-model amplitudes with zero double discontinuity, 2QE and 2PR all lie on the boundary of the allowed space and are not resolvable.</figcaption></figure>
<figure><img src="image_54.png" /><figcaption>**Figure 53.** The same plot as above in $d = 3$, only now we zoomed on the point close to $\tau_{0,0} \simeq -4.67$. On this plot we still cannot resolve the 2QE and 2PR amplitudes very well, but we see that the toy model amplitude, as well as the amplitudes with $\rho_{\text{acnode}}(s, t)$ turned on have higher values of $\tau_{1,0}$, and are thus less extremal.</figcaption></figure>

## 7.2 $c_2$ plots

Here we provide the plots for the $c_2$ observables in $d = 4$, figure 54, 55 and in $d = 3$ in figure 56, 57.

<figure><img src="image_55.png" /><figcaption>**Figure 54.** $(\frac{\lambda}{32\pi}, \frac{c_2}{32\pi})$ plane in $d = 4$. The space of allowed amplitudes, see [45, 46], (to a good approximation) lies above the one-loop result in $\phi^4$ theory (dashed line). Various amplitudes constructed in this paper all lie very close to the one-loop result. In particular, on this scale the toy-model, 2QE and 2PR amplitudes are indistinguishable.</figcaption></figure>
<figure><img src="image_56.png" /><figcaption>**Figure 55.** $(\frac{\lambda}{32\pi}, \frac{c_2}{32\pi})$ plane in $d=4$, zoomed around $\lambda = 3\pi$. In addition to the amplitudes on the previous plot we add the data points generated by turning on $\rho_{MP}$ proportional to the double spectral density of the acnode graph (rescaled in a way to have the Landau curve asymptotes to $16m^2$). We see that among the amplitudes constructed in this paper the 2QE amplitude has the minimal $c_2$ for a given $\lambda$.</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>**Figure 56.** $(\frac{\lambda}{32\pi}, \frac{c_2}{32\pi})$ plane in $d=3$. The space of allowed amplitudes, see [45], (to a good approximation) lies above the one-loop result in $\phi^4$ theory (dashed line). Various amplitudes constructed in this paper all lie very close to the one-loop result. In particular, on this scale the toy-model, 2QE and 2PR amplitudes are indistinguishable.</figcaption></figure>
<figure><img src="image_58.png" /><figcaption>Figure 57. $(\frac{\lambda}{32\pi}, \frac{c_2}{32\pi})$ plane in $d=3$, zoomed around $\lambda = -3\pi$.</figcaption></figure>

## 7.3 Grid dependence

One could ask that since our result seem very precise, what is the convergence with the number of points and can one add error bars? In this paper, as was discussed above, the main error we try to quantify is to what precision unitarity is satisfied.

With the perspective of producing an actual model for the pion-S-matrix, an actual study of the errors needs to be performed. Here, we display the result of varying the $\rho(s, t)$-grid. We investigate the effect of two different cutoffs, $10^{-10}$ and $10^{-12}$, and for these we vary the number of points.

<figure><img src="image_59.png" /><figcaption>Figure 58. Here we plot the dependence of $c_2/32\pi$ for $\lambda = 2.012434211\pi$ in $d=4$ on the number of points $N$ in the grid used for the single discontinuity $\rho(s)$. We see that a simple power-like fit agrees with the data very well. The power of $N$ agrees with the expectation value for the error made when approximating an integral with the trapezoidal rule (which is what linear interpolation used in the paper effectively does). As we varied the grid for the single discontinuity we kept the grid for the double discontinuity fixed.</figcaption></figure>
<figure><img src="image_60.png" /><figcaption>**Figure 59.** Dependence of $c_2/(32\pi)$ on the number of points $N_{\text{grid}}$ in the double spectral density $\rho(s, t)$ two-dimensional grid. The number of points has been generated by our meshing procedure taking into account some density profile function, see appendix B. These grids correspond to $\text{nx} = 70, 90, 100, 110$. The two different cutoffs in the Regge limit $10^{-12}$ and $10^{-10}$ are the points closest to zero after which we connect to 0 with a linear segment, which means that we force a $1/s$ or $1/t$ fall-off at infinity. As discussed in the text, this does not pose a problem in the numerical integrals because they are dispersive, hence have a $1/(s' - s)$ kernel that renders the difference between a linear fall-off or a true logarithmic fall-off very small for such high cut-offs.</figcaption></figure>

# 8 Mandelstam representation and quantum gravity

It is interesting to ask if methods developed in this paper can be applied to quantum gravity. In this section, we will argue that one important element for our approach, the Mandelstam representation, should not hold in quantum gravity.

The starting point of our analysis is the Mandelstam representation which is based on two a priori independent assumptions about the amplitude:

* maximal analyticity

* polynomial boundedness in $s$ for any $t$ (and vice versa)

Both assumptions are highly non-obvious and have not been proven neither in the context of nonperturbative QFT, nor for gravitational theories. It is interesting to ask a simpler question: is the Mandelstam representation consistent with other properties expected from gravitational amplitudes? Note that the Mandelstam representation trivially does not hold in perturbative string theory [70–72], where the polynomial boundedness assumption breaks down: for any $N$, there exists $t > 0$ such that the Regge growth $s^{\alpha'}t > s^N$.

In this section we review the old argument of Cerulus and Martin [49] that relates polynomial boundedness of the amplitude for unphysical values of $t$ to certain properties of the amplitudes in the Regge limit and high energy scattering at fixed angles.<fnref n="27" /> While the original argument was done for gapped theories, it is a straightforward exercise to relax this assumption.

<footnote n="27">²⁷Recently, this argument was generalized to rely on the axiomatic QFT analyticity only [73]. The resulting bound however is very weak.</footnote>
Consider the scattering amplitude as a function of fixed angle $z \equiv \cos\theta$, namely $\mathcal{T}(s, z) \equiv T(s, t = -\frac{s}{2}(1-z))$. We fix $s$ to be real and positive,<fnref n="28" /> and we consider $\mathcal{T}(s, z)$ in the complex $z$-plane. Maximal analyticity implies that $\mathcal{T}(s, z)$ is analytic, modulo the two cuts $z \in (-\infty, -1] \cup [1, \infty)$ which correspond to scattering in the $u$- and $t$-channel correspondingly.

We take three real $z$'s such that $0 < z_1 < z_2 < z_3 < 1$ and we map the $z$-plane to the $\tau$-plane. We first transform

$$w(z) = \frac{1}{z} \left( 1 - \sqrt{1 - z^2} \right), \qquad (8.1)$$

which maps the cut $z$-plane inside the unit circle in the $w$-plane. We then consider the following map

$$\tau(w) = \frac{1}{w(z_1)} \left( w + \sqrt{w^2 - w(z_1)^2} \right). \qquad (8.2)$$

This mapping maps the region $-z_1 \le z \le z_1$ in the $z$-plane to the unit circle $|\tau| = 1$.

Consider now another pair of circles in the $\tau$-plane of radii $r_2 = \tau(w(z_2))$ and $r_3 = \tau(w(z_3))$. In the original $z$-plane the circles map in the oval shape region, see figure 60.

By assumption, maximal analyticity implies that the scattering amplitude is analytic in the annulus $1 \le |\tau| \le r_3$. Let us introduce the maximal value of the amplitude on a given circle

$$M_r \equiv \max_{|\tau|=r} |\mathcal{T}(s, z)|. \qquad (8.3)$$

We then apply the Hadamard three-circle theorem, see e.g. chapter 23 in [74], that states that for $1 < r_2 < r_3$ we have the following inequality

$$M_{r_2} \le M_1^{1-\frac{\log r_2}{\log r_3}} M_{r_3}^{\frac{\log r_2}{\log r_3}}. \qquad (8.4)$$

This constraint becomes particularly interesting if we choose $z_{2,3} = 1 + \frac{2t_{2,3}}{s}$, where $t_2 < t_3 < 0$ are fixed momenta and we take $s \to \infty$.

Let us introduce, following Cerulus and Martin, the upper bound on the fixed angle scattering amplitude at high energies

$$|\mathcal{T}(s, z)| \le e^{-\phi(s)}, \quad -a \le z \le a < 1. \qquad (8.5)$$

We can then rewrite (8.4) as follows

$$|T(s, t_2)| \le e^{-C(a)\phi(s)\frac{\sqrt{-t_2-\sqrt{-t_3}}}{\sqrt{s}}} \max_{|\tau|=r(s,t_3)} |\mathcal{T}(s, z)|, \quad s \to \infty, \qquad (8.6)$$

where $C(a) > 0$ and we used that $M_{r_2} \ge |T(s, t_2)|$. This equation bounds the Regge limit of the amplitude, the l.h.s., in terms of the fixed angle scattering, and the value of the amplitude in certain sub-domain of the $z$-plane, the r.h.s.

The validity of the Mandelstam representation implies that there should exist an integer $N$, such that

$$\max_{|\tau|=r(s,t_3)} |\mathcal{T}(s, z)| \le |s|^N, \quad |s| \to \infty. \qquad (8.7)$$

<footnote n="28">²⁸We approach the real axis from above as usual.</footnote>
<figure><img src="image_61.png" /><figcaption>**Figure 60.** Here we depict regions in the z-plane which map to the concentric circles in the τ-plane. Here we consider $z_1 = \frac{1}{4}$ so that the region (magenta) $-\frac{1}{4} \le z \le \frac{1}{4}$ is mapped to the unit circle in the τ-plane. The blue region corresponds to $z_2 = 0.9$ and the red region comes from $z_3 = 0.99$ in the argument.</figcaption></figure>

Assuming that $|T(s, t_2)|$ is polynomially bounded from below (it cannot decay too fast), we then get the Cerulus-Martin bound on the fixed angle scattering for amplitudes that admit Mandelstam representation

$$ \phi(s) \le c_0 \sqrt{s} \log s, \quad s \to \infty. \qquad (8.8) $$

In other words, maximal analyticity implies the following schematic Cerulus-Martin relation

$$ \text{Regge} \le \text{Fixed angle} \times \text{Mandelstam}. \qquad (8.9) $$
The Regge limit for gravitational amplitudes is controlled by the large impact parameter scattering and its absolute value behaves polynomially in $s$, see e.g. [71, 75–77]. The fixed angle scattering, on the other hand, is believed to be entropically suppressed, see e.g. [71, 78–80],

$$ \phi_{\text{QG}} \sim (\sqrt{s})^{1+\frac{1}{d-3}}, \qquad (8.10) $$

which clearly violates (8.8).

Our conclusion is that the *expected* properties of the scattering amplitude in gravitational theories for physical $t$ (polynomial behavior in the Regge limit and exponentially faster than $\sqrt{s}$ decay for fixed angle scattering) are *not* compatible with polynomial boundedness needed for the Mandelstam representation.

This does not mean that a version of our approach cannot be used to the question of unitarizing gravity or string theory, but simply that the basic iteration presented in the present paper will not work. Another method should be found, or another regime studied, for instance the eikonal one. We believe that the general idea initiated in this paper should be eventually applicable to gravity as well.

Let us also mention, that in confining gauge theories, e.g. in QCD, we get a power-law decay at high energies for fixed angle scattering [81, 82], therefore we expect that it should be possible to use the Mandelstam representation to accurately model pion scattering.

# 9 Conclusions and open directions

We now briefly summarize the results of the paper, emphasizing some of the physics we learned on the way. We then move to listing a few open directions.

## 9.1 Summary

In this paper we have explored models of nonperturbative scattering constructed by iterations of the Mandelstam representation and unitarity. The scattering amplitudes are generated from some input data $(\lambda, \eta_{\text{MP}}(s), \rho_{\text{MP}}(s, t))$:

* the subtraction constant in the Mandelstam representation (1.1), $\lambda \equiv T(\frac{4m^2}{3}, \frac{4m^2}{3})$;

* the S-wave ($J=0$) inelasticity $1 \ge \eta_{\text{MP}}(s) \ge 0$ defined in (1.10);

* the multi-particle double spectral density $\rho_{\text{MP}}(s, t)$ defined in (1.9).

As an output we get amplitude functions that satisfy: maximal analyticity, crossing, elastic and inelastic unitarity (with the level of precision controlled by the numerical algorithm that we used). The amplitude functions that we get are also UV-complete in the sense that we have control over their behavior at all energies and scattering angles.

The amplitudes that we obtained exhibit interesting physical properties: Landau curves, particle production, nontrivial Regge and fixed angle behavior, nonperturbative scaling of various quantities with the “coupling” $\lambda$. In the present paper we focused mostly on the two classes of amplitudes: the two-particle quasi-elastic scattering amplitudes (2QE) initialized by $(\lambda, 0, 0)$; and the two-particle reducible amplitudes (2PR) given
by $(\lambda, (2.24), 0)$ for which inelasticity is analytic in spin all the way to $J = 0$. The amplitudes obtained in this way are close to saturating certain minimal coupling bounds derived by other methods [45, 46].

We now list a few notable facts that we learned from studying these amplitudes.

**Dispersiveness and sign of $\phi^4$ coupling.** Interestingly, we have found that for the amplitudes constructed in the paper the subtraction constant $\lambda$ is actually *dispersive*. We call a quantity dispersive when it is given by a dispersion relation integral of a discontinuity of the amplitude with some kernel. For the quartic coupling $\lambda$, this comes out as a surprise because a naive Cauchy argument $\lambda = \int \frac{ds'}{2\pi i} \frac{T(s', \frac{4m^2}{3})}{s' - \frac{4m^2}{3}}$ fails. Indeed, as we deform the contour to express $\lambda$ through the discontinuity of the amplitude, we find that we cannot drop the contour at infinity because $T(s, \frac{4m^2}{3}) \to \text{const}$. The way $\lambda$ turned out to be dispersive is more subtle, and comes through S-wave unitarity, where $\lambda$ contributes to $\text{Re } f_0(s)$.

In three dimensions, elastic unitarity together with the asymptotic $\text{Im } f_0(4m^2) = 0$ of our amplitudes imply that $\text{Re } f_0(4m^2) = 0$.<fnref n="29" /> In four dimensions, our amplitudes satisfy $\text{Im } f_0(\infty) = 0$, which through inelastic unitarity implies that $\text{Re } f_0(\infty) = 0$. Both $\text{Re } f_0(4m^2) = 0$ in $d=3$, and $\text{Re } f_0(\infty) = 0$ in $d=4$ provide an extra equation which expresses $\lambda$ through the discontinuity of the amplitude. In the toy model, these were the sum rules (3.17) and (3.31) correspondingly.

In the context of the toy models, for which $\rho(s,t) = 0$ and we have from unitarity $\rho(s) \ge 0$, these sum rules immediately fix the sign of $\lambda$ to the one used in the paper ($\lambda < 0$ in $d=3$, $\lambda > 0$ in $d=4$). For the opposite signs of $\lambda$ we thus conclude that the double discontinuity contribution to the sum rules (and to the amplitude itself) should be essential.

Note also that the 4d amplitudes constructed in the present paper correspond to the “wrong” sign of the potential $V(\phi) = -\frac{\lambda}{4!}\phi^4$ for which we do not expect to have a stable vacuum. Of course, the construction of the present paper concerns itself only with a particular observable, namely the $2 \to 2$ scattering amplitude, and it does not in any sense define a theory. Moreover, we can imagine a theory where negative $\phi^4$ potential emerges as a leading low-energy approximation of a healthy UV complete theory.<fnref n="30" /> Still, if we are to insist on having “pure” $\phi^4$ theory with an unbounded potential we do not expect to have a consistent unitarity S-matrix. At the level of the amplitudes constructed in the paper we did not see any inconsistencies. It would be however very interesting to understand how vacuum instability can be detected in the S-matrix language.<fnref n="31" />

**Non-trivial Regge behavior.** The amplitudes constructed in the paper exhibit a non-trivial Regge behavior. For toy models we could solve for the Regge limit analytically at a price of setting $\rho(s,t) = 0$. We then observed that the Regge limit gets modified as we consider the full back-reacted solution with $\rho(s,t) \ne 0$ for which we do not have an analytic control. It would be interesting to develop better understanding of the Regge behavior of the amplitudes constructed in the present paper.

<footnote n="29">²⁹This was argued to be true in the full $d=3 \phi^4$ theory in [58, 59].</footnote>

<footnote n="30">³⁰We thank Joan Elias Miro for discussions on this point.</footnote>

<footnote n="31">³¹We thank Nima Arkani-Hamed for discussions on this point.</footnote>
**Transition between Regge / fixed-angle regimes.** Continuing with the high-energy limits of the amplitudes, we also observed a very clear transition from Regge to the fixed-angle regime. We recall that the Regge limit corresponds to fixed $t$, large $s$, while the fixed-angle is large $s$ and large $-t$, with $s/t$ fixed. Since the scattering angle is defined by $\cos(\theta) = 1 - \frac{2t}{s-4m^2}$, these two regions correspond to $\theta \ll 1$ and $\theta \sim O(1)$ correspondingly, and in theory as $|t|$ is increased, or equivalently $\theta$, scattering amplitudes should transition from Regge to fixed-angle behavior.

We observed a very clean transition, both in $d = 3$ and $d = 4$ for

$$ \frac{s}{m^2} \sim \theta^{-2} \qquad (9.1) $$

which is illustrated in figure 16, figure 22. This regime is reached when $t/m^2 \sim O(1)$, which in impact parameter space corresponds to $b \sim 1/m$.<fnref n="32" />

## 9.2 Open directions

There are many future important directions to explore, we list some of them below:

**Are the two-particle quasi-elastic (2QE) amplitudes physical?** In this paper we have argued that there exist a natural one-parameter family of amplitudes defined by the triad $(\lambda, 0, 0)$, with no bound states below the first threshold. We have also observed that these amplitudes are located very close to the bounds on couplings derived by other methods [45, 46]. Apart from the question of whether the 2QE amplitudes are truly extremal on which we comment below, one obvious question is: are they clearly unphysical? Recall that they are defined by the requirement $\eta_{MP}(s) = 0$, which means that in their graphic, $\phi^4$-like expansion, the $J = 0$ part of any graph with a multiparticle-cut is projected out (the simplest of which being the Aks graph and the frog graph). Naively this implies a certain level of non-locality to be present in the theory. However, at the level of the constraints imposed in this paper there is no problem in doing so. In other words, this non-locality is compatible with the constraints imposed on the scattering amplitude in the present paper. We expect that the 2QE amplitudes could be excluded using multi-particle unitarity and it would be very interesting to show this explicitly.

An obvious extension of the present work is to consider 2QE amplitudes in the presence of bound states. For example, it would be interesting to see if the 2QE amplitudes with a bound state at the two-particle threshold, as first observed in [41], are close to saturating the bound on the maximal coupling found in [45, 46]. A preliminary attempt of ours to implement the fixed-point algorithm in this case appears to diverge. We observed some similar phenomena in $d = 2$ before [26], and there we restored convergence by implementing the Newton-Raphson method. It would be interesting to try this here as well.

Turning to gravity, it would be very interesting to construct quasi-elastic amplitudes in this case as well. These tentative amplitudes would be both crossing-symmetric and eikonalize correctly.

<footnote n="32">³²This relation holds for our amplitudes because they have a simple Regge limit, $T(s, t) \to \text{const}$ which leaves the saddle-point of the impact-parameter Fourier transform unchanged, $T(s, b) \sim \int d^2b e^{-i\vec{q}b} T(s, -q^2)$. For reggeized string-like particles for which $T(s, t) \sim s^{-\alpha'}t$, $t \sim m^2$ corresponds to $b \sim \log(s)$ which is the famous logarithmic transverse growth of strings at high-energies.</footnote>
**Nonperturbative unitarization via two-particle reducible resummation?** An-
other interesting case that we considered is given by the following choice of the iteration
data $(\lambda, (2.24), 0)$. In this scheme, $J = 0$ inelasticity is fixed using the Froissart-Gribov
formula.

This corresponds to re-summation of the two-particle reducible graphs in $\phi^4$ theory which are generated by iterations of two-particle unitarity and crossing, see graphs in black in figure 1. These graphs are the only graphs in $\phi^4$ theory at one and two loops, but starting from three loops there are two-particle irreducible graphs, namely the frog and the open envelope which have $\eta_{\text{MP}}(s) \neq 0$ and $\rho_{\text{MP}}(s, t) \neq 0$. As we proceed to higher loops, we generate a web of planar diagrams made of bubbles, triangles, and various connected ladders of those.

Since the planar sector of QFTs is usually expected to be summable and for $\phi^4$ is
even known to have a finite radius of convergence [52], one can already suspect that a
perturbative scheme could work in our set-up. As a matter of fact, a simple diagrammatic
counting given in section 2 and with explicit examples up to $n = 10$ in table 2, indicates
that the number of such planar, two-particle reducible graphs grows indeed like a power,
which, by standard reasoning suggests a finite radius of convergence. Indeed, this is what
we also find numerically, since we observe convergence of that scheme for finitely many
couplings, summarized in table 1.

An interesting result of our analysis is that such a truncation of the Feynman series is
self-consistent at the level of the two-to-two scattering amplitude in the following sense: it
generates the amplitude function $T(s,t)$ that satisfies both elastic and inelastic unitarity,
as well as crossing. It would be very interesting to study this unitarization scheme in other
models. It would be also interesting to understand if similar iteration techniques can be
applied to re-sum families of nPR planar graphs, e.g. melonic graphs of which the frog
diagram is the first example.

**Are extremal amplitudes physical?** An interesting problem is to find amplitudes
that satisfy the basic constraints studied in the paper and minimize/maximize certain
couplings. In section 7, we have found that all the amplitudes constructed in this paper lie
close to saturating the bounds on the low-energy observables derived in [45, 46]. The most
extremal amplitudes among the ones constructed in this paper are the 2QE amplitudes
initialized with zero multi-particle input $\eta_{\text{MP}}(s) = \rho_{\text{MP}}(s,t) = 0$. Understanding if the
QE amplitudes are truly extremal would require a search in the space of multi-particle
inelasticities encoded by $(\eta_{\text{MP}}(s), \rho_{\text{MP}}(s,t))$. We see the following possibility to perform
this search. We can start with one of the nontrivial solutions found in this paper, and then
explore infinitesimal variations $\delta\eta_{\text{MP}}(s)$ and $\delta\rho_{\text{MP}}(s,t)$ while choosing a direction in which
the coupling decreases (if there is one). At the extremum, we expect to find (7.5).

Finally, the works [41, 45, 46] observed inelasticity very different from ours. In that
case it seems to be generically pushed to the UV. For us, inelasticity comes at low/medium
energy due to graphs with simple multi-particle cuts. It would be interesting to understand
if our approach can also push inelasticity to the UV and see what constraints this poses
on $\rho_{\text{MP}}(s,t)$. Conversely, we can also imagine using the results of the present paper as an
input to the primal S-matrix bootstrap, see e.g. [83] for a related discussion.
**Is inclusion of multi-particle unitarity feasible?** The existing S-matrix bootstrap methods focus on the $2 \to 2$ scattering. In the same spirit, in the present paper, the multi-particle information, contained in $\eta_{MP}(s)$ and $\rho_{MP}(s, t)$, was used as an input. It is however interesting to notice that our approach has a natural extension in which some of the multi-particle information enters the iteration algorithm. The simplest one concerns the nontrivial support of $\rho_{MP}(s, t)$ which is very easy to impose in our approach, that is to say the choice of the Landau-curve above which $\rho_{MP}(s, t)$ is non-vanishing, see [34].

For more complicated constraints, consider, for example, the open envelope diagram in figure 61. As reviewed for example in [34], in the nonperturbative theory this diagram leads to an equation which schematically takes the following form

$$ \rho_{MP}(s,t) \stackrel{?}{=} \int K_4(T_{2 \to 2})^2 (T_{2 \to 2}^*)^2, \qquad (9.2) $$

which holds for the full nonperturbative amplitude in some finite region of the Mandelstam plane. It would be important to try implementing (9.2) and see how computationally costly it is.

In this paper we analyzed the simplest example of this type of relation in section 6, where we derived the double spectral density of the acnode graph. An interesting outcome of this analysis is that the two-to-two scattering amplitudes that enter into the relation analogous to (9.2) turned out to be evaluated on the second sheet reached through the elastic unitarity cut. As recently discussed in [66], the amplitude there in principle can also be found by iterations very similar to the one considered in the present paper. It would be very interesting to derive the nonperturbative relation of the type (9.2) explicitly (as well as the one for the acnode graph) and to implement them numerically. As for $\eta_{MP}(s)$, there are similar relations which require further exploration. For example cutting the frog graph, we get $T_{2 \to 4}$ diagram which contains, schematically, a pole $\frac{T_{2 \to 2}^2}{s_{123}-m^2}$, where $s_{123}$ is an invariant mass of the three produced particles. We can imagine adding this contribution to our iteration scheme for $\eta_{MP}(s)$.

Lastly, for $\eta_{MP}(s)$, we can imagine taking experimental data as an input given that it has a very simple physical meaning, the situation is more complicated for $\rho_{MP}(s, t)$ whose moments are related to particle production in partial waves with spin via the Froissart-Gribov formula. It would be also interesting to explore if semi-classical methods for multi-particle production, see e.g. [84–86], could shed some light on multi-particle input of the iterative algorithm.

**Iteration methods and various physical models.** Given that in our algorithm physics is naturally encoded in the multi-particle data, it would be very helpful to understand how various familiar physical models are encoded in the $(\lambda, \eta_{MP}(s), \rho_{MP}(s, t))$: amplitudes that saturate the Froissart bound [62, 87], pion scattering [88], scattering in Chern-Simons matter theories [89], stringy amplitudes [72], gravitational amplitudes.<fnref n="33" /> It would be then interesting to generalize the methods applied in this paper to these theories.

<footnote n="33">³³For the case of gravitational amplitudes, an interesting question is: do quantum gravity amplitudes satisfy the Mandelstam representation? Our conclusion in section 8 based on the lore about the high-energy fixed angle behavior of scattering amplitude in gravity is that the answer to this question is “No”.</footnote>
**Taking the continuum limit.** The algorithm we presented in this paper is based on discretizing the space of Mandelstam variables $4m^2 \le s, t < \infty$ and then interpolating between the grid points. To make sure that the amplitudes of interest exist we would like to argue that the amplitudes that we have constructed numerically admit a continuum limit. We have considered two approaches to this problem, which we have explained in the text. First, we have increased the grid resolution and checked that the results do not change in any significant way. Relatedly, we increased the grid cutoffs and checked that the amplitudes that we interpolate smoothly to the IR in $d=3$ and to the UV in $d=4$. Secondly, we looked into convergence of the iteration algorithm in the space of continuous functions directly. In the latter case we showed convergence for the toy model discussed in section 3 in $3 < d < 4$. Extending the existing functional analysis proofs to $d=3$ and $d=4$ requires developing more sophisticated methods which go beyond what has been done in the literature or in the present paper.

**Increasing the convergence range.** One feature of the amplitudes presented in this paper is that both particle production and the coupling constant, as measured by $\lambda$, are relatively small. The amplitudes that we obtained are also weakly coupled at high energies in the sense that partial waves $S_J(s) \to 1$ as $s \to \infty$ (but $T(s,t) \to \text{const}$ for fixed $t$ and $s \to \infty$). It would be very interesting to extend the method of the present paper to larger couplings and inelasticities.

To approach this class of amplitudes, we need to generalize the simple fixed point method developed in the present paper. In general, there should be nothing fundamental about the *iterative solution* to these functional equations. Any method numerical functional analysis that can work is equally good.

For instance, in $d=2$ [26], we observed that gradient methods such as Newton-Raphson allowed us to extend the range of convergence of the iterative algorithm, and to describe new-topologies of amplitudes (with extra zeros or poles), sensitive to the starting point of the algorithm. We therefore expect that such methods should also improve the performance of our algorithm in $d=3,4$. One immediate bottleneck relates to the discretization procedure of our numerical integrals. One-dimensional piecewise-linear integrals are not really integrals, they can be performed analytically, because on each segment, they amount to integrating a combination of a constant and a linear function times an integration kernel: all these integrals can be done beforehand and the whole integration of an interpolating function $\rho(x)$, $\int_0^1 K(x', x)\rho(x')dx'$ becomes the action of a matrix on the vector $\rho(0), \rho(x_1), \dots, \rho(x_{n-1}), \rho(1)$ (for more details, see [26, 57]). For two-dimensional integrals, our grids intersect the integration domains in some non-trivial way and a similar procedure is not immediately applicable. This problem should be solved, in order for a gradient descent or Newton-Raphson method to be applied in the present case of $d=3,4$. Lastly, it would be very interesting to explore the possibility of using modern machine learning methods in the present context [90].
Acknowledgments

We thank Nima Arkani-Hamed, David Atkinson, Pierre Aubert, Lucia Cordova, Miguel Correia, Joan Elias Miro, Giulia Isabella, Hofie Hannesdottir, Kelian Häring, Enrico Hermann, Andreas Juettner, Denis Karateev, Madalena Lemos, Andrew McLeod, Sebastian Mizera, Julio Parra-Martinez, Joao Penedones, Jiaxin Qiao, Balt van Rees, Matthew Schwartz, Amit Sever, Jaroslav Trnka, Pierre Vanhove, Pedro Vieira, Matt Walters, Xiang Zhang for useful discussions. We thank the authors of [45, 46] for sharing and explaining their detailed results to us. We thank Dmitrii Vorobev and Dmitri Vasiliev for exploring the d-dimensional formulae as a part of the Physics Practicum 2022. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement number 949077) and from Agence Nationale de la Recherche (ANR), project ANR-22-CE31-0017. This research was supported in part by the National Science Foundation under Grant No. NSF PHY-1748958.

A Definitions

**Partial wave expansion.** Here we summarize our conventions for the partial wave expansion following [38]:

$$T(s,t) = \sum_{J=0}^{\infty} n_J^{(d)} f_J(s) P_J^{(d)}(\cos \theta), \qquad (A.1)$$

$$P_J^{(d)}(z) = {}_2F_1\left(-J, J+d-3, \frac{d-2}{2}, \frac{1-z}{2}\right).$$

The partial wave coefficients can be found from the amplitude using the orthogonality relation

$$\frac{1}{2} \int_{-1}^{1} dz (1-z^2)^{\frac{d-4}{2}} P_J^{(d)}(z) P_J^{(d)}(z) = \frac{\delta_{J\bar{J}}}{N_d n_J^{(d)}}. \qquad (A.2)$$

The normalization coefficients are taken to be

$$N_d = \frac{(16\pi)^{\frac{2-d}{2}}}{\Gamma(\frac{d-2}{2})}, \qquad n_J^{(d)} = \frac{(4\pi)^{\frac{d}{2}}(d+2J-3)\Gamma(d+J-3)}{\pi\Gamma(\frac{d-2}{2})\Gamma(J+1)}. \qquad (A.3)$$

Using these formulas we get for scattering of identical particles

$$f_J(s) = \frac{N_d}{2} \int_{-1}^{1} dz (1-z^2)^{\frac{d-4}{2}} P_J^{(d)}(z) T(s,t(z)), \quad t(z) = -\frac{s-4m^2}{2}(1-z). \qquad (A.4)$$

The Q-functions that enter the Froissart-Gribov formula

$$f_J(s) = \frac{2N_d}{\pi} \int_{z_1}^{\infty} dz (z^2 - 1)^{\frac{d-4}{2}} Q_J^{(d)}(z) T_t(s, t(z)), \quad z_1 = 1 + \frac{4m^2}{s-4m^2}, \qquad (A.5)$$

take the following form

$$Q_J^{(d)}(z) = \frac{\sqrt{\pi}\Gamma(J+1)\Gamma(\frac{d-2}{2})}{2^{J+1}\Gamma(J+\frac{d-1}{2})} \frac{1}{z^{J+d-3}} {}_2F_1\left(\frac{J+d-3}{2}, \frac{J+d-2}{2}, J+\frac{d-1}{2}, \frac{1}{z^2}\right). \qquad (A.6)$$
**Spin zero partial wave in d dimensions.** The $J = 0$ partial wave for an amplitude given by a once-subtracted Mandelstam representation (1.1), in general $d$, can be computed to be:

$$f_0(s) = \frac{(16\pi)^{\frac{3-d}{2}}}{8\Gamma(\frac{d-1}{2})} \left( \lambda + \int_{4m^2}^{\infty} \frac{ds'}{\pi} \frac{(s-s_0)\rho(s')}{(s'-s)(s'-s_0)} \right. \\ + 2 \int_{4m^2}^{\infty} \frac{dt'}{\pi} \rho(t') \left[ \frac{1}{t_0-t'} + \frac{1}{t'} {}_2F_1 \left( 1, \frac{d-2}{2}, d-2, \frac{s-4m^2}{t'} \right) \right] \\ + 2 \int_{4m^2}^{\infty} \frac{ds' dt'}{\pi^2} \frac{(s-s_0)\rho(s',t')}{(s'-s)(s'-s_0)} \left[ \frac{1}{t_0-t'} + \frac{1}{t'} {}_2F_1 \left( 1, \frac{d-2}{2}, d-2, \frac{s-4m^2}{t'} \right) \right] \\ \left. + \int_{4m^2}^{\infty} \frac{dt' du'}{\pi^2} \frac{\rho(t',u')}{(t'-t_0)(u'-u_0)} \left[ 1 - \frac{2A_2^{(d)}(s,t',u')}{s+t'+u'-4m^2} \right] \right),$$

$$A_2^{(d)}(s, t', u') = \frac{(t_0 - t')(s_0 + t_0 - (s+t')) {}_2F_1\left(\frac{1}{2}, 1; \frac{d-1}{2}; \frac{(s-4m^2)^2}{(-4m^2 + s + 2t')^2}\right)}{-4m^2 + s + 2t'} \\ + \frac{(u_0 - u')(u_0 + s_0 - (s+u')) {}_2F_1\left(\frac{1}{2}, 1; \frac{d-1}{2}; \frac{(s-4m^2)^2}{(-4m^2 + s + 2u')^2}\right)}{-4m^2 + s + 2u'}, \quad (A.7)$$

where recall that in the text we use $s_0 = t_0 = u_0 = \frac{4m^2}{3}$.

## B Grid generation

In this section, we explain the Mathematica routines that we used to generate our meshes. They rely on the use of the Finite Element Methods package, that can be called with «NDSolve`FEM`. Mathematica documentation is provided at https://reference.wolfram.com/language/FEMDocumentation/tutorial/ElementMeshCreation.html.

This package is very handy and in particular it provides interpolating routines via the command `ElementMeshInterpolation`.

**1d meshes.** The one-dimensional meshes are straightforward to generate. A grid with Nelem elements and Gradient Ratio r is generated as follows:

```
In[1]:= <<NDSolve`FEM`
```

```
In[2]:= meshx=ToGradedMesh[Line[{0},{1}}],<|"Alignment"->"BothEnds",
"ElementCount"->Nelem,"GradingRatio"->r|>,"MeshOrder"->1];
```

The grading ratio $r$, is quantity that enforces some exponential density in the grid, near the extremities (through the `"Alignment"` - > `"BothEnds"` option, which can also be set to `Right` or `Left`, see documentation). When $r = 1$ the grid is equally spaced, when $r$ increases, more points are set to the extremities.

Note that we encountered interpolation troubles with non-linear interpolation. Mathematica struggled with the singular behavior near 0 in the $x, y$ variables, and we had to resort to enforcing linear interpolation everywhere, through the `"MeshOrder"` - > 1 command.
**2d, $\rho(s, t)$ mesh.** The most delicate grid to generate in our work is the 2d grid for $\rho(s, t)$. This grid should have an exponential density of points near the Regge region, but as few points as possible in the bulk, where not much physics happens. In particular, we empirically observed that the double spectral densities generated by the iteration algorithm vary slowly at fixed $y$ in the bulk, $x \in [0.1, 0.9]$. One could have thought that many points are needed close to the Landau curves too. However due to elastic unitarity the spectral densities decay close to the leading Landau curve in a simple known fashion, see [38], and correspondingly we have observed in our analysis that this region does not require any special care.

The FEM package has a useful command `ProductMesh` that takes two meshes and produce a tensorial product of them. The elements of this tensor-product are quadrilateral elements which are rectangles. However, we are not interested in a mesh that covers the whole range $[0; 1] \times [0; 1]$, we only need to cover the wedge below the Landau curve $y = (1-x)/4$ for the double spectral density $\rho_{el}(x, y)$, see section 2.3.

The construction starts from a product of two meshes covering $[0; 1/4]$ and $[0; 1]$. Then one removes points above the Landau curve, add the boundary Landau curve. At this point we have a grid of points, from which Mathematica can in theory generate a new mesh.

This procedure turns out to not work well, because given a set of points, the automatic meshing occurs through “most neighboring points”. Since the density of points in the $x$ and $y$ directions can be very different in some regions, this sometimes gives rise to triangle mesh elements, especially near the Landau curve, which are very unnatural.

Therefore the re-meshing has to be done by hand. One keeps all the quadrilateral elements which are not cut by the Landau curve in the original tensor product, and constructs the triangle elements that connect these quadrilaterals to the Landau curve. This gives rise to a mesh such as the one showed in figure 6, with a mixture of rectangles and triangles, which now looks uniform.

Overall, we have the following steps to generate the two-dimensional grid for the double spectral density

1. Generate x-grid `meshX`,

2. Generate y-grid `meshY`,

3. Perform tensor product,

4. Remove points above the Landau curve and add points on the Landau curve boundary.

In our model, there are a few parameters which determine the grids. Let us enter the details of each grids.

**x-grid.** This grid describes the interval $[0; 1]$. It is mostly adapted to probe the Regge limit in $d = 4$. It is parametrized by 3 parameters, `nx`, `mdx`, `maxD` which are, respectively, the number of points, the Regge cutoff, and the maximal spacing in the bulk. The grid glues three meshes: one with exponential density of points to the cutoff that spans 0 to some
functionReverse[{start_,end_}, OptionsPattern[]]:= Module[ {n,len,r, base, minD}, len=end-start; n=OptionValue["ElementCount"];
minD=OptionValue["MinimalDistance"];
r =Abs@FindRoot[ len*(r^(1/n) - 1)/(r-1)-minD==0, {r,1.0000000000001,10*1/minD}, Method -> "Brent", WorkingPrecision->100,MaxIterations->10000][[1, 2]]; base=Table[(r^(i/n) - 1)/(r - 1),{i, 0, n}]; base=end-len*base; base//Reverse

**y-grid.** The y-grid has 2 parameters, ny, mdy is more straightforward. It just has ny points with an exponential density of points to the cutoff given by mdy, which in general we set to $10^{-12}$, to be consistent with the x cut-off.

In[6]:= ymesh[ny_,mdy_]:=ToGradedMesh[Line[{0},{1/4}]], < | "Alignment" -> function, "ElementCount" -> ny, "MinimalDistance" -> mdy |>, "MeshOrder" -> 1];

**Landau-curve points.** Finally, we add points on the Landau curve which are located at the position of the points on the x-grid, i.e. they are located at $(0, 1/4)$, $(x_1, lc(x_1))$, ..., $(1, 0)$ where $lc(x) = (1 - x)/4$, is, again, the Landau curve that defines the domain in which $\rho_{el}(x, y)$ is non-vanishing.

## C Elastic amplitudes in $d=2$: product of CDD zeros

This appendix reviews some material presented in [26].

In two spacetime dimensions there is no momentum transfer, $t=0$, and we only have $s \leftrightarrow u$ crossing symmetry. In this case, eq. (3.4) describes a consistent scattering amplitude, but without particle production. Such S-matrices are well known in two spacetime dimensions and they are given by the product of the so-called CDD factors, see e.g. [57]. These S-matrices can correspond either to poles, or zeros. Given our assumption that there are no bound states, the relevant S-matrices are given by a product of the CDD zeros

$$S(s) = \prod_{i=1}^{N_{\text{tot}}} S_{\text{CDD}}^{\text{zero}}(s, m_i), \\ S_{\text{CDD}}^{\text{zero}}(s, m_z) = \frac{\sqrt{s(4m^2 - s)} - \sqrt{m_z^2(4m^2 - m_z^2)}}{\sqrt{s(4m^2 - s)} + \sqrt{m_z^2(4m^2 - m_z^2)}} \quad (C.1)$$

Following the conventions of [26], we can define the connected scattering amplitude $T(s)$ as follows

$$S(s) = 1 + i \frac{T(s)}{\sqrt{s(s - 4m^2)}}. \quad (C.2)$$