Scattering amplitudes and contour deformations

TL;DR

This paper demonstrates how contour deformations can be used to compute scattering amplitudes in Minkowski space within a scalar model, highlighting challenges in locating resonance poles and proposing direct second sheet calculations.

Contribution

It introduces a method for calculating scattering amplitudes on the second Riemann sheet using contour deformations and unitarity, addressing limitations of analytic continuation methods.

Findings

Contour deformations enable access to Minkowski space scattering regions.

The scalar model shows virtual states instead of resonance poles.

Direct second sheet calculations can overcome analytic continuation difficulties.

Abstract

We employ a scalar model to exemplify the use of contour deformations when solving Lorentz-invariant integral equations for scattering amplitudes. In particular, we calculate the onshell 2 -> 2 scattering amplitude for the scalar system. The integrals produce branch cuts in the complex plane of the integrand which prohibit a naive Euclidean integration path. By employing contour deformations, we can also access the kinematical regions associated with the scattering amplitude in Minkowski space. We show that in principle a homogeneous Bethe-Salpeter equation, together with analytic continuation methods such as the Resonances-via-Pad\'e method, is sufficient to determine the resonance pole locations on the second Riemann sheet. However, the scalar model investigated here does not produce resonance poles above threshold but instead virtual states on the real axis of the second sheet, which pose difficulties for analytic continuation methods. To address this, we calculate the scattering amplitude on the second sheet directly using the two-body unitarity relation which follows from the scattering equation.

Equations

$$i\int d^{4}k\frac{1}{k^2-m^2+iϵ}\cdots =i\int d^{3}k\underset{-\infty }{\overset{\infty }{\int }}d{k}_0\frac{1}{k_0^2-{\omega }^2+iϵ}\dots$$

$$\int d^{3}k\underset{-\infty }{\overset{\infty }{\int }}d{k}_4\frac{1}{k_4^2+{\omega }^2}\cdots =\int d^4{k}_E\frac{1}{k_E^2+m^2}\dots ,$$

$$\mathcal{I}(t)=\int d^{4}k\frac{1}{k_{+}^2+m^2}\frac{1}{k_{-}^2+m^2}.$$

$$\underset{-\infty (1+iϵ)}{\overset{\infty (1+iϵ)}{\int }}d{k}_0\iff \underset{-\infty (i-ϵ)}{\overset{\infty (i-ϵ)}{\int }}d{r}_4,$$

$$\sum _{n=0}^{\infty }{e}^{-i{E}_{n}T}|n⟩⟨n|\Omega ⟩\stackrel{T\to \infty (1-iϵ)}{\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\text{relbar}\text{joinrel}\to }{e}^{-i{E}_{0}T}|0⟩⟨0|\Omega ⟩,$$

$$\int d^{4}k=\frac{1}{2}\underset{0}{\overset{\infty }{\int }}d{k}^2{k}^2\underset{-1}{\overset{1}{\int }}dz\sqrt{1-z^2}\underset{-1}{\overset{1}{\int }}dy\underset{0}{\overset{2\pi }{\int }}d\psi .$$

$$c=\frac{g^2}{(4\pi m{)}^2},\beta =\frac{\mu }{m}.$$

$$\begin{array}{rl}D(p^2)& =\frac{1}{Z}\frac{1}{p^2+M(p^2{)}^2},\\ D^{\prime }(p^2)& =\frac{1}{Z^{\prime }}\frac{1}{p^2+M^{\prime }(p^2{)}^2}\end{array}$$

$$\Psi =K{G}_0\Psi .$$

$$\Gamma ={\Gamma }_0+K{G}_0\Gamma .$$

$$\Gamma =\frac{{\Gamma }_0}{1-K{G}_0},$$

$$T=K+K{G}_{0}T,$$

$$T=\frac{K}{1-K{G}_0}$$

$$\Psi (q,P)=\underset{k}{\int }K(q,k)G_0(k,P)\Psi (k,P),$$

$$D(k_{\pm })=\frac{1}{k_{\pm }^2+m^2}=\frac{1}{k^2+\frac{1}{4}{P}^2+m^2\pm k\cdot P},$$

$$G_0(k,P)=\frac{1}{(k^2+\frac{1}{4}{P}^2+m^2{)}^2-(k\cdot P{)}^2}.$$

$$\begin{array}{rl}K(q,k)& =\frac{g^2}{2}\left[\frac{1}{(q-k{)}^2+{\mu }^2}+\frac{1}{(q+k{)}^2+{\mu }^2}\right]\\ & =\frac{g^2(q^2+k^2+{\mu }^2)}{(q^2+k^2+{\mu }^2{)}^2-4(q\cdot k{)}^2}.\end{array}$$

$$\begin{split}P&=2m\sqrt{t}\left[\begin{array}[]{c}0\\ 0\\ 0\\ 1\end{array}\right],\quad q=m\sqrt{X}\left[\begin{array}[]{c}0\\ 0\\ \sqrt{1-Z^{2}}\\ Z\end{array}\right],\\ k&=m\sqrt{x}\left[\begin{array}[]{l}\sqrt{1-z^{2}}\,\sqrt{1-y^{2}}\,\sin\psi\\ \sqrt{1-z^{2}}\,\sqrt{1-y^{2}}\,\cos\psi\\ \sqrt{1-z^{2}}\,y\\ z\end{array}\right]\end{split}$$

$$\begin{array}[]{rl}P^{2}&=4m^{2}\,t\,,\\ q^{2}&=m^{2}X\,,\\ k^{2}&=m^{2}x\,,\end{array}\quad\begin{array}[]{rl}q\cdot P&=2m^{2}\sqrt{Xt}\,Z\,,\\ k\cdot P&=2m^{2}\sqrt{xt}\,z\,,\\ q\cdot k&=m^{2}\sqrt{Xx}\,\Omega\end{array}$$

$$\begin{array}{rl}K(X,x,\Omega )& =\frac{g^2}{m^2}\frac{X+x+{\beta }^2}{(X+x+{\beta }^2{)}^2-4Xx{\Omega }^2},\\ G_0(x,z,t)& =\frac{1}{m^4}\frac{1}{(x+t+1{)}^2-4xt{z}^2}\end{array}$$

$$\begin{array}{rl}& \Psi (X,Z,t)=\frac{m^4}{(2\pi {)}^3}\frac{1}{2}\underset{0}{\overset{\infty }{\int }}dxx\underset{-1}{\overset{1}{\int }}dz\sqrt{1-z^2}\\ & \times \underset{-1}{\overset{1}{\int }}dyK(X,x,\Omega )G_0(x,z,t)\Psi (x,z,t)\end{array}$$

$$\underset{-1}{\overset{1}{\int }}dy\frac{1}{1-(a+by{)}^2}=\frac{1}{b}\left(\mathrm{atanh}{\delta }_{+}+\mathrm{atanh}{\delta }_{-}\right),$$

$$f_{\pm }(t,\alpha ,z)=\sqrt{t}z\pm i\sqrt{t(1-z^2)+{\alpha }^2},$$

$$q_{\pm }^2=-m^2\Rightarrow Z=0,X=-(1+t)$$

$$f_{\pm }(t,0,z)=t\left(z\pm i\sqrt{1-z^2}\right)=t{e}^{\pm i\phi }$$

$$z=\frac{a^2-b^2+{\alpha }^2/2}{a^2+b^2}.$$

$$\Psi (X,Z,t)=c\int dx\int dz\mathcal{K}(X,Z,x,z,t)\Psi (x,z,t).$$

$${\Psi }_{\sigma }(t)=c{\mathcal{K}}_{\sigma \tau }(t){\Psi }_{\tau }(t),$$

$$\frac{1}{{\lambda }_i(t_i)}\stackrel{!}{=}c,i=0,1,2,\dots$$

$$\text{Re}\frac{1}{{\lambda }_i(t)}=c\text{and}\text{Im}\frac{1}{{\lambda }_i(t)}=0.$$

Full text

I Introduction

The study of resonances is a central task in the nonperturbative treatment of quantum field theories. Among the observable states of a theory are bound states but also unstable resonances, which appear as bumps in cross sections and correspond to poles in the complex momentum plane on higher Riemann sheets. Among the prominent examples in QCD are the $\sigma$ meson, whose resonance pole position is now well-established Pelaez (2016), the excitations of light baryons, where multichannel partial-wave analyses of experimental precision data have led to the addition of several new states to the PDG Tanabashi et al. (2018), or the recently observed pentaquark states, where the nature of the neighboring peaks is still under debate Aaij et al. (2015, 2019).

The theoretical investigation of scattering amplitudes and their resonance structure comes with technical challenges that appear in different guises. In a lattice formulation one calculates the finite-volume energy spectrum of the theory and extracts the resonance information through the Luescher method Luescher (1986), where current efforts focus on resonances above two- and three-particle thresholds; see e.g. Hammer et al. (2017); Briceno et al. (2018); Mai and Doring (2019); Hansen and Sharpe (2019); Jackura et al. (2019) and references therein.

In continuum approaches, scattering amplitudes and their resonance information are accessible through scattering equations or Bethe-Salpeter equations (BSEs). Here the technical difficulties concern the numerical solution of four-dimensional scattering equations in the full kinematical domain and the extraction of resonance poles on higher Riemann sheets. On the one hand, the internal poles in the loop diagrams put restrictions on the accessible kinematic regions beyond which residue calculus or contour deformations into the complex momentum plane become inevitable. For example, without contour deformations one can only access low-lying excitation spectra or matrix elements in certain kinematic windows (the ‘Euclidean region’). On the other hand, to extract the properties of resonances it is necessary to access unphysical Riemann sheets, whereas straightforward numerical calculations are restricted to the first sheet only.

