# Evolution of Efimov States

Sebastian M. Dawid<EMAIL_ADDRESS>Physics Department, University of
Washington, Seattle, WA 98195-1560, USA  Md Habib E Islam<EMAIL_ADDRESS>Department of Physics, Old Dominion University, Norfolk, Virginia 23529, USA
Thomas Jefferson National Accelerator Facility, 12000 Jefferson Avenue,
Newport News, Virginia 23606, USA  Raúl A. Briceño<EMAIL_ADDRESS>Department of Physics, University of California, Berkeley, CA 94720, USA
Nuclear Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA
94720, USA  Andrew W. Jackura<EMAIL_ADDRESS>Department of Physics,
University of California, Berkeley, CA 94720, USA  Nuclear Science Division,
Lawrence Berkeley National Laboratory, Berkeley, CA 94720, USA  Department of
Physics, William & Mary, Williamsburg, Virginia 23187, USA

(August 28, 2024)

###### Abstract

The Efimov phenomenon manifests itself as an emergent discrete scaling
symmetry in the quantum three-body problem. In the unitarity limit, it leads
to an infinite tower of three-body bound states with energies forming a
geometric sequence. In this work, we study the evolution of these so-called
Efimov states using relativistic scattering theory. We identify them as poles
of the three-particle $S$ matrix and trace their trajectories in the complex
energy plane as they evolve from virtual states through bound states to
resonances. We dial the scattering parameters toward the unitarity limit and
observe the emergence of the universal scaling of energies and couplings—a
behavior known from the non-relativistic case. Interestingly, we find that
Efimov resonances follow unusual, cyclic trajectories accumulating at the
three-body threshold and then disappear at some values of the two-body
scattering length. We propose a partial resolution to this “missing states”
problem.

_Introduction:_ The discovery of the Efimov effect in 1970 revealed the
formation of an infinite number of bound states, or trimers, in a system of
three non-relativistic bosons (Efimov, 1970, 1973). Assuming they interact via
two-body forces characterized by a large scattering length, the three-body
binding energies form a geometric series with a quotient $\lambda^{2}\approx
515$. The emergence of the phenomenon is closely tied to the scale invariance
of the quantum-mechanical $1/r^{2}$ potential (Case, 1950; de Alfaro  _et al._
, 1976; Niemann and Hammer, 2015; Dawid  _et al._ , 2018), and is the best-
known example of the renormalization group limit cycle (Wilson, 1971; Glazek
and Wilson, 2002; Bedaque  _et al._ , 1999a; Hammer and Platter, 2011).

The sequence of trimers becomes infinite in the so-called unitarity limit,
i.e., when the two-body scattering length, $a$, is made arbitrarily large,
$a\to\infty$. While such behavior has not been observed in nature, several
nuclear Adhikari and Tomio (1982); Bedaque  _et al._ (2000); Garrido  _et al._
(2006); Kievsky  _et al._ (2021) and hadronic systems (Braaten and Kusunoki,
2004; Canham  _et al._ , 2009; Wilbring  _et al._ , 2017; Valderrama, 2019,
2018) may serve as proxies due to their large scattering length. Furthermore,
Efimov physics is realized experimentally using ultracold atoms submerged in a
background magnetic field tuned to introduce a Feschbach resonance and drive
the system to the unitarity limit (Kraemer  _et al._ , 2006; Zaccanti  _et
al._ , 2009; Pollack  _et al._ , 2009; Huckans  _et al._ , 2009; Williams  _et
al._ , 2009; Gross  _et al._ , 2009; Ferlaino  _et al._ , 2009; Lompe  _et
al._ , 2010; Nakajima  _et al._ , 2010; Gross  _et al._ , 2010; Berninger  _et
al._ , 2011; Wild  _et al._ , 2012; Huang  _et al._ , 2014). Given the
generality of the result, this phenomenon has ignited a rich line of research
into universality across different subfields (Fedorov and Jensen, 1993;
Delfino and Frederico, 1996; Hanna and Blume, 2006; Braaten and Hammer, 2006;
Pricoupenko, 2010; von Stecher, 2010; Wang  _et al._ , 2012; Naidon  _et al._
, 2014; Horinouchi and Ueda, 2015; Naidon and Endo, 2017; Greene  _et al._ ,
2017; Deltuva, 2021; Rosa  _et al._ , 2022; Frederico and Gattobigio, 2023;
Pricoupenko, 2023).

Figure 1: Trajectories of the first three trimer poles in the
$(\kappa/m,1/ma)$ plane where $\kappa={\rm sign}(\operatorname{Re}\Delta
E)\sqrt{|m\operatorname{Re}\Delta E|}$. The $3\varphi$ and $\varphi b$
thresholds are shown explicitly as grey and orange lines. Solid red, blue, and
green lines denote physical bound states, while dashed ones denote either
virtual bound states on the unphysical $\varphi b$ sheet or resonances on the
nearest $3\varphi$ sheet. Stars denote the emergence of a virtual state from
the logarithmic cut on the second $\varphi b$ sheet. Circles denote the
evolution of this virtual state onto a real bound state. Squares denote the
further evolution of the state to three-body resonance. Insets show behavior
of trimers near the three-body threshold.

Although the unitarity limit does not seem to exist in nature, we can expose
the universal scaling behavior by exploring the evolution of Efimov states in
vicinity of this limit. We investigate this evolution using relativistic
scattering theory, which has been derived as part of ongoing efforts to
develop a model-independent framework for studying three-body systems Hansen
and Sharpe (2014, 2015); Mai  _et al._ (2017); Blanton and Sharpe (2020a);
Briceño  _et al._ (2017, 2019a); Hansen  _et al._ (2020); Blanton and Sharpe
(2021a, 2020b); Jackura  _et al._ (2019a); Briceño  _et al._ (2019b); Blanton
and Sharpe (2021b); Jackura (2023). Building on previous work Dawid  _et al._
(2023), we identify the trimers as poles of the $S$ matrix in the complex
energy variable and study their behavior for various values of $a$, including
the $a\to\infty$ limit. We provide evidence of the discrete scaling
relationship between the binding energies of the three-body spectrum,

$\displaystyle\Delta E_{n}(a)=Q_{a}^{2}\,\Delta E_{n+1}(Q_{a}a)\,,$ (1)

where $\Delta E_{n}$ is the binding energy of the $n^{\mathrm{th}}$ bound
state, and $Q_{a}$ is a scaling quotient that asymptotes to Efimov’s $\lambda$
in the unitarity limit. This scaling relationship holds as the states evolve
from bound states to unstable resonances, verifying that the relativistic
framework recovers the known non-relativistic results.

Furthermore, by studying the analytic structure of the scattering amplitude,
we find a much richer picture of the trimers’ behavior than previously
identified. We discuss intriguing properties of their evolution across various
unphysical Riemann sheets of the complex energy plane, such as the formation
of cyclic trajectories of the three-body poles and the emergent scaling
property of the associated residues. The behavior of the Efimov resonances is
sufficiently puzzling that it motivates us to conjecture about the structure
of the three-boson amplitudes and to call for further investigation of these
states. Before presenting our findings, we briefly review the framework needed
to obtain them.

_Relativistic scattering theory:_ We consider the scattering of three
identical spinless bosons of mass $m$, which we label as “$\varphi$”, in their
c.m. frame. We fix the total angular momentum of the system to $J=0$, as well
as neglect contributions from the two-particle subsystems of angular momenta
higher than zero. The $3\varphi\to 3\varphi$ scattering amplitude depends on
the total relativistic energy $E$ and two more variables. We describe the
system by splitting the scattering states into a _spectator_ particle and a
_pair_ , formed from the two other bosons, and use initial and final spectator
momenta, $k$ and $p$, as the remaining kinematic parameters. In what follows,
we use a notation with an implicit energy dependence.

Physical states are associated with poles of the amplitude, with a residue
corresponding to the coupling of the state to the open scattering channel.
Lehmann–Symanzik–Zimmerman (Lehmann  _et al._ , 1955; Zimmermann, 1958; Fried,
1960; Duncan, 2012) reduction implies that this identification also holds for
poles off the real energy axis. Causality assures that a complex-valued pole
can not reside on the “physical” energy plane and must instead appear in
unphysical Riemann sheets generated by square-root and logarithmic branch cuts
of the scattering amplitude. Depending on the location in the complex plane
and the sheet, these poles are associated with bound states (real-valued,
physical sheet), virtual states (real-valued, unphysical sheet), or resonances
(complex-valued, unphysical sheet).

The relativistic three-body amplitude, $\mathcal{M}_{3}(p,k)$, exhibits poles
associated with trimers in the $E^{2}$ plane. Near the $n^{\text{th}}$ pole,
it behaves like

$\displaystyle\mathcal{M}_{3}(p,k)$
$\displaystyle=-\frac{\Gamma_{n}(p)\,\Gamma_{n}(k)}{E^{2}-E_{n}^{2}}+\mathcal{O}\\!\left(E^{0}\right)\,,$
(2)

where $E_{n}$ is the trimer energy. Bound or virtual states have
$\operatorname{Im}E_{n}=0$, while resonances $\operatorname{Im}E_{n}\neq 0$.
The residue, i.e., the coupling of the $n^{\text{th}}$ trimer to the
$3\varphi$ state, factorizes into momentum-dependent vertex factors
$\Gamma_{n}(k)$ that are closely related to the Faddeev wave functions in the
non-relativistic limit.

As implied by the unitarity of the $S$ matrix, the amplitude is described by a
set of integral equations Hansen and Sharpe (2015); Mai  _et al._ (2017);
Jackura  _et al._ (2019b); Jackura (2023). They depend on two dynamical
inputs, the $2\varphi\to 2\varphi$ scattering amplitude, $\mathcal{M}_{2}$,
and the three-body $K$ matrix, $\mathcal{K}_{3}$, which describes short-
distance dynamics of three particles. Given these two objects, one can solve
integral equations to obtain the scattering amplitude Hansen and Sharpe
(2015); Jackura  _et al._ (2021); Dawid  _et al._ (2023). As argued in the
supplemental material, the universal scaling behavior is independent of
$\mathcal{K}_{3}$, and we set it to zero in the remainder of this letter. Most
of the techniques we use have been developed in Ref. Dawid  _et al._ (2023)
and references within. We discuss some new details in the supplemental
material.

The three-body scattering amplitude is given by,

$\displaystyle\mathcal{M}_{3}(p,k)$
$\displaystyle=-\mathcal{M}_{2}(p)\,G(p,k)\,\mathcal{M}_{2}(k)$
$\displaystyle\qquad-\mathcal{M}_{2}(p)\int_{k^{\prime}}G(p,k^{\prime})\,\mathcal{M}_{3}(k^{\prime},k)\,,$
(3)

where $G$ is the $S$-wave-projected propagator describing particle exchange
between the pairs. It is a kinematic function with a logarithmic branch cut
(Hansen and Sharpe, 2015; Jackura  _et al._ , 2019b; Dawid  _et al._ , 2023).
Finally, the integral measure is $\int_{k}\equiv\int_{0}^{k_{\rm max}}{\rm
d}k\,k^{2}/(2\pi)^{2}\omega_{k}$, where $k_{\rm max}$ is the maximum allowed
value of the momentum and $\omega_{k}=\sqrt{m^{2}+k^{2}}$ is the spectator
energy. The cutoff momentum is fixed by our choice of $\mathcal{K}_{3}=0$.
Changes in the regularization lead to a different three-body $K$ matrix,
assuring that the resultant amplitude is independent of the cutoff.

For $\mathcal{M}_{2}$, we use the leading-order effective range expansion,

$\displaystyle\mathcal{M}_{2}(k)=\frac{16\pi\varepsilon_{k}}{-1/a-iq_{k}}\,,$
(4)

where $\varepsilon_{k}=\sqrt{(E-\omega_{k})^{2}-k^{2}}$ is the pair’s energy
in its c.m. frame, and $q_{k}=\sqrt{\varepsilon_{k}^{2}/4-m^{2}}$ is the
relative momentum between the particles in the pair. Due to the square root in
the definition of the relative momentum, the amplitude is defined on two
Riemann branches in the complex $\varepsilon_{k}$ variable; the first being
the physical sheet with $\operatorname{Im}q_{k}>0$, while the second
(unphysical) sheet corresponding to $\operatorname{Im}q_{k}<0$.

Regardless of the value of $|ma|\geq 1$, the $\mathcal{M}_{2}$ amplitude has a
pole in the $\varepsilon_{k}$ variable, corresponding to a state with mass
$m_{b}=2\sqrt{m^{2}-1/a^{2}}$. It resides on the real axis below the two-body
threshold, $\varepsilon_{k}<2m$. If $a>0$, the pole is on the first sheet and
is associated with a two-body bound state. Otherwise, it is a virtual state on
the second sheet.

Figure 2: Trajectories of the first three resonances on the nearest unphysical
Riemann sheet of the complex $\Delta E$ plane. Energies of the 2nd and 3rd
trimers are rescaled by $\lambda^{2}$ and $\lambda^{4}$, respectively. At
large $|ma|$ and close to the threshold, all trajectories exhibit discrete
scaling symmetry. As $|ma|$ decreases, the scaling symmetry breaks down,
although, for the second and third resonant states, the discrepancy between
the trajectories remains small.

