one may analyze it using the final and initial invariant mass squared of
pairs, $\sigma_{p},\sigma_{k}$. In practical applications, we find that the
momentum representation proves more useful in the study of analytic
continuation. It is because the OPE cuts take a simpler shape in this form.
They wrap around the origin of the complex $q$ plane and have associated
parity copies allowing for less problematic choices of the deformed
integration contours. On the other hand, the invariants-space OPE cuts follow
the movable upper integration limit and have a more complicated, fishing-hook-
like shape.

However, in some cases, the invariant-space equations are simpler to
manipulate. One such case is a derivation of the positions of the OPE branch
points. Ultimately, it is desirable to have two representations since one can
prove more useful than the other in analyses concerned with different physical
systems and the LQCD data. In particular, bound-state and resonance poles
occur in $\mathcal{M}_{2}$ at fixed values of the two-body invariant mass,
making it an intuitively better variable to consider. Moreover, the variables
$\sigma_{p},\sigma_{k}$ are Lorentz invariants and do not change for different
values of the total invariant mass $s$ in contrary to momenta $(p,k)$. In this
appendix, we concisely present the invariants representation of the ladder
equation focusing on the analytical structure of the building blocks of the
equation.

The $S$-wave projected ladder equation, Eq. (8), is written in terms of
Lorentz invariants $\sigma_{p},\sigma_{k}$, as,

$\displaystyle
d_{S}(\sigma_{p},\sigma_{k})=-G_{S}(\sigma_{p},\sigma_{k})-\int\limits_{0}^{\sigma_{\rm
max}}d\sigma_{q}\,K(\sigma_{p},\sigma_{q})\,d_{S}(\sigma_{q},\sigma_{k})\,.$
(61)

Variable $\sigma_{q}$ is the invariant mass squared of the intermediate pair
in the OPE process. The integration kernel is,

$\displaystyle
K(\sigma_{p},\sigma_{q})=\frac{1}{2\pi}G_{S}(\sigma_{p},\sigma_{q})\tau(s,\sigma_{q})\mathcal{M}_{2}(\sigma_{q})\,,$
(62)

where the implicit $s$ dependence is assumed. The integration is performed in
the interval $[0,\sigma_{\rm max}]$, where $\sigma_{\rm
max}=(\sqrt{s}-m)^{2}$. It corresponds to $q=0$, while $\sigma_{q}=0$
corresponds to $q=q_{\rm max}$ in the integral of Eq. (8).

The integration kernel contains three objects. The three-body phase space is,

$\displaystyle\tau(\sigma_{q})=\frac{\lambda^{1/2}(s,\sigma_{q},m^{2})}{8\pi
s}\,.$ (63)

It has an explicit pole at $s=0$ and the branch points at $\sigma_{\tau
1}=(\sqrt{s}-m)^{2}$ and $\sigma_{\tau 2}=(\sqrt{s}+m)^{2}$. We orient both
associated cuts to the right. In particular, for real $s$, this results in a
single branch cut running between the two branch points. We note that the
upper integration limit coincides with the former branch point. The two-body
amplitude is given in Eq. (13) as a function of $\sigma_{q}$. As can be seen,
it has a left-hand cut at $\sigma_{q}=0$ and a right-hand cut at
$\sigma_{q}=4m^{2}$ required by the unitarity. It also develops a pole on the
first complex sheet at $\sigma_{b}$. The $\mathcal{M}_{2}$ amplitude can be
rewritten in a “propagator” form that makes the presence of the pole explicit,

$\displaystyle\mathcal{M}_{2}(\sigma_{q})=\frac{R(\sigma_{q})}{\sigma_{q}-\sigma_{b}-i\epsilon_{b}}\,,$
(64)

where the residue,

$\displaystyle
R(\sigma_{q})=-(32\pi)^{2}\,\sigma_{q}\left(\mathcal{K}_{2}^{-1}+i\rho\right)\,.$
(65)

In Eq. (64), we included infinitesimal $i\epsilon_{b}$ in the denominator
(different than the $i\epsilon$ in the OPE amplitude) to shift the pole
position above the real $\sigma_{q}$ axis. It is necessary when solving for
the physical amplitude, as discussed in Ref. Jackura  _et al._ (2021), and is
equivalent to the integration contour deformation. For $s=s_{\varphi b}$, the
upper limit of the integration coincides with the pole, leading to the
unitarity branch point in the ladder solution $d_{S}(\sigma_{p},\sigma_{k})$.
For $\operatorname{Re}s<s_{\varphi b}$, the integration interval does not
coincide with the singularities of $\mathcal{M}_{2}$.

Figure 18: The OPE amplitude $G_{S}(\sigma_{p},\sigma_{q}+i\delta)$ in units
of $1/m^{2}$, evaluated slightly above the real $\sigma_{q}$ axis. The
kinematic parameters are $s/m^{2}=8.3$ and $\sigma_{p}=\sigma_{b}$ for
$ma=16$. Infinitesimal $\delta=10^{-4}$. Singularities described in the text
and defined in Eqs. (71), (75), (75) are highlighted with dashed lines. Figure
19: Cut structure of the imaginary part of the OPE,
$\operatorname{Im}G_{S}(\sigma_{q},\sigma_{k})$, in the complex $\sigma_{q}$
plane for typical value of kinematical variables. We use hard cut-off and plot
the amplitude in units of $1/m^{2}$. We set $ma=16$ and consider
$\operatorname{Re}s<s_{\varphi b}$. Two-body invariant mass,
$\sigma_{k}=\sigma_{b}\approx 3.984m^{2}$. The red circle corresponds to
$\sigma_{-}$ and the blue to the $\sigma_{+}$ branch point. The OPE cut (black
curve) runs between these two points. We indicate the upper integration limit
$\sigma_{\rm max}=(\sqrt{s}-m)^{2}$ by a yellow point. It is enclosed by the
cut and contour deformation is required. We present four cases: (a) Purely
real kinematical parameters, for which the cut takes a shape of a circle
crossed by a line; (b) Purely real kinematical parameters, with $s$ closer to
$s_{\varphi b}\approx 8.976m^{2}$, for which the cut moves the right and
shrinks; (c) The non-zero, positive imaginary part of $s$, for which the cut
opens, moves upward and resembles a fishing hook; (d) Same as previously, but
for larger $\operatorname{Re}s$, the structure shrinks and $\sigma_{+}$
approaches $\sigma_{\rm max}$. In particular, when $\sigma_{k}=m(m+\sqrt{s})$,
the OPE branch point coincides with the upper integration limit.

The $S$-wave projection of the OPE amplitude is,

$\displaystyle
G_{S}(\sigma_{p},\sigma_{k})=-\frac{s\,H(\sigma_{p},\sigma_{k})}{\lambda^{1/2}(s,\sigma_{p},m^{2})\lambda^{1/2}(s,\sigma_{k},m^{2})}\log\left(\frac{z(\sigma_{p},\sigma_{k})-\lambda^{1/2}(s,\sigma_{p},m^{2})\lambda^{1/2}(s,\sigma_{k},m^{2})}{z(\sigma_{p},\sigma_{k})+\lambda^{1/2}(s,\sigma_{p},m^{2})\lambda^{1/2}(s,\sigma_{k},m^{2})}\right)\,,$
(66)

where the $z(\sigma_{p},\sigma_{k})$ function is defined as,

$\displaystyle
z(\sigma_{p},\sigma_{k})=2\,s\,(\sigma_{k}+i\epsilon)-(s+\sigma_{k}-m^{2})(s+m^{2}-\sigma_{p})\,.$
(67)

This representation follows from Eq. (II.1), where one expresses the external
spectator’s momenta through the relation given in Eq. (4). The function
$H(\sigma_{p^{\prime}},\sigma_{p})$ is the smooth/ hard regularization scheme,
as described in the paragraph containing Eq. (32). In Fig. 18, we present
$G_{S}(\sigma_{p},\sigma_{q})$ as a function of real-valued $\sigma_{q}$,
choosing the smooth cut-off, defined in Eq. (32).

The OPE amplitude, Eq. (66), considered as a function of $\sigma_{p}$ for
fixed $s$ and $\sigma_{k}$, has two logarithmic branch points connected with a
cut. Its parametrization is obtained from the condition,

$\displaystyle
z(\sigma_{p},\sigma_{k})+x\,\lambda^{1/2}(s,\sigma_{p},m^{2})\lambda^{1/2}(s,\sigma_{k},m^{2})=0\,,$
(68)

which is an equation satisfied by the pole positions of the integrand in the
right-hand side of Eq. (II.1). Solving for $\sigma_{p}$ yields,