Both issues are typical obstacles in the functional approach of Dyson-Schwinger equations (DSEs) and BSEs, where one determines quark and gluon correlation functions and solves BSEs to arrive at hadronic observables. Progress has been made in the calculation of hadron spectra and matrix elements, see e.g. Cloet and Roberts (2014); Eichmann et al. (2016a); Sanchis-Alepuz and Williams (2018); Burkert and Roberts (2019) and Refs. therein, but the treatment of resonances is still in its early stages. Beyond technical aspects there are also conceptual challenges: when constructing hadron matrix elements from quarks and gluons, the bound states (such as $\pi$ and $\rho$ mesons, $N$ and $\Delta$ baryons) and decay mechanisms that turn these bound states into resonances (e.g. $\rho\to\pi\pi$ or $\Delta\to N\pi$) must both emerge from the underlying quark-gluon interactions. The resonance mechanism corresponds to the dynamical emergence of internal hadronic poles in matrix elements ($\pi\pi$, $N\pi$, etc.), which are also responsible for the so-called ‘meson cloud’ effects. Such properties have been studied by resumming quark and gluon topologies to intermediate meson propagators Fischer and Williams (2008); Sanchis-Alepuz et al. (2014); Williams (2018); Miramontes and Sanchis-Alepuz (2019) or by going to multiquark systems where these poles are generated dynamically Heupel et al. (2012); Eichmann et al. (2016b); Wallbott et al. (2019). In addition to internal hadronic poles, however, also the dynamical singularities encoded in the elementary quark and gluon correlation functions restrict the kinematical domains of matrix elements and must be taken into account when calculating observables in the far spacelike, timelike or lightlike regions.

In this work we focus on the methodological aspects, namely how to calculate scattering amplitudes in the kinematical domains where contour deformations become necessary (which is usually referred to as ‘going to Minkowski space’), and how that information can be used to extract the resonance information on higher Riemann sheets. Numerical contour deformations have been employed in the literature in the calculation of two- and three-point functions, see e.g. Refs. Maris (1995); Eichmann (2009); Strauss et al. (2012); Windisch et al. (2013a, b); Pawlowski et al. (2018); Williams (2018); Miramontes and Sanchis-Alepuz (2019); Weil et al. (2017); here we apply them to a four-point function in the form of a two-body scattering amplitude.

To illustrate the generic features, we exemplify the procedure for the simplest example, namely a scalar theory with a scalar exchange, which for a massless exchange particle becomes the well-known Wick-Cutkosky model Wick (1954); Cutkosky (1954); Nakanishi (1969). We calculate the $2\to 2$ scattering amplitude of the theory as well as its homogeneous Bethe-Salpeter (BS) amplitude and inhomogeneous BS vertex. We will see that in principle already the eigenvalues of the homogeneous BSE are sufficient to extract the resonance information. It turns out, however, that the model does not produce resonances above threshold but virtual states on the second Riemann sheet, which will pose difficulties for standard analytic continuation methods. Instead, we solve the scattering equation directly and employ the two-body unitarity property which allows us to calculate the scattering amplitude also on the second sheet.

Recent progress has also been made in calculating propagators and BS amplitudes in Minkowski space, see Carbonell and Karmanov (2010); Frederico et al. (2014); Jia and Pennington (2017); Leitão et al. (2017a, b); de Paula et al. (2017); Biernat et al. (2018); Ydrefors et al. (2019); Frederico et al. (2019); Solis et al. (2019) and references therein. Here we want to point out that there is no intrinsic difference between Euclidean and Minkowski space approaches: To obtain scattering amplitudes in the complex plane, contour deformations are necessary both in a Euclidean and Minkowski metric. When implemented properly, the resulting amplitude obtained with a Euclidean path deformation is identical to the Minkowski space amplitude. We discuss this in Sec. II and use a Euclidean metric for the remainder of this work; Euclidean conventions can be found in Appendix A of Ref. Eichmann and Ramalho (2018).

The paper is organized as follows. After illustrating contour deformations for simple examples in Sec. II, we establish the scalar model that we employ in Sec. III. In Sec. IV we solve its homogeneous and inhomogeneous BSEs, explain the contour deformation procedure and the analytic continuation to the second Riemann sheet. In Sec. V we solve the $2\to 2$ scattering equation of the theory, discuss the two-body unitary relation and present our results. We conclude in Sec. VI. Technical details on contour deformations are relegated to Appendix A.

II Contour deformations

To motivate the idea of contour deformations, we first illustrate the problem with simple examples. To begin with, consider an integral with only one propagator pole in the loop such as e.g. in a Fourier transform:

[TABLE]

where $\omega=\sqrt{\bm{k}^{2}+m^{2}}$. The integrand has poles on the real $k_{0}$ axis, whose locations depend on $\bm{k}^{2}$ and start at $k_{0}=\pm m$ (see Fig. 1). In the standard Minkowski treatment one exploits the $i\epsilon$ term to shift the poles away from the real axis, performs the $k_{0}$ integration by closing the integration contour at complex infinity, picks up the appropriate residues, and finally integrates over the $\bm{k}^{2}$ dependence of those residues.

In Euclidean conventions one defines $k_{4}=ik_{0}$, but because real and momentum space rotate in opposite directions one has $d^{4}k_{E}=-id^{4}k$. The integral then becomes

[TABLE]

where the integration proceeds from left to right in the complex $k_{4}$ plane and the poles lie on the imaginary axis.

Now suppose we interchange the $d^{3}k$ and $dk_{4}$ integrations and integrate over $d^{3}k$ first. Along the former pole positions we now find branch cuts in the complex $k_{4}$ plane starting at $k_{4}=\pm im$, illustrated in the right panel of Fig. 1, and instead of closing the contour analytically we integrate numerically by avoiding the cuts. Clearly, the two strategies are equivalent: if the integration path crossed a cut, one would not pick up all the poles and obtain a wrong result.

The Euclidean and Minkowski paths give identical results in this example. One can rotate the Minkowski path counterclockwise because there are no further poles in its way and the opposite integration direction is canceled by the sign in the integration measure — a Wick rotation is possible. In the following we drop the subscript ‘E’ and continue to work with Euclidean conventions.

Now consider an integral with two poles in it, such as the loop diagram in Fig. 2:

[TABLE]

The total momentum is $P^{\mu}$, the internal momenta are $k_{\pm}^{\mu}=k^{\mu}\pm P^{\mu}/2$, and we defined111Note that $t$ is related to the usual definition of the Mandelstam variable $\tilde{s}$ through $\tilde{s}=-4m^{2}t$. We adapted our notation to the Compton scattering kinematics in Sec. V.3; in the following we therefore refer to the resonant channel as the $t$ channel and to the crossed channels as $s$ and $u$ channels. the dimensionless variable $t=P^{2}/(4m^{2})$. The integral is Lorentz invariant and thus only depends on $t$. This is the simplest example of a two-point correlation function like a self-energy or vacuum polarization, which in principle can produce the singularity structure shown in Fig. 2. Bound states appear on the negative real axis of $t$ and resonances above the threshold $t<-1$ on higher Riemann sheets. The perturbative integral (3) can at best produce a two-particle cut but if the internal propagators and vertices were dressed and non-perturbative, the integral could also generate bound-state and resonance poles.

Because $\mathcal{I}(t)=\mathcal{I}^{\ast}(t^{\ast})$ is an analytic function, it is sufficient to consider the upper half plane in $t$ only: The real part is symmetric around the real axis and the imaginary part is antisymmetric. It is then more convenient to plot the function in the complex $\sqrt{t}$ plane, which confines it to the upper right quadrant (Fig. 3). In this case the bound states appear on the imaginary axis below threshold ($\text{Im}\sqrt{t}<1$), the cut starts at the threshold and the resonances lie above threshold on a higher Riemann sheet. In this way one can directly read off the real and imaginary parts of the masses $M_{i}$, which appear at $\text{Im}\sqrt{t}=\text{Re}\,M_{i}/(2m)$.

Suppose we want to calculate $\mathcal{I}(t)$ for some $t\in\mathds{C}$. Fig. 4 shows the resulting cuts in the complex plane of $r_{4}=k_{4}/m$. There are four vertical cuts centered around the external point $\pm\sqrt{t}$. Since we divided out the mass, the vertical distance between $\sqrt{t}$ and the onset of the cuts is equal to 1. As before, the Euclidean integration path proceeds from left to right.

If $\text{Im}\sqrt{t}>1$, however, the cuts cross the real axis and the straight Euclidean path (we denote it by E’) would cross the cuts. Hence we must deform the contour to avoid the cuts: The correct Euclidean path is E. As a consequence, E = E’ only below the threshold $\text{Im}\sqrt{t}<1$, i.e., in the colored region in Fig. 3, where a naive Euclidean integration is sufficient and gives the correct result. Above threshold, one has to deform the contour to obtain the correct value of the integral. The situation can be generalized to unequal masses or complex propagator poles, but the principle is the same: a straight Euclidean integration path is only valid in a limited domain of complex $t$.

What would be the corresponding Minkowski path? Apparently it cannot proceed along the vertical axis as in Fig. 1: It does not matter whether we start slightly on the right and end up slightly on the left because there are no singularities on the imaginary axis. In fact, the $i\epsilon$ prescription entails

[TABLE]

since it originates from the need to isolate the interacting vacuum $|\Omega\rangle$ in a correlation function,

[TABLE]

and thereby remove the higher energy contributions $E_{n}$ of the free $n-$particle states $|n\rangle$. The integration path between $T=\pm\infty(1-i\epsilon)$ in the action of the quantum field theory thus corresponds to $k_{0}\to\pm\infty(1+i\epsilon)$ and $r_{4}\to\pm\infty\,(i-\epsilon)$. Therefore, the proper Minkowski path is the diagonal line from bottom right to top left, which must also be deformed to avoid the cuts, cf. Fig. 5.