_Analytic continuation:_ Following our previous work in Refs. Jackura  _et
al._ (2021); Dawid  _et al._ (2023), we numerically solve Eq. (Evolution of
Efimov States) to obtain $\mathcal{M}_{3}$ in the complex $E^{2}$ plane on the
physical and the nearest unphysical sheets. Analytic continuation to the
complex plane depends on the nature of the singularities of the three-body
scattering amplitude encoded in Eq. (Evolution of Efimov States).

In addition to potential poles, the $\mathcal{M}_{3}$ amplitude has a
logarithmic branch cut inherited from the partial-wave projected propagator,
$G$. Furthermore, it has two possible physical thresholds manifesting as
corresponding branch points. These are the square-root bound-state-spectator
threshold at $E_{\mathrm{thr.}}^{(\varphi b)}=m+m_{b}$ and the logarithmic
three-body threshold at $E_{\mathrm{thr.}}^{(3\varphi)}=3m$ (Landau, 1959;
Eden  _et al._ , 1966). Unphysical sheets are associated with these two
singularities.

The emergence of these thresholds in Eq. (Evolution of Efimov States) has a
non-perturbative origin. The three-body amplitude inherits the singularities
of $\mathcal{M}_{2}$ in the external momentum variables, $p$ and $k$. The
threshold branch points of $\mathcal{M}_{3}(p,k)$ in the $E^{2}$ plane emerge
from the second term of the integral equation when these energy-dependent
singularities in the $k^{\prime}$ variable coincide with the origin of the
integration interval, $k^{\prime}=0$ (Dawid  _et al._ , 2023; Eden  _et al._ ,
1966; Burkhardt, 1969). The branch cut at $E=E_{\mathrm{thr.}}^{(\varphi b)}$
arises from the collision of point $k^{\prime}=0$ with the two-body bound-
state pole. The $E_{\mathrm{thr.}}^{(3\varphi)}$ branch point appears from the
collision with the square-root branch point of $\mathcal{M}_{2}$.

This observation suggests a procedure for the analytic continuation of the
amplitude defined by the integral equation. Namely, to extend
$\mathcal{M}_{3}$ to the unphysical Riemann sheets of the $E^{2}$ plane,
either through the $\varphi b$ or the $3\varphi$ cut, one needs to avoid
integrating over the discontinuities associated with the above-mentioned
collisions, i.e., avoid coincidence of the integration interval with the pole
or threshold cut in $k^{\prime}$. We accomplish this by deforming the
integration contour into the complex $k^{\prime}$ momentum plane in Eq.
(Evolution of Efimov States). In doing this, we ensure that the deformed
integration path avoids logarithmic branch points of $G(p,k^{\prime})$ and all
other singularities induced by the non-perturbative nature of the equation
Dawid  _et al._ (2018). We give a more detailed description of this procedure
in the supplemental material.

Figure 3: Riemann surfaces of the three-body amplitude in $\Delta E$. The
trajectories of the trimers are shown, along with three bound states positions
for some $a$. On the nearest unphysical sheet, the 2nd and 3rd trimer exhibit
the cyclic behavior as shown in Fig. 2. Mirror poles are found by continuing
up to sheet $-1$. We postulate that the higher trimers come from the further
unphysical sheets ($\geq+2$), where their cyclic behavior repeats.

_Efimov trimers:_ Close to the unitarity limit, i.e., for $|ma|\gg 1$, we find
that $\mathcal{M}_{3}$ develops multiple bound state poles. We observe that
their binding energies, $\Delta E_{n}=E_{n}-E_{\text{thr.}}^{(3\varphi)}$,
obey the discrete scaling symmetry given in Eq. (1), which is characteristic
of the Efimov phenomenon. The quotient $Q_{a}\to\lambda$ as $ma\to\infty$,
confirming that the relativistic framework recovers the expected Efimov
scaling. In addition, in the supplement, we show that the vertex functions
$\Gamma_{n}(p)$ agree with the prediction of Ref. (Hansen and Sharpe, 2017)
for non-relativistic momenta, $p\ll m$.

By dialing $ma$ to smaller values, we trace the trimers on their trajectories
that span across multiple Riemann sheets associated with the dimer-particle
and three-particle cuts. The trimers evolve from the virtual states (small,
positive $a$) through bound states (large $a$ of both signs) to resonances
(small, negative $a$). In Fig. 1, we present their trajectories on the so-
called Efimov plot, i.e., in the $(\kappa,1/a)$ plane, where $\kappa={\rm
sign}(\operatorname{Re}\Delta E)\sqrt{|m\operatorname{Re}\Delta E|}$.

At small values of $a>0$, the ground state exhibits a noticeable deviation
from the scaling behavior. Nevertheless, its qualitative features (e.g., pole
trajectory, residues) remain analogous to shallow trimers of non-relativistic
binding energies.

All excited states follow similar trajectories. They emerge as virtual states
on the unphysical $\varphi b$ sheet from the logarithmic cut inherited by
$\mathcal{M}_{3}$ from the one-particle exchange amplitude, $G$. They approach
the dimer-particle threshold and move to the first sheet, becoming bound
states. They remain bound states for large negative values of $a$ and evolve
to become resonances on the nearest unphysical sheet associated with the
logarithmic $3\varphi$ threshold cut.

We trace their motion on this sheet, in the complex $\Delta E$ variable, and
present it in Fig. 2. At a given finite value of $ma$, there is only a single
resonance pole in the unphysical sheet. As we decrease $ma$ from zero to
$ma=-8.71$, the “ground state” resonance moves to the three-body threshold on
an arc from complex infinity. It is natural since $ma\to 0$ corresponds to no
dynamics and the removal of all but the free states from the spectrum.

By contrast, the excited three-body resonances follow cyclic trajectories,
which start and end at the three-body threshold and accumulate near this
point. By rescaling them by an appropriate power of $\lambda^{2}$, we observe
that they nearly overlap, providing additional evidence of the discrete scale
invariance in the three-boson system. The loop-like trajectories of the Efimov
resonances were previously noticed in Ref. (Hyodo  _et al._ , 2014), where a
non-relativistic approach was used.

Moreover, we discover an interesting pattern as these excited trimers move
between the physical and unphysical Riemann sheets. Namely, the first excited
resonance of energy $\Delta E_{2}$ emerges from the threshold on the
unphysical sheet at the same value, $ma\approx-8.71$, at which the ground
state resonance reaches this point and becomes a bound state on the physical
energy plane. It leads to a “missing poles” problem: one pole reaches the
threshold, and two emerge. This behavior is repeated for all states, i.e.,
whenever the $n^{\text{th}}$ resonance enters the threshold, the
$(n+1)^{\text{th}}$ resonance appears on the unphysical sheet and the
$n^{\text{th}}$ bound state on the physical one. Moreover, we find that the
residues of all three poles converge to the same value when they approach the
$3\varphi$ branch point.

This puzzling behavior violates our expectation that the number of poles,
equivalent to the number of physical states, must be conserved when one varies
the theory parameters. The only exception is the instance of lifted spectrum
degeneracy, which we verify does not happen in our system by studying the
order of trimer poles.

We propose a possible resolution to this puzzle by noting that the three-body
scattering amplitude has infinitely many unphysical sheets; see Fig. 3. We
label the two nearest sheets $\pm 1$ while denoting the physical sheet by $0$.
States in the $-1$ plane are the complex-conjugate or “ _mirror_ ” poles of
those in the $+1$ sheet. This is unlike the two-body case, where the Schwarz
reflection principle ensures that a resonance pole has its mirror image on the
same sheet. For the three-body amplitude, similarly to a complex logarithm, it
implies a reflection between the $\pm n$ sheets.

We conjecture that the missing poles come from the higher Riemann branches,
one from each. The $n^{\text{th}}$ state approaches the threshold from complex
infinity on the $n^{\text{th}}$ Riemann sheet and moves to the
$(n-1)^{\text{th}}$ one, where it starts evolving on a cyclic trajectory.
Eventually, in the unitarity limit, it travels to the physical energy plane,
contributing to the geometric series of bound states. We depict this idea in
Fig. 3, where the dashed lines represent trajectories of $n\geq 4$ states.
Reference (Hyodo  _et al._ , 2014) did not address the “missing poles” issue.

One could verify this conjecture by analytically continuing the amplitude to
the higher Riemann sheets. Although we are currently unable to extend our
solution to the other sheets, we have performed numerical extrapolations
presented in the supplemental material that further support this conjecture.
Our proposal only partially resolves the puzzle. Whenever a resonance
approaches the threshold, its “mirror” image does the same. Yet, we observe
only one bound state emerging on the physical energy plane. The “mirror” poles
seem to vanish when meeting their complex-conjugate partners at the threshold,
which violates our expectation about the conservation of the number of states.

_Discussion:_ To summarize, we found and presented the emergence of the Efimov
effect from the relativistic three-body scattering equations. In particular,
we discovered evidence of the discrete scaling symmetry in the trajectories of
resonances in the nearby unphysical sheet of complex energy. By studying the
evolution of the spectrum onto unphysical sheets, we make several observations
suggesting that the Efimov phenomenon is closely related to the logarithmic
nature of the three-body unitarity cut, i.e., the presence of infinitely many
branches. At the same time, our conjecture about the behavior of the trimer
trajectories cannot be the end of the story because of the mirror fashion in
which $+$ and $-$ sheets contribute trimer poles to the physical sheet.

The “missing pole” problem is not just a mathematical curiosity but points to
a deficiency in our knowledge about the analytic structure and properties of
three-particle scattering amplitudes. This, in turn, affects our understanding
of the nature of particles that couple strongly to three-particle states
(Adams  _et al._ , 1998; Aghasyan  _et al._ , 2018; Antoniazzi  _et al._ ,
1994; Choi  _et al._ , 2003; Aaij  _et al._ , 2022a, b, 2014; Cheng  _et al._
, 2016; Chang  _et al._ , 2017; LHC, 2023). Having relativistic scattering
amplitudes that satisfy unitarity and whose analytic structure we can fully
control will impact a broad set of experimental, phenomenological, and lattice
QCD studies. As a result, we close by encouraging further investigations along
these lines.

_Acknowledgements:_ The authors thank T. Hyodo, S. Sharpe, M. Baker, and S.
Mizera for valuable discussions. SMD is supported by U.S. Department of Energy
Contract no. DE-SC0011637. RAB and MHI acknowledge the support of the USDOE
Early Career award, contract DE-SC0019229. MHI acknowledges the support from
the Jefferson Science Associates/Jefferson Lab graduate fellowship program.
AWJ acknowledges the support of the USDOE ExoHad Topical Collaboration,
contract DE-SC0023598.