$\displaystyle\sigma_{\text{cut},\pm}(s,\sigma_{k},x)=(s+m^{2})+\frac{2s\sigma_{k}(s+\sigma_{k}-m^{2})\pm\sqrt{sx^{2}\lambda(s,\sigma_{k},m^{2})}\sqrt{4m^{2}B_{x}-\sigma_{k}B_{1}}}{B_{x}}\,,$
(69)

where function

$\displaystyle B_{x}\equiv
B_{x}(s,\sigma_{k})=(x^{2}-1)\lambda(s,\sigma_{k},m^{2})-4s\sigma_{k}\,.$ (70)

We have set $\epsilon=0$. The above formula is analogous to the momentum-space
parametrization of Eq. (35). Equation (69) is symmetric with respect to
$x\to-x$ change; thus, we can take $x$ in the $[0,1]$ interval. Two solutions
labeled “$\pm$” do not describe two “parity copies” of the cut, but two
smoothly connected halves of the same cut attached to a different branch
point. We call the branch points $\sigma_{\pm}$, and obtain them from the
above parametrization by setting $x=1$,

$\displaystyle\sigma_{\pm}=\sigma_{\text{cut},\pm}(s,\sigma_{k},1)=\frac{1}{2}(s-\sigma_{k}+3m^{2})\pm
16\pi\,\rho(\sigma_{k})\,\lambda^{1/2}(s,\sigma_{k},m^{2})\,.$ (71)

The expression for $\sigma_{\pm}$ has singularities in $\sigma_{k}$ and $s$
since the second term of Eq. (71) contains both the triangle function and two-
body phase space. They have practical consequences for the implementation of
the integral equation solution. For example, the ordering between
$\operatorname{Re}\sigma_{\pm}$ (i.e., which point is on the left and which on
the right in the complex plane) depends on the relative value of
$\operatorname{Im}s$ and $\operatorname{Im}\sigma_{k}$. It affects the choice
of the integration contour; considering only real, positive values of
$\sigma_{k}<4m^{2}$ and $\operatorname{Re}s<s_{\varphi b}$, the
$\lambda^{1/2}$ has a cut in $s$ below $(\sqrt{\sigma_{k}}+m)^{2}$. For $s\to
s^{*}$, the real parts of branch points transform into each other,
$\operatorname{Re}\sigma_{-}\leftrightarrow\operatorname{Re}\sigma_{+}$.

Similarly to the momentum-representation OPE amplitude, for
$\operatorname{Im}s=0$, the branch cut wraps around the real $\sigma_{p}$
axis, resulting in the circular cut, as seen in Fig. 19. It occurs when $s$ is
decreased below the value of $s$ given in Eq. (45), at which the branch point
$\sigma_{+}$ collides with $\sigma_{\rm max}=(\sqrt{s}-m)^{2}$. The cut
encloses the upper integration limit $\sigma_{\rm max}$. For a non-zero
imaginary part of $s$ or $\sigma_{k}$ (or non-zero $\epsilon$), the circle
“opens”.

One finds the point where the cut passes the real axis by looking for the
solution of condition (68) with a vanishing imaginary part. We can rewrite it
as,

$\displaystyle\frac{z(\sigma_{p},\sigma_{k})}{\lambda^{1/2}(s,\sigma_{p},m^{2})\lambda^{1/2}(s,\sigma_{k},m^{2})}+x=0\,.$
(72)

Using the fact that $x$ is purely real, the crossing in the real axis
satisfies,

$\displaystyle\operatorname{Im}\left[z(\sigma_{p},\sigma_{k})\left(\lambda^{1/2}(s,\sigma_{p},m^{2})\lambda^{1/2}(s,\sigma_{k},m^{2})\right)^{*}\right]=0\,.$
(73)

Assuming real $s<(\sqrt{\sigma_{k}}+m)^{2}$, $\sigma_{p}<\sigma_{\rm max}$ and
$\sigma_{k}$, we simplify it to,

$\displaystyle\operatorname{Re}z(\sigma_{p},\sigma_{k})=0\,,$ (74)

by noticing that the
$\lambda^{1/2}(s,\sigma_{p},m^{2})\lambda^{1/2}(s,\sigma_{k},m^{2})$ factor is
purely imaginary. Thus, we find that the circular cut crosses the real
$\sigma_{p}$ axis at,

$\displaystyle\sigma_{c1}=\frac{(s-m^{2})(s+m^{2}-\sigma_{k})}{(s-m^{2}+\sigma_{k})}\,.$
(75)

For complex $s$, as $\operatorname{Im}s\to 0$, the cut approaches the real
$\sigma_{p}$ axis at another point, which we call $\sigma_{c2}$. Referring to
panel (a) of Fig. 19 for illustration, it is the point where the line and the
circle cross each other. To express it in terms of $s$ and $\sigma_{k}$, we
write $\sigma_{p}=\sigma_{c2}+i\delta$, where $i\delta$ is a positive,
infinitesimal imaginary part. It constitutes a parametrization of the line
tangent to the cut near the real $\sigma_{p}$ axis. Again, we start from the
condition (73). For real $s<(\sqrt{\sigma_{k}}+m)^{2}$, this becomes,

$\displaystyle\delta\,(s+\sigma_{k}-m^{2})\,\operatorname{Im}\lambda^{1/2}(s,\sigma_{p},m^{2})+(2s\sigma_{k}-(s+\sigma_{k}-m^{2})(s+m^{2}-\sigma_{c2}))\operatorname{Re}\lambda^{1/2}(s,\sigma_{p},m^{2})=0\,,$
(76)

where this time $\lambda^{1/2}(s,\sigma_{p},m^{2})$ has in general non-zero
real and imaginary parts. We expand the triangle function around $\delta=0$,

$\displaystyle\lambda^{1/2}(s,\sigma_{p},m^{2})=\lambda^{1/2}(s,\sigma_{c2},m^{2})-\frac{i\delta\left(m^{2}+s-\sigma_{c2}\right)}{\lambda^{1/2}(s,\sigma_{c2},m^{2})}+\mathcal{O}(\delta^{2})\,,$
(77)

which, assuming $\lambda^{1/2}(s,\sigma_{c2},m^{2})$ is purely imaginary,
leads to,

$\displaystyle-2s\left(-2m^{4}+m^{2}(2s+\sigma_{k})+\sigma_{k}(\sigma_{c2}-s)\right)+\mathcal{O}(\delta)=0\,.$
(78)

Neglecting terms of order $\delta$, the solution of the equation becomes,

$\displaystyle\sigma_{c2}=\frac{\left(m^{2}-s\right)\left(2m^{2}-\sigma_{p}\right)}{\sigma_{p}}\,.$
(79)

We note that $\sigma_{c1},\sigma_{c2}$ correspond to $q_{c1}$ and $q_{c2}$
given in Eqs. (39), (40) The generalized values, $q^{\prime}_{c1,2}$, provided
in App. C, correspond to the points where the OPE branch cut, considered in
the complex $\sigma_{q}$ plane, crosses a line
$\operatorname{Im}(\sqrt{s}-m)^{2}=$ const. Control over the functional form
of those points is essential when preparing the deformed integration contour,
which enters the closed circle for $\operatorname{Im}s=0$ through the
$\sigma_{c2}$ or $q_{c2}$.

##  Appendix B Short introduction to analytic continuation

This appendix should serve as a pedagogical summary of concepts used in Sec.
IV, where we discuss an analytic continuation of the ladder equation. It is
based on Refs. Cohen (2007); Eden  _et al._ (1966); Burkhardt (1969); Lang
(1985) which may be consulted for more details.

To understand our treatment of the integral equation, it is beneficial to
consider a simpler case of an analytic continuation of a complex integral. We
define a generic,

$\displaystyle
I(z)=\\!\\!\\!\int\limits_{\mathcal{C}(w_{1},w_{2})}\\!\\!\\!f(w,z)\,dw\,,$
(80)

where the integrand $f(w,z)$ is a complex function of argument $w$ and depends
on a complex parameter $z$. Integration is performed over a path
$\mathcal{C}(w_{1},w_{2})$ which starts at $w_{1}=w_{1}(z)$ and ends at
$w_{2}=w_{2}(z)$. As indicated, these two points can also depend on $z$. The
homogeneous term of the ladder equation, Eq. (11), has an analogous form;
however, we do not know the equivalent of $f(w,z)$ beforehand, since $d_{S}$
is an unknown of the integral equation.