That this is indeed the correct path can also be seen by putting the point $\sqrt{t}$ back onto the imaginary axis (right panel in Fig. 5). In that case all cuts also lie on the imaginary axis and one can use the $i\epsilon$ term to displace the cuts, while the Minkowski path is the straight line from bottom to top. As a consequence, all cuts in the upper half plane of $r_{4}$ are shifted to the right of the path and all cuts in the lower half plane to its left; closing the contour on either side gives the same result. This is also what happens in the left panel of Fig. 5, where the upper cuts appear on the right of the Minkowski path and the lower cuts on its left — but that path is just the same as the Euclidean contour.

In general there is no intrinsic difference between quantities in ‘Euclidean’ or ‘Minkowski space’. In both cases one has to deform integration contours to avoid cuts and the final result is the same. A naive Euclidean integration path would give the wrong result above threshold, whereas a naive Minkowski integration (in the sense of a straight vertical path) becomes meaningless once $t$ is complex. The $i\epsilon$ prescription in the action only tells us where the integration starts and ends; that one additionally has to deform contours at the level of correlation functions is implicit in their definition. Therefore, what we mean by ‘Euclidean conventions’ is merely a Euclidean metric with metric tensor $\delta^{\mu\nu}$ as opposed to a Minkowski metric $g^{\mu\nu}$. A collection of Euclidean conventions can be found e.g. in Appendix A of Ref Eichmann and Ramalho (2018).

As an aside, there is also no inherent problem with perturbation theory which is usually done in Euclidean space. The ‘Wick rotation’ only amounts to rewriting the integrals in a Euclidean metric, which is always possible as long as the singularities of their integrands are kept in mind and there is no contribution from the contour at complex infinity. If one employs Feynman parametrizations, then the integrals in momentum space become integrals over Feynman parameters and one has to analyze the singularity structure of the integrands in Feynman parameter space instead of momentum space.

In the remainder of this paper we will need contour deformations for integrals with three and four internal poles, which are also not just integrals but appear in integral equations. Because all expressions are Lorentz invariant, it is convenient to preserve manifest Lorentz invariance by splitting the loop integral into $k^{2}$ and a three-dimensional solid angle via hyperspherical variables:

[TABLE]

Equivalently, one could use the radial variable $k=\sqrt{k^{2}}$ such that $dk^{2}\,k^{2}/2=dk\,k^{3}$. As a consequence, the cuts are no longer vertical lines in the complex $k^{2}$ plane but pick up more complicated shapes which we discuss in Sec. IV.2. In any case, the strategy is the same: Depending on the external kinematics, one must deform the straight Euclidean integration path from $k^{2}=0$ to $k^{2}\to\infty$ to avoid the cuts in the complex $k^{2}$ plane.

III Scalar model

We consider the simplest scalar model that is capable of producing resonances: Two scalar particles $\phi$ and $\varphi$ with masses $m$ and $\mu$, respectively, and a three-point interaction $\sim g\phi\phi\varphi$. This leaves two parameters:

[TABLE]

$\beta$ is the mass ratio and the coupling $g$ is dimensionful, so we defined a dimensionless coupling constant $c$. As a consequence, the mass $m$ drops out from all equations and only sets the scale.

In principle one should dress the propagators by solving their coupled Dyson-Schwinger equations (DSEs) shown in Fig. 6. However, tree-level propagators are good enough for our purposes because the dressing effects in the scalar theory are relatively small. If we define a mass function for each propagator, $M(p^{2})$ and $M^{\prime}(p^{2})$, via

[TABLE]

where $Z$ and $Z^{\prime}$ are the respective renormalization constants, then an exemplary solution of the coupled DSEs for a given value of $c$ (and using tree-level vertices only) is shown in Fig. 7. The three dots on each curve correspond to three different renormalization points at which the renormalized masses $m$ and $\mu$ were specified as input; the renormalization constants $Z$ and $Z^{\prime}$ are determined in the solution process. One can see that both mass functions are essentially flat over a large momentum range and only begin to drop in the ultraviolet.

If the coupling $c$ becomes large enough, both curves for the squared mass functions eventually cross zero and become negative in the ultraviolet. This reflects the vacuum instability of the $\phi\phi\varphi$ theory and implies that the model is only physically acceptable for small couplings, which has also consequences for the so-called anomalous BSE solutions Ahlig and Alkofer (1999). Since we employ the scalar theory only as a toy model for calculating resonance properties, we will not further consider this point and restrict ourselves to tree-level propagators in what follows.

To extract the bound state and resonance properties of the model, we consider the three equations depicted in Fig. 8. The first is the homogeneous BSE for the BS amplitude $\Psi$, which in a compact form reads

[TABLE]

$K$ stands for the BS kernel and $G_{0}$ for the product of the two propagators. We employ a simple ladder exchange, which for a vanishing exchange-particle mass $\mu=0$ becomes the Wick-Cutkosky model Wick (1954); Cutkosky (1954); Nakanishi (1969). The homogeneous BSE yields the eigenvalues of $KG_{0}$ for given quantum numbers $J^{PC}$, which depend on the total momentum variable $t\in\mathds{C}$. If an eigenvalue satisfies $\lambda_{i}(t)=1$, this corresponds to a pole in the scattering amplitude. Hence, in principle one can extract both the bound-state and resonance information from the homogeneous BSE.

The second equation in Fig. 8 is the inhomogeneous BSE for the vertex $\Gamma$:

[TABLE]

The three-point vertex is essentially the scattering amplitude in a given $J^{PC}$ channel: If $G$ denotes the full four-point function satisfying $G=G_{0}+G_{0}KG$, then the vertex is defined as $G_{0}\Gamma=G\,\Gamma_{0}$, i.e., modulo external propagators it is the contraction of $G$ with the tree-level vertex $\Gamma_{0}$ carrying quantum numbers $J^{PC}$. Thus, it contains the bound-state and resonance poles in that channel. In a symbolic form one can write the solution of Eq. (9) as

[TABLE]

which shows explicitly that whenever an eigenvalue of $KG_{0}$ becomes 1, one has found a pole in the vertex.

Finally, the third equation in Fig. 8 is the scattering equation for the full scattering amplitude $T$:

[TABLE]

where $T$ is the connected and amputated part of $G$ defined via $G=G_{0}+G_{0}TG_{0}$. Since we denoted the total momentum variable by $t$, we will refer to the ‘horizontal’ channel as the $t$ channel in the following. To ensure crossing symmetry in the $s$ and $u$ channels, we have symmetrized the kernel in the scattering equation; for the two preceding equations this is not necessary because both terms yield the same result. Symbolically, the solution of the scattering equation has the form

[TABLE]

and therefore $T$ contains all singularities, including the exchange-particle poles encoded in the kernel $K$.

In the following we investigate these equations in detail: the homogeneous and inhomogeneous BSEs in Sec. IV and the scattering equation in Sec. V.

IV Resonances from the BSE

IV.1 Explicit form of the BSE

The homogeneous BSE reads explicitly:

[TABLE]

where $\int_{k}=d^{4}k/(2\pi)^{4}$ is the integral measure, $q$ is the relative and $P$ the total momentum, $k$ is the relative momentum in the loop, and the propagator momenta are $q_{\pm}=q\pm P/2$ and $k_{\pm}=k\pm P/2$. We restrict ourselves to scalar bound states with quantum numbers $J^{PC}=0^{++}$; in that case $\Psi(q,P)$ is a scalar as well.

The internal propagators are given by

[TABLE]

so their product $G_{0}(k,P)=D(k_{+})\,D(k_{-})$ becomes

[TABLE]

Symmetrizing the kernel is not necessary at this stage but we do it nevertheless for later convenience. The kernel is then the sum of $s$ and $u$-channel exchanges:

[TABLE]

Using the hyperspherical variables from Eq. (5), we choose the rest frame

[TABLE]

so that the radial integration variable is $k^{2}=m^{2}x$ and its external counterpart is $q^{2}=m^{2}X$. In principle, however, we never need to specify a frame if we define those variables in a Lorentz-invariant way:

[TABLE]

where $\Omega=zZ+\sqrt{1-z^{2}}\sqrt{1-Z^{2}}\,y$ . As a result, the kernel and two-body propagator become

[TABLE]

and the BSE takes the form

[TABLE]

where the mass $m$ drops out. The integration over the angle $\psi$ is trivial because no Lorentz invariant depends on it and the integral over $y$ can be performed analytically,

[TABLE]

with $\operatorname{atanh}\delta=[\ln(1+\delta)-\ln(1-\delta)]/2$ and $\delta_{\pm}=b/(1\pm a)$.

The inhomogeneous BSE has the same form as Eq. (20) if we replace $\Psi(X,Z,t)\to\Gamma(X,Z,t)$ and add an inhomogeneous term $\Gamma_{0}=1$ on the right-hand side corresponding to quantum numbers $J^{PC}=0^{++}$.

IV.2 Contour deformation

To proceed, we analyze the singularity structure of the integrand in Eq. (20). The equation is solved at fixed $t$, which remains an external parameter. After performing the $y$ and $z$ integrations, the poles in $K$ and $G_{0}$ produce cuts in the complex plane of the radial variable $x$, and if those cuts cross the positive real axis the ‘naive’ integration path $0<x<\infty$ is forbidden. Since we discuss the complex $\sqrt{t}$ plane instead of $t$, we also analyze the cut structure in the complex $\sqrt{x}$ plane; and we do not change the integration domain of the variable $z$ which remains in the interval $-1<z<1$.

The singularities encoded in the quantities (19) can be written in a closed form by defining the function

[TABLE]

where $t\in\mathds{C}$ and $\alpha>0$ are parameters and the resulting cut is parametrized by $-1<z<1$:

• (1)

The propagator poles encoded in $G_{0}$ generate a cut $\sqrt{x}=\mathcal{C}_{1}=f_{\pm}(t,1,z)$.

• (2)

The exchange particle poles in $K$ produce a cut $\sqrt{x}=\mathcal{C}_{3}=f_{\pm}(X,\beta,\Omega)$, with $\Omega$ defined below (18).

• (3)

We may need the vertex $\Gamma(X,Z,t)$ at the kinematical point where also the constituent particles are onshell. In that case

[TABLE]

so that the kernel produces another cut at $\sqrt{x}=\mathcal{C}_{2}=f_{\pm}(-(1+t),\beta,\sqrt{1-z^{2}}\,y)$. The combination of $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ will become especially relevant for the onshell scattering amplitude in Sec. V.

• (4)

By means of the inhomogeneous BSE, the vertex $\Gamma(X,Z,t)$ develops bound-state poles and multi-particle cuts which are functions of $t$ only. Since the equation is solved at a given $t$, the only constraints that these singularities provide is that one cannot obtain a solution directly along the cut or at a pole location. In practice we avoid the imaginary $\sqrt{t}$ axis and solve the equation for $\text{Re}\sqrt{t}>0$, which circumvents both problems.

• (5)

In principle the BS amplitude or vertex can dynamically develop singularities in the relative momentum variables $X$ and $Z$. The variable $Z$ is protected because we keep it in the interval $-1<Z<1$, and as long as the singularities in $X$ only appear on the timelike $X$ axis they do not pose any restriction.

Let us analyze the function $\sqrt{x}=f_{\pm}(t,\alpha,z)$ in the complex $\sqrt{x}$ plane depending on the values of $t\in\mathds{C}$ and $\alpha>0$. We denote $\sqrt{t}=a+ib$, where $a,b\geq 0$ because $\sqrt{t}$ is in the upper right quadrant. Because we can identify $-\sqrt{x}$ with $+\sqrt{x}$, it is sufficient to consider only one of the cuts, for example $f_{-}(t,\alpha,z)$. The endpoints of the cut correspond to $z=\pm 1$ and thus $\sqrt{x}=\pm\sqrt{t}-i\alpha$. The cut starts on the vertical line $\sqrt{x}=a$ with distance $\alpha$ below the point $\sqrt{t}$ and ends on the opposite line $\sqrt{x}=-a$ with distance $\alpha$ below $-\sqrt{t}$.

For $\alpha=0$ the cut becomes a half-circle

[TABLE]

with $z=\cos\varphi$. For $\alpha>0$ it produces more complicated shapes as shown in Fig. 9. When lowering $z$ from $+1$ to $-1$, the modulus of $f_{-}(t,\alpha,z)$ always grows; and as long as $\alpha<2b$, any cut passes through the complex conjugate point $\sqrt{t^{\ast}}=a-ib$ when

[TABLE]

In other words, all possible cuts remain inside the colored region in Fig. 9, which is defined by the vertical lines $\sqrt{x}=\pm a$ and the circle with radius $|\sqrt{t}|=\sqrt{a^{2}+b^{2}}$. If $\alpha>b$, a contour deformation is not necessary because the cut does not cross the real axis.

Let us apply our findings to Eq. (20). The propagator cut $\mathcal{C}_{1}$ has the form discussed above with $\alpha=1$. If $b=\text{Im}\sqrt{t}<1$, it does not cross the real axis and no action is required — a ‘naive’ Euclidean integration is possible. If $b>1$, we must deform the integration contour; at this point any possible path that avoids the cut is sufficient.

The kernel cut $\mathcal{C}_{3}$ is analogous but with $t$ replaced by $X$ and $\alpha=\beta$. If Eq. (20) were just an integral and not an integral equation, then for $\text{Im}\sqrt{X}<1$ no further steps would be necessary. However, the fact that we feed the amplitude $\Psi(X,Z,t)$ back into the integral on the r.h.s. entails that the paths for $X$ and $x$ must match, i.e., we solve the equation along the same path for $X$ and $x$. Therefore, at every point $\sqrt{x}$ along the path the kernel generates another cut $\mathcal{C}_{3}$ that must be avoided.

All those cuts are still confined inside regions analogous to Fig. 9, so that each point $\sqrt{x}$ defines a corresponding cut region. To ensure that we never cross any cut along the path, we must stay outside of the regions belonging to the previous points on the path. This means once we have reached the point $\sqrt{x}$, the allowed region to proceed is bounded by a vertical line and a circle. For $\sqrt{x}$ in the upper right quadrant we may turn left only up to a vertical line and we may turn right and return to the real axis not faster than on a circle. In other words, both $\text{Re}\sqrt{x}$ and $|\sqrt{x}|$ must never decrease along the path.

A possible path satisfying these constraints is drawn in Fig. 9: The first section is a straight line connecting the origin with the point $\sqrt{t}$; the second is an arc that starts at $\sqrt{t}$ and returns to the real axis with increasing modulus; and the third is a straight line along the real axis up to infinity. In this way we can cover the entire complex plane and solve the homogeneous BSE (20) as well as the inhomogeneous BSE for any $\sqrt{t}\in\mathds{C}$.

If the BS amplitude or vertex are taken fully onshell, we must also consider the cut $\mathcal{C}_{2}$. Its form is analogous except that the point $t$ is replaced by $-(1+t)$, so that for fixed $t$ one must circumvent two cuts in the complex $\sqrt{x}$ plane. This situation is discussed in Appendix A.

IV.3 Eigenvalues of the homogeneous equation

To solve the homogeneous BSE (20), we write it as

[TABLE]

We pulled out the coupling $c$ since in this model it is just a constant overall factor in the kernel and the propagators do not depend on it because they remain at tree level. The contour deformation discussed above only enters in the integral measure $dx$. Eq. (26) can be written as a matrix-vector equation

[TABLE]

where the matrix indices $\sigma$, $\tau$ absorb the momentum dependence in the variables $X$ and $Z$. In practice we calculate the eigenvalues $\lambda_{i}(t)$ of the matrix $\mathcal{K}(t)$, which aside from $t$ only depend on the mass ratio $\beta$ but no longer on the coupling $c$. A solution of the equation and thus a pole in the scattering amplitude corresponds to the points $t_{i}$ where the condition

[TABLE]

is satisfied. The respective eigenvector $\Psi_{\sigma}(t_{i})$ is the onshell BS amplitude of that state.

The situation is illustrated in Fig. 10. The left panel shows the first three inverse eigenvalues of $\mathcal{K}(t)$ as a function of $\text{Im}\sqrt{t}$, i.e., along the imaginary axis up to the threshold in Fig. 3. The right panel magnifies the eigenvalue for the ground state. We chose a mass ratio $\beta=4$ where the eigenvalues are well separated. For $c\lesssim 6$ there is no bound state. For $6\lesssim c\lesssim 11$, the ground state is bound since the condition $1/\lambda_{0}=c$ can be satisfied; its mass corresponds to $\text{Im}\sqrt{t_{0}}=M_{0}/(2m)$.

For larger couplings, the condition is still satisfied for spacelike values of $t$: The eigenvalues continue to drop for $P^{2}>0$ and the inverse eigenvalues grow for $\sqrt{t}>0$ on the real axis. Therefore, with increasing coupling the pole in the scattering amplitude slides down on the imaginary axis of $\sqrt{t}$ in Fig. 3, becomes tachyonic and continues to move along the positive real axis towards infinity. Recalling Sec. III, this signals again that for physically acceptable solutions the coupling of the model must be restricted to small enough values.

If we increase the coupling further, then from Fig. 10 eventually also the first excited state becomes bound, until its mass drops to zero and becomes tachyonic; followed by the second excited state, etc. Hence the question: what is the nature of a state before it becomes bound? In particular, what happens for small couplings where all states, including the ground state, are above threshold (here for $c\lesssim 6$)?

One should note that the onshell BS wave function $\chi(q,P)=G_{0}(q,P)\,\Psi(q,P)$, which in quantum field theory is defined as the Fourier transform of the matrix element of two field operators between the vacuum and the asymptotic one-particle state $|P\rangle$, technically only refers to the solution of Eq. (27) at a particular bound-state pole for real $t_{i}$. Eq. (27), on the other hand, provides an analytic continuation of $\Psi(t)$ for general $t\in\mathds{C}$ which follows from the scattering equation, so that $\Psi(t)$ is the residue of the scattering amplitude also for resonance poles in the complex plane. In particular, the matrix $\mathcal{K}(t)$ contains the scattering information in the $t$ channel for given quantum numbers and Eq. (28) is the general condition for a pole in the scattering matrix, irrespective of whether it is real or complex.

Fig. 11 shows the same three eigenvalues as before (the lower panels zoom in on the ground state), but now including the results above threshold obtained with the contour deformation. The 15 different curves correspond to evenly spaced vertical lines in the complex $\sqrt{t}$ plane between $0.01<\text{Re}\sqrt{t}<0.65$, where the smallest curve in each plot is the result closest to the imaginary axis. Very close to the axis one requires better and better numerics to get stable results above threshold; the curve for the ground state in Fig. 11 is stable but artifacts in the excited states can already be seen and they grow for higher-lying states.

In Fig. 11 one clearly sees the non-analyticity at the threshold, where the real and imaginary parts of all three eigenvalues have a kink and a branch cut opens in the imaginary parts. The condition $1/\lambda_{i}(t)=c$ is never satisfied in the complex plane and thus the first Riemann sheet is free of singularities as expected. In principle, however, it can be met on the second Riemann sheet, which is the analytic continuation of the first sheet above threshold. Because the coupling is real, the pole position is the intersection of

[TABLE]

IV.4 Continuation to the second sheet