## References

  * Efimov (1970) V. Efimov, Phys. Lett. B 33, 563 (1970).
  * Efimov (1973) V. Efimov, Nucl. Phys. A 210, 157 (1973).
  * Case (1950) K. M. Case, Phys. Rev. 80, 797 (1950).
  * de Alfaro  _et al._ (1976) V. de Alfaro, S. Fubini, and G. Furlan, Nuovo Cim. A 34, 569 (1976).
  * Niemann and Hammer (2015) P. Niemann and H. W. Hammer, Few Body Syst. 56, 869 (2015), arXiv:1504.04511 [nucl-th] .
  * Dawid  _et al._ (2018) S. M. Dawid, R. Gonsior, J. Kwapisz, K. Serafin, M. Tobolski, and S. D. Głazek, Phys. Lett. B 777, 260 (2018), arXiv:1704.08206 [quant-ph] .
  * Wilson (1971) K. G. Wilson, Phys. Rev. D 3, 1818 (1971).
  * Glazek and Wilson (2002) S. D. Glazek and K. G. Wilson, Phys. Rev. Lett. 89, 230401 (2002), [Erratum: Phys.Rev.Lett. 92, 139901 (2004)], arXiv:hep-th/0203088 .
  * Bedaque  _et al._ (1999a) P. F. Bedaque, H. W. Hammer, and U. van Kolck, Phys. Rev. Lett. 82, 463 (1999a), arXiv:nucl-th/9809025 .
  * Hammer and Platter (2011) H.-W. Hammer and L. Platter, Phil. Trans. Roy. Soc. Lond. A 369, 2679 (2011), arXiv:1102.3789 [nucl-th] .
  * Adhikari and Tomio (1982) S. K. Adhikari and L. Tomio, Phys. Rev. C 26, 83 (1982).
  * Bedaque  _et al._ (2000) P. F. Bedaque, H. W. Hammer, and U. van Kolck, Nucl. Phys. A 676, 357 (2000), arXiv:nucl-th/9906032 .
  * Garrido  _et al._ (2006) E. Garrido, D. V. Fedorov, and A. S. Jensen, Phys. Rev. Lett. 96, 112501 (2006).
  * Kievsky  _et al._ (2021) A. Kievsky, L. Girlanda, M. Gattobigio, and M. Viviani, Ann. Rev. Nucl. Part. Sci. 71, 465 (2021), arXiv:2102.13504 [nucl-th] .
  * Braaten and Kusunoki (2004) E. Braaten and M. Kusunoki, Phys. Rev. D 69, 074005 (2004), arXiv:hep-ph/0311147 .
  * Canham  _et al._ (2009) D. L. Canham, H. W. Hammer, and R. P. Springer, Phys. Rev. D 80, 014009 (2009), arXiv:0906.1263 [hep-ph] .
  * Wilbring  _et al._ (2017) E. Wilbring, H. W. Hammer, and U.-G. Meißner, (2017), arXiv:1705.06176 [hep-ph] .
  * Valderrama (2019) M. P. Valderrama, Phys. Rev. D 99, 034010 (2019), arXiv:1811.10173 [hep-ph] .
  * Valderrama (2018) M. P. Valderrama, Phys. Rev. D 98, 034017 (2018), arXiv:1805.10584 [hep-ph] .
  * Kraemer  _et al._ (2006) T. Kraemer, M. Mark, P. Waldburger, J. G. Danzl, C. Chin, B. Engeser, A. D. Lange, K. Pilch, A. Jaakkola, H.-C. Nägerl, and R. Grimm, Nature 440, 315 (2006).
  * Zaccanti  _et al._ (2009) M. Zaccanti, B. Deissler, C. D’Errico, M. Fattori, M. Jona-Lasinio, S. Müller, G. Roati, M. Inguscio, and G. Modugno, Nature Physics 5, 586 (2009).
  * Pollack  _et al._ (2009) S. E. Pollack, D. Dries, and R. G. Hulet, Science 326, 1683 (2009).
  * Huckans  _et al._ (2009) J. H. Huckans, J. R. Williams, E. L. Hazlett, R. W. Stites, and K. M. O’Hara, Phys. Rev. Lett. 102, 165302 (2009).
  * Williams  _et al._ (2009) J. R. Williams, E. L. Hazlett, J. H. Huckans, R. W. Stites, Y. Zhang, and K. M. O’Hara, Phys. Rev. Lett. 103, 130404 (2009).
  * Gross  _et al._ (2009) N. Gross, Z. Shotan, S. Kokkelmans, and L. Khaykovich, Phys. Rev. Lett. 103, 163202 (2009).
  * Ferlaino  _et al._ (2009) F. Ferlaino, S. Knoop, M. Berninger, W. Harm, J. P. D’Incao, H.-C. Nägerl, and R. Grimm, Phys. Rev. Lett. 102, 140401 (2009).
  * Lompe  _et al._ (2010) T. Lompe, T. B. Ottenstein, F. Serwane, K. Viering, A. N. Wenz, G. Zürn, and S. Jochim, Phys. Rev. Lett. 105, 103201 (2010).
  * Nakajima  _et al._ (2010) S. Nakajima, M. Horikoshi, T. Mukaiyama, P. Naidon, and M. Ueda, Phys. Rev. Lett. 105, 023201 (2010).
  * Gross  _et al._ (2010) N. Gross, Z. Shotan, S. Kokkelmans, and L. Khaykovich, Phys. Rev. Lett. 105, 103203 (2010).
  * Berninger  _et al._ (2011) M. Berninger, A. Zenesini, B. Huang, W. Harm, H.-C. Nägerl, F. Ferlaino, R. Grimm, P. S. Julienne, and J. M. Hutson, Phys. Rev. Lett. 107, 120401 (2011).
  * Wild  _et al._ (2012) R. J. Wild, P. Makotyn, J. M. Pino, E. A. Cornell, and D. S. Jin, Phys. Rev. Lett. 108, 145305 (2012).
  * Huang  _et al._ (2014) B. Huang, L. A. Sidorenkov, R. Grimm, and J. M. Hutson, Phys. Rev. Lett. 112, 190401 (2014).
  * Fedorov and Jensen (1993) D. V. Fedorov and A. S. Jensen, Phys. Rev. Lett. 71, 4103 (1993).
  * Delfino and Frederico (1996) A. Delfino and T. Frederico, Phys. Rev. C 53, 62 (1996).
  * Hanna and Blume (2006) G. J. Hanna and D. Blume, Phys. Rev. A 74, 063604 (2006).
  * Braaten and Hammer (2006) E. Braaten and H. W. Hammer, Phys. Rept. 428, 259 (2006), arXiv:cond-mat/0410417 .
  * Pricoupenko (2010) L. Pricoupenko, Phys. Rev. A 82, 043633 (2010).
  * von Stecher (2010) J. von Stecher, Journal of Physics B: Atomic, Molecular and Optical Physics 43, 101002 (2010).
  * Wang  _et al._ (2012) J. Wang, J. P. D’Incao, B. D. Esry, and C. H. Greene, Phys. Rev. Lett. 108, 263001 (2012).
  * Naidon  _et al._ (2014) P. Naidon, S. Endo, and M. Ueda, Phys. Rev. Lett. 112, 105301 (2014).
  * Horinouchi and Ueda (2015) Y. Horinouchi and M. Ueda, Phys. Rev. Lett. 114, 025301 (2015).
  * Naidon and Endo (2017) P. Naidon and S. Endo, Rept. Prog. Phys. 80, 056001 (2017), arXiv:1610.09805 [quant-ph] .
  * Greene  _et al._ (2017) C. H. Greene, P. Giannakeas, and J. Perez-Rios, Rev. Mod. Phys. 89, 035006 (2017), arXiv:1704.02029 [cond-mat.quant-gas] .
  * Deltuva (2021) A. Deltuva, Phys. Rev. C 103, 064001 (2021), arXiv:2106.01458 [nucl-th] .
  * Rosa  _et al._ (2022) D. S. Rosa, T. Frederico, G. Krein, and M. T. Yamashita, Phys. Rev. A 106, 023311 (2022), arXiv:2209.00575 [physics.atom-ph] .
  * Frederico and Gattobigio (2023) T. Frederico and M. Gattobigio, (2023), arXiv:2303.14952 [physics.atm-clus] .
  * Pricoupenko (2023) L. Pricoupenko, (2023), arXiv:2303.12637 [cond-mat.quant-gas] .
  * 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] .
  * 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] .
  * 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] .
  * Jackura  _et al._ (2019a) 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 (2019a), 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 (2023) A. W. Jackura, Phys. Rev. D 108, 034505 (2023).
  * Dawid  _et al._ (2023) S. M. Dawid, M. H. E. Islam, and R. A. Briceño, Phys. Rev. D 108, 034016 (2023).
  * Lehmann  _et al._ (1955) H. Lehmann, K. Symanzik, and W. Zimmermann, Nuovo Cim. 1, 205 (1955).
  * Zimmermann (1958) W. Zimmermann, Nuovo Cim. 10, 597 (1958).
  * Fried (1960) Z. Fried,  _Bound states in the Lehmann, Symanzik and Zimmermann formulation of field theory_ , Ph.D. thesis (1960).
  * Duncan (2012) A. Duncan,  _The Conceptual Framework of Quantum Field Theory_ (Oxford University Press, 2012).
  * Jackura  _et al._ (2019b) 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 (2019b), arXiv:1809.10523 [hep-ph] .
  * 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] .
  * Landau (1959) L. Landau, Nuclear Physics 13, 181 (1959).
  * 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).
  * Hansen and Sharpe (2017) M. T. Hansen and S. R. Sharpe, Phys. Rev. D 95, 034501 (2017), arXiv:1609.04317 [hep-lat] .
  * Hyodo  _et al._ (2014) T. Hyodo, T. Hatsuda, and Y. Nishida, Phys. Rev. C 89, 032201 (2014), arXiv:1311.6289 [hep-ph] .
  * Adams  _et al._ (1998) G. S. Adams _et al._ (E852 Collaboration), 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._ , 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._ (2022a) R. Aaij _et al._ (LHCb), Nature Phys. 18, 751 (2022a), arXiv:2109.01038 [hep-ex] .
  * Aaij  _et al._ (2022b) R. Aaij _et al._ (LHCb), Nature Commun. 13, 3351 (2022b), arXiv:2109.01056 [hep-ex] .
  * Aaij  _et al._ (2014) R. Aaij _et al._ (LHCb), Phys. Rev. D 90, 112004 (2014), arXiv:1408.5373 [hep-ex] .
  * Cheng  _et al._ (2016) H.-Y. Cheng, C.-K. Chua, and Z.-Q. Zhang, Phys. Rev. D 94, 094015 (2016), arXiv:1607.08313 [hep-ph] .
  * Chang  _et al._ (2017) P. Chang, K.-F. Chen, and W.-S. Hou, Prog. Part. Nucl. Phys. 97, 261 (2017), arXiv:1708.03793 [hep-ph] .
  * LHC (2023) (2023), arXiv:2303.04062 [hep-ex] .
  * Brayshaw (1968a) D. D. Brayshaw, Phys. Rev. 176, 1855 (1968a).
  * Brayshaw (1968b) D. D. Brayshaw, Phys. Rev. 167, 1505 (1968b).
  * Glöckle (1978) W. Glöckle, Phys. Rev. C 18, 564 (1978).
  * Matsuyama and Yazaki (1991) A. Matsuyama and K. Yazaki, Nucl. Phys. A 534, 620 (1991).
  * Konishi  _et al._ (2017) A. Konishi, O. Morimatsu, and S. Yasui, (2017), arXiv:1705.02569 [hep-ph] .
  * Bedaque  _et al._ (1999b) P. F. Bedaque, H. W. Hammer, and U. van Kolck, Nucl. Phys. A 646, 444 (1999b), arXiv:nucl-th/9811046 .
  * Dietz  _et al._ (2022) S. Dietz, H.-W. Hammer, S. König, and A. Schwenk, Phys. Rev. C 105, 064002 (2022), arXiv:2109.11356 [nucl-th] .
  * Lindesay and Noyes (1980) J. V. Lindesay and H. P. Noyes, SLAC-PUB-2515 (1980).
  * Lindesay and Noyes (1986) J. V. Lindesay and H. P. Noyes, SLAC-PUB-2932 (1986).
  * Frederico (1992) T. Frederico, Phys. Lett. B 282, 409 (1992).
  * Carbonell and Karmanov (2003) J. Carbonell and V. A. Karmanov, Phys. Rev. C 67, 037001 (2003).
  * Ydrefors  _et al._ (2017) E. Ydrefors, J. H. Alvarenga Nogueira, V. Gigante, T. Frederico, and V. A. Karmanov, Phys. Lett. B 770, 131 (2017), arXiv:1703.07981 [nucl-th] .
  * Ydrefors  _et al._ (2019) E. Ydrefors, J. H. Alvarenga Nogueira, V. A. Karmanov, and T. Frederico, Phys. Lett. B 791, 276 (2019), arXiv:1903.01741 [hep-ph] .
  * Mohseni  _et al._ (2021) K. Mohseni, A. J. Chaves, D. R. da Costa, T. Frederico, and M. R. Hadizadeh, Phys. Lett. B 823, 136773 (2021), arXiv:2111.02015 [nucl-th] .
  * Frederico and Ydrefors (2021) T. Frederico and E. Ydrefors, Few Body Syst. 62, 8 (2021).
  * 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] .
  * Hansen and Sharpe (2016) M. T. Hansen and S. R. Sharpe, Phys. Rev. D 93, 096006 (2016), [Erratum: Phys.Rev.D 96, 039901 (2017)], arXiv:1602.00324 [hep-lat] .
  * Blanton  _et al._ (2019) T. D. Blanton, F. Romero-López, and S. R. Sharpe, JHEP 03, 106 (2019), arXiv:1901.07095 [hep-lat] .
  * Briceño  _et al._ (2018) R. A. Briceño, M. T. Hansen, and S. R. Sharpe, Phys. Rev. D 98, 014506 (2018), arXiv:1803.04169 [hep-lat] .

## Supplemental Material

###  .1 Three-body integral equations

Here, we review the key expressions needed to solve for the three-body
scattering amplitude. We employ the relativistic integral equations defined in
Ref. Hansen and Sharpe (2015) for the $\mathcal{M}_{3}$. We note that similar
three-body models were studied using non-relativistic Faddeev equations
(Brayshaw, 1968a, b; Glöckle, 1978; Matsuyama and Yazaki, 1991; Hyodo  _et
al._ , 2014; Konishi  _et al._ , 2017) and non-relativistic EFTs (Bedaque  _et
al._ , 1999a, b; Dietz  _et al._ , 2022). First attempts to understand the
system of three bosons in the relativistic framework have been made in Refs.
(Lindesay and Noyes, 1980, 1986). However, their result disagreed with the
subsequent Light-Front calculation (Frederico, 1992) obtained in an equivalent
regularization scheme. Its redefinition (Carbonell and Karmanov, 2003) led to
an approximate agreement between two approaches but also implied the
unphysical behavior of the trimer mass becoming imaginary at finite values of
$a$. We do not observe any indication of the Thomas collapse reported in this
work. One can find additional studies in Refs. (Ydrefors  _et al._ , 2017,
2019; Mohseni  _et al._ , 2021; Frederico and Ydrefors, 2021) and a summary of
the situation in Ref. (Naidon and Endo, 2017).

Strictly speaking, in this article, we work with the pair-spectator amplitude
$\mathcal{M}^{(u,u)}_{3}$, where specific spectators are chosen. The
superscript emphasizes that the initial and final states are “unsymmetrized”.
It becomes the genuine $3\varphi\to 3\varphi$ amplitude $\mathcal{M}_{3}$ only
after symmetrization with respect to different spectator choices (Hansen and
Sharpe, 2015). The amplitude is expressed as a sum of two terms,