If we know the analytic structure of $f(w,z)$, we can infer the analytic
structure of $I(z)$. In general, singularities of $I(z)$ appear for those
values of $z$ for which: a) $f(w,z)$ has explicit, $w$-independent
singularities in $z$; b) $z$-dependent singularity in $w$ coincides with the
lower limit of integration, $w_{1}$; c) $z$-dependent singularity in $w$
coincides with the upper limit of integration, $w_{2}$; d) two movable
singularities of $f(w,z)$, pinch the integration contour; e) movable
singularities of $f(w,z)$ require contour deformation to complex infinity. We
note it is sufficient to know singularities of $f(w,z)$ to establish
singularities of $I(z)$ and not the value of $f(w,z)$ at every point of the
complex plane.

We illustrate this with a typical example of a real integral,

$\displaystyle I(x)=\int_{-1}^{1}\frac{dw}{w-x}\,,$ (81)

where the integration variable lies on the real axis between $w_{1}=-1$ and
$w_{2}=+1$. For $x$ outside of the integration range, we can easily evaluate
the integral and obtain,

$\displaystyle I(x)=\log\left(\frac{x-1}{x+1}\right)\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\
x\in(-\infty,-1)\cup(1,\infty)\,.$ (82)

The integral is not defined for $x\in[w_{1},w_{2}]$ due to the pole
singularity at $w=x$. However, having the explicit functional form, given in
Eq. (82), it is possible to assign a meaning to this function in this range.
Namely, we promote the real $I(x)$ to a function $I(z)$ of a complex variable
$z$ which is equal to $x$ on the real axis. Due to the multi-valued nature of
the complex logarithm, $I(z)$ has two branch points, at $w_{1}$ and $w_{2}$,
and two associated cuts. These can be chosen arbitrarily, corresponding to
different definitions of the function on the first Riemann sheet. For example,
we can align both cuts with the real axis and orient them to the right, which
results in a single branch cut in the interval $[w_{1},w_{2}]$. For this
choice, the function is undefined on this short segment of the real axis,
which is clear since the original integral in (81) was ill-defined there.

However, the cuts can be oriented in other directions, e.g., to cover
$(-\infty,-1)\cup(1,\infty)$, such that the function has a well-defined value
$I(x)$ for $-1<x<1$. To establish a relation between the complex function
$I(z)$ with its cuts moved away from $[-1,+1]$ and the original defining
integral, we can promote it to a complex integral along a general complex
contour $\mathcal{C}(-1,+1)$. The chosen integration path determines the cut
structure of the resulting $I(z)$. The pole of the integrand leads to
singularity only if it coincides with the integration path; thus, if the
integration contour avoids the interval $[-1,+1]$, the integral is well-
defined there. The function varies continuously as we cross the $[-1,+1]$
interval vertically and becomes equal to its value on the nearest Riemann
sheet of the “principal” definition, according to the Cauchy theorem, see Fig.
20.

Figure 20: Left: Imaginary part of the integral $I(x)$ in Eq. (81), defined as
a complex function. Right: Corresponding integration path
$\mathcal{C}(-1,+1)$. Continuation of the function to the real interval
$[-1,+1]$ requires contour deformation. Function (b) is given there by the
function (a) evaluated on the Riemann sheet connected to the bottom half-
plane.

It is easy to see that in the example of Eq. (81), we encounter cases (b) and
(c). The $x$-dependent pole of $f(w,x)=1/(w-x)$ coincides with $w_{1}$ and
$w_{2}$, resulting in the two branch points of the logarithm at these points.
These branch points have a fixed position that cannot be altered by contour
deformation, as every deformed path must begin and end at the same locations.
Consequently, the shape of the new branch cut is determined by the deformed
integration path.

In addition to illustrating branch cuts emergence, the above example suggests
how a contour deformation allows one to extend the definition of $I(z)$ to a
point $z_{0}$ where it was originally ill-defined. In other words, contour
deformation defines the analytic continuation of $I(z)$ to a new region of the
complex plane. This fact is well-known in the general S-matrix theory Eden
_et al._ (1966) and has been widely used in the phenomenological studies of
the scattering processes, e.g., see Ref. Mai (2021); Albaladejo  _et al._
(2022); Doring  _et al._ (2009).

In Secs. III and IV, we explain how singularities of the solution
$d_{S}(p,k)$, considered as a function of $s$ emerge from conditions (a)-(e)
applied to the integration kernel $K(p,q)$ and the unknown function
$d_{S}(q,k)$. For example, $K$ contains the two-body amplitude
$\mathcal{M}_{2}$ that exhibits an $s$-dependent, complex pole, $q=q_{b}$.
Collision of $q_{b}$ with the lower integration endpoint, $q=0$, leads to the
unitarity branch point at $s=s_{\varphi b}$ [condition (b)]. Collision of the
OPE branch point $p_{+}$ with $q=0$ leads to a branch point at
$s=(m^{2}-\sigma_{k})^{2}/m^{2}$ [condition (b)]. Moreover, the final
amplitude inherits an explicit “short” cut from the OPE amplitude in the
inhomogeneous term, considered as a function of $s$ for fixed $(p,k)$
[condition (a)].

To summarize, although we do not know the solution $d_{S}(p,k)$ the defining
integral equation can be used to infer singularities of the amplitude even
without solving it. Moreover, the ladder amplitude can be analytically
continued to the kinematic regions of interest via the contour deformation.

##  Appendix C Numerical methods

In this section, we describe numerical methods that were used to obtain
analytically continued solutions presented in Sec. V. Similarly to Ref.
Jackura  _et al._ (2021) we employ the Nyström method Nyström (1930); Delves
and Mohamed (1988), i.e., we discretize momentum variables and rewrite the
problem as an algebraic equation. The methods presented below are applicable
in more general studies of the three-body integral equations. They are
relatively well-known but we describe them here to increase the reader’s ease
in reproducing the results presented in this study.

###  C.1 Definition of the deformed contour

For reader’s convenience, below we reproduce the partial-wave projected,
amputated ladder equation, Eq. (8),

$\displaystyle
d_{S}(p,k)=-G_{S}(p,k)-\int\limits_{0}^{q_{\text{max}}}dq\,K(p,q)\,d_{S}(q,k)\,.$
(83)

where the integration kernel, $K(p,q)$ is defined in Eq. (10). We indicated
the finite range of the integration with the upper limit $q_{\text{max}}$,
which is defined by the cut-off function, Eq. (32). The motion of the OPE cuts
in the complex $p$ and $q$ planes necessitates contour deformation in Eq. (83)
for a large range of values of $s$ and $p$. In the following, we consider
values of kinematic variables $(p,k;s)$ for which the real-$q$ axis is crossed
by a cut. To compute the solution, we deform the integration path,

$\displaystyle\int\limits_{0}^{q_{\text{max}}}dq\longrightarrow\int_{\mathcal{C}}=\int_{0}^{1}dt\,\gamma^{\prime}(t)\,.$
(84)

The complex contour $\mathcal{C}$ is defined by a parametrization
$q=\gamma(t)$ where real parameter $t\in[0,1]$. For a given set of kinematic
variables, the curve has fixed endpoints, $\gamma(0)=0$, and
$\gamma(1)=q_{\text{max}}$. The integral equation becomes,

$\displaystyle
d_{S}(p,k)=-G_{S}(p,k)-\int_{0}^{1}dt\,\gamma^{\prime}(t)\,K(p,\gamma(t))\,d_{S}(\gamma(t),k)\,.$
(85)

The Nyström method is applied to the ladder equation in the above form.

Node | Im$(s)\leq$0 | Im$(s)>$0
---|---|---
$q_{0}$ | 0 | 0
$q_{1}$ | $-\frac{2}{3\sqrt{2}}(1+i)|q_{c2}|$ | $-\frac{2}{3\sqrt{2}}(1+i)|q_{c2}|$
$q_{2}$ | $\frac{1}{2}\operatorname{Re}\left(p_{+}-p_{-}\right)-|q_{c2}|i$ | $-p_{\text{cut},+}(s,q_{b},-x_{0})$
$q_{3}$ | $\frac{2}{3}(1-2i)|q_{c2}|$ | $p_{\text{cut},+}(s,q_{b},x_{0})$
$q_{4}$ | $\frac{3}{2}|q_{c2}|$ | $\frac{2}{3}(1-2i)|q_{c2}|$
$q_{5}$ | $q_{\text{max}}$ | $\frac{3}{2}|q_{c2}|$
$q_{6}$ | $-$ | $q_{\text{max}}$
Table 4: Example nodes defining piece-wise linear contour $\mathcal{C}$ for
positive and negative values of Im$(s)$. Momentum $q_{c2}$ is defined in Eq.
(40) and momentum $p_{\pm}$ in Eq. (37). Both are evaluated at
$\sigma_{k}=\sigma_{b}$. Moreover,
$x_{0}=|\operatorname{Re}\left(z(q_{c2},q_{b})/2q_{c2}q_{b}\right)|$.