To analytically continue the eigenvalues to the second Riemann sheet, we employ the Schlessinger point method or Resonances-via-Padé (RVP) method Schlessinger (1968), which has been advocated recently as a tool for locating resonance positions in the complex plane Haritan and Moiseyev (2017); Tripolt et al. (2017, 2019). It amounts to a continued fraction

[TABLE]

which is simple to implement by an iterative algorithm. Given $n$ input points $z_{i}$ with $i=1\dots n$ and a function whose values $f(z_{i})$ are known, one determines the $n$ coefficients $c_{i}$ and thereby obtains an analytic continuation of the original function for arbitrary values $z\in\mathds{C}$. The continued fraction can be recast into a standard Padé form in terms of a division of two polynomials.

The situation in the complex $\sqrt{t}$ plane is shown in the upper panel of Fig. 12. If the two Riemann sheets are plotted next to each other, then crossing over from the first to the second sheet is analytic above the threshold $\text{Im}\sqrt{t}=1$. The alignment of the two sheets (left or right) is not important because on each sheet the eigenvalues are analytic functions whose real parts are symmetric around the imaginary $\sqrt{t}$ axis and whose imaginary parts are antisymmetric. Bound states can appear on the first sheet along the imaginary $\sqrt{t}$ axis below the threshold and resonances in the complex plane of the second sheet above threshold. In addition, typical features of $s$-wave amplitudes are virtual states Glöckle (1983); Hanhart et al. (2014); Pelaez (2016), which are poles on the imaginary $\sqrt{t}$ (or real $t$) axis on the second sheet below threshold.

Note that with the two sheets aligned side by side there is now a cut below threshold, for $0<\text{Im}\sqrt{t}<1$, which separates the two sheets — crossing over to the second sheet is smooth above threshold whereas below it is not. This poses a difficulty for the RVP method because a Padé ansatz cannot reproduce branch cuts but only poles. We therefore adapt the strategy of Ref. Hanhart et al. (2014) and unfold the two sheets by considering the $\sqrt{1+t}$ plane shown in the bottom panel of Fig. 12. Bound states now appear on the positive real axis, which eventually turns into the spacelike axis, whereas virtual states would appear as poles on the negative real axis. The lower half plane is again mirror-symmetric and does not provide new information. Importantly, there is no longer a cut in the complex $\sqrt{1+t}$ plane, so one can analytically continue to the second sheet also below threshold.

Our setup is shown in the left panel of Fig. 13. The right half plane is the first Riemann sheet in $\sqrt{t}$, where we calculate the BSE eigenvalues either below threshold without or above threshold with contour deformations. The left half plane is the second sheet. To determine the sensitivity of the RVP algorithm to the input region where the function is defined, we choose five different domains which are shown by the colored rectangles. Inside each rectangle we pick $n=20,21,\dots 80$ input points randomly. For each rectangle and each $n$ we then determine the resulting pole positions from Eq. (29). Finally we average those pole positions over $n$ to get an estimate on the sensitivity to $n$. For each of the five input regions, the resulting pole averages with their standard deviations are given by the circles in Fig. 13.

The center panel in Fig. 13 shows the result for $\beta=4$ and $c=8$. This merely serves as a check because in this case we already know the result; from Fig. 10 and Eq. (28) one can read off the bound-state pole position $\sqrt{t}\approx 0.84\,i$. Indeed, the RVP method accurately reproduces this value for either of the three input regions below threshold, whereas the two upper regions lead to scattered poles around that value and thus larger circles.

In the right panel of Fig. 13 we show the result for a weak coupling $c=1$ where the bound state has presumably turned into a resonance (cf. Fig. 10). In this case, however, we do not find resonance poles above threshold but rather below the threshold on the second sheet. For a given $n$ the method typically finds one pole in the $\sqrt{1+t}$ plane slightly above or below the negative real axis but sometimes also two. Transformed back to $\sqrt{t}$, the poles can dive under the cut and thus appear on the right half plane but they still lie on the second sheet. This is of course an artifact since the method is not aware of the analyticity property where the resulting function on the lower half plane of $\sqrt{1+t}$ must be the complex conjugate of that on the upper half plane; it merely performs an analytic continuation. Keeping this in mind, the results indeed point towards the occurrence of virtual states instead of resonances above the threshold.

As an independent check we also solve the inhomogeneous BSE (9) for the vertex $\Gamma(X,Z,t)$. Its practical form is analogous to Eqs. (26–27) with the addition of an inhomogeneity

[TABLE]

and instead of the eigenvalues of $\mathcal{K}(t)$ the inhomogeneous BSE determines the vertex $\Gamma_{\sigma}(t)$ directly. Its singularity structure coincides with the singularities of the full scattering amplitude for quantum numbers $J^{PC}=0^{++}$. Eq. (31) can be solved by iteration, except in the vicinity of the poles where we employ matrix inversion to avoid convergence problems:

[TABLE]

Once $\Gamma(X,Z,t)$ is known, we solve the equation one more time to obtain the onshell vertex $\Gamma(-(1+t),0,t)$ which is a function of $t$ only, cf. Eq. (23).

Fig. 14 shows the resulting onshell vertex for $\beta=4$ and five different values of the coupling. The different curves correspond to the same grid as for the eigenvalues in Fig. 11. For $c=7$ and $c=6$ there is a clear bound-state pole. For $c\leq 5$ the coupling is below the value where bound states can occur and the ground state becomes unbound (cf. Fig. 10). Instead of producing a resonance bump, however, the pole simply disappears and what remains is just the threshold cusp. The cusp is the onset of the cut below threshold: On the $\text{Im}\sqrt{t}$ axis, the continuation of $\text{Re}\,\Gamma$ from above threshold becomes the function on the second sheet, just like for the eigenvalues in Fig. 11. Combined with the RVP method, the resulting pole positions from the inhomogeneous BSE are compatible with those in Fig. 13 obtained from the ground-state eigenvalue of the homogeneous BSE.

We therefore conclude that the scalar model does not produce resonances above threshold. Our results so far yield poles on the second sheet below threshold, which suggests virtual states. This poses an obstacle for analytic continuation methods since one has to bridge a considerable distance in the complex $\sqrt{1+t}$ plane. That the RVP method is well suited for finding resonance poles above thresholds has been demonstrated in Ref. Tripolt et al. (2017). Thus, if the model did produce clear resonance bumps, the homogeneous (or inhomogeneous) BSE in combination with contour deformations and the RVP method would form an adequate toolbox to determine the pole locations in the complex plane.

On the other hand, there is a method that provides direct access to the second sheet: two-body unitarity, which follows from the scattering equation and allows one to calculate the scattering amplitude on the second sheet if it is known on the first sheet. We discuss it in more detail in Sec. V.2; to utilize it, we must first solve the full scattering equation in Eq. (11).

V Scattering amplitude

V.1 Onshell scattering amplitude

The scattering equation for the $2\to 2$ scattering amplitude $T(q,p,P)$ is shown in Fig. 15 and reads explicitly:

[TABLE]

The ingredients are the same as before, Eqs. (15–16), except that the amplitude depends on a second relative momentum $p$ with $p_{\pm}=p\pm P/2$.

Ultimately we are only interested in the onshell scattering amplitude where all particles are on their mass shells, $p_{\pm}^{2}=q_{\pm}^{2}=-m^{2}$, which entails $p\cdot P=q\cdot P=0$ and $p^{2}=q^{2}=-m^{2}\,(1+t)$. This leaves two independent variables: the total momentum transfer $t=P^{2}/(4m^{2})$ and the crossing variable $\lambda=p\cdot q/m^{2}$. The latter can be written as $\lambda=(1+t)\,Y$, where the hyperspherical variable $Y=\cos\theta$ is the cosine of the scattering angle $\theta$ in the CM frame. In the following we denote the onshell scattering amplitude by $T_{\text{on}}(t,Y)$.

In Fig. 16 we illustrate the singularity structure of $T_{\text{on}}(t,Y)$ in the Mandelstam plane of the variables $t$ and $\lambda$. The bound states, resonance poles and $t$-channel cuts generated in the nonperturbative solution of the scattering equation are identical to those obtained with the (in-)homogeneous BSEs. They appear at fixed $t$ and are independent of $\lambda$, so they form horizontal lines in the Mandelstam plane. The physical $t$ channel opens above threshold $t<-1$.

By contrast, the singularities that depend on $\lambda$ are purely perturbative. On the one hand, the exchange kernel $K(q,p)$ has poles in the $s$ and $u$ channels. The Mandelstam variables $s$ and $u$ are given by

[TABLE]

and thus $\lambda=(s-u)/(4m^{2})$. The exchange poles appear at $s=\mu^{2}$ or $u=\mu^{2}$, which corresponds to

[TABLE]

as shown in Fig. 16. On the other hand, by expanding the scattering equation into a ladder series, each perturbative loop diagram produces further cuts if $s>4\mu^{2}$ or $u>4\mu^{2}$ and therefore

[TABLE]

The boundaries of the three channels are the lines with $s=0$ and $u=0$ corresponding to $\lambda=\mp(1+t)$ and therefore $Y=\mp 1$.

Because $T_{\text{on}}(t,Y)$ is invariant under the transformations $p\to-p$ or $q\to-q$, the scattering amplitude is symmetric in the crossing variable $\lambda$ and thus can depend on $Y$ only quadratically. The $Y$ dependence can be absorbed in a partial-wave expansion:

[TABLE]

where $P_{l}(Y)$ are the Legendre polynomials and inside the $t$-channel region one has $-1<Y<1$. Because the dependence on $Y$ is quadratic, only even partial waves $l=0,2,4\dots$ survive.