$\displaystyle\mathcal{M}^{(u,u)}_{3}(p,k)=\mathcal{D}^{(u,u)}(p,k)+\mathcal{M}_{\text{df},3}^{(u,u)}(p,k)\,.$
(5)

The first term is called the _ladder_ amplitude and describes the scattering
driven by one-particle exchanges, with the three-body couplings “turned off”.
The second term depends on the three-body $K$ matrix, $\mathcal{K}_{3}$,
introduced in the main text. Compared to Refs. (Hansen and Sharpe, 2014, 2015)
we drop the “df” label from this object to simplify the notation. We also drop
the ${(u,u)}$ superscript from the $\mathcal{M}_{3}$, $\mathcal{D}$, and
$\mathcal{M}_{\text{df},3}$ remembering these are “unsymmetrized” objects.

Following Refs. Jackura  _et al._ (2021); Hansen  _et al._ (2021); Dawid  _et
al._ (2023), we perform the partial-wave projection, setting all relevant
angular momenta to zero. The equation for the partial-wave projected version
of $\mathcal{D}$ is,

$\displaystyle\mathcal{D}(p,k)=-\mathcal{M}_{2}(p)\,G(p,k)\,\mathcal{M}_{2}(k)-\mathcal{M}_{2}(p)\int\limits_{0}^{k_{\rm
max}}\\!\frac{{\rm d}k^{\prime}\,k^{\prime
2}}{(2\pi)^{2}\omega_{k^{\prime}}}\,G(p,k^{\prime})\,\mathcal{D}(k^{\prime},k)\,.$
(6)

All objects in the above equation implicitly depend on the total energy
variable, $E$. The one-particle exchange (OPE) amplitude $G$ is the boson
exchange propagator,

$\displaystyle
G(\bm{p},\bm{k})=\frac{H(p,k)}{(E-\omega_{p}-\omega_{k})^{2}-(\bm{p}+\bm{k})^{2}-m^{2}+i\epsilon}\,,$
(7)

which, after the partial-wave projection, becomes,

$\displaystyle G(p,k)$
$\displaystyle=-\frac{H(p,k)}{4pk}\,\log\left(\frac{z(p,k)+i\epsilon-2pk}{z(p,k)+i\epsilon+2pk}\right)\,.$
(8)

Here $z(p,k)=(E-\omega_{k}-\omega_{p})^{2}-k^{2}-p^{2}-m^{2}$ and $H(p,k)$ is
a cutoff function. In previous work Jackura  _et al._ (2021); Dawid  _et al._
(2023), we employed smooth regularization choice. Here, we exclusively
consider hard cutoff functions that are defined to be equal to $1$ in the
range of integration and $0$ otherwise. The maximal value the momentum takes
is given by $k_{\rm max}=\sqrt{E^{2}-m^{2}}/{2E}$ which corresponds to minimal
$\varepsilon_{k^{\prime}}=0$ of the pair in the intermediate three-body state.

The second term of Eq. (5), amplitude $\mathcal{M}_{\text{df},3}$, can be
written as a function of $\mathcal{D}$. Here we only consider the simplest
scenario (the so-called isotropic approximation) where the three-body $K$
matrix depends solely on the total energy of the system,
$\mathcal{K}_{3}(p,k)\equiv\mathcal{K}_{3}(E)$. The three-body matrix
$\mathcal{K}_{3}$ is practically an unconstrained function of $E$ that can
consist of a polynomial of arbitrary order and a sum of poles. In this case,
$\mathcal{M}_{\text{df},3}$ takes a form,

$\displaystyle\mathcal{M}_{\text{df},3}(p,k)$
$\displaystyle=\mathcal{L}(p)\,\frac{1}{\mathcal{K}_{3}^{-1}+F_{3}^{\infty}}\,\mathcal{L}(k)\,,$
(9)

where the “endcap” functions describing two-body rescattering contributions
are,

$\displaystyle\mathcal{L}(p)$
$\displaystyle=\frac{1}{3}-\mathcal{M}_{2}(p)\rho(p)-\int_{k^{\prime}}\mathcal{D}(p,k^{\prime})\rho(k^{\prime})\,.$
(10)

The three-body kinematic function, $F_{3}^{\infty}$, describes the effects of
propagation of on-shell three-body states interacting either within two-body
subchannels or via one-particle exchanges. It is expressed in terms of
$\mathcal{M}_{2}$ and $\mathcal{D}$,

$\displaystyle
F_{3}^{\infty}=\int_{k^{\prime}}\rho(k^{\prime})\,\mathcal{L}(k^{\prime})=\int_{k^{\prime}}\rho(k^{\prime})\left(\frac{1}{3}-\mathcal{M}_{2}(k^{\prime})\rho(k^{\prime})\right)-\int_{p^{\prime}}\int_{k^{\prime}}\rho(p^{\prime})\mathcal{D}(p^{\prime},k^{\prime})\rho(k^{\prime})\,.$
(11)

Finally, the phase space function, $\rho$, is defined as,

$\rho(k)=\,\frac{-i}{16\pi\varepsilon_{k}}\sqrt{\varepsilon_{k}^{2}/4-m^{2}}\,.$
(12)

Note that this differs slightly from what was presented in, for example, Refs.
(Hansen and Sharpe, 2014, 2015), where the authors used the smooth cutoff
function. As a result, the phase space was proportional to the cutoff
function. We also dropped the “3” subscript from the definition of the phase
space to simplify the notation.

###  .2 Cancellation of $\mathcal{K}_{3}=0$ poles

We note that separation into long- and short-range forces is ambiguous since
$\mathcal{K}_{3}$ is a regularization-scheme dependent object (Jackura, 2023).
Cutoff dependence of $\mathcal{D}$ is compensated by the changes of
$\mathcal{M}_{\text{df},3}$ under variations of $\mathcal{K}_{3}$, assuring
the scattering amplitude is independent of the regularization choice. In the
main text we consider the case $\mathcal{K}_{3}=0$, for which
$\mathcal{M}_{\text{df},3}=0$ and $\mathcal{M}_{3}=\mathcal{D}$. The purpose
of this and the subsequent section is to argue that this particular choice
does not affect our main conclusions. In particular, we would like to verify
that near the unitarity, $\mathcal{M}_{3}$ develops Efimov bound states for a
non-zero three-body $K$ matrix.

From Eq. (9), we see that $\mathcal{M}_{\text{df},3}$ contributes poles to the
full three-body amplitude whenever,

$\displaystyle\mathcal{K}_{3}^{-1}+F_{3}^{\infty}\Big{|}_{E=E_{n}}=0\,,$ (13)

As a first step, we would like to show that all the poles of $\mathcal{M}_{3}$
are those described by the above equation. Although it is evident that Eq.
(13) does indeed describe some states, it is not immediately clear that it
governs all of them. In particular, the first term in Eq. (5), $\mathcal{D}$,
can contribute additional poles to the full three-body amplitude.

Assume $\mathcal{D}$ has a pole at position $E_{n,0}$, i.e., it can be
expanded as,

$\displaystyle\mathcal{D}(p,k)$
$\displaystyle=-\frac{\Gamma_{n,0}(p)\Gamma_{n,0}(k)}{E^{2}-E_{n,0}^{2}}+\mathcal{O}\\!\left(E^{0}\right)\,,$
(14)

where the subscript “$0$” is meant to emphasize that this corresponds to the
$\mathcal{K}_{3}=0$ limit of $\mathcal{M}_{3}$. We note that when
$\mathcal{D}$ has a pole in $E^{2}$, from Eqs. (10), (11) it is clear that so
do $\mathcal{L}$ and $F_{3}^{\infty}$. In its vicinity, these functions take
the form,

$\displaystyle\mathcal{L}(p)$
$\displaystyle=\frac{\Gamma_{n,0}(p)}{E^{2}-E_{n,0}^{2}}\int_{k^{\prime}}\Gamma_{n,0}(k^{\prime})\rho(k^{\prime})\,+\mathcal{O}\\!\left(E^{0}\right),$
(15) $\displaystyle F_{3}^{\infty}$
$\displaystyle=\frac{1}{E^{2}-E_{n,0}^{2}}\left[\int_{k^{\prime}}\Gamma_{n,0}(k^{\prime})\rho(k^{\prime})\right]^{2}\,+\mathcal{O}\\!\left(E^{0}\right)\,.$
(16)

Using these expansions in Eq. (9) we find,

$\displaystyle\mathcal{M}_{{\rm df},3}(p,k)$
$\displaystyle=\frac{\Gamma_{n,0}(p)\Gamma_{n,0}(k)}{E^{2}-E_{n,0}^{2}}+\mathcal{O}\\!\left(E^{0}\right),$
(17)

which exactly cancels the pole contributions from $\mathcal{D}$ in Eq. (5). As
a result, the only poles present in $\mathcal{M}_{3}$ are those given by Eq.
(13).

###  .3 Behavior of $F_{3}^{\infty}$

Figure 4: Inverse $-1/F_{3}^{\infty}$ as a function of $\Delta E$, for
$ma=10^{6}$. Pole positions of $\mathcal{M}_{3}$ are determined by finding
values of $\Delta E$ at which $\mathcal{K}_{3}$ crosses the red curves. Two
typical forms of $\mathcal{K}_{3}$ are shown: a constant
$\mathcal{K}_{3}^{(0)}/m=2.5\cdot 10^{5}$ in brown, and a NLO threshold-
expanded isotropic version in orange (Hansen and Sharpe, 2016; Blanton  _et
al._ , 2019). For the latter: $\mathcal{K}_{3}^{(0)}/m=10^{5}$,
$\mathcal{K}_{3}^{(1)}/m=4.5\cdot 10^{9}$, and
$\Delta=(E^{2}-(E_{\text{thr.}}^{(3\varphi)})^{2})/(E_{\text{thr.}}^{(3\varphi)})^{2}$.
Dots represent binding energies in the limit $\mathcal{K}_{3}=0$ for which
$\mathcal{M}_{3}=\mathcal{D}$.

Equation (13) describes states appearing in the $\mathcal{K}_{3}=0$ limit
since the $F_{3}^{\infty}$ function has the same poles as $\mathcal{D}$, i.e.,
it diverges whenever $E=E_{n}$. Close to unitarity, by solving Eq. (Evolution
of Efimov States) numerically, we find that $F_{3}^{\infty}$ becomes a log-
periodic function of binding energy, $\Delta E$. Moreover, the sign of the
$F_{3}^{\infty}$’s first derivative is constrained below the threshold. These
features are shown in Fig. 4. Thus, in the $\mathcal{K}_{3}=0$ case, the
bounds states obey the discrete scaling symmetry characteristic of Efimov’s
phenomenon. As we go to higher $a$, poles of $F_{3}^{\infty}$ accumulate near
the energy threshold, and the energy gap between them becomes increasingly
small.

In that region, regardless of its functional form, $\mathcal{K}_{3}$ is to a
good approximation constant111Except for instances when it has a pole exactly
at the threshold., see an orange line in Fig. 4 for an example. Therefore,
near the threshold, we recover the Efimov scaling for any model of
$\mathcal{K}_{3}$—an effect ensuring the universality of our result. It
confirms the cutoff independence of this behavior at unitarity since various
models of the three-body $K$ matrix correspond to different regularization
choices.

###  .4 Analytic continuation

Here, we briefly describe the analytic continuation of the integral equation
defining the $\mathcal{D}(p,k)$ amplitude. We wish to extend it from the real
axis to the complex energy, which is achieved by the method of the integration
contour deformation. Generalization of the integration interval from a
straight real line to a complex curved path is motivated by the presence of
movable (energy-dependent) singularities of $G(p,k^{\prime})$,
$\mathcal{M}_{2}(k^{\prime})$, and $\mathcal{D}(k^{\prime},k)$ that can cross
real $k^{\prime}$ momentum axis when $E^{2}$ becomes complex-valued. These
singularities must be avoided at each value of $E^{2}$ and $p$; otherwise, the
final result becomes contour-dependent, invalidating the uniqueness of the
analytic continuation.

The main focus of this article, i.e., identification of the trimer states as
poles of the amplitude, requires a solution of the integral equation for real
energies below the ${\varphi b}$ and ${3\varphi}$ thresholds and for complex
energy in the unphysical Riemann sheets associated with these two branch cuts.
We first briefly summarize methods employed in Refs. (Jackura  _et al._ ,
2021; Dawid  _et al._ , 2023) to analyze the $\varphi b$ system and only then
describe the continuation of the amplitude through the three-body cut. We
focus on the $\mathcal{K}_{3}=0$ case for clarity of the presentation. Once we
compute the ladder amplitude from Eq. (6), we can obtain the full
$\mathcal{M}_{3}$ amplitude by utilizing formulas provided in Sec. .1 of this
supplement.

####  .4.1 The dimer-particle threshold

In Ref. (Jackura, 2023), we explained how the scattering amplitude between a
two-body bound state and a spectator, $\mathcal{M}_{\varphi b}$, is obtained
from the three-body scattering amplitude using the LSZ reduction formula. The
amplitude $\mathcal{D}(p,k)$ has poles at values of the pairs’ energies
squared, $\varepsilon_{p}^{2},\varepsilon_{k}^{2}$, equal to that of the bound
state, $m_{b}^{2}$. The residue of the three-body amplitude at these poles is
proportional to $\mathcal{M}_{\varphi b}$. Explicitly, expanding
$\mathcal{D}(p,k)$ in their vicinity, we find,