We note that every self-consistent contour that avoids singularities of the
OPE and the integration kernel is a legitimate choice. In practice, the
contour used in the solution routine must evolve with values of $s$, $k$, and
the scattering length $a$, since the position of the OPE cuts depends on these
parameters. Due to the complicated shapes of the cuts we employ contours
defined in a piece-wise linear manner, which allows for more control than
explicitly given, fixed functions. A contour is defined by a set of $n+1$
nodes $\\{q_{i}\\}_{i\in[0,n]}$, which connect lines constituting the
integration path. The $i$-th line is defined as,

$\displaystyle\gamma_{i}(t)=\left(\frac{q_{i+1}-q_{i}}{t_{i+1}-t_{i}}\right)\,t+\frac{t_{i+1}q_{i}-t_{i}q_{i+1}}{t_{i+1}-t_{i}}\,,\leavevmode\nobreak\
\leavevmode\nobreak\ t\in[t_{i},t_{i+1}]\,,\leavevmode\nobreak\
\leavevmode\nobreak\ i=0,1,\dots,n-1\,,$ (86)

where $t_{0}=0<t_{1}<\dots<t_{n-1}<t_{n}=1$. Two example sets of nodes for two
different cases of Im$(s)$ are given in Tab. 4. They are suitable for
$k=q_{b}$, and a relatively large range of complex $s$ and positive $a$.
Example contours created using these nodes are shown in Fig. 9. Note that for
$\operatorname{Im}s>0$ the $G_{S}$ amplitude is evaluated on the second sheet
between points $q_{2}$ and $q_{3}$. One can use contours that have a different
number of nodes depending on the shape of the cut and other practical
considerations.

We note that points $q_{c1}$ and $q_{c2}$, derived in Eqs. (39), (40) are used
in the definition of both contours. They roughly describe the size of the
“circle” and thus are useful in devising an integration path that avoids the
OPE amplitude cuts. Although we are satisfied with this prescription, one can
also generalize those points to a case when $\sigma_{k}$ and $s$ are complex.
This describes the “open” circle scenario. The generalized points are called
$q^{\prime}_{c1}$ and $q^{\prime}_{c2}$. They are derived from the condition
$\operatorname{Im}[z(p,k)/2pk]=0$. Below, we show an example derivation of
$q^{\prime}_{c1}$; the other point is obtained analogously. First, we observe
that the above condition implies,

$\displaystyle\operatorname{Im}[z(p,k)p^{*}k^{*}]=0\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \Rightarrow\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname{Im}z(p,k)\operatorname{Re}k^{*}+\operatorname{Re}z(p,k)\operatorname{Im}[k^{*}]=0\,.$
(87)

since $p=q^{\prime}_{c1}$ is real. (For $q_{c2}$ we assume purely imaginary
$p=q^{\prime}_{c2}$.) We observe that,

$\displaystyle z(p,k)=\sigma_{k}-2(\sqrt{s}-\omega_{k})\omega_{p}\,.$ (88)

Thus,

$\displaystyle\operatorname{Im}z(p,k)=\operatorname{Im}\sigma_{k}-2\omega_{p}\operatorname{Im}[\sqrt{s}-\omega_{k}]\,,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\
\operatorname{Re}z(p,k)=\operatorname{Re}\sigma_{k}-2\omega_{p}\operatorname{Re}[\sqrt{s}-\omega_{k}]\,.$
(89)

This leads to a linear equation for $\omega_{p}$, which can be solved,

$\displaystyle\omega_{p}=\frac{1}{2}\frac{\operatorname{Im}[\sigma_{k}k^{*}]}{\operatorname{Im}[(\sqrt{s}-\omega_{k})k^{*}]}\,.$
(90)

Thus position where the OPE cut crosses the real axis is,

$\displaystyle
q^{\prime}_{c1}=\sqrt{\frac{1}{4}\left(\frac{\operatorname{Im}[\sigma_{k}\,k^{*}]}{\operatorname{Im}[(\sqrt{s}-\omega_{k})\,k^{*}]}\right)^{2}-m^{2}}\,.$
(91)

Similarly, we can obtain a point where it crosses the imaginary axis,

$\displaystyle
q^{\prime}_{c2}=\sqrt{\frac{1}{4}\left(\frac{\operatorname{Re}[\sigma_{k}\,k^{*}]}{\operatorname{Re}[(\sqrt{s}-\omega_{k})\,k^{*}]}\right)^{2}-m^{2}}\,.$
(92)

Reflection of these points with respect to the origin of the complex momentum
plane gives the remaining crossover points of the OPE. We note that for real
$s$, $q_{c1}=q^{\prime}_{c1}$, but $q_{c2}\neq q^{\prime}_{c2}$, since
$q_{c2}$ is not a point of crossover.

Figure 21:  Example interpolation of the piece-wise linear contour with a
smooth one for various values of parameter $\alpha$. Black lines represent
cuts of $G_{S}(q_{b},q)$ for $ma=16$. Left panel: $s/m^{2}=8.6+0.05i$. Right
panel: $s/m^{2}=8.6-0.05i$.

In certain cases, we find that “smoothing” the integration contour leads to a
better numerical convergence of the amplitudes. Derivative
$\gamma^{\prime}(t)$ in Eq. (83) is discontinuous for the piece-wise linear
path, which might prevent one from using certain types of quadratures when
discretizing the integral equation. To smooth out the $\gamma(t)$ function
around points $q_{i}=\gamma(t_{i})$ one may, for example, perform an
interpolation of the contour using cardinal Hermite splines (Schoenberg,
1973). To achieve continuity of $\gamma^{\prime\prime}(t)$ we use the 5th
order polynomials, defined as,

$\displaystyle\left[\begin{array}[]{c}p_{1}(t)\\\ p_{2}(t)\\\ p_{3}(t)\\\
p_{4}(t)\\\ p_{5}(t)\\\
p_{6}(t)\end{array}\right]=\left[\begin{array}[]{cccccc}1&0&0&-10&15&-6\\\
0&0&0&10&-15&6\\\ 0&1&0&-6&8&-3\\\ 0&0&0&-4&7&-3\\\ 0&0&1/2&-3/2&+3/2&-1\\\
0&0&0&1/2&-1&1/2\end{array}\right]\left[\begin{array}[]{c}1\\\ t\\\ t^{2}\\\
t^{3}\\\ t^{4}\\\ t^{5}\end{array}\right]\,.$ (111)

The integration contour is given by,

$\displaystyle\gamma_{i}(t)=p_{1}(y_{i})\,q_{i}+p_{2}(y_{i})\,q_{i+1}+p_{3}(y_{i})\,q^{\prime}_{i}+p_{4}(y_{i})\,q^{\prime}_{i+1}+p_{5}(y_{i})\,q^{\prime\prime}_{i}+p_{6}(y_{i})\,q^{\prime\prime}_{i+1}\,,$
(112)

for $t\in[t_{i},t_{i+1}]$ and $i=0,1,\dots,n-1$. Here
$y_{i}=(t-t_{i})/(t_{i+1}-t_{i})$. The first and second tangents are defined
as,

$\displaystyle q^{\prime}_{i}$ $\displaystyle=$
$\displaystyle\alpha\,\left(\frac{q_{i+1}-q_{i}}{t_{i+1}-t_{i}}+\frac{q_{i}-q_{i-1}}{t_{i}-t_{i-1}}\right)\,,$
(113) $\displaystyle q^{\prime\prime}_{i}$ $\displaystyle=$
$\displaystyle\alpha\,\left(\frac{q^{\prime}_{i+1}-q^{\prime}_{i}}{t_{i+1}-t_{i}}+\frac{q^{\prime}_{i}-q^{\prime}_{i-1}}{t_{i}-t_{i-1}}\right)\,\,.$
(114)

For $i=0$ we use only the first and for $i=n$ second term in (113), and define
$q^{\prime\prime}_{0}=q^{\prime\prime}_{n}=0$. Real parameter $\alpha$ defines
the “tension” of the interpolating curve. For $\alpha=0$ one obtains the
piece-wise linear contour. An example of smooth contours is given in Fig. 21.

###  C.2 Discretization of spectator momenta

To apply Nyström method to the Eq. (85), at fixed $s$ and $k$, one has to fix
the integration contour $\mathcal{C}$, and evaluate the $p$ variable on
$\mathcal{C}$, i.e., rewrite $p=\gamma(u)$, $u\in[0,1]$. One then discretizes
both real variables, $t$ and $u$, to rewrite the integral equation as a matrix
one. In the simplest numerical approach, we use a uniform mesh of $N+1$
points: $p_{i}=\gamma(u_{i})$ and $q_{i}=\gamma(t_{i})$, where
$i\in\\{0,...,N\\}$ and $t_{i}=u_{i}=i/N$. Thus, each linear path of a contour
contains the number of discrete points proportional to its length. This
represents a simple extension of the “brute force” method from Ref. Jackura
_et al._ (2021). The integral is replaced with a sum,