If $\beta$ is large enough and the exchange-particle poles are far away from the $t$-channel region, the $\lambda$ dependence is small and the partial-wave expansion converges rapidly. Moreover, the exchange poles only appear in the inhomogeneous term of Eq. (33), so they drop out in the difference $T_{\text{on}}-K_{\text{on}}$ whose remaining $\lambda$-dependent singularities are the cuts in the $s$ and $u$ channels. Thus, the only nonperturbative singularities of the scattering amplitude are those in the variable $t$.

V.2 Two-body unitarity

The scattering equation provides a direct way to access the second Riemann sheet via two-body unitarity, see e.g. Ref. Gribov (2012) for a detailed discussion. If we take the inverse of the scattering equation, $T^{-1}=K^{-1}-G_{0}$, and subtract it for two different kinematical configurations we obtain

[TABLE]

Multiplying with $T_{+}$ from the left and $T_{-}$ from the right yields

[TABLE]

Using $G=G_{0}+G_{0}\,K\,G=G_{0}+G_{0}\,T\,G_{0}$, the second term on the r.h.s. above can be rearranged as

[TABLE]

If $T_{+}$ and $T_{-}$ are the scattering amplitudes along the $t$-channel cut ($t<-1$) with $t_{\pm}=t\pm i\epsilon$ slightly displaced on either side, then the difference of the kernels drops out because the ladder kernel $K$ does not depend on $t$. The difference of the tree-level propagators forces the scattering amplitudes inside the loop integral onto the mass shell and results in the unitarity relation

[TABLE]

where

[TABLE]

From the knowledge of the discontinuity along the cut one therefore arrives at a relation between the onshell scattering amplitude on the first and second sheet:

[TABLE]

which is an inhomogeneous integral equation that determines the amplitude $T_{\text{II}}$ on the second sheet once $T_{\text{I}}$ is known. Inserting the partial-wave decomposition (37), it becomes the algebraic relation

[TABLE]

where $f_{l}(t)_{\text{I}}$ and $f_{l}(t)_{\text{II}}$ denote the partial waves on the first and second sheet. Thus, the poles on the second sheet correspond to the zeros of the denominator on the first sheet:

[TABLE]

Once the scattering amplitude is solved from its scattering equation (33), unitarity is therefore automatic and one has direct access to the second Riemann sheet.

V.3 Half-offshell amplitude

There is one complication that remains to be addressed. Because the internal momenta $k_{\pm}$ in Eq. (33) are sampled in offshell kinematics, the scattering equation does not provide a self-consistent relation for the onshell amplitude but only for its half-offshell counterpart. Therefore, we must first solve the equation for the half-offshell amplitude; once completed, we perform ‘one more iteration’ to arrive at the scattering amplitude for onshell external kinematics.

For the half-offshell amplitude we relax the conditions $q_{\pm}^{2}=-m^{2}$ so that $q\cdot P$ and $q^{2}$ remain general. As a consequence, the amplitude now depends on four independent variables. On kinematical grounds the half-offshell amplitude is similar to onshell Compton scattering with two virtual photons, so we can use the same Lorentz-invariant momentum variables to analyze it Eichmann and Ramalho (2018):222The kinematical analogue of Fig. 15 is the annihilation process $N\bar{N}\to\gamma^{\ast}\gamma^{\ast}$. To compare with the notation in Ref. Eichmann and Ramalho (2018), replace the momenta $\{q,p,P\}$ with $\{\Sigma,-p,-\Delta\}$.

[TABLE]

The amplitude then depends on four variables $\eta_{+}$, $\eta_{-}$, $\omega$ and $\lambda$, where $\omega$ and $\lambda$ can again only enter quadratically.

In analogy to (17), we alternatively use the hyperspherical variables $t$, $X$, $Y$, $Z$ defined by

[TABLE]

where we attached minus signs to $P$ and $p$ to comply with the Compton scattering kinematics. This corresponds to the Lorentz-invariant definitions

[TABLE]

and likewise for the internal variables $x$, $z$ and $\psi$ by taking dot products with the loop momentum $k^{\mu}$.

The left panel of Fig. 17 shows the resulting kinematic domain in $\eta_{+}$ and $\eta_{-}$ for $\omega=0$. As before, bound states and resonances appear at fixed $t$ below and above threshold $t=-1$, respectively. At a given $t$, the domain of the scattering equation where it is solved self-consistently is

[TABLE]

The onshell scattering amplitude $T_{\text{on}}(t,Y)$ corresponds to $X=-(1+t)$ and $Z=0$ and thus

[TABLE]

The onshell Mandelstam plane is then the line $\eta_{+}=-1$ in Fig. 17. In the half-offshell case the Mandelstam variables $s$ and $u$ become

[TABLE]

and the Mandelstam plane is expressed through the variables $\eta_{-}$ and $\lambda$. The exchange-particle poles correspond to $\eta_{-}=1-\beta^{2}\pm 2\lambda$; for $\lambda=0$, the intersection of the two poles at $s=u=\mu^{2}$ (cf. Fig. 16) would form a vertical line $\eta_{-}=1-\beta^{2}$ in Fig. 17.

If we also switch on the remaining variable $\omega$, then for $t>0$ the integration domain forms a cone around the $\eta_{+}$ axis, which is shown in the right panel of Fig. 17. This is so because from Eq. (49) one has

[TABLE]

For $t<0$ the integration domain forms a similar cone around the $\eta_{-}$ axis and with $\omega$ imaginary. Hence, to obtain a self-consistent solution one must first solve the equation inside these cones, and afterwards solve the equation once more by setting $X=-(1+t)$ and $Z=0$, where the half-offshell amplitude is integrated over to obtain the onshell scattering amplitude $T_{\text{on}}(t,\lambda)$.

To cast the scattering equation (33) into a form analogous to Eq. (20), we express it in hyperspherical variables:

[TABLE]

where $K$ and $G_{0}$ are the same as before,

[TABLE]

except that the argument $\widetilde{\Omega}$ also depends on the angle $\psi$:

[TABLE]

The innermost $\psi$ integration is no longer trivial but it can still be performed analytically:

[TABLE]

where

[TABLE]

The external kernel in the first line of Eq. (53) follows from replacing the appropriate momentum arguments.

The analysis of branch cuts is analogous to Sec. IV.2. We keep the integration domains of $z$ and $y$ inside their intervals in Eq. (49) and deform the integration contour in $\sqrt{x}$ only. In this way we cover the area $|\lambda|<|1+t|\Leftrightarrow|Y|<1$, i.e., the region between the lines $s=0$ and $u=0$ in the Mandelstam plane. The cuts in the integrand are as follows:

• (1)

The propagator poles encoded in $G_{0}(k,P)$ generate a cut $\sqrt{x}=\mathcal{C}_{1}=f_{\pm}(t,1,z)$; the function $f_{\pm}$ has been defined in Eq. (22).

• (2)

The exchange particle poles in $K(q,k)$ produce a cut $\sqrt{x}=\mathcal{C}_{3}=f_{\pm}(X,\beta,\widetilde{\Omega})$.

• (3)

The exchange-particle poles in the $s$ and $u$ channels, which are generated in the internal scattering amplitude $T(x,y,z,t)$ itself, produce cuts in the same place as those in $K(k,p)$, cf. Eq. (16), namely at $\sqrt{x}=\mathcal{C}_{2}=f_{\pm}(-(1+t),\beta,\sqrt{1-z^{2}}\,y)$.

• (4)

Finally, if we solve the equation once more in the final step to obtain the onshell scattering amplitude, the exchange-particle poles in $K(q,k)$ produce a cut $\sqrt{x}=\mathcal{C}_{4}=f_{\pm}(-(1+t),\beta,\sqrt{1-z^{2}}\,\Upsilon)$.

Since the third argument of $f_{\pm}$ always takes values in the interval $(-1,1)$, the cuts $\mathcal{C}_{2}$ and $\mathcal{C}_{4}$ lead to the same condition and do not need to be discussed separately. As before, the bound-state poles and cuts which are dynamically generated in the equation depend on $t$ only and we avoid them by discarding the imaginary axis $\text{Re}\sqrt{t}=0$.

The essential complication here is the appearance of two cuts, $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$, which both depend on the external point $t$ and must be avoided. As a consequence, the kinematically safe region where no contour deformation is necessary shrinks to the intersection of the two conditions

[TABLE]

Moreover, it is no longer possible to cover the entire complex $\sqrt{t}$ plane using a contour deformation in $\sqrt{x}$ only; to do so, one would have to deform contours in the remaining variables $z$ and/or $y$ as well. The cut $\mathcal{C}_{3}$, on the other hand, leads to the same condition as before, namely that the deformed path in $\sqrt{x}$ avoiding the cuts $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ must be chosen such that $\text{Re}\sqrt{x}$ and $|\sqrt{x}|$ never decrease along it. The details of the contour deformation procedure are given in Appendix A.

In practice the half-offshell scattering equation (53) turns again into a matrix-vector equation analogous to Eqs. (27) and (31):

[TABLE]

where the multi-indices $\sigma$, $\tau$ encode the dependence on the variables $X$, $Y$ and $Z$. To avoid convergence problems in the iterative solution in the vicinity of the poles, we solve the equation through matrix inversion:

[TABLE]

V.4 Results for the scattering amplitude

Following the steps above, we first solve the half-offshell amplitude $T(X,Y,Z,t)$ from its scattering equation (53) and obtain the onshell amplitude $T_{\text{on}}(t,Y)$ by solving the equation once more for $X=-(1+t)$ and $Z=0$. Finally we extract the partial-wave amplitudes $f_{l}(t)$ through

[TABLE]

and determine the amplitude on the second Riemann sheet from Eq. (44).

The integration domain $Y\in[-1,1]$ is the area enclosed by the lines $s=0$ and $u=0$ in the Mandelstam plane of Fig. 16. In the $t$-channel region ($t<-1$) the exchange-particle poles do not enter in the integration domain, whereas for $t>-1$ they do appear above a certain value of $t$,