$\displaystyle\mathcal{D}(p,k)=\frac{g^{2}\,\mathcal{M}_{\varphi
b}(E)}{(\varepsilon_{p}^{2}-m_{b}^{2})(\varepsilon_{k}^{2}-m_{b}^{2})}+\cdots\,,$
(18)

where $g^{2}=128\pi m_{b}/a$. Continuing the external spectator momenta, $p$
and $k$, to the value corresponding to the two-body bound state pole, $q_{b}$,

$\displaystyle q_{b}=\frac{\lambda^{1/2}(E^{2},m_{b}^{2},m^{2})}{2E}\,,$ (19)

brings $\varepsilon_{p}^{2},\varepsilon_{k}^{2}\to m_{b}^{2}$, allowing to
extract the residue from the above formula. Here
$\lambda(x,y,z)=x^{2}+y^{2}+z^{2}-2xy-2xz-2yz$ is the Källén triangle
function.

Figure 5: Imaginary part of $\mathcal{M}_{2}(k^{\prime})$ in the complex
$k^{\prime}$ plane for $ma=-6$ and $s/m^{2}=9.01-0.01i$. Black lines represent
branch cuts, and blue dots are branch points. Three example integration
contours are shown. On the left panel, $\mathcal{C}_{1}$ (yellow, straight
line) is the original integration interval of Eq. (Evolution of Efimov States)
defining the solution, $\mathcal{M}_{3}(p,k)$, on the physical Riemann sheet
of the complex $E^{2}$ plane. To continue the three-body amplitude through the
three-body cut, one deforms $\mathcal{C}_{1}\to\mathcal{C}_{2}$. Amplitude
$\mathcal{M}_{2}$ in the integration kernel must be evaluated on its second
sheet when $k^{\prime}$ belongs to the green, dashed piece of
$\mathcal{C}_{2}$. On the right panel, we show an example contour,
$\mathcal{C}_{3}$, defining the ladder amplitude on the +2 unphysical Riemann
sheet of the $E^{2}$ plane.

The trimer bound and virtual states are poles of the $\mathcal{M}_{\varphi b}$
amplitude below the $\varphi b$ threshold that lie on the first and second
Riemann sheets, respectively. From the two-body $S$ matrix unitarity, the
analytic continuation of the amplitude to the unphysical Riemann sheet,
$\mathcal{M}_{\varphi b}^{\text{II}}$, is known explicitly,

$\displaystyle\mathcal{M}_{\varphi
b}^{\text{II}}(E)=\frac{\mathcal{M}_{\varphi b}(E)}{1+2i\rho_{\varphi
b}(E)\mathcal{M}_{\varphi b}(E)}\,,$ (20)

where $\rho_{\varphi b}(E)$ is a $\varphi b$ phase-space factor,

$\displaystyle\rho_{\varphi b}(E)=\frac{q_{b}}{8\pi E}\,.$ (21)

Thus, the knowledge of $\mathcal{M}_{\varphi b}$ on the first Riemann sheet is
sufficient to recover both the positions of three-body bound and virtual
states.

The main complication in extending the three-body amplitude below the dimer-
particle threshold $E_{\text{thr}}^{(\varphi b)}$ is the logarithmic branch
cut of the homogeneous and inhomogeneous terms of the integral equation in the
$p$ variable. It originates from the OPE partial-wave projected propagator,
Eq. (8), present in both of these terms and thus is an “inherited” singularity
of the $\mathcal{D}$ amplitude, appearing in its left-hand argument Dawid  _et
al._ (2023). At a fixed value of momentum $k$ and total energy, the
inhomogeneous term of the equation contributes a fixed cut that has to be
avoided by the deformed integration contour. In particular, for external
momentum $k=q_{b}$ and $E<E_{\text{thr}}^{(\varphi b)}$, the cut may take a
complicated shape resembling a circle, see Fig. 4 in Ref. (Dawid  _et al._ ,
2023).

The OPE amplitude in the homogeneous part of the equation is evaluated at
complex momenta $k^{\prime}$ along the contour chosen to avoid this fixed
quasi-circular shape. Since the right-hand argument of $G(p,k^{\prime})$ takes
infinitely many values, this amplitude contributes a whole set of
singularities to the left-argument dependence of the ladder amplitude that we
call a domain of non-analyticity. Note that this area of the complex
$k^{\prime}$ plane is defined by the specific integration path chosen. The
deformed contour must detour this region of the complex plane, a condition we
call self-consistency of the integration contour. Only self-consistent
integration paths define proper analytic continuation of the integral equation
solution. An example domain of non-analyticity is found in Fig. 8 of Ref.
(Dawid  _et al._ , 2023).

We refer the reader interested in more details to this work. Once we verified
that all singularities of the integration kernel and of the ladder amplitude
$\mathcal{D}(k^{\prime},k)$ are bypassed, we may attempt to solve the integral
equation numerically along the deformed contour, as described in App. C
therein.

####  .4.2 Three-particle threshold

Analytic continuation of the ladder amplitude through the three-body
(logarithmic) cut to the associated unphysical Riemann sheet requires further
discussion that is new to this work.

Figure 6: Imaginary part of the amplitude $m^{2}\,d(p,k)$ in the complex
$E^{2}$ plane for $ma=-8.1$ and $\varepsilon_{p}=\varepsilon_{k}=2m^{2}$. The
left panel represents the solution on sheet 0, while the right panel is
obtained via the analytic continuation to the nearest unphysical sheet.

To simplify the discussion, let us amputate the external two-body amplitudes
from the integral equation, obtaining,

$\displaystyle d(p,k)=-G(p,k)-\int\limits_{0}^{k_{\rm max}}\\!\frac{{\rm
d}k^{\prime}\,k^{\prime
2}}{(2\pi)^{2}\omega_{k^{\prime}}}\,G(p,k^{\prime})\,\mathcal{M}_{2}(k^{\prime})\,d(k^{\prime},k)\,,$
(22)

where $\mathcal{D}(p,k)=\mathcal{M}_{2}(p)\,d(p,k)\,\mathcal{M}_{2}(k)$. In
addition to deciding whether the behavior of the singularities of the OPE
amplitude demands integration contour deformation, one has to consider the
singularities of the two-body amplitude, $\mathcal{M}_{2}$, in the integration
kernel. As described in the main text, the singularities of the amplitude in
the $E$ variable emerge from the collision of the integration end point with
the singularities of $\mathcal{M}_{2}$ in the integration kernel. The three-
body right-hand cut of the three-body amplitude arises when the right-hand,
two-body branch point of $\mathcal{M}_{2}$ at

$\displaystyle k^{\prime}=\frac{\lambda^{1/2}(E^{2},m^{2},(2m)^{2})}{2E}$ (23)

coincides with the lower limit of integration, $k^{\prime}=0$. A particular
choice of the integration contour behavior around this branch point determines
the Riemann sheet of the obtained integral equation solutions; see Fig. 5.

Namely, when the integration contour circumvents the branch point from the
top, and the $\mathcal{M}_{2}$ is evaluated exclusively on its first Riemann
sheet, the integral equation defines the amplitude on the first (physical)
Riemann sheet above the three-particle threshold. It is not the only available
option. The integration contour can pass the branch point from the bottom. In
this case, the three-body is computed on its first unphysical (+1) Riemann
sheet of the complex $E$ plane. Since there is a cut extending from the
$\varepsilon_{k^{\prime}}=(2m)^{2}$ point, to properly define the analytic
continuation of the three-body amplitude, one must add a relevant
discontinuity function to $\mathcal{M}_{2}$, when evaluating the integration
kernel in this case, as indicated by the dashed green line in Fig. 5.

Finally, when the integration contour encircles the branch point $n$ times, as
shown in the right panel of Fig. 5, the amplitude is computed on the
$(1+n)^{\rm th}$ Riemann sheet, using the notation of Fig. 3. We remind the
reader that the three-body cut is logarithmic, i.e., contrary to algebraic
discontinuities such as in the square root function, it is associated with
infinitely many Riemann sheets. Winding around the branch point generates them
in the integral-equation representation of the three-body amplitude.

In Fig. 6 we present example solution of Eq. (22) for $ma=-8.1$ and external
momenta $p=k$ chosen such that
$\varepsilon_{k}^{2}=\varepsilon_{p}^{2}=2m^{2}$. The solution on the physical
Riemann sheet (left panel) exhibits a clear resonance “bump” above the real
axis, hinting at a pole on the nearest unphysical sheet. The right panel
presents the result of the analytic continuation to the +1 branch of the
amplitude, which reveals a ground state pole at position
$E_{1}^{2}/m^{2}\approx 9.0050-0.0012i$. It is a pole indicated by the red
color in Figs. 1 and 2 in the main text. Another pole, which we discuss in the
next section, is found in the upper half plane, at position $E/m^{2}\approx
8.9944+0.0087i$. We also observe another “bump” in the amplitude above the
real axis on sheet +1, a hint that another pole could be expected on sheet +2.

####  .4.3 Three-particle threshold — sheets $\geq+2$

Unfortunately, we found it practically impossible to define integration
contours of the type presented in the right panel of Fig. 5, that would be
self-consistent, i.e., that would avoid additional singularities generated by
the one-particle exchange amplitude. We note that the domain of non-
analyticity always seems to contain the branch point of
$\mathcal{M}_{2}(k^{\prime})$, meaning that it is unattainable to surround
encircle it with a loop-like integration contour without crossing it. A
rigorous solution to the problem of non-consistency of the integration paths
required to extend the three-body amplitude to sheets $\geq+2$ needs
additional study and is beyond the scope of this paper. As we note in the main
text, continuation to the higher Riemann sheets is most probably the solution
to the “missing pole” problem and thus will be of practical interest in the
future study of the three-body problem.

Figure 7: Evolution of the “quirky” pole in the upper half-plane of the +1
Riemann sheet. The pole approaches the threshold from complex infinity. As
$|ma|$ decreases, its trajectory forms a fractal pattern with rescaling
quotient converging to $\lambda^{2}$.

###  .5 Additional pole on sheet +1

In addition to the regular Efimov resonances in the +1 Riemann sheet of the
complex $E^{2}$ plane, we also find a new, “quirky” state in the upper half-
plane, $\operatorname{Im}E^{2}>0$. This trimer pole is visible on the right
panel of Fig. 6, in the upper left corner of the frame. We trace its
trajectory with changing two-body scattering length and present it in Fig. 7,
denoting its binding energy as $\Delta E_{q}$. We observe as it approaches the
three-body threshold from a deep region of the complex plane, supposedly the
complex infinity. Interestingly, as the magnitude of $ma$ increases, it
bypasses the three-body branch point and turns back. This pattern keeps
repeating for arbitrarily high values of $ma$. The fractal trajectory it forms
has a self-similarity factor converging to the Efimov’s constant,
$\lambda^{2}$.

We find this state puzzling, and at this point we can not draw any conclusions
on it nature. The pole approaches the three-body branch point in the unitarity
limit and could play a role in the accumulation of bound states poles at the
threshold. It could be a physically interesting phenomenon or simply an
artifact of the relativistic integral equations. At this time, neither option
can be ruled out, and further investigations are needed.

Figure 8: The $ma$ dependence of the normalized trimer residues. Real (left)
and imaginary (right) parts of the bound state residue is shown with a red,
solid line. Residues of resonance poles are presented as dashed lines. The
$+2$ sheet pole, obtained through the extrapolation procedure, is shown with
an error band estimate. Around $ma=-8.71$ the poles approach the threshold and
all residues converge to the same value, around $-0.061$.

###  .6 Matching of pole residues at the three-body threshold

As explained in the main text, by studying evolution of the trimer poles, we
found that the $n^{\text{th}}$ resonance on the +1 sheet reaches the
$3\varphi$ threshold at the same value of the two-body scattering length at
which the $n^{\text{th}}$ bound state emerges on the physical energy plane,
and the $(n+1)^{\text{th}}$ resonance appears on the +1 sheet. For instance,
for $n=1$ this occurs around $ma=-8.71$, while for $n=2$ around $ma=-203$.
Moreover, even though the residues behave in a seemingly unrelated manner away
from the $E_{\text{thr}}^{(3\varphi)}$, they converge to the same number at
the value of $ma$ where the transition between Riemann sheets occurs.

In Fig. 8, we present the numerical evidence of this surprising behavior for
the $n=1$ case. We fix $p=k$ by choosing
$\varepsilon_{p}=\varepsilon_{k}=2m^{2}$ and compute the residues of the 1st
bound state and the 1st and 2nd resonances. We normalize them by
$\mathcal{M}_{2}(p)^{2}$ to cancel enhancement from the two-body rescatterings
in the final and initial pairs.

Although we can not analytically continue the amplitude to the Riemann sheets
higher than the first unphysical one, for some values of the scattering
length, we see its enhancement above the real axis on the $+1$ sheet. It is
reminiscent of the narrow resonant “bump” and suggests the existence of a pole
in the lower half-plane of the $+2$ sheet. Working under the assumption that
this pole is located near the threshold at $ma=-8.71$, we perform a Breit-
Wigner fit to the $\mathcal{D}(p,p)$ amplitude $i\epsilon$ above the real
energy axis to extract its position and residue. We use the model,