$\displaystyle\bm{d}_{i}=-\bm{G}_{i}-\sum_{j=0}^{N-1}\bm{K}_{ij}\,\bm{d}_{j}\,,$
(115)

where,

$\displaystyle\bm{d}_{i}$ $\displaystyle=$ $\displaystyle
d_{S}(\gamma(u_{i}),p)\,,$ (116) $\displaystyle\bm{G}_{i}$ $\displaystyle=$
$\displaystyle G_{S}(\gamma(u_{i}),p)\,,$ (117) $\displaystyle\bm{K}_{ij}(s)$
$\displaystyle=$
$\displaystyle\left[\gamma(t_{j+1})-\gamma(t_{j})\right]\,K(\gamma(u_{i}),\gamma(t_{j}))\,.$
(118)

We used bold font to indicate that $d_{S}$, $G_{S}$, and kernel $K$ became
vectors and a matrix in the discrete $(u_{i},t_{j})$ space. In Eq. (118), we
employed the simplest rectangular rule with a forward derivative. One can also
apply other methods (e.g. trapezoidal, Simpson, etc.) and use the exact value
of $\gamma^{\prime}(t)$ at a discrete point $t_{j}$. The solution of the
algebraic equation is,

$\displaystyle\bm{d}_{i}^{(\text{sol})}=-\sum_{j=0}^{N-1}\,[\mathbbm{1}+\bm{K}]^{-1}_{ij}\,\bm{G}_{j}\,.$
(119)

Assuming we know $\bm{d}^{(\text{sol})}$, the final amplitude is obtained by
extrapolating the solution to the momentum of interest, e.g., $p=q_{b}$,

$\displaystyle\mathcal{M}_{\varphi
b}(s)=-g^{2}\,G_{S}(q_{b},q_{b})-g^{2}\,\sum_{j=0}^{N-1}\,\left[\gamma(t_{j+1})-\gamma(t_{j})\right]\,K(q_{b},\gamma(t_{j}))\,\bm{d}^{(\text{sol})}_{j}\,.$
(120)

The conceptually simple rectangular rule is an elementary numerical technique
that yields improving results with larger $N$. However, its convergence with
the matrix sizes is relatively slow and can be accelerated with alternative,
more sophisticated discretization techniques. For instance, one can use
Gaussian quadratures Golberg (1979); Atkinson (1997), or spline-based method
Glöckle  _et al._ (1982); Horacek and Malina (1977); Jackura  _et al._ (2021).
We find that Gauus-Chebyshev (GC) and Gauss-Legendre (GL) quadratures can be
easily employed and offer a great improvement in the convergence of the
solutions. Conceptually, implementation of a Gauss quadrature is achieved by
replacing,

$\displaystyle\int_{0}^{1}\,dt\,g(t)=\sum_{n=0}^{N-1}w_{i}\,g(t_{i})\,,$ (121)

for a function $g(t)$. Here, $t_{i}$ and $w_{i}$ are pre-computed mesh points
and corresponding weights, respectively. In practice, this amounts to the
replacement of Eq. (118) with,

$\displaystyle\bm{K}_{ij}(s)$ $\displaystyle=$
$\displaystyle\frac{1}{2}\,w_{j}\,\gamma^{\prime}(t^{\prime}_{j})\,K(\gamma(u^{\prime}_{i}),\gamma(t^{\prime}_{j}))\,.$
(122)

Since the GC and GL quadratures are defined for the integration interval
$t\in[-1,1]$, we map linearly $[-1,1]\to[0,1]$, hence the $1/2$ factor in the
equation above. Primed variables are obtained from the Gauss points $t_{j}$ as
$t^{\prime}_{j}=(1+t_{j})/2$. In our C++ implementation of the ladder
equation, we use available GL quadratures (weights and points) from Ref.
Bogaert (2014). The solution in the GL method is still given by Eq. (119)
while extrapolation $p\to q_{b}$ and $\mathcal{M}_{\varphi b}$ is achieved
through,

$\displaystyle\mathcal{M}_{\varphi
b}(s)=-g^{2}\,G_{S}(q_{b},q_{b})-\frac{1}{2}\,g^{2}\sum_{j=0}^{N-1}\,w_{i}\,\gamma^{\prime}(t^{\prime}_{j})\,K(q_{b},\gamma(t^{\prime}_{j}))\,\bm{d}^{(\text{sol})}_{j}\,.$
(123)

Discretization procedures described in this subsection are also applied to the
homogeneous version of the ladder equation, Eq. (25). Position of the three-
body bound state pole in $s$ is obtained either from the determinant
condition, Eq. (26) or identification of zeros of
$1/\operatorname{Re}(d_{S}(p,p))$ for some choice of external momenta. In both
cases, we accomplish it numerically by using the secant method with precision
$\Delta s=10^{-12}$.

###  C.3 Analysis of the systematic effects

In Ref. Jackura  _et al._ (2021), the authors studied systematic effects of
the numerical approaches by considering two limits: matrix size $N$ going to
infinity, and the two-body pole position shift, $i\epsilon$, going to zero.
Here, we do not deal with poles coinciding with the integration contour, which
usually cause numerical instabilities. Thus, the analysis of systematic
effects is greatly simplified and the precision of the solutions is improved
compared to the previously studied, more demanding case.

Figure 22: Mesh size dependence of the binding energy of the deepest three-
body bound state, $E_{1}$, obtained from the BF and GL methods. On the top-
left panel, we present convergence for two cases, $ma=16$, and $ma=10^{4}$. On
the bottom-left panel, we show the corresponding value of the relative error,
$\Delta E_{1}$ as a function of the $N_{\text{min}}$, as described in the
text. On the right panel, we present a contour used to obtain these results,
with an example of Gaussian nodes for $N=15$. The depicted function is
$\operatorname{Im}G(p,q)$ for $\sigma_{p}=2m^{2}$, with cuts represented by
black lines. It is given in units of $1/m^{2}$. The smooth cut-off was used.
Figure 23: Mesh size dependence of the $\mathcal{M}_{\varphi b}(s)$ for
$s/m^{2}=8.7-0.1i$, obtained from the BF and GL methods. On the top-left
panel, we present the convergence of the real and imaginary parts of the
amplitude. On the bottom-left panel, we show the corresponding value of
$\Delta\mathcal{M}_{\varphi b}$ as a function of the $N_{\text{min}}$. On the
right panel, we present a smooth contour ($\alpha=0.45$) used to obtain these
results, highlighting example Gaussian nodes for $N=30$. The depicted function
is $\operatorname{Im}G(q_{b},q)$, with cuts represented by black lines. It is
given in units of $1/m^{2}$. The smooth cut-off was used.

We find that the GL method leads to a fast convergence of the results when
smooth integration contours are used. Typically, a mesh of $N\approx 100$
points is sufficient to obtain results that cease to depend on the matrix size
within desired precision. Piece-wise linear contours may cause unwanted
oscillations of $\bm{d}^{(\text{sol})}$ considered as a function of $N$. In
this case, the GL method amplitudes still converge faster than the one
obtained from the BF results, but it is harder to analyze them systematically.
We find, the BF method always leads to a smooth controllable dependence of the
$N$-dependent amplitudes, regardless of the type of contour. However, it
requires the implementation of extrapolation to continuum, $N\to\infty$.
Practically, this means one has to calculate the amplitude
$\bm{d}^{(\text{sol})}$ using a set of few matrix sizes, usually of the order
$N\approx 10^{3}$, and then fit the result with the polynomial formula,

$\displaystyle\bm{d}^{(\text{sol})}(N)=A+\frac{B}{N}\,,$ (124)

where the asymptote $A=\bm{d}^{(\text{sol})}(\infty)$ is taken as the
continuum result Jackura  _et al._ (2021). Higher orders of $1/N$ can be
included to improve convergence. From this perspective, the GL method is much
more effective, since it allows one to use a single, relatively small value of
$N$ to obtain the desired outcome with high confidence. We note, that
convergence of the BF method can be improved via different means, e.g.,
Richardson extrapolation Delves and Mohamed (1988); however, we do not
implement any acceleration techniques in this work.