[TABLE]

which follows from Eq. (34) as the intersection of $s=\mu^{2}$ and $u=0$. As discussed above, the exchange-particle poles only appear in the inhomogeneous term of the scattering equation. If we split the partial waves into

[TABLE]

then the exchange poles only contribute to the first term, whose integrand from Eqs. (53–54) takes the form

[TABLE]

with $B=\beta^{2}/\left[2(1+t)\right]-1$. The $f_{l}^{(K)}(t)$ can thus be obtained analytically:

[TABLE]

and so on for higher partial waves. These expressions have poles for $t=t_{P}$ and branch cuts for $t>t_{P}$. The dynamically generated bound-state and resonance poles, on the other hand, only appear in the second term $f_{l}^{(T-K)}(t)$ whose integrand $T_{\text{on}}-K_{\text{on}}$ is almost independent of $Y$ because the only singularities in that variable are the $s$- and $u$-channel cuts beginning at $t>\beta^{2}-1$.

Fig. 18 shows the leading partial wave $f_{0}(t)$ in the complex $\sqrt{t}$ plane for $\beta=4$ and four values of the coupling $c=\{1,3,5,7\}$. In this case the exchange poles are relatively far away from the displayed region ($t_{P}=3$) so that only the dynamical poles are visible. As before we plot the first Riemann sheet on the right half plane and the second sheet on the left. Recalling Fig. 10, bound states can only appear in the interval $6\lesssim c\lesssim 11$. One can clearly see that below $c\sim 6$ the ground state of the system is a virtual state with its pole on the imaginary $\sqrt{t}$ axis (or real $t$ axis) on the second sheet. With increasing coupling strength $c$ it moves up towards the threshold. At $c\sim 6$ the pole crosses over to the first sheet, where it becomes a bound state, and slides down along the axis towards the origin. For $c\gtrsim 11$, it becomes tachyonic and continues on the real $\sqrt{t}$ axis.

The resulting pole trajectory is shown in Fig. 18. For $c=1$ the location of the virtual pole is $\sqrt{t}\approx 0.40\,i$, which is compatible with the RVP estimate in Fig. 13 for the two input regions above threshold. Thus, the homogeneous BSE in combination with the RVP method did indeed predict the poles of scattering amplitude on the second sheet in the right ballpark. To determine the precise locations, however, we had to solve the scattering equation to have access to the unitarity relation (43).

The pattern in Fig. 19 repeats itself for the excited states. If the coupling increases further, eventually the first excited state appears as a virtual state and follows the same trajectory, until it becomes bound (from Fig. 10 at $c\sim 66$) and then again tachyonic ($c\sim 93$), followed by the second excited state and so on.

We should note that tachyonic bound states above a critical value of the coupling are not necessarily tied to the vacuum instability of the cubic interaction but they can also appear as truncation artifacts. The conditions for tachyonic poles are that (a) the eigenvalues of the system are monotonically decreasing functions of $t$ when passing from the timelike ($t<0$) to the spacelike side ($t>0$), and (b) that the coupling, when pulled out of the kernel, can be tuned independently without affecting the eigenvalues $\lambda_{i}(t)$. The first condition is a generic feature of BSEs since the eigenvalues of $\mathcal{K}(t)$ always inherit the falloff with $t$ encoded in the propagators; it means that in Fig. 10 the inverse eigenvalues continue to grow when extending the plots to the spacelike region on the left. Then, for a large enough coupling the condition (28) can always be met at some value $t>0$. The second feature, on the other hand, is special to our model in the sense that there is only one overall coupling that can be pulled out of the kernel. In general, when solving the DSEs for the propagators they also become functions of $c$ and so does the kernel beyond simple truncations. Increasing the coupling typically increases the self-energies and thus the inverse BSE eigenvalues, so that Eq. (28) may no longer have a solution for spacelike $t$. This is effectively what happens in QCD, where the coupling $\alpha_{s}$ in the BSE kernel cannot be dialled independently of the quark propagator; instead, the strong dependence of the propagator on $\alpha_{s}$ is manifest in dynamical chiral symmetry breaking and also changes its singularity structure in the process.

The results discussed so far for $\beta=4$ are generic and qualitatively also hold for different values of $\beta$. This can be inferred from Fig. 20, where the bound-state regions in the $(\beta,c)$ plane are plotted for the lowest three states of the system. Inside the lowest band the ground state is bound; for smaller values of $c$ it is a virtual state and for larger values it becomes tachyonic. The same pattern repeats itself for the excited states.

For smaller values of $\beta$ eventually also the exchange-particle poles from Eq. (65) become important. Fig. 21 exemplifies $f_{0}(t)$ for $\beta=2$ and $c=12$, which lies between the first and second band in Fig. 20. In this case the ground state has become tachyonic and is no longer visible in the plot; instead, the large structure on the first sheet is the exchange-particle pole at $t_{P}=\beta^{2}/4-1=0$ from Eq. (62). One can also see the first excited state on the second sheet, which has not yet become bound and is still a virtual state below threshold.

In the limit $\beta=0$ of a massless exchange particle, which is the Wick-Cutkosky model, all ‘normal’ eigenvalues in Fig. 10 would end at the threshold so that $1/\lambda_{i}(t=-1)=0$. In this case there are only bound states (and tachyons for large enough couplings) but no virtual states. The states that do not satisfy this property, i.e. $1/\lambda_{i}(t=-1)\neq 0$, are the so-called anomalous states Wick (1954); Cutkosky (1954); Nakanishi (1969); Ahlig and Alkofer (1999). To investigate what becomes of them when they cross the threshold remains the subject of future work.

VI Conclusions and outlook

We have investigated contour deformations as a tool to extract resonance properties from bound-state and scattering equations. As an example we have solved the homogeneous and inhomogeneous Bethe-Salpeter equation for a scalar model with a ladder-exchange interaction and we calculated the $2\to 2$ scattering amplitude from its scattering equation. We find that the model does not produce resonances above threshold but rather virtual states on the second Riemann sheet.

Our analysis was carried out in several steps. First, we employed contour deformations to access the kinematic regions associated with amplitudes in ‘Minkowski space’. We pointed out that there is no intrinsic difference between Minkowski and Euclidean space approaches; in both cases one needs to deform integration contours and the result must be the same. We used a Euclidean metric, which allowed us to formulate the four-dimensional integrals in a manifestly Lorentz-invariant way, and performed the contour deformations in the radial variable.

We demonstrated that the resonance information is in principle already contained in the homogeneous Bethe-Salpeter equation, because it determines an analytic continuation of the Bethe-Salpeter amplitude which in quantum field theory is only defined at a particular (real) bound-state pole. To extract that information, however, one needs analytic continuation methods. We employed the Schlessinger point or Resonances-via-Padé method, which is well-suited for situations with resonances above threshold. However, the scalar model we investigated here produces virtual states below threshold whose pole locations are more difficult to determine accurately by analytic continuation than nearby resonances.

To address this, we solved the half-offshell scattering equations for the $2\to 2$ scattering amplitude, from where we extracted the onshell scattering amplitude and employed the two-body unitarity relation which provides direct access to the second Riemann sheet. This allows us to determine the pole positions precisely and confirms that the resonance poles of the model are indeed virtual states below the threshold.

Our analysis leaves room for several future investigations. One can study scattering amplitudes with complex propagator poles, the nature of anomalous states, or extend the approach to three-body systems. Our determination of the cut structure does not depend on the scalar nature of the system and can be applied without changes to the scattering of particles with any spin, such as $NN$ or $N\pi$ scattering, as long as self-energies can be neglected. Moreover, contour deformations provide a general way to overcome singularity restrictions for example in QCD, which can help to access properties of resonances but also highly excited states, form factors in the far spacelike and timelike regions, or general matrix elements in kinematic regions where a naive Euclidean integration contour is no longer applicable.

Acknowledgements

This work was funded by the Fundação para a Ciência e a Tecnologia (FCT) under Grants IF/00898/2015 and UID/FIS/00777/2013.

Appendix A Contour deformation with two cuts

In this appendix we provide details on our contour deformation procedure in the calculation of the scattering amplitude. As discussed in Sec. V.3, the complication in this case arises from the appearance of the two cuts

[TABLE]

in the complex $\sqrt{x}$ plane which must be circumvented, where the function $f_{\pm}$ is given by

[TABLE]

With $\sqrt{t}$ in the upper right quadrant, $\text{Im}\,[-(1+t)]<0$ and therefore $\sqrt{-(1+t)}$ is in the lower right quadrant. Hence we can write

[TABLE]

where $a,b,c,d>0$ and

[TABLE]

with $A>0$, $|B|<A$.

The typical situation in the complex $\sqrt{x}$ plane is illustrated in Fig. 22. The point $\sqrt{t}$ in the upper right quadrant defines the vertical line $\mathcal{V}_{1}$ and a circle with radius $R_{1}=\sqrt{a^{2}+b^{2}}$. The point $\sqrt{-(1+t)}$ defines the line $\mathcal{V}_{2}$ and a circle with radius $R_{2}=\sqrt{c^{2}+d^{2}}=\sqrt{A}$. The two relevant cuts are then given by

[TABLE]

where $z\in(-1,1)$ is a generic variable parametrizing the cuts and the opposite branches are not relevant for the further discussion. At $z=1$, the cuts start on $\mathcal{V}_{1}$ and $\mathcal{V}_{2}$ at the respective points

[TABLE]