$\displaystyle\mathcal{D}(p,p)=-\frac{\Gamma_{m}^{2}}{E^{2}+i\epsilon-
E_{m}^{2}}+f(E)\,,$ (24)

where we again fix $p=k$ such that $\varepsilon_{k}=\varepsilon_{p}=2m^{2}$.
The $\Gamma_{m}^{2}$ is a complex residue, $E_{m}$ is the pole position of the
modeled trimer, and $f(E)$ is a background quadratic function. We perform
several fits, differing by the value of $\epsilon$ and the form of $f(E)$. Our
model depends, in total, on four real parameters when $f(E)=0$ up to ten for
$f(E)=c_{0}+c_{1}E+c_{3}E^{2}$, where all $c_{i}$ are complex numbers. At each
$ma$, we find little variability in values of $\Gamma_{m}^{2}$ and $E_{m}^{2}$
between different fits and take the largest difference between any two of them
as the error estimate.

The pole position extracted in this way approaches the three-body threshold as
$ma\to-8.71$. In Fig. 8, we also show evidence that the residue of this pole
converges to the same value as the residues of the remaining ones, supporting
our conjecture about the missing poles. It indicates that the sheet $+2$ pole
is the 2nd Efimov state, as suggested in Fig. 3. However, we note that this
result is obtained in an approximate, model-dependent manner and is less
reliable than the exact analytic continuation of the amplitude.

###  .7 Scaling of bound-state energies and vertex functions

Here, we present additional evidence of the recovery of the discrete scaling
symmetry from the three-body integral equations. In agreement with Efimov’s
prediction, it emerges not only in the unitarity limit but for finite values
of $a$ as well. All characteristic features of the three-body states (e.g.,
energies, trajectories, pole residues) are with good precision described by
the same functions of $a$, modulo rescaling of dimensional quantities by
appropriate powers of $Q_{a}$, e.g., $a\to Q_{a}a$ or $E\to E/Q_{a}^{2}$.

Figure 9: Left panel: ratios of subsequent binding energies as the functions
of $ma$. The scattering length for the $\Delta E_{2}/\Delta E_{3}$ is rescaled
by $\lambda$, as indicated by the label on the horizontal axis. Efimov’s
ratio, $\lambda^{2}$, is presented as a gray solid line. Right panel:
normalized residues of the first three bound state poles of the
$\mathcal{M}_{\varphi b}$ amplitude as the functions of $ma$. The scattering
lengths for the subsequent residues are rescaled by appropriate powers of
$\lambda$, as indicated by the labels on the horizontal axis.

As expected, the bound states in the relativistic system behave like Efimov
states when their binding energies become small compared to the scale-
invariance breaking quantities like $m$ or the cutoff $k_{\text{max}}$. As
shown on the left panel of Fig. 9, the ratio of subsequent binding energies
approaches the square of Efimov’s constant in the unitary limit. Extrapolated
values are $0.7\%$ away from $\lambda^{2}$ for the $\Delta E_{1}/\Delta E_{2}$
ratio, and $10^{-3}$% for the $\Delta E_{2}/\Delta E_{3}$ one, which we
consider an excellent agreement. For finite values of $ma$, both ratios show
similar functional dependence on the scattering length. After we rescale it by
a factor of $\lambda$ for the second one, they overlap as shown in Fig. 10.
Although the energies exhibit dependence on the $\mathcal{K}_{3}$, their
ratios are mostly independent of its value, in agreement with the argument
made in Sec. .3 of this supplement.

On the right panel of Fig. 8, we show that the analogous property holds for
the residues associated with the trimer poles. At positive values of $a$, we
compute the $\mathcal{M}_{\varphi b}$ amplitude introduced in Eq. (18), and
extract the residues from the expansion,

$\displaystyle\mathcal{M}_{\varphi b}=-\frac{|\Gamma_{\varphi
b}|^{2}}{E^{2}-E_{n}^{2}}+\mathcal{O}(E)\,.$ (25)

Residue $|\Gamma_{\varphi b}|^{2}$ describes the coupling strength between the
trimer and the $\varphi b$ state. We normalize it with the two-body coupling
$g^{2}$ and plot it as a function of $ma$. As we can see, rescaling the
characteristic length scales brings all residues close to each other, with
almost no difference between the second and the third trimer and a
$\mathcal{O}(10\%)$ discrepancy between these two and the first one. Despite
the relatively small binding energy of the deep bound state, we understand it
as a remnant of relativistic effects in the formation of this state.

Finally, in Fig. 10, we reproduce numerically the known analytic form for the
non-relativistic vertex function, which was derived in Ref. Hansen and Sharpe
(2017) and confirmed numerically in Ref. Briceño  _et al._ (2018),

$\displaystyle|\Gamma(k)|^{2}=|c||A|^{2}\frac{256\pi^{5/2}}{3^{1/4}}\frac{m^{2}\kappa^{2}}{k^{2}(\kappa^{2}+3k^{2}/4)}\frac{\sin^{2}\Big{(}s_{0}\sinh^{-1}\big{(}\sqrt{3}k/2\kappa\big{)}\Big{)}}{\sinh^{2}(\pi
s_{0}/2)}\,.$ (26)

Here, $\kappa=-\sqrt{|m\Delta E|}$ and $s_{0}\approx 1.006$. Normalization
constant $|c|=96.351$, while $A$ is close to $1$ when $ma\to\infty$. The plot
in Fig. 10 was obtained for the first three trimers at $ma=-10^{6}$, which
corresponds to binding energies: $\Delta E_{1}/m=-2.8\cdot 10^{-2}$, $\Delta
E_{2}/m=-5.5\cdot 10^{-5}$, and $\Delta E_{3}/m=-1.1\cdot 10^{-7}$. The
presented agreement of the non-relativistic result with our finding is another
confirmation that the observed states are undoubtedly Efimov in nature. As
expected, we can observe a discrepancy between our numerical result and the
analytic formula for momenta that can be considered relativistic,
$p/m=\mathcal{O}(1)$.

Figure 10: Residues of the $\mathcal{M}_{3}$ amplitude at the first three
trimer poles for $ma=-10^{6}$. Dashed black lines represent the analytic
prediction of Ref. Hansen and Sharpe (2017).


# Buffering variability in cell regulation motifs close to criticality

Daniele Proverbio1,2<EMAIL_ADDRESS>Arthur N. Montanari1 Alexander
Skupin1,3,4 Jorge Gonçalves1,5 1- Luxembourg Centre for Systems Biomedicine,
University of Luxembourg, 6 Avenue du Swing, 4367, Belvaux, Luxembourg  2-
College of Engineering, Mathematics and Physical Sciences, University of
Exeter, EX4 4QL, Exeter, UK  3- Department of Physics and Material Science,
University of Luxembourg, 162a Avenue de la Faiencerie, 1511 Luxembourg,
Luxembourg  4- Department of Neuroscience, University of California San Diego,
9500 Gilman Drive, La Jolla, CA, United States  5- Department of Plant
Sciences, University of Cambridge, CB2 3EA, Cambridge, UK

###### Abstract

Bistable biological regulatory systems need to cope with stochastic noise to
fine-tune their function close to bifurcation points. Here, we study stability
properties of this regime in generic systems to demonstrate that cooperative
interactions buffer system variability, hampering noise-induced regime shifts.
Our analysis also shows that, in the considered cooperativity range, impending
regime shifts can be generically detected by statistical early warning signals
from distributional data. Our generic framework, based on minimal models, can
be used to extract robustness and variability properties of more complex
models and empirical data close to criticality.

Many biological systems self-regulate their functions through bistable
circuits, which have been associated to genetic Angeli  _et al._ (2004); Kheir
Gouda  _et al._ (2019) as well as growth feedbacks Deris  _et al._ (2013). In
particular, positive feedback loops have long been studied in systems and
synthetic biology De Mot  _et al._ (2016); Acar  _et al._ (2005); Guinn  _et
al._ (2020); they regulate crucial functions like enzymatic activity or gene
transcriptional changes during cell fate decisions Huang  _et al._ (2007);
Fiorentino  _et al._ (2020). Autoactivating positive feedback loops, simple
circuit motifs promoting bistability and fine regulation of dynamical states
close to self-organised criticality, are of particular importance Alon (2019);
Tripathi  _et al._ (2020). Cellular heterogeneity, i.e., random cell-to-cell
variations Komin and Skupin (2017), can further direct transitions Kaern  _et
al._ (2005); Weber and Buceta (2013) and induce regime shifts between
alternative stable states of gene expression or of protein concentrations
Thomas  _et al._ (2014). Positive feedback loops with stochastic fluctuations
have been observed in a variety of system including the transcription network
of E. coli Milo  _et al._ (2002) or in the regulation of $\beta$-galactosidase
Ozbudak  _et al._ (2004), which results from a sudden transition from low
(“off”) to high (“on”) level states of the lac operon at a critical point of
an inducer concentration.

There are mainly two ways in which bistable systems can switch between
alternative steady states Ashwin  _et al._ (2012): transitions driven by
bifurcations (which, due to loss of system resilience Scheffer  _et al._
(2009), may be anticipated by small random fluctuations) and transitions
driven by large random jumps. Cells and other biological systems are
hypothesized to live close to criticality to quickly respond to changing
environmental conditions Mora and Bialek (2011), but they should not respond
to random environmental changes (noise) in order to maintain their
evolutionary fitness. Close to criticality, the dynamical motifs have reduced
resilience and the system can exhibit increasing variability in response to
noise Mojtahedi  _et al._ (2016); Sharma  _et al._ (2016). This is typical of
nonlinear systems approaching a critical bifurcation and corresponds to
augmented sensitivity to random perturbations and diverging response time, a
phenomenon known as critical slowing down (CSD) Scholz  _et al._ (1987); Byrd
_et al._ (2019).

Mechanisms to buffer variability while maintaining the critical state are thus
necessary to finely regulate desired transitions Ozbudak  _et al._ (2004) or
to better cope with undesired shifts Scheffer  _et al._ (2009). To this end,
two main strategies can cooperate: moving the system state away from the
bifurcation point, or deepening the basin of attraction to avoid random
fluctuations pushing the system state to undesired attractors. These
strategies correspond to changes in different environmental or regulatory
conditions in cellular systems Dai  _et al._ (2015), allowing organisms to
exploit different mechanisms to buffer variability close to criticality. In
mono-stable systems, noise can be bound by the action of molecular compounds
like microRNAs Siciliano  _et al._ (2013) as well as by temporal relays of
signalling molecules Lestas  _et al._ (2010); Del Vecchio  _et al._ (2016).
For critical regimes in bistable processes, a key mechanism to buffer
variability is identified here: the cooperative interactions tuning the
activation function of positive feedback loops.

This buffering mechanism can be analysed by considering the simple and well-
known adimensional model for stochastic autoactivating positive feedbacks
Santillán (2008); Strogatz (2018):

$\dot{x}=f(x)+\eta(t)=K+c\frac{x^{n}}{1+x^{n}}-x+\eta(t)\;.$ (1)

This model describes the Michaelis-Menten kinetics of a transcriptional factor
activator (denoted by $x$) Frigola  _et al._ (2012); Weber and Buceta (2013),
which activates its own transcriptions when bound to a responsive element
(Fig. S1 sup ). System (1) arises from a model reduction of a two-variables
genetic toggle switch, under the assumption of slow-fast timescale separation
between the two variables Strogatz (2018). Here, $f(x)$ groups the
deterministic terms, with steady state $\tilde{x}$ ($\dot{x}|_{\tilde{x}}=0$),
$K$ is the basal expression rate, and $c$ is the maximum production rate with
critical value $c_{0}$ marking a saddle-node bifurcation (Fig. 1). The
dissociation constant in the denominator of the Hill function was normalised
to 1 without loss of generality Smolen  _et al._ (1998). The noise term
$\eta(t)$ accounts for intrinsic stochasticity of biological processes Hasty
_et al._ (2000). We consider additive Gaussian white noise to approximate the
fast degrees of freedom associated to a mean field regime Berglund and Gentz
(2006); Sharma  _et al._ (2016), with the statistical properties
$\langle\eta\rangle=0$,
$\langle\eta(t)\eta(t^{\prime})\rangle=2\sigma\delta(t-t^{\prime})$, where
$\sigma$ represents its intensity.

The Hill coefficient $n$, which describes the nonlinear cooperative binding
mechanisms, is usually interpreted as the number of transcription factors that
cooperatively promote transcription Santillán (2008). The smallest value
inducing bistability in the circuit is $n=2$, while $n\to\infty$ yields the
logic approximation for the activating Hill function,

$\lim_{n\to\infty}\frac{x^{n}}{1+x^{n}}=\Theta(x-1)\,,$ (2)

where $\Theta(\cdot)$ is the Heaviside step function, making system (1) a
discrete switch without bistability nor CSD.