The analysis of the numerical uncertainty of $\bm{d}^{(\text{sol})}$ and its
extrapolations can be performed as described in Chap. 4 of Ref. Delves and
Mohamed (1988). Since the estimated error of our results is satisfactorily
small, we use simpler, rough estimates. We note that the convergence of the
results typically depends on the distance between the singularities of the
kernel and the integration path. For instance, three-body pole positions are
obtained from the ladder equation at external momenta $p^{\prime},p$ for which
the OPE cuts are far from the integration range. In the GL method, this leads
to a relative difference of the order $10^{-4}$% between the $N=15$ and
$N=1000$ results, and virtually no difference between $N=50$ and $N=1000$
values. Thus, for the bound-state pole positions, we take the finite-$N$ GL
result with the error given by the precision of the root-finding algorithm,
which we set to $\Delta s=10^{-12}$. We find that the extrapolated BF result
converges to the GL one when large matrices are used for the fit.

In Fig. 22, we show convergence of the binding energy of the ground-state
timer, $E_{1}=\sqrt{s_{\varphi b}}-\sqrt{s_{b}}$, with matrix size $N$. The
bottom left panel shows the relative difference $\Delta
E_{1}=100\times\left|E_{1,\text{BF}}(N_{\text{min}})-E_{1,\text{GL}}\right|/E_{1,\text{GL}}$,
where $E_{1,\text{BF}}(N_{\text{min}})$ is the extrapolated BF result obtained
from fitting the Eq. (124) in the interval $[N_{\text{min}},1000]$. We see
that, as the larger matrices are used in the fit, the extrapolated BF result
converges to the GL one, reaching an acceptable relative difference of
$10^{-2}$% at $N_{\text{min}}\approx 500$. Since the BF method requires
computation at several values of $N$ to achieve this level of agreement, we
point to a significant advantage of the GL over the BF method.

When the OPE cuts approach the origin of the complex $q$ plane and enclose the
lower limit of the integration, the convergence of the results becomes slower.
The GL method amplitudes exhibit oscillatory behavior with $N$ and do not
stabilize entirely at any finite value of the matrix size. However, despite
this behavior, they still converge very quickly with oscillations damped by
orders of magnitude within a relatively small range of $N$. Due to the
oscillations, one can not easily extrapolate the GL values, e.g., by using a
version of Eq. (124). Instead, for a given value of $s,k$, one computes
$\bm{d}^{(\text{sol})}$ at a few close values of $N$, and takes their average
as the final result, with the largest difference between the two of used
values as an error estimate. Applying this procedure to different values of
complex $s$ reveals that for sufficiently large $N$, the error estimate is
much smaller than $10^{-2}$%, allowing one to use a finite-$N$ result as a
sufficient approximation of the continuum one.

For illustration, in Fig. 23 we present example results for the amplitude
$\mathcal{M}_{\varphi b}$ computed at $s/m^{2}=8.7-0.1i$ and $ma=16$. In the
top left panel, for $N<100$ we see large oscillations of the GL amplitude,
that are quickly damped and hardly noticeable for larger values of $N$. The
bottom panel shows the “quality measure” of the solution,
$\Delta\mathcal{M}_{\varphi b}=100\times\left|(\mathcal{M}_{\varphi
b}(N_{\text{min}})-\mathcal{M}_{\varphi b})/\mathcal{M}_{\varphi b}\right|$,
for both methods. We assume that the “correct solution”, $\mathcal{M}_{\varphi
b}$, is well approximated by an average of GL results obtained for
$N=950,955,...,1000$. For the BF method, the $\mathcal{M}_{\varphi
b}(N_{\text{min}})$ is the extrapolated result obtained from fitting Eq. (124)
in the interval $[N_{\text{min}},1000]$. For the GL method
$\mathcal{M}_{\varphi b}(N_{\text{min}})$ is an average of three values of the
amplitude computed at matrix sizes $N=N_{\text{min}},N_{\text{min}}+5,$ and
$N_{\text{min}}+10$. We see that the GL method offers a reduction of such
defined error by several orders of magnitude compared to the BF approach at a
given matrix size. The actual improvement depends on the value of $s$ and the
contour smoothness parameter $\alpha$.

In this study, we consider a sub-percent precision of our results as entirely
satisfactory. Such uncertainty is much smaller than anticipated errors from
the lattice data that would enter our integral equations through the inclusion
of non-zero $\mathcal{K}_{\text{df},3}$. In most applications, we choose to
use the GL method with matrix size $N=500$, which should result in a relative
error of at most $10^{-2}$%. As discussed above, when the OPE cuts are far
from the integration contour, like in the case of
$\sigma_{p}=\sigma_{k}=2m^{2}$ which we used to extract bound-state pole
positions, the error is expected to be many orders of magnitude smaller.