As $z$ decreases, $|\mathcal{C}_{1}|$ and $|\mathcal{C}_{2}|$ increase. At $z=0$, $\mathcal{C}_{1}$ passes through $\sqrt{-(1+t)}$ at $\mathcal{V}_{2}$. At the point $\sqrt{t^{\ast}}$, again at $\mathcal{V}_{1}$, the cut $\mathcal{C}_{1}$ leaves the circle with radius $R_{1}$ and continues until it reaches its endpoint $-a-i(b+1)$. The cut $\mathcal{C}_{2}$ in Fig. 22 already starts outside of its circle with radius $R_{2}$ and progresses until its endpoint $-c+i(\beta+d)$.

A proper integration path $\sqrt{x}=\mathcal{I}$ in Fig. 22 would start at the origin, pass through the region between $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ and return to the real axis afterwards. In addition, it must satisfy the constraints arising from the exchange-particle cut $\mathcal{C}_{3}$, namely that $\text{Re}\,\mathcal{I}$ and $|\mathcal{I}|$ must never decrease along the path. For example, once $\mathcal{I}$ has picked up a nonvanishing real part it is no longer possible to return to the imaginary axis.

Depending on the alignment of the cuts, different path deformations may be necessary. To this end we identify the following regions:

• •

$\mathcal{C}_{1}$ never crosses the positive real axis if $|\text{Im}\sqrt{t}|<1$, cf. Eq. (58), which entails $b<1$.

• •

$\mathcal{C}_{2}$ never crosses the real axis if $|\text{Im}\sqrt{-(1+t)}|<\beta$, which entails

[TABLE]

• •

The line $\mathcal{V}_{1}$ is on the left of $\mathcal{V}_{2}$ if

[TABLE]

which corresponds to

[TABLE]

and implies $R_{1}>R_{2}$. If $\mathcal{V}_{1}$ is on the right, the situation is reversed and $a>c$.

• •

If $a<c$, a contour deformation is possible as long as $\mathcal{P}_{2}$ does not touch $\mathcal{C}_{1}$, which leads to the condition

[TABLE]

with $G=bc+3ad$.

• •

If $a>c$, a contour deformation is possible as long as $\mathcal{P}_{1}$ does not touch $\mathcal{C}_{2}$, which corresponds to

[TABLE]

Based on these constraints, we divide up the complex $\sqrt{t}$ plane into regions which are distinguished by a different contour deformation procedure. They are plotted in Fig. 23 for three values of the mass ratio $\beta$, where the constraints $b=1$ and $d=\beta$ correspond to the dashed lines and $a=c$ to the solid line. Fig. 24 illustrates exemplary configurations for each region:

Region I ($b<1$, $d<\beta$): No contour deformation is necessary because none of the cuts cross the real axis, as exemplified in the upper left panel of Fig. 24.

Region Ia: The simplest extension is to test whether the cuts $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$ can be passed by a straight line starting at the origin, which is true if $\arg\mathcal{P}_{2}>\arg\mathcal{P}_{1}$. Because both points are on the right half plane, we can drop the $\arctan$ and the condition becomes

[TABLE]

An exemplary configuration is plotted in the bottom left panel of Fig. 24. The optimal integration path is the average of the two lines connecting the origin with $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$, respectively, where we integrate up to a radius which is larger than both $R_{1}$ and $R_{2}$ before going back to the real axis. The resulting region in the complex $\sqrt{t}$ plane is shown by the dashed opaque areas in Fig. 23. They enclose region I but also extend into the other regions discussed below.

Region II ($b<1$, $d>\beta$, $a>c$): Only $\mathcal{C}_{2}$ crosses the real axis, as illustrated in the top center panel of Fig. 24. Here the path deformation is still simple because we only need to circumvent $\mathcal{C}_{2}$ and can return to the real axis afterwards.

To find the optimal integration path, consider the function $f_{+}((c-id)^{2},\alpha,z)$. For $\alpha=\beta$ this is just the cut $\mathcal{C}_{2}$. As we decrease $\alpha$, we deform its shape; $\mathcal{P}_{2}$ slides down along the line $\mathcal{V}_{2}$ until the contour eventually touches the cut $\mathcal{C}_{1}$ at the point $\mathcal{P}_{1}$. The corresponding value of $\alpha$ is given by the r.h.s. of Eq. (76). Decreasing $\alpha$ further, the contour crosses the cut $\mathcal{C}_{1}$ and for $\alpha=0$ it becomes the half-circle starting at $c-id$. Therefore, the choice

[TABLE]

which averages over the cut $\mathcal{C}_{2}$ and the contour touching $\mathcal{P}_{1}$, is guaranteed to lie between the cuts $\mathcal{C}_{1}$ and $\mathcal{C}_{2}$.

Our integration path thus proceeds from the origin to the point

[TABLE]

on $\mathcal{V}_{2}$, goes along $\mathcal{I}_{+}(z)$ until it reaches the real axis at

[TABLE]

and from there continues to $\sqrt{x}\to\infty$. Because the real part and modulus of the contour always increase along that path, the constraints from the exchange-particle cut $\mathcal{C}_{3}$ are satisfied. In this way we can cover the entire region II without crossing any branch cut.

Region III ($b>1$, $d<\beta$, $a<c$): This is the opposite case where only $\mathcal{C}_{1}$ crosses the real axis (Fig. 22 and top right panel in Fig. 24). Here we employ the function $f_{-}((a+ib)^{2},\alpha,z)$. For $\alpha=1$ this is the cut $\mathcal{C}_{1}$ and for $\alpha=0$ it becomes the circle passing through $a+ib$. As we decrease $\alpha$, $\mathcal{P}_{1}$ will slide up along the line $\mathcal{V}_{1}$ until the contour touches the point $\mathcal{P}_{2}$ of the cut $\mathcal{C}_{2}$. The resulting value for $\alpha$, however, is not the one in Eq. (75), which was obtained by equating $\mathcal{P}_{2}$ and $\mathcal{C}_{1}$ with $\alpha=1$. Instead, we need to keep $\alpha$ general which results in

[TABLE]

where $G=bc+3ad$ and

[TABLE]

Here we took care of another subtlety: If $|\mathcal{P}_{2}|>R_{1}$ (which is the case in Fig. 24), then the cut $\mathcal{C}_{2}$ no longer reaches $f_{-}((a+ib)^{2},\alpha,z)$ which remains within the circle with radius $R_{1}$. In that case we replace $\mathcal{P}_{2}$ by the intersection of that circle with $\mathcal{V}_{2}$, which from $|\mathcal{P}_{2}|=R_{1}$ corresponds to $\beta=d+\sqrt{a^{2}+b^{2}-c^{2}}$ and leads to the above definition. The analogous situation for region II can never happen because $|\mathcal{P}_{1}|<R_{2}$.

The integration path then proceeds from the origin to

[TABLE]

on $\mathcal{V}_{1}$, goes along $\mathcal{I}_{-}(z)$ until it reaches the real axis at

[TABLE]

and from there continues to $\sqrt{x}\to\infty$. The constraints from the exchange-particle cut $\mathcal{C}_{3}$ are again satisfied, so we can cover the entire region III without crossing any branch cut.

Region IV ($b>1$, $a>c$): This case is more complicated and illustrated in the bottom center of Fig. 24. Here $\mathcal{C}_{1}$ crosses the real axis; whether $\mathcal{C}_{2}$ also crosses or not is not relevant. The integration path is similar as in Region II, except that one must still circumvent the cut $\mathcal{C}_{1}$ after passing the real axis. Because the real part of the contour must never decrease, one can stay on $\mathcal{I}_{+}(z)$ only as long as it does not become vertical. At that point we switch to another contour that brings us back to the real axis. In Fig. 23 one can see that Region IV only appears for $\beta>1$.

Region V ($d>\beta$, $a<c$) is the opposite case shown in the bottom right of Fig. 24. Here $\mathcal{C}_{2}$ crosses the real axis, whereas $\mathcal{C}_{1}$ may or may not intersect with the real axis. The integration path is similar to Region III; once again, one must switch contours if $\mathcal{I}_{-}(z)$ becomes vertical. One can see in Fig. 23 that Region V only covers a small portion of the complex $\sqrt{t}$ plane and it also only appears for $\beta<1$.

The union of all Regions I–V covers the entire accessible area in the complex $\sqrt{t}$ plane defined by Eqs. (75–76). Region Ia also extends into the other areas and provides a useful cross check since in the overlap regions one can test different contour deformation strategies.

In principle there are many possible ways to generalize or even automatize the contour deformation. For example, outside of Region Ia one could use the two rays connecting the origin with $\mathcal{P}_{1}$ and $\mathcal{P}_{2}$ to define inner and outer zones and construct integration paths accordingly. Alternatively, one could connect the tangent on one cut with the endpoint of the other and then integrate along straight lines. Such techniques can be useful in more general situations, for example when the ingredients of the equation are nonperturbative and contain not only tree-level poles but also poles or cuts in the complex plane.

Finally, what turned out very useful was a pole analysis in the complex $z$ plane, which is the angular integration variable in the BSE (20) or scattering equation (53). The $z$ integration goes from $z=-1$ to $z=+1$ and for a given path deformation in $\sqrt{x}$ the poles of the integrands never cross the integration path in $z$. Still, the closer the path in $\sqrt{x}$ comes to one of those singularities (e.g. in region Ia of Fig. 24, when it passes through the cuts), the closer the poles in the complex $z$ plane come to the real axis and therefore the integrand as a function of $z$ varies strongly in the vicinity of those poles. Here it is useful to perform an adaptive integration by splitting the $z$ integration into intervals and accumulate the grid points around the nearest singularities as sketched in Fig. 25. For example, for the propagator poles corresponding to $\mathcal{C}_{1}$ their locations are determined by solving $\sqrt{x}=f_{\pm}(t,1,z)$ for $z$. In that way we typically gain a factor $\sim 10^{2}\dots 10^{3}$ in CPU time while maintaining the same numerical accuracy.