This Letter investigates the dependence of resilience properties Dai  _et al._
(2015) on the cooperativity index $n$. This way, we test how cells can keep
production rates c close to their critical values and nonetheless increase
resilience and buffer variability using other regulation mechanisms. As shown
in Fig. 1, increasing $n$ from 2 yields different bifurcation diagrams, where
critical points shift to the left and the distance between the upper stable
manifold and the unstable manifold decreases when the system gets close to
criticality. We focus on systems residing on the upper branch (to be
consistent with the mean field assumption) and moving left towards the saddle-
node bifurcation point. This way, we investigate how biological circuits can
buffer variability close to critical states by exploiting dynamical
mechanisms. We also assess the parameter range where CSD-based early warning
signals correctly indicate impending regime shifts.

Figure 1: Bifurcation diagrams $(\tilde{x},c)$ of Eq. (1) for different Hill
coefficients $n$. Stable and unstable branches are represented by blue and red
colors, respectively. Black dots identify the bifurcation (saddle-node) points
$c_{0}$. Each value of $n$ is displayed close to its corresponding diagram.

To characterise the system stability properties, we analyse the stationary
potentials and probability density functions (PDF) depending on $n$, in
analogy to previous works Friedman  _et al._ (2006); Kumar  _et al._ (2014).
Consider the forward Fokker-Plank equation for the probability density
function $P(x,t)$ associated with Eq. (1):

$\frac{\partial P(x,t)}{\partial t}=-\frac{\partial}{\partial
x}\left[f(x)P(x,t)\right]+\frac{\partial^{2}}{\partial x^{2}}\left[\sigma
P(x,t)\right]\,,$ (3)

where $f(x)$ lumps the deterministic terms of Eq. (1). The stationary solution
$P_{s}(x)$ takes the form Gardiner (1985)

$\displaystyle P_{s}(x)$ $\displaystyle=N_{c}e^{-\phi(x)}\,,$ (4)
$\displaystyle\phi(x)$
$\displaystyle=\frac{1}{2}\ln\sigma-\frac{1}{\sigma}\int^{x}f(x^{\prime})dx^{\prime}\,,$
(5)

where $\phi(x)$ describes the adjoint stochastic potential whose depth is
related to system resilience, i.e., its ability to recover after a
perturbation. $N_{c}$ is a normalization constant such that
$\int_{\Omega}P_{s}(x)=1$ ($\Omega$ is the domain).

Figure 2: Effect of the Hill coefficient $n$ on the stochastic potential
$\phi(x)$ (left column) and on the stationary probability density function
$P_{s}(x)$ (right column), when only additive white noise is present. (a, b)
$c=0.8$ (“off” state); (c, d) $c=c_{0}(n)+0.05$ (multistable region); (e, f)
$c=2.7$ (“on” state favoured). In all cases, $K=0.1$ and $\sigma=0.05$.

Fig. 2 shows the dependency of $\phi(x)$ and $P_{s}(x)$ on $n$ when the system
is either in an “off” state far away from criticality (Fig. 2a,b), close to
the criticality (Fig. 2c,d) or beyond it, where the “on” state is favoured
(Fig. 2e,f). Increasing the cooperativity index $n$ does not alter the
underlying bistability, but modifies the depth of the potential and increases
the separation of alternative states (Fig. S2 sup ). For the “off” and “on”
states, the corresponding equilibria exhibit significantly deep attractor
basins with only minor dependence on $n$, as also indicated in the bifurcation
diagrams (Fig. 1) and by $P_{s}(x)$. Close to critical points, the picture
changes. The potential $\phi(x)$ displays two commensurable wells, which are
more evident and symmetric for larger $n$, suggesting that both states become
equally occupied in noisy environments. For increasing $n$, $P_{s}(x)$
displays sharper peak separation between the bistable states: the system
diffuses less to intermediate states and is more constrained around single
equilibrium values, as anticipated due the steeper potential barriers in
$\phi(x)$. Random deviations are thus suppressed faster and transitions from
one state to another are sharper and thus more robust against noise.

We now focus on how $n$ influences variability measures, like variance and
autocorrelation, close to criticality. Obtaining globally analytic expressions
is challenging, in particular for high values of $n$. Hence, we focus on a
local analysis close to the bifurcation points and employ a geometrical
methodology. To derive generic results for critical manifolds, we use their
local topological equivalence to bifurcation normal forms Kuehn and Bick
(2021). The normal forms associated to dynamical systems are simplified
minimal-order forms to which all systems exhibiting a certain type of
bifurcation are, around the equilibrium, topologically equivalent Kuznetsov
(2013). Supplementary Material sup  provides a brief background to normal
forms and terminology used in this work. For saddle-node bifurcations like in
Fig. 1, the associated normal form is Kuznetsov (2013)

$\dot{x}=p-x^{2}\,,$ (6)

with two equilibrium manifolds $\tilde{x}_{1,2}=\pm\sqrt{p}$, one stable ($+$)
and the other unstable ($-$). Note that the normal form corresponds to a
parabola. To study the behaviour of stochastic solutions near the stable
manifold, consider the evolution of its first-order perturbation, $y=\delta
x|_{\tilde{x}_{1}}$ exposed to the same additive white noise $\eta(t)$ as in
Eq. (1) Kuehn (2011). Since $k=2\sqrt{p}$ is the distance of the control
parameter value from its critical value $p_{0}=0$, note that $k$ is
proportional to $c-c_{0}$ from the original system, following normal form
properties. The corresponding Langevin equation accounting for mean field
fluctuations around the stable equilibrium is then given by

$\dot{y}=-ky+\eta(t)\,.$ (7)

Eq. (7) is a typical Ornstein–Uhlenbeck (OU) process with exact solutions for
statistical moments Gardiner (1985).

To connect the quantitative effects of $n$ with the more qualitative
topological form Eq. (6), recall that $n$ widens or narrows the local
parabolic shape of the original bifurcation diagram for Eq. (1) (Fig. 1). Eq.
(6) thus needs to be augmented with a term $\rho$ to modify the focal width of
its parabolic stable manifold, which corresponds to the width of the parabola
at the focal point. This leads to

$\dot{x}=p-\rho x^{2}\,.$ (8)

In this formulation, $\rho$ corresponds to the focal width of the normal form.
Supplementary Material sup  contains analytical derivations for the
approximation of system (1) to the normal form (8), and its relationships with
the geometrical results. Propagating $\rho$ into Eq. (7) adds a tuning term to
the bifurcation parameter, $k\to\sqrt{\rho}k$. Hence, the corresponding OU
process for a semi-quantitative saddle-node normal form is

$\dot{y}=-\sqrt{\rho}ky+\eta(t)\,.$ (9)

Among its statistical moments and power spectral properties, we are primarily
interested in quasi-steady-state variance (Var) and lag-1 autocorrelation
(AC1), measures of system variability close to criticality. They have been
proposed as proxies for system resilience and early warning signals (EWS) of
impending bifurcation points Scheffer  _et al._ (2009); Trefois  _et al._
(2015). Based on our mapping to the OU process (Eq. (9)), the analytical
solutions for Var and AC1 take the form Gardiner (1985):

$\text{Var}=\frac{\sigma}{\sqrt{\rho}k},\quad\text{AC1}=e^{-\sqrt{\rho}k}\,.$
(10)

Figure 3: Theoretical dependence of EWS measures (a) Var and (b) AC1 on $\rho$
(related to focal width) and $k$ (distance measure from critical parameter
values). Vertical lines represent slices for fixed values of $\hat{\rho}$ (Eq.
(13)) corresponding to the mean $F_{W}$ shown in panel c. The associated $n$
values increase from left to right. (c) Relationship between $F_{W}$ (Eq.
(11)) and corresponding values of $n$, with best fit. Error bars correspond to
one standard deviation (Eq. (12)). (d) Generic escape rate $\tau$ as a
function of noise level $\sigma$ and $k$ (Eq. (15)).

Eqs. (10) are generic for noisy saddle-node bifurcations. Fig. 3a,b shows them
as functions of $\rho$ and $k$. To connect with the original autoactivating
feedback system, we estimate the focal width of the bifurcation diagrams for
each $n$ by fitting a parabolic form
$c=\alpha\tilde{x}^{2}+\beta\tilde{x}+\gamma$ to the data points of each
bifurcation diagram in the vicinity of the saddle point. Using Matlab Curve
Fitting toolbox also provides uncertainties over
$\vec{\theta}=[\alpha,\beta,\gamma]$, resulting from small deviations from a
perfect parabolic shape. By definition, the fitted focal width is

$F_{W}=2|\tilde{x}(c_{F})-\tilde{x}_{F}|\,,$ (11)

where ($\tilde{x}_{F}$, $c_{F}$) are the coordinates of the parabolic focus.
To get a reasonable estimate of the corresponding uncertainties, the
associated standard deviation is derived from the fitted parameter
uncertainties $\text{std}(\bar{\theta}_{i})$ using a first-order approximated
propagation method Taylor (1997):

$\text{std}(F_{W})=\left[\sum_{i}\left(\frac{\partial
F_{W}}{\partial\theta_{i}}\text{std}(\theta_{i})\right)^{2}\right]^{\frac{1}{2}}\,.$
(12)

The relationship between $F_{W}$ and $n$ is plotted in Fig. 3c, with the
corresponding $\text{std}(F_{W})$. The pattern decreases quadratically,
thereby marking a rapid decrease followed by almost plateauing. Hence, a
bounded and relatively small cooperativity index is in principle sufficient to
effectively buffer variability close to criticality.

The estimated $\hat{\rho}$ values from fitted focal widths, for $n=2$ to
$n=8$, are obtained as

$\hat{\rho}=\xi(F_{W})^{-1}\,,$ (13)

where $\xi$ is a tuning parameter proportional to the Hill function
(Supplementary Material sup ). Mean $\hat{\rho}$ values are marked in Fig.
3a,b with solid vertical lines. Consistently with the trend observed in Fig.
3c, the mean values spread as $n$ increases (from left to right). Low $n$
values yield higher sensitivity to noise, as both Var and AC1 show
substantially higher values for small cooperativity indices $n$, even when $k$
is large (i.e., further away from the critical point, but still within the
bistable region, cf. Fig. 1). Thus, values of $\rho$ can belong to two
regions: one, where the values for both metrics are high for all $k$ (left
side of Fig. 3a,b), or another one where both metrics maintain low values for
most $k$ and increase rapidly close to criticality (right side of Fig. 3a,b).
The region $\hat{\rho}\to\infty$ corresponds to the logic approximation (2)
with $n\to\infty$, where Var and AC1 also change abruptly in a step-wise
manner. The ultra-sensitive region $\hat{\rho}\to 0$ is spanned by increasing
dissociation constants (Fig. S3 sup ), and potentially by changing other
parameters, here not explicitly considered, or by different activation
functions describing, for example, wild-type vs mutant organisms Hasty  _et
al._ (2000). Other pathways like growth feedbacks Deris  _et al._ (2013) will
likely correspond to additional regions in the parameter space. These
investigations are left to future studies.

We finally investigate the performance of EWS against impending bifurcation
points. The motivation is the following: consider complex systems lacking
validated mechanistic models; in our case, this would translate to a scenario
where $n$—or even the precise activation function—of an eukaryotic cell is
poorly identifiable Hasty  _et al._ (2000). This consideration leads to
questioning if we can identify statistical signals, computed on empirical
data, that provide reliable information about the system’s loss of resilience.
Increasing trends of Var and AC1 have been widely suggested to work as EWS
Scheffer  _et al._ (2009); Trefois  _et al._ (2015) but their robustness
remains elusive. To study how generic they are in the identified parameter
range and to account for mean trends and uncertainties, we numerically
integrate the original stochastic system (1) using the Euler-Maruyama method.
To mimic cell populations slowly evolving close to equilibrium, we sample
$10^{4}$ time points over 200 repeated experiments in dependence of $c$. This
leads to a distribution of statistical indicators (e.g., see Fig. 4a, inset).

To distinguish between bifurcation-induced transitions, anticipated by loss of
resilience, and noise-induced transitions, we measure the scale between the
distance to the bifurcation point and the noise level by the Kramers escape
rate
$\tau=2\pi(\sqrt{|U^{\prime\prime}(\tilde{x}_{1})U^{\prime\prime}(\tilde{x}_{2})|})^{-1}\exp[(U(\tilde{x}_{2})-U(\tilde{x}_{1}))/\sigma]$
Gardiner (1985). For any saddle-node bifurcation manifold (8) equipped with
additive noise, $U(\tilde{x}_{2})-U(\tilde{x}_{1})=(4/3)(k^{3}/\sqrt{\rho})$
and $|U^{\prime\prime}(\tilde{x}_{1,2})|=2\sqrt{p\rho}$. Hence,

$\tau=\frac{\pi}{k\sqrt{\rho}}\exp\left[\frac{4}{3}\frac{k^{3}}{\sqrt{\rho}\sigma}\right]\,.$
(14)

Values lying at the exponential boundary of

$\tau\simeq\mathcal{O}(\exp[k^{3}/\sigma])\,$ (15)

(Fig. 3d) provide comparable ranges of control parameters and noise levels for
all simulations with different $n$. They distinguish two regimes, one where
few noise-induced transitions might occur ($\tau\lesssim 2$, Fig. 4a,b) and
another regime primarily determined by bifurcation-driven resilience loss
($\tau\gtrsim 2$, Fig. 4c,d). For the considered $\sigma=0.02$, the system is
very close to critical points.