## References

  * Ketzer  _et al._ (2020) B. Ketzer, B. Grube, and D. Ryabchikov, Prog. Part. Nucl. Phys. 113, 103755 (2020), arXiv:1909.06366 [hep-ex] .
  * Aaij  _et al._ (2022a) R. Aaij _et al._ (LHCb), JHEP 01, 065 (2022a), arXiv:2111.03001 [hep-ex] .
  * Davier  _et al._ (2014) M. Davier, A. Höcker, B. Malaescu, C.-Z. Yuan, and Z. Zhang, Eur. Phys. J. C 74, 2803 (2014), arXiv:1312.1501 [hep-ex] .
  * Garzia (2018) I. Garzia (BESIII), Nuovo Cim. C 41, 102 (2018).
  * Roper (1964) L. D. Roper, Phys. Rev. Lett. 12, 340 (1964).
  * Arndt  _et al._ (2006) R. A. Arndt, W. J. Briscoe, I. I. Strakovsky, and R. L. Workman, Phys. Rev. C 74, 045205 (2006), arXiv:nucl-th/0605082 .
  * Adams  _et al._ (1998) G. S. Adams _et al._ (E852), Phys. Rev. Lett. 81, 5760 (1998).
  * Aghasyan  _et al._ (2018) M. Aghasyan _et al._ (COMPASS), Phys. Rev. D 98, 092003 (2018), arXiv:1802.05913 [hep-ex] .
  * Antoniazzi  _et al._ (1994) L. Antoniazzi _et al._ (E705), Phys. Rev. D 50, 4258 (1994).
  * Choi  _et al._ (2003) S. K. Choi _et al._ (Belle), Phys. Rev. Lett. 91, 262001 (2003), arXiv:hep-ex/0309032 .
  * Aaij  _et al._ (2020) R. Aaij _et al._ (LHCb), JHEP 08, 123 (2020), arXiv:2005.13422 [hep-ex] .
  * Aaij  _et al._ (2022b) R. Aaij _et al._ (LHCb), JHEP 01, 131 (2022b), arXiv:2109.07360 [hep-ex] .
  * Aaij  _et al._ (2022c) R. Aaij _et al._ (LHCb), Nature Phys. 18, 751 (2022c), arXiv:2109.01038 [hep-ex] .
  * Aaij  _et al._ (2022d) R. Aaij _et al._ (LHCb), Nature Commun. 13, 3351 (2022d), arXiv:2109.01056 [hep-ex] .
  * Briceño  _et al._ (2021) R. A. Briceño, J. V. Guerrero, M. T. Hansen, and A. M. Sturzu, Phys. Rev. D 103, 014506 (2021), arXiv:2007.01155 [hep-lat] .
  * Lüscher (1991) M. Lüscher, Nuclear Physics B 354, 531 (1991).
  * Luscher (1986a) M. Luscher, Commun. Math. Phys. 104, 177 (1986a).
  * Luscher (1986b) M. Luscher, Commun. Math. Phys. 105, 153 (1986b).
  * Rummukainen and Gottlieb (1995) K. Rummukainen and S. A. Gottlieb, Nucl. Phys. B 450, 397 (1995), arXiv:hep-lat/9503028 .
  * Kim  _et al._ (2005) C. h. Kim, C. T. Sachrajda, and S. R. Sharpe, Nucl. Phys. B 727, 218 (2005), arXiv:hep-lat/0507006 .
  * Briceno and Davoudi (2013) R. A. Briceno and Z. Davoudi, Phys. Rev. D 88, 094507 (2013), arXiv:1204.1110 [hep-lat] .
  * Hansen and Sharpe (2012) M. T. Hansen and S. R. Sharpe, Phys. Rev. D 86, 016007 (2012), arXiv:1204.0826 [hep-lat] .
  * Briceno (2014) R. A. Briceno, Phys. Rev. D 89, 074507 (2014), arXiv:1401.3312 [hep-lat] .
  * Dudek  _et al._ (2014) J. J. Dudek, R. G. Edwards, C. E. Thomas, and D. J. Wilson (Hadron Spectrum), Phys. Rev. Lett. 113, 182001 (2014), arXiv:1406.4158 [hep-ph] .
  * Alexandrou  _et al._ (2017) C. Alexandrou, L. Leskovec, S. Meinel, J. Negele, S. Paul, M. Petschlies, A. Pochinsky, G. Rendon, and S. Syritsyn, Phys. Rev. D 96, 034525 (2017), arXiv:1704.05439 [hep-lat] .
  * Prelovsek  _et al._ (2021) S. Prelovsek, S. Collins, D. Mohler, M. Padmanath, and S. Piemonte, JHEP 06, 035 (2021), arXiv:2011.02542 [hep-lat] .
  * Brett  _et al._ (2018) R. Brett, J. Bulava, J. Fallica, A. Hanlon, B. Hörz, and C. Morningstar, Nucl. Phys. B 932, 29 (2018), arXiv:1802.03100 [hep-lat] .
  * Woss  _et al._ (2019) A. J. Woss, C. E. Thomas, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. D 100, 054506 (2019), arXiv:1904.04136 [hep-lat] .
  * Woss  _et al._ (2021) A. J. Woss, J. J. Dudek, R. G. Edwards, C. E. Thomas, and D. J. Wilson (Hadron Spectrum), Phys. Rev. D 103, 054502 (2021), arXiv:2009.10034 [hep-lat] .
  * Wilson  _et al._ (2019) D. J. Wilson, R. A. Briceno, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. Lett. 123, 042002 (2019), arXiv:1904.03188 [hep-lat] .
  * Wilson  _et al._ (2015a) D. J. Wilson, R. A. Briceno, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. D 92, 094502 (2015a), arXiv:1507.02599 [hep-ph] .
  * Briceno  _et al._ (2018) R. A. Briceno, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. D 97, 054513 (2018), arXiv:1708.06667 [hep-lat] .
  * Andersen  _et al._ (2018) C. W. Andersen, J. Bulava, B. Hörz, and C. Morningstar, Phys. Rev. D 97, 014506 (2018), arXiv:1710.01557 [hep-lat] .
  * Wilson  _et al._ (2015b) D. J. Wilson, J. J. Dudek, R. G. Edwards, and C. E. Thomas, Phys. Rev. D 91, 054008 (2015b), arXiv:1411.2004 [hep-ph] .
  * Briceno  _et al._ (2017) R. A. Briceno, J. J. Dudek, R. G. Edwards, and D. J. Wilson, Phys. Rev. Lett. 118, 022002 (2017), arXiv:1607.05900 [hep-ph] .
  * Gayer  _et al._ (2021) L. Gayer, N. Lang, S. M. Ryan, D. Tims, C. E. Thomas, and D. J. Wilson (Hadron Spectrum), JHEP 07, 123 (2021), arXiv:2102.04973 [hep-lat] .
  * Dudek  _et al._ (2016) J. J. Dudek, R. G. Edwards, and D. J. Wilson (Hadron Spectrum), Phys. Rev. D 93, 094506 (2016), arXiv:1602.05122 [hep-ph] .
  * Moir  _et al._ (2016) G. Moir, M. Peardon, S. M. Ryan, C. E. Thomas, and D. J. Wilson, JHEP 10, 011 (2016), arXiv:1607.07093 [hep-lat] .
  * Rendon  _et al._ (2020) G. Rendon, L. Leskovec, S. Meinel, J. Negele, S. Paul, M. Petschlies, A. Pochinsky, G. Silvi, and S. Syritsyn, Phys. Rev. D 102, 114520 (2020), arXiv:2006.14035 [hep-lat] .
  * Silvi  _et al._ (2021) G. Silvi _et al._ , Phys. Rev. D 103, 094508 (2021), arXiv:2101.00689 [hep-lat] .
  * Briceño  _et al._ (2018a) R. A. Briceño, J. J. Dudek, and R. D. Young, Rev. Mod. Phys. 90, 025001 (2018a).
  * Hansen and Sharpe (2014) M. T. Hansen and S. R. Sharpe, Phys. Rev. D 90, 116003 (2014), arXiv:1408.5933 [hep-lat] .
  * Hansen and Sharpe (2015) M. T. Hansen and S. R. Sharpe, Phys. Rev. D 92, 114509 (2015), arXiv:1504.04248 [hep-lat] .
  * Blanton and Sharpe (2020a) T. D. Blanton and S. R. Sharpe, Phys. Rev. D 102, 054520 (2020a), arXiv:2007.16188 [hep-lat] .
  * Briceño  _et al._ (2017) R. A. Briceño, M. T. Hansen, and S. R. Sharpe, Phys. Rev. D 95, 074510 (2017), arXiv:1701.07465 [hep-lat] .
  * Briceño  _et al._ (2019a) R. A. Briceño, M. T. Hansen, and S. R. Sharpe, Phys. Rev. D 99, 014516 (2019a), arXiv:1810.01429 [hep-lat] .
  * Hansen  _et al._ (2020) M. T. Hansen, F. Romero-López, and S. R. Sharpe, JHEP 07, 047 (2020), [Erratum: JHEP 02, 014 (2021)], arXiv:2003.10974 [hep-lat] .
  * Blanton and Sharpe (2021a) T. D. Blanton and S. R. Sharpe, Phys. Rev. D 104, 034509 (2021a), arXiv:2105.12094 [hep-lat] .
  * Blanton and Sharpe (2020b) T. D. Blanton and S. R. Sharpe, Phys. Rev. D 102, 054515 (2020b), arXiv:2007.16190 [hep-lat] .
  * Mai  _et al._ (2017) M. Mai, B. Hu, M. Döring, A. Pilloni, and A. Szczepaniak, Eur.Phys.J. A53, 177 (2017), arXiv:1706.06118 [nucl-th] .
  * Mai and Döring (2017) M. Mai and M. Döring, Eur. Phys. J. A 53, 240 (2017), arXiv:1709.08222 [hep-lat] .
  * Hammer  _et al._ (2017a) H.-W. Hammer, J.-Y. Pang, and A. Rusetsky, JHEP 09, 109 (2017a), arXiv:1706.07700 [hep-lat] .
  * Hammer  _et al._ (2017b) H. W. Hammer, J. Y. Pang, and A. Rusetsky, JHEP 10, 115 (2017b), arXiv:1707.02176 [hep-lat] .
  * Döring  _et al._ (2018) M. Döring, H. W. Hammer, M. Mai, J. Y. Pang, t. A. Rusetsky, and J. Wu, Phys. Rev. D 97, 114508 (2018), arXiv:1802.03362 [hep-lat] .
  * Jackura  _et al._ (2019a) A. Jackura, C. Fernández-Ramírez, V. Mathieu, M. Mikhasenko, J. Nys, A. Pilloni, K. Saldaña, N. Sherrill, and A. P. Szczepaniak (JPAC), Eur. Phys. J. C 79, 56 (2019a), arXiv:1809.10523 [hep-ph] .
  * Dawid and Szczepaniak (2021) S. M. Dawid and A. P. Szczepaniak, Phys. Rev. D 103, 014009 (2021), arXiv:2010.08084 [nucl-th] .
  * Müller  _et al._ (2022) F. Müller, J.-Y. Pang, A. Rusetsky, and J.-J. Wu, JHEP 02, 158 (2022), arXiv:2110.09351 [hep-lat] .
  * Jackura  _et al._ (2019b) A. W. Jackura, S. M. Dawid, C. Fernández-Ramírez, V. Mathieu, M. Mikhasenko, A. Pilloni, S. R. Sharpe, and A. P. Szczepaniak, Phys. Rev. D 100, 034508 (2019b), arXiv:1905.12007 [hep-ph] .
  * Briceño  _et al._ (2019b) R. A. Briceño, M. T. Hansen, S. R. Sharpe, and A. P. Szczepaniak, Phys. Rev. D 100, 054508 (2019b), arXiv:1905.11188 [hep-lat] .
  * Blanton and Sharpe (2021b) T. D. Blanton and S. R. Sharpe, Phys. Rev. D 103, 054503 (2021b), arXiv:2011.05520 [hep-lat] .
  * Jackura (2022) A. W. Jackura, (2022), arXiv:2208.10587 [hep-lat] .
  * Hörz and Hanlon (2019) B. Hörz and A. Hanlon, Phys. Rev. Lett. 123, 142002 (2019), arXiv:1905.04277 [hep-lat] .
  * Blanton  _et al._ (2020) T. D. Blanton, F. Romero-López, and S. R. Sharpe, Phys. Rev. Lett. 124, 032001 (2020), arXiv:1909.02973 [hep-lat] .
  * Mai  _et al._ (2020) M. Mai, M. Döring, C. Culver, and A. Alexandru, Phys. Rev. D 101, 054510 (2020), arXiv:1909.05749 [hep-lat] .
  * Culver  _et al._ (2020) C. Culver, M. Mai, R. Brett, A. Alexandru, and M. Döring, Phys. Rev. D 101, 114507 (2020), arXiv:1911.09047 [hep-lat] .
  * Fischer  _et al._ (2021) M. Fischer, B. Kostrzewa, L. Liu, F. Romero-López, M. Ueding, and C. Urbach, Eur. Phys. J. C 81, 436 (2021), arXiv:2008.03035 [hep-lat] .
  * Brett  _et al._ (2021) R. Brett, C. Culver, M. Mai, A. Alexandru, M. Döring, and F. X. Lee, Phys. Rev. D 104, 014501 (2021), arXiv:2101.06144 [hep-lat] .
  * Blanton  _et al._ (2021) T. D. Blanton, A. D. Hanlon, B. Hörz, C. Morningstar, F. Romero-López, and S. R. Sharpe, JHEP 10, 023 (2021), arXiv:2106.05590 [hep-lat] .
  * Alexandru  _et al._ (2020) A. Alexandru, R. Brett, C. Culver, M. Döring, D. Guo, F. X. Lee, and M. Mai, Phys. Rev. D 102, 114523 (2020), arXiv:2009.12358 [hep-lat] .
  * Draper  _et al._ (2023) Z. T. Draper, A. D. Hanlon, B. Hörz, C. Morningstar, F. Romero-López, and S. R. Sharpe, (2023), arXiv:2302.13587 [hep-lat] .
  * Hansen  _et al._ (2021) M. T. Hansen, R. A. Briceño, R. G. Edwards, C. E. Thomas, and D. J. Wilson (Hadron Spectrum), Phys. Rev. Lett. 126, 012001 (2021), arXiv:2009.04931 [hep-lat] .
  * Garofalo  _et al._ (2023) M. Garofalo, M. Mai, F. Romero-López, A. Rusetsky, and C. Urbach, JHEP 02, 252 (2023), arXiv:2211.05605 [hep-lat] .
  * Jackura  _et al._ (2021) A. W. Jackura, R. A. Briceño, S. M. Dawid, M. H. E. Islam, and C. McCarty, Phys. Rev. D 104, 014507 (2021), arXiv:2010.09820 [hep-lat] .
  * Romero-López  _et al._ (2019) F. Romero-López, S. R. Sharpe, T. D. Blanton, R. A. Briceño, and M. T. Hansen, JHEP 10, 007 (2019), arXiv:1908.02411 [hep-lat] .
  * Raposo and Hansen (2023) A. B. a. Raposo and M. T. Hansen, in  _39th International Symposium on Lattice Field Theory_ (2023) arXiv:2301.03981 [hep-lat] .
  * Nyström (1930) E. J. Nyström, Acta Mathematica 54, 185 (1930).
  * Delves and Mohamed (1988) L. Delves and J. Mohamed,  _Computational Methods for Integral Equations_ (Cambridge University Press, 1988).
  * Hwa (1964) R. C. Hwa, Phys.Rev. 134, B1086 (1964).
  * Holman (1965) W. J. Holman, Phys.Rev. 138, B1286 (1965).
  * Grisaru (1966) M. T. Grisaru, Phys.Rev. 146, 1098 (1966).
  * Rubin  _et al._ (1966) M. Rubin, R. Sugar, and G. Tiktopoulos, Phys. Rev. 146, 1130 (1966).
  * Rubin  _et al._ (1967) M. Rubin, R. Sugar, and G. Tiktopoulos, Phys. Rev. 159, 1348 (1967).
  * Brayshaw (1968a) D. D. Brayshaw, Phys. Rev. 176, 1855 (1968a).
  * Glockle (1978) W. Glockle, Phys. Rev. C 18, 564 (1978).
  * Orlov and Turovtsev (1984) Y. Orlov and V. Turovtsev, Zh. Eksp. Teor. Fiz. 86, 1600 (1984).
  * Eichmann  _et al._ (2019) G. Eichmann, P. Duarte, M. Peña, and A. Stadler, Phys. Rev. D 100, 094001 (2019), arXiv:1907.05402 [hep-ph] .
  * Sadasivan  _et al._ (2020) D. Sadasivan, M. Mai, H. Akdag, and M. Döring, Phys.Rev. D101, 094018 (2020), [Erratum: Phys.Rev.D 103, 019901 (2021)], arXiv:2002.12431 [nucl-th] .
  * Brayshaw (1968b) D. D. Brayshaw, Phys. Rev. 167, 1505 (1968b).
  * Huber  _et al._ (2022) M. Q. Huber, W. Kern, and R. Alkofer, (2022), arXiv:2212.02515 [hep-ph] .
  * Ydrefors  _et al._ (2020) E. Ydrefors, J. H. Alvarenga Nogueira, V. A. Karmanov, and T. Frederico, Phys. Rev. D 101, 096018 (2020), arXiv:2005.07943 [hep-ph] .
  * Sadasivan  _et al._ (2022) D. Sadasivan, A. Alexandru, H. Akdag, F. Amorim, R. Brett, C. Culver, M. Döring, F. X. Lee, and M. Mai, Phys. Rev. D 105, 054020 (2022), arXiv:2112.03355 [hep-ph] .
  * Efimov (1970) V. Efimov, Phys. Lett. B 33, 563 (1970).
  * Naidon and Endo (2017) P. Naidon and S. Endo, Rept. Prog. Phys. 80, 056001 (2017), arXiv:1610.09805 [quant-ph] .
  * Green  _et al._ (2021) J. R. Green, A. D. Hanlon, P. M. Junnarkar, and H. Wittig, Phys. Rev. Lett. 127, 242003 (2021), arXiv:2103.01054 [hep-lat] .
  * Hansen and Sharpe (2017) M. T. Hansen and S. R. Sharpe, Phys. Rev. D 95, 034501 (2017), arXiv:1609.04317 [hep-lat] .
  * Briceño  _et al._ (2018b) R. A. Briceño, M. T. Hansen, and S. R. Sharpe, Phys. Rev. D 98, 014506 (2018b), arXiv:1803.04169 [hep-lat] .
  * Cohen (2007) H. Cohen,  _Complex Analysis with Applications in Science and Engineering_ (Springer US, 2007).
  * Eden  _et al._ (1966) R. J. Eden, P. V. Landshoff, D. I. Olive, and J. C. Polkinghorne,  _The analytic S-matrix_ (Cambridge Univ. Press, Cambridge, 1966).
  * Burkhardt (1969) H. Burkhardt,  _Dispersion Relation Dynamics: A Phenomenological Introduction to S-matrix Theory_ (North-Holland Publishing Company, 1969).
  * Lang (1985) S. Lang,  _Complex Analysis_, Graduate Texts in Mathematics (Springer New York, 1985).
  * Mai (2021) M. Mai, Eur. Phys. J. ST 230, 1593 (2021), arXiv:2010.00056 [nucl-th] .
  * Albaladejo  _et al._ (2022) M. Albaladejo, Łukasz Bibrzycki, S. M. Dawid, C. Fernández-Ramírez, S. Gonzàlez-Solís, A. N. H. Blin, A. W. Jackura, V. Mathieu, M. Mikhasenko, V. I. Mokeev, E. Passemar, A. Pilloni, A. Rodas, J. A. Silva-Castro, W. A. Smith, A. P. Szczepaniak, and D. Winney, Progress in Particle and Nuclear Physics 127, 103981 (2022).
  * Doring  _et al._ (2009) M. Doring, C. Hanhart, F. Huang, S. Krewald, and U. G. Meissner, Nucl. Phys. A 829, 170 (2009), arXiv:0903.4337 [nucl-th] .
  * Schoenberg (1973) I. J. Schoenberg,  _Cardinal Spline Interpolation_ (Society for Industrial and Applied Mathematics, 1973) https://epubs.siam.org/doi/pdf/10.1137/1.9781611970555 .
  * Golberg (1979) M. A. Golberg,  _Solution Methods for Integral Equations_, Mathematical Concepts and Methods in Science and Engineering (Springer New York, 1979).
  * Atkinson (1997) K. E. Atkinson,  _The Numerical Solution of Integral Equations of the Second Kind_, Cambridge Monographs on Applied and Computational Mathematics (Cambridge University Press, 1997).
  * Glöckle  _et al._ (1982) W. Glöckle, G. Hasberg, and A. R. Neghabian, Z. Phys. A 305, 217 (1982).
  * Horacek and Malina (1977) J. Horacek and L. Malina, Czech. J. Phys. B 27, 1 (1977).
  * Bogaert (2014) I. Bogaert, SIAM Journal on Scientific Computing 36, A1008 (2014), https://doi.org/10.1137/140954969 .