Figure 4: Average trends of (a) Var and (b) AC1 as a function of $c$ close to
bifurcation $c_{0}$, in the regime where noise-induced transitions might occur
(blue area of Fig. 3d). Simulations are presented over 200 realizations for
different $n$. For each $c$, the indicators spread into distributions, as
exemplified for two $c$ values in inset (a). (c,d) Average trends of Var and
AC1 farther from bifurcation $c_{0}$. (e,f) Evolution of p-values between Var
and AC1 distributions at each $c-c_{0}$ and the “reference” distribution. The
“reference” distribution corresponds to (e) $c-c_{0}=0.2$ (starting “closer to
bifurcation”) and (f) $c-c_{0}=0.35$ (“farther from bifurcation”). Dashed
lines represent typically used p-values in biological experiments.

Fig. 4a–d display average values for Var and AC1 from numerical simulations.
When the dynamics is mostly characterised by the bifurcation (Fig. 4c,d), both
measures display patterns consistent with those predicted in Fig. 3a,b and
increasing $n$ better buffers variability. When the noise level becomes
comparable to the potential depths, Vars for different $n$ become very close
to each other due to the more prominent role of noise-induced uncertainties
(Fig. 4a). By contrast, AC1s (Fig. 4b) remain separated due to their lower
sensitivity to noise (cf. Eq. (10)), but with less marked—and, therefore,
harder-to-detect—trends, similarly to those observed in real-world data
Proverbio  _et al._ (2022).

For online applications (i.e., as new data come in and without future
knowledge of the system evolution), it is necessary to quantify whether an
observed increasing trend is statistically significant, assessing whether it
corresponds to a EWS or some spurious fluctuation Boettiger and Hastings
(2012). To do so, we look for significant p-values between the computed
distributions close to criticality and those far from the bifurcation
(“reference”) (inset in Fig. 4a). Fig. 4e,f show patterns of p-values from Var
and AC1, averaged over all $n$. The p-values are computed closer to the
bifurcation in Fig. 4e ($c-c_{0}$ corresponding to the parameters in Fig.
4a,b) and farther from the bifurcation in Fig. 4f ($c-c_{0}$ corresponding to
the parameters in Fig. 4c,d). The p-values cross their significance levels
(either 0.1 or 0.05 Andrade (2019), dashed lines in Fig. 4e,f) before the
bifurcation point. This assesses that significant increasing trends of
statistical indicators can be detected prior to the transitions, thus
constituting reliable early warning signals. This analysis thus certifies the
potential use of proposed EWS to detect approaching bifurcation points in
biological motifs, providing a quantification of how much in advance the EWS
become significant depending on the reference and on the p-values threshold.

Overall, our study characterised fundamental dynamical mechanisms to buffer
systems’ variability in critical regimes. We determined parameter ranges,
corresponding to plausible cooperativity values for the positive feedback loop
motif, where both variance and autocorrelation display low relative
sensitivity to additive noise. In other ranges, however, the system poorly
buffers its variability. Investigating whether these ranges could correspond
to other dynamical mechanisms is demanded for future studies. Moreover, state-
dependent noise can be further incorporated in model (1) to make it closer to
biological reality Hasty  _et al._ (2000). Although it was not explicitly
considered in this paper, primarily to focus on the effects of a single
parameter on the system’s stability properties, in Supplementary Material sup
(also including Refs Bruggeman and Teusink (2018); Pan  _et al._ (2010); Zahri
(2014)) we investigate its influence on the variability metrics considered
above. Notably, it does not alter trends of AC(1) but affect those of Var
(Fig. S4 sup ), as expected from dependencies on $\sigma$ in Eq. (10). This is
relevant in further buffering fluctuation amplitudes near criticality and
calls for caution when using Var as an EWS indicator in systems strongly
characterised by such type of noise. A deep investigation of interplays
between noise types and resilience properties is left for future works. These
may unravel alternative ways by which cells regulate their states or support
the hypothesis of a self-organised fine-tuning in “safe” parameter spaces.
Overall, our analysis contributes with quantitative insights to analytical and
experimental studies of bistable systems’ resilience and connects general and
system-specific predictions Dai  _et al._ (2015).

We also assessed the sensitivity of proposed EWS to an additional regulation
mechanism and suggested how to on-line quantify significant increasing trends
from distributional data. We showed that extra parameters in the dynamical
system do not alter the warning capabilities of indicators associated with
CSD. In the considered parameters’ range, they are sufficiently generic to
detect resilience loss. As several indicators have been developed upon, to
detect cell-fate decisions Mojtahedi  _et al._ (2016) and to possibly
anticipate undesired shifts to, e.g., cancerous states Yang  _et al._ (2018);
Aihara  _et al._ (2022), our results constitute an important step to interpret
and apply them correctly. However, their use should be treated carefully if
other quantitative mechanisms between noise and bifurcations might be at play,
like certain types of state-dependent noise, possibly shadowing theoretical
trends. Following our methodology, future studies might inquire other
indicators and their behaviour under changing $n$ and additional conditions.
Our framework can also be easily extended to inquire the performance of
buffers and EWS in different dynamical models and experimental setups.

Acknowledgements The authors thank J. Fuentes for valuable discussions. D.P.
is supported by the Luxembourg National Research Fund (FNR) PRIDE DTU CriTiCS
(10907093) and A.S. by the FNR (C14/BM/7975668/CaSCAD) and by the NIH NBCR
(NIH P41 GM103426).

## References

  * Angeli  _et al._ (2004) D. Angeli, J. E. Ferrell, and E. D. Sontag, P. Natl. Acad. Sci. USA 101, 1822 (2004).
  * Kheir Gouda  _et al._ (2019) M. Kheir Gouda, M. Manhart, and G. Balázsi, Proceedings of the National Academy of Sciences 116, 25162 (2019).
  * Deris  _et al._ (2013) J. B. Deris, M. Kim, Z. Zhang, H. Okano, R. Hermsen, A. Groisman, and T. Hwa, Science 342, 1237435 (2013).
  * De Mot  _et al._ (2016) L. De Mot, D. Gonze, S. Bessonnard, C. Chazaud, A. Goldbeter, and G. Dupont, Biophys. J. 110, 710 (2016).
  * Acar  _et al._ (2005) M. Acar, A. Becskei, and A. Van Oudenaarden, Nature 435, 228 (2005).
  * Guinn  _et al._ (2020) M. T. Guinn, Y. Wan, S. Levovitz, D. Yang, M. R. Rosner, and G. Balázsi, Front. Genet. 11, 1292 (2020).
  * Huang  _et al._ (2007) S. Huang, Y.-P. Guo, G. May, and T. Enver, Develop. Biol. 305, 695 (2007).
  * Fiorentino  _et al._ (2020) J. Fiorentino, M.-E. Torres-Padilla, and A. Scialdone, Ann. Rev. Gen. 54, 167 (2020).
  * Alon (2019) U. Alon,  _An introduction to systems biology: design principles of biological circuits_ (CRC press, 2019).
  * Tripathi  _et al._ (2020) S. Tripathi, D. A. Kessler, and H. Levine, Phys. Rev. Lett. 125, 088101 (2020).
  * Komin and Skupin (2017) N. Komin and A. Skupin, Current Opinion in Systems Biology 3, 154 (2017).
  * Kaern  _et al._ (2005) M. Kaern, T. C. Elston, W. J. Blake, and J. J. Collins, Nat. Rev. Gen. 6, 451 (2005).
  * Weber and Buceta (2013) M. Weber and J. Buceta, PloS One 8, e73487 (2013).
  * Thomas  _et al._ (2014) P. Thomas, N. Popović, and R. Grima, P. Natl. Acad. Sci. USA 111, 6994 (2014).
  * Milo  _et al._ (2002) R. Milo, S. Shen-Orr, S. Itzkovitz, N. Kashtan, D. Chklovskii, and U. Alon, Science 298, 824 (2002).
  * Ozbudak  _et al._ (2004) E. M. Ozbudak, M. Thattai, H. N. Lim, B. I. Shraiman, and A. Van Oudenaarden, Nature 427, 737 (2004).
  * Ashwin  _et al._ (2012) P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox, Philos. T. R. Soc. A 370, 1166 (2012).
  * Scheffer  _et al._ (2009) M. Scheffer, J. Bascompte, W. A. Brock, V. Brovkin, and S. R. e. a. Carpenter, Nature 461, 53 (2009).
  * Mora and Bialek (2011) T. Mora and W. Bialek, Journal of Statistical Physics 144, 268 (2011).
  * Mojtahedi  _et al._ (2016) M. Mojtahedi, A. Skupin, J. Zhou, I. G. Castaño, R. Y. Leong-Quong, H. Chang, K. Trachana, A. Giuliani, and S. Huang, PLoS Biol. 14, e2000640 (2016).
  * Sharma  _et al._ (2016) Y. Sharma, P. S. Dutta, and A. Gupta, Phys. Rev. E 93, 032404 (2016).
  * Scholz  _et al._ (1987) J. Scholz, J. Kelso, and G. Schöner, Phys. Lett. A 123, 390 (1987).
  * Byrd  _et al._ (2019) T. A. Byrd, A. Erez, R. M. Vogel, C. Peterson, M. Vennettilli, G. Altan-Bonnet, and A. Mugler, Phys. Rev. E 100, 022415 (2019).
  * Dai  _et al._ (2015) L. Dai, K. S. Korolev, and J. Gore, P. Natl. Acad. Sci. USA 112, 10056 (2015).
  * Siciliano  _et al._ (2013) V. Siciliano, I. Garzilli, C. Fracassi, S. Criscuolo, S. Ventre, and D. Di Bernardo, Nat. Commun. 4, 1 (2013).
  * Lestas  _et al._ (2010) I. Lestas, G. Vinnicombe, and J. Paulsson, Nature 467, 174 (2010).
  * Del Vecchio  _et al._ (2016) D. Del Vecchio, A. J. Dy, and Y. Qian, Journal of The Royal Society Interface 13, 20160380 (2016).
  * Santillán (2008) M. Santillán, Math. Mod. Nat. Phen. 3, 85 (2008).
  * Strogatz (2018) S. H. Strogatz,  _Nonlinear dynamics and chaos_ (CRC press, 2018).
  * Frigola  _et al._ (2012) D. Frigola, L. Casanellas, J. M. Sancho, and M. Ibañes, PloS One 7, e31407 (2012).
  * (31) “See Supplementary Material at [URL] for supplemental background, calculations and figures.” .
  * Smolen  _et al._ (1998) P. Smolen, D. A. Baxter, and J. H. Byrne, Am. J. Physiol.-Cell Ph. 274, C531 (1998).
  * Hasty  _et al._ (2000) J. Hasty, J. Pradines, M. Dolnik, and J. J. Collins, P. Natl. Acad. Sci. USA 97, 2075 (2000).
  * Berglund and Gentz (2006) N. Berglund and B. Gentz,  _Noise-induced phenomena in slow-fast dynamical systems: a sample-paths approach_ (Springer Science & Business Media, 2006).
  * Friedman  _et al._ (2006) N. Friedman, L. Cai, and X. S. Xie, Phys. Rev. Lett. 97, 168302 (2006).
  * Kumar  _et al._ (2014) N. Kumar, T. Platini, and R. V. Kulkarni, Phys. Rev. Lett. 113, 268105 (2014).
  * Gardiner (1985) C. W. Gardiner,  _Handbook of Stochastic Methods_ (Springer, Boca Raton, 1985).
  * Kuehn and Bick (2021) C. Kuehn and C. Bick, Sci. Adv. 7, eabe3824 (2021).
  * Kuznetsov (2013) Y. A. Kuznetsov,  _Elements of applied bifurcation theory_, Vol. 112 (Springer Science & Business Media, 2013) p. 591.
  * Kuehn (2011) C. Kuehn, Physica D 240, 1020 (2011).
  * Trefois  _et al._ (2015) C. Trefois, P. M. Antony, J. Goncalves, A. Skupin, and R. Balling, Curr. Opin. Biotech. 34, 48 (2015).
  * Taylor (1997) J. R. Taylor,  _An Introduction to Error Analysis_ (University Science Books, Mill Valley, California, 1997).
  * Proverbio  _et al._ (2022) D. Proverbio, F. Kemp, S. Magni, and J. Gonçalves, PLOS Computational Biology 18, e1009958 (2022).
  * Boettiger and Hastings (2012) C. Boettiger and A. Hastings, Proceedings of the Royal Society B 279, 4734 (2012).
  * Andrade (2019) C. Andrade, Indian journal of psychological medicine 41, 210 (2019).
  * Bruggeman and Teusink (2018) F. J. Bruggeman and B. Teusink, Curr. Opin. Sys. Biol. 8, 144 (2018).
  * Pan  _et al._ (2010) W. Pan, Z. Wanf, H. Gao, Y. Li, and M. Du, Int. J. Robust Nonlin. Contr. 18, 557 (2010).
  * Zahri (2014) M. Zahri, J. Taibah Univ. Sci. 8, 186 (2014).
  * Yang  _et al._ (2018) B. Yang, M. Li, W. Tang, W. Liu, S. Zhang, L. Chen, and J. Xia, Nat. Commun. 9, 1 (2018).
  * Aihara  _et al._ (2022) K. Aihara, R. Liu, K. Koizumi, X. Liu, and L. Chen, Gene 808, 145997 (2022).