††thanks: joint first author††thanks: joint first author††thanks:
corresponding author<EMAIL_ADDRESS>
# Runaway Transition in Irreversible Polymer Condensation with Cyclisation:
Supplementary Information

Maria Panoukidou  School of Physics and Astronomy, University of Edinburgh,
Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK  Simon Weir  School of Physics
and Astronomy, University of Edinburgh, Peter Guthrie Tait Road, Edinburgh,
EH9 3FD, UK  Valerio Sorichetti  Laboratoire de Physique Théorique et Modèles
Statistiques (LPTMS), CNRS, Université Paris-Saclay, F-91405 Orsay, France
Yair Gutierrez Fosado  School of Physics and Astronomy, University of
Edinburgh, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK  Martin Lenz
Laboratoire de Physique Théorique et Modèles Statistiques (LPTMS), CNRS,
Université Paris-Saclay, F-91405 Orsay, France  PMMH, CNRS, ESPCI Paris, PSL
University, Sorbonne Université, Université de Paris, F-75005, Paris, France
Davide Michieletto  School of Physics and Astronomy, University of Edinburgh,
Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK  MRC Human Genetics Unit,
Institute of Genetics and Cancer, University of Edinburgh, Edinburgh EH4 2XU,
UK

## Methods

### Molecular Dynamics Simulations

#### Force-Field and equilibration process

We model a 6,500 bp-long linear DNA molecule as a bead-spring polymer made of
$l_{0}=174$ beads. The total number of polymer chains is $N_{m}=200$. The
polymers are modelled via the Kremer-Grest model [1]. Each bead has a diameter
$\sigma=13$ nm (or $\sim 38$ bp), modelled as a truncated and shifted Lennard-
Jones potential (WCA)

$U_{\text{LJ}}(r)=4\epsilon\left[(\sigma/r)^{12}-(\sigma/r)^{6}+1/4\right],$
(1)

for $r<r_{c}=2^{1/6}\sigma$ and 0 otherwise. Here $r$ represents the distance
between beads and $\epsilon=1.0$ parametrises the strength of the potential.
The radius of the bead, $\sigma$, defines the length units in our system.
Consecutive beads are connected through a permanent Finite Extensible Non-
linear Elastic (FENE) bond

$U_{\text{FENE}}(r)=-0.5KR_{0}^{2}\log{\left[1-\left(r/R_{0}\right)^{2}\right]}$
(2)

with $K=30\epsilon/\sigma^{2}$ and $R_{0}=1.5\sigma$. The bending stiffness of
the polymer is controlled by a Kratky–Porod interaction

$U_{b}(r)=\dfrac{k_{B}Tl_{p}}{\sigma}(1-\cos{\theta}),$ (3)

which regulates the angle ($\theta$) defined by the two tangent vectors
connecting three consecutive beads along the polymer. Here,
$l_{p}=4\sigma=150$ bp is the persistence length of DNA. The solvent is
simulated implicitly using a Langevin thermostat, so that the time evolution
of our system is governed by the stochastic partial differential equations

$m\ddot{\bm{r}}=-\xi\dot{\bm{r}}-\bm{\nabla}U+\sqrt{2k_{B}T\xi}\bm{\eta}$ (4)

where $\bm{r}$ is the position of a particle, $\xi$ its friction, $m$ its
mass, $U$ the sum of the interaction potentials discussed above and
$\bm{\eta}$ white noise with unit variance. The diffusion timescale is
$\tau_{B}=\xi\sigma_{b}^{2}/k_{B}T$. The integration of the Langevin equation
is done with a velocity-Verlet algorithm, using a time step $\Delta
t=0.012\tau_{B}$ in LAMMPS [2].

Various monomer densities were considered, ranging from $10^{-2}c^{*}$ to
$1c^{*}$, where $c^{*}=0.012\sigma^{-3}$ is the monomer concentration at which
the polymers start to overlap. The overlap concentration $c^{*}$ was measured
by computing the radius of gyration $R_{g}$ of the polymers in equilibrium at
very large dilution ($c^{*}=l_{0}/(4/3\pi R_{g}^{3}$)). All the systems were
equilibrated for a sufficient amount of time to ensure that the polymer chains
have moved at least a distance equal to $R_{g}$.

#### Topology reconstruction

After the equilibration step, 40 replicas of production runs were started for
each monomer number density considered. The ligation is performed
stochastically, and is attempted every 100 time steps between two polymer end
beads that are closer than $R_{c}=2^{1/6}\sigma$ using the `fix bond/create`
LAMMPS command. The probability of successful ligation (i.e., bond formation)
is set to $f=0.1$. Once ligated, the bond formed between the polymers is
irreversible and cannot be broken, therefore accounting for the formation of a
covalent bond between the DNA fragments. During the ligation process snapshots
of the system are taken every $10^{6}$ time steps on both the 3D coordinates
of the beads and the bond list at those time steps. From the bond list we can,
later on, reconstruct the topology of the individual polymers, i.e. if fused
with others to form linear chains or if circularised.

Once the production runs finished, the trajectories were analysed and the
topology was reconstructed using our python code
(https://git.ecdf.ed.ac.uk/taplab/dna-ligation.git). We will refer to this
code from now on as TR code. The code takes as input the instantaneous
trajectory and bond list from LAMMPS and checks for newly formed linear and
ring chains. The output of the python code is a file containing the number and
length of linear chains that have formed in a given simulation time step.
Similar files are produced for the ring chains. These files are then used to
calculate the average length and the number of linear/ring chains figures.

The starting point of the topology reconstruction algorithm is an array
$\mathbf{b}$ of size $N_{b}\times 2$; each row $b_{i}=(id_{1},id_{2})$
represents the IDs of atoms that are bonded within the system:

$\mathbf{b}=\begin{bmatrix}id_{1}&id_{2}\\\ id_{3}&id_{4}\\\ \vdots&\vdots\\\
id_{N_{n-1}}&id_{N_{n}}\end{bmatrix}\,,$

Since not all particles are linked together, some do not appear in the array
$\mathbf{b}$. In order to avoid operations with large sparse arrays, the
matrix $\mathbf{b}$ is mapped to $\mathbf{\tilde{b}}$ that contains only
indexes from $1$ to the maximum number of atoms connected
$M:\mathbf{b}\to\mathbf{\tilde{b}}$.

The next step of the TR algorithm is to create a connectivity matrix
$\mathbf{C}$ based on the list $\mathbf{\tilde{b}}$. Each row of $\mathbf{C}$
represents an atom index and consists of three components
$\mathbf{C}(id_{i},:)=(id_{i-1},id_{i+1},flag)$. Since in our case a particle
can be linked with two more particles the first two components of each row
$id_{i-1},id_{i+1}$ represent the connections of particle $id_{i}$ (note that
$id_{i-1}$ and $id+{i+1}$ are not necessarily consecutive in 1D but can be any
other particle bonded to particle $i$). The third component, $flag$, takes
only the values ${0,1}$ and accounts for the particles $id_{i}$ that already
belong to a polymer chain. The flag column of $\mathbf{C}$ is initialised to
zeros. During the reading process of the connectivity matrix the algorithm
switches the flags to 1 of the particles that are already considered to belong
in a chain. Rings are extracted in the same manner and a ring is found if the
current atom index is the same as the starting atom index. This reading
process outputs $N_{c}$ arrays which have different lengths and each of them
contains the particle (mapped) ids that are connected in a polymer chain. The
final step of the TR algorithm is to map the atom indexed back to the original
ones: $M^{-1}:\mathbf{\tilde{b}}\to\mathbf{b}$. This algorithm is very generic
and can be applied also in cases where the atoms are not initially in polymers
as in our case, but rather individual atoms that can connect with each other
during the simulation.

### Experiments

We consider a 6500 bp-long plasmid (referred to as 1288 here) which is
converted into linear form by using a restriction enzyme (XhoI). This
linearisation step is checked on a gel electrophoresis. The equilibrium radius
of gyration of this linear DNA molecule is about $R_{g}\simeq
l_{p}\sqrt{L/3l_{p}}\simeq 0.2$ $\mu$m (in agreement with diffusion data from
Ref. [3]). This yields an overlap concentration $c^{*}=3M/(4N_{A}\pi
R_{g}^{3})\simeq 0.2$ mg/ml where $N_{A}=6.023\,10^{23}$ 1/mol is the Avogadro
number and $M=L\times 650$ g/mol is the molecular weight of the plasmid. For
the low DNA concentration experiments we set ourselves at $0.01c^{*}$, i.e.
$c=2$ $\mu$g/ml. To perform ligation we use T4 ligase (NEB, M0202L, 1U
corresponds to 0.5 ng or 0.00735 pmoles of protein according to Ref. [4]), and
work at 1x T4 ligase reaction buffer concentration, which contains 1 mM ATP.

To classify the topology of the DNA under ligation, we perform time-resolved
gel electrophoresis. We prepare a master solution of DNA at desired
concentration, 1x ligase buffer and 2 U/$\mu$l T4 ligase. Aliquots are taken
at different elapsed times, the T4 ligase is heat-inactivated at 65∘ for 15
minutes. Aliquots are treated with Nb.BbvCI Nickase (NEB, R0631L) to relax the
supercoiled population. A subaliquot is created and treated with exonuclease
(PlasdmiSafe, Lucigen) to remove the linear population. The resulting aliquots
are run on a gel: we load 20ng of DNA from each aliquot onto a 1% agarose gel
prepared using 1x TAE buffer. A standard $\lambda$DNA - HindIII digest
(NEB,N3012S) marker is also loaded. The gel is run at $\sim$ 2.5V/cm for 5
hours and post-stained with SybrGold (ThermoFisher) for 30 minutes. A Syngene
G-box and Genesys software is used to image the gels.

The combination of nickase (relaxing the DNA supercoiling) and exonuclease
(fully digesting linear DNA molecules) allowed the topology of the DNA in each
band to be unambiguously identified. Further, the $\lambda$DNA - HindIII
digest marker confirmed the bands were of the correct size for monomer and
dimer lengths. To extract the relative amount of molecules in each lane we
compute, using ImageJ, the intensity of each lane and account for the fact
that the band with dimers has chains that are twice as long. We then normalise
against the sum of the three bands to obtain the relative fraction of chains
in each population.

The dynamics of the systems are measured using microrheology. Solutions are
made by mixing 8 $\mu$l of 1288 linearised plasmid at different concentrations
to a final concentration in the range 2ng/$\mu$l-500ng/$\mu$l with 1 $\mu$l of
40 U/ul T4 ligase and 1 $\mu$l of T4 ligase reaction buffer. Control solutions
are prepared at the same time and in the same manner substituting additional
TE for the T4 ligase. The samples are kept at room temperature on a roller for
several days. The samples are then spiked with $a=800$ nm PVP-coated
polystyrene beads, pipetted and sealed onto a slide and imaged using an
inverted microscope. We take a 30-minute movie and we analyse the movies using
a particle tracking algorithm (trackpy [5]) and extract the trajectories and
mean squared displacements(MSD) of the tracers $\langle\Delta
r^{2}(t)\rangle=\langle\left[\bm{r}(t+\tau)-\bm{r}(t)\right]^{2}\rangle$.
Diffusion coefficients are extracted by fitting to the MSD’s via MSD$=2Dt$.
The viscosity is obtained using the Stokes-Einstein relation [6],
$\eta=k_{B}T/(3\pi Da)$.

## Generic considerations on the proposed Smoluchowski equation

In this section we explain in detail how we solve the modified Smoluchowski
equation reported in the main text. The linear polymers undergo irreversible
ligation with rate $k_{1}(i,j)$ ($i,j$=lengths of the filaments reacting) and
cyclisation with rate $k_{0}(l)$. The concentrations of filaments of length
$l$ at time $t$, $\rho(l,t)$, and of rings, $\rho_{r}(l,t)$, are thus governed
by the following equations:

$\displaystyle\dot{\rho}(l,t)$
$\displaystyle=\frac{1}{2}\sum_{ij;i+j=l}k_{1}(i,j)\rho(i,t)\rho(j,t)+$
$\displaystyle-\rho(l,t)\sum_{i=1}^{\infty}k_{1}(l,i)\rho(i,t)-k_{0}(l)\rho(l,t)$
(5) $\displaystyle\dot{\rho}_{r}(l,t)$ $\displaystyle=k_{0}(l)\rho(l,t)\,.$
(6)

Once a chain undergoes cyclisation, it cannot react anymore and thus is
removed from the ensemble of linear chains and transferred to the ensemble of
rings. The dynamics is also constrained by the requirement that no monomers
are created or destroyed and the total mass is conserved:

$\sum_{l=1}^{\infty}l[\rho(l,t)+\rho_{r}(l,t)]=N/V=n\quad\forall t\,$ (7)

where $N$ is the total number of monomers and $V$ is the system’s volume.
Assuming that the reaction takes place on a time scale larger than the Rouse
relaxation time, the length-dependence of the annealing rate is [7]

$\displaystyle k_{1}(i,j)$
$\displaystyle=\kappa_{1}(D_{i}+D_{j})(R_{i}+R_{j})$ (8)
$\displaystyle=\kappa_{1}\left(l_{i}^{-\alpha}+l_{j}^{-\alpha}\right)\left(l_{i}^{\nu}+l_{j}^{\nu}\right)\,,$
(9)

where $\kappa_{1}$ is a constant coefficient and $l_{i}=il_{0}$ is the length
of a polymer with degree of polymerisation $i$, with $l_{0}$ the monomer size.
The cyclisation rate is taken to be $k_{0}(l)=\kappa_{0}l^{\mu}$, where
$\mu=-4\nu$. Note that we take $\mu=-4\nu$ instead of the Shimada-Yamakawa
exponent $\mu=-3\nu$ because we are dealing with a non-equilibrium, sub-
diffusive search process (see below). As we describe in detail later (see
“Looping probability and the exponent $\mu$” section), both theory and
simulations suggest that the _irreversible_ looping probability scales with
the length of the chain as $l^{-4\nu}$ and not as $l^{-3\nu}$.

### Time scaling of the mean length - concentrated regime

We can solve the Smoluchovski equation in the concentrated limit, where we
expect the ring formation to be negligible. We can thus use a mean-field
expression for the mean length $l$ and estimate (at scaling level) the rate as

$k\sim DR\sim l^{\nu-\alpha}\,.$ (10)

The timescale of chain growth is related to this rate and, in a mean-field
picture, we can write $d\rho(t)/dt=-k\rho(t)^{2}$ or

$k\sim\dfrac{1}{\rho(t)t}$ (11)

according to dimensionality arguments. In these equations, $\rho(t)$ is the
number density of chains with mean length $l$. This also means that

$t\sim\dfrac{1}{\rho(t)k}\sim\dfrac{l^{-\nu+\alpha}}{\rho(t)}\,.$ (12)

We now point out that the average number of fragments is related to the
average length of the fragments as $\rho(t)=l_{0}\rho(0)/l$. This is because
mass is conserved so $\rho(0)l_{0}=\sum_{i}\rho_{i}(t)(il_{0})$, or in mean
field, $\rho(0)l_{0}=\rho l$ at any time. Hence we find

$t\sim\dfrac{1}{\rho(t)k}\sim\dfrac{l^{-\nu+\alpha}}{\rho(t)}\sim\dfrac{l^{-\nu+\alpha+1}}{\rho(0)l_{0}}\sim
l^{1+\alpha-\nu}\,.$ (13)

To find how the mean length grows in time we can just invert this equation to
find

$l(t)\sim t^{1/(1+\alpha-\nu)}\sim t^{1/(1-\lambda)}\equiv t^{\gamma}\,,$ (14)

with $\lambda=\nu-\alpha$. Notice that $\alpha=1$ for Rouse dynamics and
$\alpha=2$ for reptation. Recall that here we are considering the diffusion in
3D of the chain, not its curvilinear diffusion, hence why $\alpha=2$ for
reptating chains. At the same time, we consider $\nu=1/2$ for ideal chains and
$\nu=3/5$ for self-avoiding chains. For instance, by plugging in these numbers
we obtain $l\sim t^{\gamma}$ with $\gamma\simeq 0.66-0.7$ for Rouse (so above
overlap but not entangled) and $\gamma\simeq 0.4$ for reptation (so above
overlap and entangled).

### Time-dependence of the mean length - dilute regime

In the very dilute regime, we can assume that rings and short n-mers should be
more favourable. We thus now solve the Smoluchovski equation analytically
assuming that only monomers, dimers and monomer rings are present. Let us
assume that the concentration of rings, monomers and dimers is denoted by
$n_{r},n_{1}$ and $n_{2}$, respectively. The dynamical system that describes
the evolution of those populations is the following:

$\displaystyle\dfrac{dn_{r}(t)}{dt}=k_{0}(1)n_{1}(t)$ (15)
$\displaystyle\dfrac{dn_{1}(t)}{dt}=-k_{1}(1,1)n_{1}^{2}(t)-k_{0}(1)n_{1}(t)$
(16) $\displaystyle\dfrac{dn_{2}(t)}{dt}=\dfrac{1}{2}k_{1}(1,1)n_{1}^{2}(t)$
(17)

We solve Eq. (16) neglecting the second term as $n_{1}^{2}\ll 1$ in the
infinite dilution limit:

$n_{1}(t)=n_{1}(0)e^{-k_{0}(1)t}$ (18)

The concentration of monomer rings is thus

$\dfrac{dn_{r}(t)}{dt}=k_{0}n_{1}(0)e^{-k_{0}(1)t}$ (19)

which yields

$n_{r}(t)=n_{1}(0)(1-e^{-k_{0}(l)t})\,.$ (20)

Substituting in Eq. (17), we get

$\dfrac{dn_{2}(t)}{dt}=\dfrac{1}{2}k_{11}n_{1}^{2}(t)=\dfrac{1}{2}k_{11}\left(n_{1}(0)e^{-k_{0}(1)t}\right)^{2}$
(21)

which yields

$n_{2}(t)=\dfrac{1}{4k_{0}}k_{12}n_{1}^{2}(0)\left(1-e^{-2k_{0}(1)t}\right)\,.$
(22)

Assuming these three are the only contributions to the system, then the mean
length is given by the following relation

$\displaystyle\langle l(t)\rangle$
$\displaystyle=\dfrac{l_{0}n_{1}(t)+l_{0}n_{o}(t)+2l_{0}n_{2}(t)}{n_{1}(t)+n_{o}(t)+n_{2}(t)}=$
$\displaystyle=l_{0}\dfrac{n_{1}(t)+n_{o}(t)+2n_{2}(t)}{n_{1}(t)+n_{o}(t)+n_{2}(t)}\,.$
(23)

In denser solutions, where the population is more polydisperse, the
Smoluchowski equation cannot be solved analytically and we refer to the
previous section for a scaling prediction and in the section “Looping
probability and the exponent $\mu$” for a perturbative approach in the limit
of small cyclisation. In the next section we explain the numerical method that
we implement to solve it.

### Time-dependence of the mean length - concentrated regime

To numerically solve the Smoluchowski equation Eq. (6) we wrote a code in
MATLAB using the explicit Euler scheme. The numerical evaluation of Eq. (6) is
then iterated to find the best free parameters $\kappa_{0}$ and $\kappa_{1}$
that fit the MD or experimental mean length versus time data. We recall the
expression of the ligation rate:
$k_{1}(i,j)=\kappa_{1}[(il_{0})^{\nu}+(jl_{0})^{\nu}][(il_{0})^{-\alpha}+(jl_{0})^{-\alpha}]$
and of the cyclisation rate $k_{0}(l)=\kappa_{0}l^{\mu}$. The fit is done
using the nonlinear least squares MATLAB embedded function (i.e. lsqcurvefit).
The rate of ring formation $\kappa_{0}$ and the rate of linear chains
formation $\kappa_{1}$ are extracted from this fit by considering four
independent sets of 10 independent simulations replicas so that we obtain an
error on the fits.

Solving the Smoluchowski equation to fit the data from MD simulations consists
of two main parts:

  1. 1.

We create an objective function for the lsqcurvefit (called Obj_smoluchowski)
that takes as input the array of initial coefficient guess
$K_{0}=(\kappa_{1},\kappa_{0})$ and the time data array $xdata$. It returns
the average length as a function of time, array $ydata$. In the objective
function:

    1. (a)

An array $\mathbf{L}=\mathbf{n}\cdot N_{0}$ is initialised where
$n=\\{1,2,\dots,N_{m}=200\\}$ and $N_{0}=174$. This represents the set of
lengths that can be found in the system (recall that we initialise our MD
simulations with 200 chains of 174 beads each). Also, the arrays with the
number density of linear and ring chains are initialised as follows,
$\mathbf{n_{l_{0}}}=(N_{m}/vol,0,\dots,0)_{1\times N_{m}}$ and
$\mathbf{n_{r_{0}}}=(0,\dots,0)_{1\times N_{m}}$ since initially all the
molecules are linear chains.

    2. (b)

_for_ $t$ = {1 to simulation final step time} do _call_
$(\mathbf{n_{L_{new}}},\mathbf{n_{R_{new}}})=exEuler\\_smoluchowski\left(\mathbf{n_{L}},\mathbf{n_{R}},K\right)$
function (see point 2 below)

    3. (c)

update arrays $\mathbf{n_{L}}=\mathbf{n_{L_{new}}}$ and
$\mathbf{n_{R}}=\mathbf{n_{R_{new}}}$. Calculate the total average length
$\mathbf{l_{total}}$ as

$\mathbf{l_{total}(t)}=\frac{\mathbf{n_{L_{new}}}\cdot\mathbf{L}+\mathbf{n_{R_{new}}}\cdot\mathbf{L}}{\sum_{i}^{N_{m}}{n_{L_{new}}^{i}}+\sum_{i}^{N_{m}}{n_{R_{new}}^{i}}}$

    4. (d)

_exit for loop_ and parse $\mathbf{l_{total}(t)}$ to $ydata$

  2. 2.

The exEuler_smoluchowski function takes as input the initial number densities
of linear and ring chains and the reaction rates
$\mathbf{n_{L}},\mathbf{n_{R}},K$. Based on the given rates K, it outputs
final number density arrays $\mathbf{n_{L_{new}}},\mathbf{n_{R_{new}}}$, after
the reactions have taken place. When this function is called, the number
density of linear and ring chains of each population are updated according to
Eq. (6). The monomer, dimer, and so on populations are increased according to
the first two terms of Eq. (6) while the number of them that is converted into
rings is subtracted by the $\mathbf{n_{L_{new}}}$ and added to the
$\mathbf{n_{R_{new}}}$ array.
In the first two terms of Eq. (6) the rate $k_{1}(i,j)$ is not a scalar
quantity by rather a matrix that follows the relation Eq. (9). The extracted
coefficient against which the fitting is optimised is the scalar $\kappa_{1}$.
Similarly, for the sink term of Eq. (6), the equation
$k_{0}(l)=\kappa_{0}l^{-4\nu}$ is used and the fitting coefficient exported is
the scalar $\kappa_{0}$.

The coefficients $K=(\kappa_{1},\kappa_{0})$ are updated iteratively by the
lsqcurvefit algorithm to best fit the data. Once the optimum values are
obtained the algorithm terminates.

### Looping probability and the exponent $\mu$

In equilibrium the looping probability of a chain of length $L$ is given by
the Shimada-Yamakawa formula [8, 9]. For $L\gg l_{p}$ (with $l_{p}$ the
persistence length) the looping probability decays as $P(L)\sim L^{\mu}$ with
$\mu=-3\nu$. This looping probability also holds for an irreversible, non-
equilibrium scenario if the process is reaction-limited. This is because the
chain ends would have the time to explore many conformations and to diffuse
the whole volume of the chain, $V\sim L^{3\nu}$, before reacting (as it would
happen in equilibrium). In a diffusion-limited, irreversible ligation process,
one should instead compute the time it takes for an end to diffuse over a
certain length distance $\xi$. The dynamics of the end is described by the
Rouse model [10], so that $\xi=b(D_{1}t/b^{2})^{1/4}$, where $b$ is the size
of a monomer and $D_{1}$ its free diffusion coefficient. Then, setting $\xi=R$
(the size of the chain) one obtains $(R/b)^{4}=(D_{1}t/b^{2})$, which implies
$k_{0}\sim t^{-1}\sim R^{-4}\sim L^{-4\nu}$. So considering $\mu=-4\nu$
effectively takes into account the fact that the chain ends are performing a
sub-diffusive search process within the polymer coil, as expected for Rouse
dynamics.

Figure 1: Computing the rate of loop formation as a function of polymer
length. (Left) $N_{\text{rings}}$ as a function of time for a system with
$N_{0}=400$ chains with $M=87$ monomers each. The black line represents the
average of the numerical derivatives of the curves in the limit $t\to 0$.
(Right) Numerical derivative $dN_{\text{rings}}/dt$ showing a power law decay
as a function of initial polymer length. The exponent of the power law is
found to be $-2.6$, close to the values $-4\nu$ which matches our calculation
assuming the irreversible loop formation to be a diffusion-limited and sub-
diffusive search process (see text).

We have used MD simulations to check the dependence on the polymer length of
our loop formation. To do this we ran short dilute simulations and starting
from polymers with different lengths, monitored the number of rings
$N_{\text{rings}}$ over time and then taken the derivative
$dN_{\text{rings}}/dt$ in the limit $t\to 0$. As shown in Fig. 1, the rate of
ring formation $dN_{\text{rings}}/dt=k_{0}(L)=\kappa_{0}L^{\mu}$ follows a
power law which we can fit with an exponent that is close to $\mu=-4\nu$, with
$\nu=0.588$ the exponent expected for self-avoiding chains. In light of this
evidence, we decided to use $\mu=-4\nu$ in the main paper.

## Calculation of mass converted into rings at infinite time

To estimate the amount of mass that is converted into rings at long times, we
do a perturbative calculation valid in the limit of small $\kappa$. We start
from the continuum Smoluchovski equation:

$\displaystyle\dfrac{d\rho(l,t)}{dt}=\dfrac{1}{2}\int_{0}^{l}K(i,l-i)\rho(i,t)\rho(l-i,t)di+$
$\displaystyle-\int_{0}^{\infty}K(i,l)\rho(l,t)\rho(i,t)di-\dfrac{1}{2}\kappa
l^{\mu}\rho(l,t)\,,$ (24)

where we define $K\equiv k_{1}/\kappa_{1}$ and recall that $\kappa\equiv
2\kappa_{0}/(n_{0}\kappa_{1})$ which is thus a scaling function such that
$K(ai,al)\sim a^{\lambda}k_{1}(i,l)$ where $\lambda=\nu-\alpha$ [11].

We now treat $\kappa_{0}$ perturbatively, starting with $\kappa_{0}=0$. In
this case, there is no mass lost into rings and we can thus write a
conservation law

$\int_{0}^{\infty}l\rho(l,t)dl=1\quad\forall t\,.$ (25)

Even for $\kappa$ non zero, we assume the loss of mass to cyclisation remains
finite and of order $\kappa$. We will check the self-consistency of this
assumption below. Using the mass conservation and Eq. (24) we can write the
following scaling relations: $l^{2}\rho=1$, $\rho
t^{-1}=l^{1+\lambda}\rho^{2}$. We therefore obtain:

$l\sim t^{1/(1-\lambda)}\,.$ (26)

which is the same as Eq. (14). Note that we must have $\lambda<1$ for the
average length of polymers to increase over time. We can also write the
density distribution as

$\rho\sim t^{-2/(1-\lambda)}$ (27)

which in the limit of long times or large lengths may be written as

$\rho(l,t)\sim
t^{-2/(1-\lambda)}\mathcal{G}\left(\dfrac{l}{t^{1/(1-\lambda)}}\right)$ (28)

where $\mathcal{G}$ is a scaling function that only depends on the ratio
$l/t^{1/(1-\lambda)}$.

We now introduce the ring length distribution $\rho_{r}(l,t)$ and its
evolution equation as

$\dfrac{d\rho_{r}(l,t)}{dt}=2\kappa_{0}l^{\mu}\rho(l,t)\,.$ (29)

Since at time $t=0$ there are no rings, we can then write

$\rho_{r}(l,t\to\infty)=2\kappa_{0}l^{\mu}\int_{0}^{\infty}\rho(l,t)dt\,.$
(30)

We can plug in the result we obtained for the distribution of length of linear
chains Eq. (28) to yield

$\displaystyle\rho_{r}(l,t\to\infty)$
$\displaystyle=2\kappa_{0}l^{\mu}\int_{0}^{\infty}t^{-2/(1-\lambda)}\mathcal{G}\left(\dfrac{l}{t^{1/(1-\lambda)}}\right)dt$
$\displaystyle=2\kappa_{0}l^{\mu}(1-\lambda)l^{-(1+\lambda)}\int_{0}^{\infty}x^{-\lambda}\mathcal{G}(x)dx$
(31)

where we defined $x=l/t^{1/(1-\lambda)}$. Thus, the number density of polymers
that are converted into rings over infinite time is

$\rho_{r}^{\infty}(l)=2\kappa_{0}(1-\lambda)l^{\mu-1-\lambda}\int_{0}^{\infty}x^{\lambda}\mathcal{G}(x)dx\,.$
(32)

Since $\lambda<1$ and assuming the $\mathcal{G}(x)=\mathcal{O}(1)$ when $x\to
0$, the integral converges in 0. For convergence of this integral at $\infty$
we also require that the scaling function decays faster than $x^{\lambda-1}$.

Assuming this functional form for the distribution of ring lengths at infinite
time, we now compute the total average mass transformed into rings at infinite
time as

$\displaystyle M_{\text{rings}}^{\infty}$
$\displaystyle=\int_{1}^{\infty}l\rho_{r}^{\infty}(l,t)dl$
$\displaystyle=2\kappa_{0}(1-\lambda)\int_{0}^{\infty}x^{-\lambda}\mathcal{G}(x)dx\int_{1}^{\infty}l^{\mu-\lambda}dl\,.$
(33)

The convergence of this integral requires that $\lambda-\mu>1$ and in this
case we get

$M_{\text{rings}}^{\infty}=2\kappa_{0}\dfrac{1-\lambda}{\lambda-\mu-1}\int_{0}^{\infty}x^{-\lambda}\mathcal{G}(x)dx\,.$
(34)

From this equation we see that the fraction of mass in rings at infinite time
$M_{\text{rings}}^{\infty}/M_{0}$ converges to a finite value proportional to
$\kappa_{0}$ (and hence $<1$ at small $\kappa_{0}$).

With this calculation we have thus shown that at small enough but non-zero
$\kappa_{0}$, the fraction of mass turning into rings is finite if
$\lambda-\mu=\nu-\alpha+4\nu>1$ or $\nu>\alpha/5$ which is valid for any type
of polymer in the non-entangled ($\alpha=1$) regime, and for ideal and self-
avoiding polymers in the entangled regime ($\alpha=2$). This implies that in
these regimes we expect the cyclisation probability to decay fast enough and
cannot prevent the runaway of the $M_{0}-M_{\text{rings}}^{\infty}$ mass into
linear chains that keep growing in time.

Consistently with this, in both MD and MC simulations we never observe the
formation of rings larger than 10 initial monomers.

## Description of the DSMC algorithm

The modified Smoluchowski equation proposed here (6) can only be solved
analytically for certain forms of the ligation rate $k_{1}(i,j)$ and of the
cyclisation rate $k_{0}(l)$. Moreover, the computational cost of directly
solving the system of ordinary differential equations as described in Sec.
Time-dependence of the mean length - concentrated regime increases steeply
with the system size. Thus, in order to find solutions to the Smoluchowski
equation for larger system sizes, we employ here Direct Simulation Monte Carlo
(DSMC) [12, 13, 14]. DSMC is a powerful stochastic method to solve
differential equations such as Eq. (6), and which samples the correct ligation
kinetics in the limit of large system sizes. The starting point for the Monte
Carlo algorithm is an array $\mathbf{m}$ of length $N$, each element $i$ of
which contains a number $m_{i}$ which represents the mass/length of the chains
$i$:

$\mathbf{m}=(m_{1},m_{2},\dots,m_{N})\,.$

A value of $0$ corresponds to absence of a certain chain. Moreover, to satisfy
mass conservation we ensure that $\sum_{i=1}^{N}m_{i}=N$ (total number of
monomers) is true for any time during the simulation. We will also consider an
analogous array of length $N$ denoted $\mathbf{r}$ (initially empty), where we
save the masses of the rings.

For an initially monodisperse condition, we set
$\mathbf{m}_{0}=(1,1,\dots,1)$. After the array m is initialized, we run the
DSMC simulation, which consists in repeating a large number of times a Monte
Carlo step (described below). The execution is terminated when the the system
has reached a state in which there is a single linear chain and a number of
non-reactive rings.

#### Description of the Monte Carlo step

In this paragraph, we will describe the single Monte Carlo (MC) step, which is
repeated a large number of times during the numerical resolution of Eq. (6)
performed using the DSMC algorithm. With reference to Eq. (6), we define
$n_{f}\equiv V\sum_{i=1}^{\infty}n_{i}$ (total number of chains). Before the
start of the simulation, we give an estimate of the maximum annealing rate
$k_{\text{max}}$ and of the maximum cyclisation rate $k_{0}^{\text{max}}$. The
exactness of the algorithm does not depend on this initial choice, however,
choosing values that are too far from the actual maximum rates can lead to a
reduced efficiency [12].

During every MC step, we either attempt to perform an ligation reaction (with
probability $p$) or a cyclisation one (with probability $1-p$). The value of
$p$ is calculated initially and then updated during the course of the
simulation in such a way that the average number densities $n(l)$ satisfy (6).
At the beginning of each MC step, $p$ is evaluated as

$p^{-1}=1+\frac{2Nk_{o}^{\text{max}}}{(n_{f}-1)nk_{\text{max}}}\,.$ (35)

We will show below that this choice also guarantees that the simulation
samples the correct number of cyclisation and ligation events per unit volume
and unit time as required by Eq. (6).

We define a waiting time variable that is set to zero at the beginning of the
simulation. After each reaction, a waiting time increment is generated. These
increments are also chosen in order to guarantee the correct number of
ligation and fragmentation reactions per unit time/volume, as detailed below.
We can now describe the MC step, during which the following actions are
performed:

  1. 1.

We evaluate the probability of annealing $p$ according to Eq. (35). The
explicit form of $p$, Eq. (35), will be discussed in detail below.

  2. 2.

We pick a random number $0\leq r\leq 1$ from a uniform distribution. If $r\leq
p$, we attempt a ligation event:

    1. (a)

We pick a pair of elements of the array $\mathbf{m}$, denoted $\alpha,\beta$
at random. Since there are $n_{f}(n_{f}-1)$ ordered pairs of chains, the
probability to pick a specific pair is $[n_{f}(n_{f}-1)]^{-1}$. Let the length
associated with these elements be $m_{\alpha}=i$ and $m_{\beta}=j$.

    2. (b)

We evaluate the ligation rate $k_{1}(i,j)$ for the two chains. If
$k_{1}(i,j)>k_{\text{max}}$, we set $k_{1}^{\text{max}}=k_{1}(i,j)$ and return
to (1). Otherwise, we continue.

    3. (c)

We pick another random number $r^{\prime}$, and perform the ligation if
$r^{\prime}\leq k_{1}(i,j)/k_{1}^{\text{max}}$. If ligation is unsuccessful,
we return to (1). Otherwise, we continue.

    4. (d)

We increment the waiting time by $\Delta
t^{\text{lig}}_{i,j}=\frac{2AN}{n_{f}(n_{f}-1)nk_{ij}}$. Here $A$ is a
parameter, the only condition on which is that it must be between $0$ and $1$,
as we will discuss in more detail below.

    5. (e)

After incrementing the waiting time, we update $\mathbf{f}$ by setting
$m_{\alpha}=0$ and $m_{\beta}=i+j$.

  3. 3.

If $r>p$, we attempt a cyclisation event:

    1. (a)

We pick a chain $\gamma$ at random with probability $n_{f}^{-1}$. Let
$m_{\gamma}=l$.

    2. (b)

We evaluate the cyclisation rate $k_{0}(l)$. If $k_{0}(l)>k_{0}^{\text{max}}$,
set $k_{0}^{\text{max}}=k_{0}(l)$ and return to (1). Otherwise, we continue.

    3. (c)

We extract another random number $0\leq r^{\prime}\leq 1$ from a uniform
distribution, and perform cyclisation if $r^{\prime}\leq
k_{0}(l)/k_{0}^{\text{max}}$. If cyclisation is unsuccessful, we return to
(1). Otherwise, we continue.

    4. (d)

We increment the waiting time by $\Delta
t^{\text{cyc}}_{l}=\frac{1-A}{n_{f}k_{o}(l)}$, with $A$ defined above in step
(2).

    5. (e)

We record the value of $l$ in $\mathbf{r}$ and set $m_{\gamma}=0$.

We now prove that the definitions of $p$ (Eq. (35)), the waiting time
increments $\Delta t^{\text{lig}}_{i,j}$ (for ligation) and $\Delta
t^{\text{cyc}}_{i,k-i}$ (for cyclisation) give a number of ligation and
cyclisation events per unit time which is consistent with the Smoluchowski
equation Eq. (6). Over a single MC step, the mean number of ligation events
involving the ordered pair of filaments $(\alpha,\beta)$ is

$\langle\\#L_{\alpha,\beta}\rangle\equiv\frac{p}{n_{f}(n_{f}-1)}\frac{k_{1}(m_{\alpha},m_{\beta})}{k_{1}^{\text{max}}}\,.$
(36)

We note that in the algorithm we consider $(m_{\alpha},m_{\beta})$ as an
ordered pair, and thus in (36) we consider the reaction $(i,j)\to l$ as
distinct from $(j,i)\to l$. The mean number of ligation events involving any
two chains with lengths $i,j$ can be obtained by multiplying the above
quantity by $2(1-\delta_{ij}/2)V^{2}n_{i}n_{j}$. The factor
$2(1-\delta_{ij}/2)$ takes into account the fact that, as mentioned above, for
$i\neq j$, there are two ways to perform the ligation, whereas for $i=j$ there
is only one. The factor $V^{2}n_{i}n_{j}$ is the product of the volume
fractions of filaments of lengths $i$ and $j$. We thus have

$\displaystyle
2V^{2}n_{i}n_{j}\left(1-\frac{\delta_{ij}}{2}\right)\times\frac{p}{n_{f}(n_{f}-1)}\frac{k_{1}(i,j)}{k_{1}^{\text{max}}}=$
$\displaystyle Vk_{1}(i,j)n_{i}n_{j}\left(1-\frac{\delta_{ij}}{2}\right)\Delta
t\,,$ (37)

where we have equated the mean number of ligation events involving any two
chains with lengths $i,j$ to the value required by the Smoluchowski equation.
Recalling that $n=N/V$, we thus find

$\Delta t=\frac{2pN}{n_{f}(n_{f}-1)nk_{1}^{\text{max}}}\,.$ (38)

Eq. (38) relates the time interval $\Delta t$ to the probability of ligation.
We will now obtain a second equality involving $p$ and $\Delta t$, which will
allow us to prove that the expression Eq. (35) for $p$ guarantees the correct
number of ligation and cyclisation events per unit time.

The mean number of cyclisation events involving chains $\gamma$ is

$\langle\\#C_{\gamma}\rangle\equiv\frac{(1-p)k_{0}(m_{\gamma})}{n_{f}k_{0}^{\text{max}}}\,.$
(39)

To obtain the mean number of cyclisations of a generic $l-$mer we need to
multiply this quantity by $Vn_{l}$, i.e., the volume fraction of filaments of
length $l$. Equating this quantity to the expected number of rings formed in a
time interval $\Delta t$ we obtain

$Vn_{l}\times\frac{(1-p)k_{0}(l)}{n_{f}k_{0}^{\text{max}}}=k_{0}(l)n_{l}V\Delta
t\,,$ (40)

and hence

$\Delta t\equiv\frac{1-p}{n_{f}k_{o}^{\text{max}}}\,.$ (41)

By equating the two expressions for $\Delta t$, Eq. (38) and Eq. (41), we find
Eq. (35). We have thus proven that the latter is the correct expression of
$p$, which gives the correct number of cyclisation and ligation events per
unit time and unit volume, as required by the Smoluchowski equation.

Finally, we will prove below that the constants $A$ and $1-A$ introduced when
calculating the waiting time increments are consistent with Eq. (38) and Eq.
(41). To show this, it is sufficient to observe that the total time increment
during an MC step is:

$\displaystyle\Delta t$ $\displaystyle=\sum_{0\leq\alpha<\beta\leq
n_{f}-1}\langle\\#L_{\alpha,\beta}\rangle\Delta
t^{\text{lig}}_{m_{\alpha},m_{\beta}}+\sum_{i=1}^{m_{\gamma}-1}\langle\\#C_{\gamma}\rangle\Delta
t^{\text{cyc}}_{m_{\gamma}}$ $\displaystyle=\sum_{0\leq\alpha<\beta\leq
n_{f}-1}\left[\frac{pk_{m_{\alpha},m_{\beta}}}{n_{f}(n_{f}-1)k_{\text{max}}}\right]\left[\frac{2AN}{n_{f}(n_{f}-1)nk_{m_{\alpha},m_{\beta}}}\right]$
$\displaystyle\ \ \
+\sum_{i=1}^{m_{\gamma}-1}\left[\frac{(1-p)k_{o}(l)}{nk_{o}^{\text{max}}}\right]\frac{1-A}{n_{f}k_{o}(l)}$
$\displaystyle=\frac{2ApN}{n_{f}(n_{f}-1)nk_{\text{max}}}+\frac{(1-A)(1-p)}{n_{f}k_{o}^{\text{max}}}$

One can see that this equality is consistent with Eq.(38) and Eq.(41). We note
that the algorithm samples on average the correct kinetics independently of
the value of $A$, as long as $0\leq A\leq 1$. Here we take $A=1$, meaning that
the waiting time increment is calculated only after a successful ligation
reaction, but not after a successful cyclisation reaction.

## References

  * Kremer and Grest [1990] K. Kremer and G. S. Grest, The Journal of Chemical Physics 92, 5057 (1990).
  * Plimpton [1995] S. Plimpton, J. Comp. Phys. 117, 1 (1995).
  * Robertson  _et al._ [2006] R. M. Robertson, S. Laib, and D. E. Smith, Proc. Natl. Acad. Sci. USA 103, 7310 (2006).
  * Taylor and Hagerman [1990] W. H. Taylor and P. J. Hagerman, Journal of Molecular Biology 212, 363 (1990).
  * Crocker  _et al._ [2000] J. C. Crocker, M. T. Valentine, E. R. Weeks, T. Gisler, P. D. Kaplan, A. G. Yodh, and D. A. Weitz, Phys. Rev. Lett. 85, 888 (2000).
  * Hansen and McDonald [2013] J.-P. Hansen and I. R. McDonald,  _Theory of simple liquids: with applications to soft matter_ (Academic press, 2013).
  * De Gennes [1982] P. G. De Gennes, The Journal of Chemical Physics 76, 3316 (1982).
  * Shimada and Yamakawa [1984] J. Shimada and H. Yamakawa, Macromolecules 17, 689 (1984).
  * Rosa and Everaers [2008] A. Rosa and R. Everaers, PLoS computational biology 4, 1 (2008).
  * Doi and Edwards [1988] M. Doi and S. Edwards,  _The theory of polymer dynamics_ (Oxford University Press, 1988).
  * Meakin and Ernst [1988] P. Meakin and M. H. Ernst, Phys. Rev. Lett. 60, 2503 (1988).
  * Garcia  _et al._ [1987] A. L. Garcia, C. Van Den Broeck, M. Aertsens, and R. Serneels, Physica A 143, 535 (1987).
  * Liffman [1992] K. Liffman, J. Comput. Phys. 100, 116 (1992).
  * Kruis  _et al._ [2000] F. E. Kruis, A. Maisels, and H. Fissan, AIChE J. 46, 1735 (2000).


# Effects of Efimov states on quench dynamics in a three-boson trapped system

A. D. Kerin and A. M. Martin  School of Physics, University of Melbourne,
Parkville, VIC 3010, Australia

###### Abstract

We investigate the effects of Efimov states on the post-quench dynamics of a
system of three identical bosons with contact interactions, in a spherically-
symmetric three-dimensional harmonic trap. The quench we consider is in the
s-wave contact interaction and we focus on quenches from the non-interacting
to strongly interacting regimes and vice-versa. The calculations use the
hyperspherical solutions of the three-body problem enable us to evaluate the
semi-analytical results of the Ramsey and particle separation, post quench. In
the case where the interactions are quenched from the non-interacting to
strongly interacting regime we find convergent aperiodic solutions for both
the Ramsey signal and the particle separation. In contrast for quenches from
the strongly interacting regime to the non interacting regime both the Ramsey
signal and particle separation are periodic functions. However, in this case
we find that the solutions for the particle separation diverges indicating
that in such a system large oscillations may be observable.

##  I Introduction

Efimov states are a unique type of many-body quantum state where short-range
interactions create effective long-range forces due to the exchange
interaction. They were first predicted by Efimov Efimov (1971) and first
observed by Kraemer et al. Kraemer  _et al._ (2006). Efimov states appear in
systems of as few as three bodies (Kartavtsev and Malykh, 2007, 2008; Endo
_et al._ , 2011; Petrov, 2012; Werner, 2008; Jonsell  _et al._ , 2002; Kerin
and Martin, 2022a) and while the circumstances under which Efimov states
appear is well understood the specifics of those states, e.g. their energies,
can vary depending on the specifics of the system. The influence of Efimov
states is highly relevant to a number of topics of cold-gas research including
their effects on, among other quantities, the two and three-body contacts or
three-body decay rates after a quench Fletcher  _et al._ (2017); Colussi  _et
al._ (2018); Eigen  _et al._ (2018); D’Incao  _et al._ (2018); Colussi  _et
al._ (2019); Musolino  _et al._ (2019, 2022).

In this work we consider the dynamics of a system of three identical bosons,
interacting via a contact interaction, in an isotropic harmonic trap. Such
systems can be constructed in experiment Serwane  _et al._ (2011); Murmann
_et al._ (2015); Zürn  _et al._ (2013, 2012); Stöferle  _et al._ (2006) in the
form of dilute ultracold gases. Specifically, we consider the dynamics of the
system after the a quench in the contact s-wave interactions. We consider two
quench pathways, from the non-interacting regime to the strongly interacting
(unitary) regime and vice-versa. We utilise known solutions of the static case
Cui (2012); D’Incao  _et al._ (2018); Jonsell  _et al._ (2002); Blume and
Daily (2010); Kerin and Martin (2022a); Kestner and Duan (2007) to calculate
the Ramsey signal and particle separation as functions of time following the
quench and investigate the effects of different Efimov energy spectra upon the
dynamics. Such solutions have been used to calculate thermodynamics quantities
such as virial coefficients or Tan contacts Liu  _et al._ (2009, 2010);
Stöferle  _et al._ (2006); Cui (2012); Rakshit  _et al._ (2012); Kaplan and
Sun (2011); Mulkerin  _et al._ (2012a, b); Nascimbène  _et al._ (2010); Ku
_et al._ (2012); Levinsen  _et al._ (2017); Daily and Blume (2010); Bougas
_et al._ (2021); Colussi (2019); Colussi  _et al._ (2019); Enss  _et al._
(2022). In this work we consider a three-dimensional system complementing
previous investigations into quench dynamics in two-dimensional Bougas  _et
al._ (2022) and one-dimensional Pecak  _et al._ (2016); Volosniev (2017);
Kehrberger  _et al._ (2018); Sowiński and García-March (2019) systems.

We note that the predictions in this paper are experimentally testable with
current techniques. Notably Ref. Guan  _et al._ (2019) prepared a harmonically
trapped system of two 6Li atoms, quenched in trap geometry and measured the
particle separation. It is also possible to experimentally obtain systems of
three harmonically trapped atoms Serwane  _et al._ (2011); Murmann  _et al._
(2015); Zürn  _et al._ (2013, 2012); Stöferle  _et al._ (2006) and the quench
in s-wave scattering length is possible using tools such as Feshbach resonance
Fano (1935); Feshbach (1958); Tiesinga  _et al._ (1993); Chin  _et al._
(2010). Additionally experiments measuring the Ramsey signal of trapped cold
gases after a quench have been performed Cetina  _et al._ (2016).

This paper is structured in the following way. Sec. II provides an overview of
the hyperspherical solution to the problem of three identical bosons in a
spherical harmonic trap interacting via a contact interaction, including a
review of Efimov states. In Sec. III we use the static solutions to calculate
observables of the post-quench system. We consider the non-interacting to
unitary (forwards) and vice versa (backwards) quenches. In these two cases we
calculate the Ramsey signal, the overlap of the pre- and post-quench states,
and the expectation of the particle separation. For the forwards quench we
find that both quantities can be calculated semi-analytically and in the
reverse case the Ramsey signal is still calculable but the particle separation
diverges.

##  II Overview of the Three-Body Problem

To begin, the Hamiltonian of three identical non-interacting bodies in an
isotropic three-dimensional harmonic trap is

$\displaystyle\hat{H}=\sum_{k=1}^{3}\left[\frac{-\hbar^{2}}{2m}\nabla_{k}^{2}+\frac{m\omega^{2}r_{k}^{2}}{2}\right],$
(1)

where $\vec{r}_{k}$ is the position of the $k^{\rm th}$ particle, $m$ is the
particle mass, and $\omega$ is the trapping frequency. For convenience we
define the length scales

$\displaystyle a_{\mu}=\sqrt{\frac{\hbar}{\mu\omega}},\quad
a_{M}=\sqrt{\frac{\hbar}{M\omega}},$ (2)

where $\mu=m/2$, and $M=3m$.

We use the Bethe-Peierls boundary condition to model the contact interactions
Bethe and Peierls (1935)

$\displaystyle\lim_{r_{ij}\rightarrow
0}\left[\frac{d(r_{ij}\Psi)}{dr_{ij}}\frac{1}{r_{ij}\Psi}\right]_{r_{ij}\rightarrow
0}=\frac{-1}{a_{\rm s}},$ (3)

where $\Psi$ is the total three-body wavefunction,
$r_{ij}=|\vec{r}_{i}-\vec{r}_{j}|$, and $a_{\rm s}$ is the s-wave scattering
length.

The wavefunction of three identical harmonically trapped atoms subject to Eq.
(3) is known Werner and Castin (2006a); Liu  _et al._ (2009). In particular
the hyperspherical formulation Werner and Castin (2006a) is a closed form
description of the wavefunction in the strongly interacting (unitary) and non-
interacting regimes.

We define the hyperradius $R$ and hyperangle $\alpha$

$\displaystyle R^{2}=\sqrt{r^{2}+\rho^{2}},\quad\alpha=\arctan{(r/\rho)},$ (4)

where

$\displaystyle\vec{r}$ $\displaystyle=$
$\displaystyle\vec{r}_{2}-\vec{r}_{1},$ (5) $\displaystyle\rho$
$\displaystyle=$
$\displaystyle\frac{2}{\sqrt{3}}(\vec{r}_{3}-\frac{\vec{r}_{1}+\vec{r}_{2}}{2}),$
(6)

and

$\displaystyle\vec{C}$ $\displaystyle=$
$\displaystyle\frac{\vec{r}_{1}+\vec{r}_{2}+\vec{r}_{3}}{3}$ (7)

is the centre-of-mass (COM) coordiante. The COM Hamiltonian is a simple
harmonic oscillator (SHO) Hamiltonian. The COM wavefunction is unaffected by
Eq. (3) and is a SHO wavefunction of argument $\vec{C}$ and lengthscale
$a_{M}$. In hyperspherical coordinates the relative Hamiltonian is given

$\displaystyle\hat{H}_{\rm rel}$ $\displaystyle=$
$\displaystyle\frac{-\hbar^{2}}{2\mu}\Bigg{(}\frac{\partial^{2}}{\partial
R}+\frac{1}{R^{2}\sin(\alpha)\cos(\alpha)}\frac{\partial^{2}}{\partial\alpha^{2}}\left(\cos(\alpha)\sin(\alpha)\right)$
$\displaystyle+$ $\displaystyle\frac{5}{R}\frac{\partial}{\partial
R}-\frac{4}{R^{2}}-\frac{\hat{\Lambda}_{r}^{2}}{R^{2}\sin(\alpha)}-\frac{\hat{\Lambda}_{\rho}^{2}}{R^{2}\cos(\alpha)}\Bigg{)}+\frac{\mu\omega^{2}R^{2}}{2}.$

We define an ansatz wavefunction of the form

$\displaystyle\psi_{\rm 3brel}$ $\displaystyle=$ $\displaystyle
N_{qls}\frac{F_{qs}(R)}{R^{2}}(1+\hat{P}_{13}+\hat{P}_{23})\frac{\varphi_{ls}(\alpha)}{\sin(2\alpha)}Y_{lm}(\hat{\rho}),\qquad$
(9)

where $N_{qls}$ is the normalisation constant, $F_{qs}$ is the hyperradial
wavefunction and
$\phi_{ls}=(1+\hat{P}_{13}+\hat{P}_{23})\varphi_{ls}(\alpha)Y_{lm}(\hat{\rho})/\sin(2\alpha)$
is the hyperangular wavefunction. The exchange operators $\hat{P}_{13}$ and
$\hat{P}_{23}$ exchange the positions of particles one and three and particles
two and three respectively.

Requiring the ansatz to be an eigenfunction of the Hamiltonian leads to the
hyperangular and hyperradial equations

$\displaystyle s^{2}\varphi_{ls}(\alpha)$ $\displaystyle=$
$\displaystyle-\varphi_{ls}^{\prime\prime}(\alpha)+\frac{l(l+1)}{\cos^{2}(\alpha)}\varphi_{ls}(\alpha),$
(10) $\displaystyle E_{\rm rel}F_{qs}(R)$ $\displaystyle=$
$\displaystyle\frac{-\hbar^{2}}{4\mu}\left(F_{qs}^{\prime\prime}(R)+\frac{F_{qs}^{\prime}(R)}{R}\right)$
(11) $\displaystyle+$ $\displaystyle\left(\frac{\hbar^{2}s^{2}}{4\mu
R^{2}}+\mu\omega^{2}R^{2}\right)F_{qs}(R),$

and noting that a divergence at $\alpha=\pi/2$ is non-physical gives the
condition

$\displaystyle\varphi_{ls}\left(\frac{\pi}{2}\right)$ $\displaystyle=$
$\displaystyle 0.$ (12)

Eqs. (10)-(12) determine the form of $F_{qs}(R)$ and $\varphi_{ls}(\alpha)$
Werner and Castin (2006a, b); Liu  _et al._ (2010)

$\displaystyle F_{qs}(R)$ $\displaystyle=$
$\displaystyle\begin{cases}\left(\tilde{R}\right)^{s}e^{-\tilde{R}^{2}/2}L_{q}^{s}\left(\tilde{R}^{2}\right),&s^{2}>0\\\
&\\\ \dfrac{1}{\tilde{R}}W_{\dfrac{E_{\rm
rel}}{2\hbar\omega},\dfrac{s}{2}}(\tilde{R}^{2}),&s^{2}<0\end{cases},$ (13)
$\displaystyle\varphi_{ls}(\alpha)$ $\displaystyle=$
$\displaystyle\cos^{l+1}(\alpha)$ $\displaystyle\times$
$\displaystyle{}_{2}F_{1}\left(\frac{l+1-s}{2},\frac{l+1+s}{2};l+\frac{3}{2};\cos^{2}(\alpha)\right),$

where $L_{q}^{s}$ is the associated Laguerre polynomial, $W_{E_{\rm
rel}/2\hbar\omega,s/2}$ is the Whittaker function, ${}_{2}F_{1}$ is the
hypergeometric function, $\tilde{R}=R/a_{\mu}$, $q\in\mathbb{Z}_{\geq 0}$, and
$s^{2}\in\mathbb{R}$ are the energy quantum numbers and $l\in\mathbb{Z}_{\geq
0}$ is the angular momentum quantum number. In the rest of this work we
consider only the $l=0$ case for reasons that are elucidated in the appendix.
As such we omit angular momentum indices ($l$) in subsequent notation.

In this framework the $s$-eigenvalues can only be fully specified in the non-
interacting and unitary regime. In the unitary limit applying Eq. (3) to Eq.
(9) gives the transcendental equation

$\displaystyle
0=\frac{d\varphi_{s}^{\prime}(\alpha)}{d\alpha}\Big{|}_{\alpha=0}+\frac{8}{\sqrt{3}}\varphi_{s}\left(\frac{\pi}{3}\right),$
(15)

which determines the $s$-eigenvalues, some solutions are presented below in
Table 1. In the non-interacting limit applying Eq. (3) to Eq. (9) gives $s$ as

$\displaystyle s=\begin{cases}2\\\ 2n+6\\\ \end{cases},$ (16)

for $l=0$, where $n\in\mathbb{Z}_{\geq 0}$.

$n$ | $s$
---|---
0 |  $i\cdot$ 1.006…
1 | 4.465…
2 | 6.818…
3 | 9.324…
Table 1: The three-body $s$-eigenvalues at unitarity for the 3 boson case for
$l=0$ to three decimal places.

For $s^{2}>0$, the universal case, the energy of the wavefunction is $E_{\rm
rel}=(2q+s+1)\hbar\omega$, which is implicitly determined by Eq. (11), recall
$q\in\mathbb{Z}_{\geq 0}$. For $s^{2}<0$, the Efimov case, the energy is not
uniquely determined by requiring the wavefunction, Eq. (9), be an
eigenfunction of the Hamiltonian, Eq. (LABEL:eq:HypersphericalHamiltonian), it
is instead a free parameter. Hence we require an additional condition to fix
the energy. The Efimov hyperradial wavefunction oscillates increasingly
rapidly as $R\rightarrow 0$ and we set a condition to fix the phase of the
oscillation Werner and Castin (2006a); Jonsell  _et al._ (2002)

$\displaystyle\arg\Gamma$ $\displaystyle\left[\frac{1+s-E_{\rm
rel}/\hbar\omega}{2}\right]$ (17)
$\displaystyle=-|s|\ln(\frac{R_{t}}{a_{\mu}})+\arg\Gamma(1+s)\mod\pi,$

where $R_{t}$ is the three-body parameter, an arbitrary parameter with units
of distance. $R_{t}$ determines the energies of the Efimov states. Physically
speaking $R_{t}$ is required because, in the Efimov case, Eq. (11) has an
attractive potential term proportional to $1/R^{2}$ which allows for
arbitrarily small interparticle distances. At small distances the contact
interaction assumption breaks down and the short range nature of the
interaction become significant. The Efimov energies are plotted as a function
of $R_{t}$ in Fig. 1. The energy spectrum is unbounded from below and above,
we label the states with $q\in\mathbb{Z}$, defining the $q=0$ state to be the
lowest energy state with $E_{\rm rel}>0$ at $R_{t}=\exp(\pi/|s|)a_{\mu}$. Note
that the energy of the $q=N$ state evaluated at $R_{t}=a_{\mu}$ is equal to
the energy of the $q=N-1$ state evaluated at $R_{t}=\exp(\pi/|s|)a_{\mu}$.

Figure 1: The energy spectrum for Efimov states as defined by Eq. (17).
Calculated using $s=i\cdot 1.006\dots$. The upper limit on the horizontal axis
is $R_{t}/a_{\mu}=e^{\pi/|s|}\approx 22.7$, and the vertical black line is
$R_{t}=a_{\mu}$.

##  III Quench Dynamics

In this paper we calculate the Ramsey signal and particle separation after a
quench in $a_{s}$. As part of this we need to calculate various integrals of
the hyperspherical wavefunction and the details of these are presented in the
appendix. The COM wavefunction is unaffected by Eq. (3) and so is unaffected
by a quench in $a_{\rm s}$. As such only the relative motion impacts the
system behaviour.

The time-dependent post-quench relative wavefunction is given

$\displaystyle\ket{\psi(t)}$ $\displaystyle=$ $\displaystyle e^{-i\hat{H}_{\rm
rel}t/\hbar}\ket{F_{q_{\rm i}s_{\rm i}}\phi_{s_{\rm i}}}$ (18)
$\displaystyle=$ $\displaystyle\sum_{q,s}\bra{F_{qs}\phi_{s}}\ket{F_{q_{\rm
i}s_{\rm i}}\phi_{s_{\rm i}}}e^{-iE_{qs}t/\hbar}\ket{F_{qs}\phi_{s}},\quad$

where $\hat{H}_{\rm rel}$ is the post-quench relative Hamiltonian, quantum
numbers with subscript i refer to the initial state, quantum numbers with no
subscripts are the post-quench eigenvalues and $E_{qs}$ are the post-quench
eigenenergies.

###  III.1 Ramsey signal

The Ramsey signal Cetina  _et al._ (2016) is defined as the wavefunction
overlap of the initial and final states,

$\displaystyle S(t)$ $\displaystyle=$ $\displaystyle\bra{\Psi_{\rm
i}(t)}\ket{\Psi^{\prime}(t)}=\sum_{j=0}^{\infty}|\bra{\Psi_{\rm
i}(0)}\ket{\Psi^{\prime}_{j}}|^{2}e^{-i(E^{\prime}_{j}-E_{\rm
i})t/\hbar},\qquad$ (19)

where $\Psi_{\rm i}$ is the initial state with energy $E_{\rm i}$,
$\Psi^{\prime}$ is the post-quench state and the $\Psi^{\prime}_{j}$s are the
eigenstates of the post-quench system with energy $E^{\prime}_{j}$, $j$ is
summing over all post-quench eigenstates.

As mentioned above the COM wavefunction is unaffected by the quench and
integrates to one. Hence the Ramsey signal depends only on the relative
wavefunction,

$\displaystyle S(t)$ $\displaystyle=$ $\displaystyle\sum_{q,s}|\bra{F_{q_{\rm
i}s_{\rm i}}\phi_{s_{\rm i}}}\ket{F_{qs}\phi_{s}}|^{2}e^{-i(E_{qs}-E_{q_{\rm
i}s_{\rm i}})t/\hbar}.\quad$ (20)

To evaluate the Ramsey signal we need to evaluate the hyperradial integral
$\bra{F_{qs}}\ket{F_{q_{\rm i}s_{\rm i}}}$ and the hyperangular integral
$\bra{\phi_{s}}\ket{\phi_{s_{\rm i}}}$. The appendix contains the details of
the evaluation of these integrals. With the integrals known we can then
calculate the Ramsey signal for the forwards and backwards quenches. There is
a degree of freedom in the choice of the Efimov energy spectrum, determined by
the value of $R_{t}$. Whatever the value of $R_{t}$ the normalisation is
preserved but the post-quench behaviour is nonetheless affected by the choice
of $R_{t}$.

Figure 2: Ramsey signal of the system quenched from non-interacting to
unitarity. In each panel $s_{\rm i}=2$ and the solid red line corresponds to
$q_{i}=0$, the dashed green line to $q_{i}=1$, and the dotted blue line to
$q_{i}=2$. The upper panel uses $R_{t}=a_{\mu}$ to calculate the Efimov energy
spectrum, the middle panel uses $R_{t}=5a_{\mu}$ and the lower
$R_{t}=10a_{\mu}$. These Ramsey signals are evaluated using Eq. (19) with 40
terms in each of the sums, 1600 terms total. We find that the summation is
convergent.

The Ramsey signal is the weighted sum of oscillators,
$S(t)=Ae^{-iat}+Be^{-ibt}+Ce^{-ict}+\dots$, where the weights are the square
overlaps between initial state and post-quench eigenstates and the angular
frequencies are the differences between the initial energy and post-quench
eigenenergies. The magnitude is similarly a weighted sum of oscillators but
the angular frequencies of the oscillatory terms are the differences between
post-quench eigenenergies, $(a-b),(b-c),(a-c),\dots$. The phase of the Ramsey
signal is dominated by the phase of the most heavily weighted terms.

In Fig. 2 we plot the Ramsey signal of the forwards quench for a number of
initial states and values of $R_{t}$. The calculations of the Ramsey signal
for the forwards quench are performed including only the $q\geq-1$ Efimov
energies except for the $R_{t}=a_{\mu}$ calculation which includes only the
$q\geq 0$ Efimov energies. The neglected energies are significantly lower
(e.g. for $R_{t}=a_{\mu}$ the $q=-1$ Efimov energy is
$\approx-566\hbar\omega$) and do not contribute meaningfully. Unlike with the
two-body case Kerin and Martin (2020) the magnitude of the Ramsey signal of
the forwards quench is aperiodic. This is because the post-quench
eigenenergies are irrational because the unitary $s$-eigenvalues are
irrational as are the Efimov energies in general. This means the angular
frequencies in Eq. (20) ($a$, $b$, $c\dots$ from the previous paragraph) are
irrational as are the differences between them, hence the magnitude and phase
of the Ramsey signal are aperiodic.

In the $R_{t}=a_{\mu}$, $(q_{\rm i},s_{\rm i})=(0,2)$ case (solid red line in
the upper panel of Fig. 2) the post quench states with the largest overlaps
are $(E_{q=1},s)\approx(2.27\hbar\omega,i\cdot 1.006)$, with square overlap of
$\approx 0.666$, $(E_{q=0},s)\approx(-0.85\hbar\omega,i\cdot 1.006)$, with
square overlap of $\approx 0.14$ and $(q,s)=(0,4.465\dots)$, with square
overlap $\approx 0.105$. The two largest modes in the magnitude have periods
of $\approx 2\pi/3\omega$. The phase is dominated by a period of $\approx
2.7\pi/\omega$.

In the $R_{t}=5a_{\mu}$, $(q_{\rm i},s_{\rm i})=(0,2)$ case (solid red line in
the middle panel of Fig. 2) the most significant terms are
$(E_{q=0},s)\approx(1.077\hbar\omega,i\cdot 1.006)$, with a square overlap of
$\approx 0.583$, $(E_{q=1},s)\approx(3.37\hbar\omega,i\cdot 1.006)$ with a
square overlap of $\approx 0.24$ and $(q,s)=(0,4.465\dots)$ with a square
overlap of $\approx 0.105$. The two largest modes in the magnitude have
periods of $\approx 0.87\pi/\omega$ and $\approx 0.45\pi/\omega$. The phase is
dominated by the period $\approx\pi/\omega$.

In the $R_{t}=10a_{\mu}$, $(q_{\rm i},s_{\rm i})=(0,2)$ case (solid red line
in the lower panel of Fig. 2) the terms with the largest overlaps are
$(E_{q=0},s)\approx(1.603\hbar\omega,i\cdot 1.006)$ with a square overlap of
$\approx 0.72$, $(q,s)=(0,4.465\dots)$ with a square overlap of $\approx
0.105$ and $(E_{q=1},s)\approx(3.875\hbar\omega,i\cdot 1.006)$ with a square
overlap of $\approx 0.092$. This leads to two main modes in the magnitude with
periods of $\approx 0.5\pi/\omega$ and $\approx 0.9\pi/\omega$ and the phase
is dominated by the period of $\approx 1.4\pi/\omega$.

Figure 3: Ramsey signal of the system quenched from unitarity to non-
interacting. In each panel $s_{\rm i}=i\cdot 1.006\dots$ and the solid red
line corresponds to $q_{i}=0$, the dashed green line to $q_{i}=1$, and the
dotted blue line to $q_{i}=2$. The upper panel uses $R_{t}=a_{\mu}$ to
calculate the Efimov energy spectrum, and thus the energy of the initial
state. The middle panel uses $R_{t}=5a_{\mu}$ and the lower $R_{t}=10a_{\mu}$.
These Ramsey signals are evaluated using Eq. (19) with 40 terms in each of the
sums, 1600 terms total. We find that the summation is convergent.

In Fig. 3 we plot the Ramsey signal of the backwards quench for a system
initially in an Efimov state for a variety of Efimov energies. Unlike in the
forwards quench the magnitude of the Ramsey signal of the backwards quench is
periodic. This is because the non-interacting eigenenergies are all odd
integer multiples of $\hbar\omega$. The difference between the post-quench
eigenenergies are even integers leading to the magnitude having period
$\pi/\omega$. However the phase is dominated by the largest term in Eq. (20)
and the angular frequencies of each term are irrational because the initial
state is an Efimov state which, in general, has an irrational energy. This
leads to the irregularity in the phase.

For $R_{t}=a_{\mu}$, $(E_{q_{\rm i}=0},s_{\rm i})\approx(-0.850,i\cdot 1.006)$
(solid red line of the upper panel of Fig. 3) the largest terms are the
overlaps with $(q,s)=(0,2)$ with square overlap $\approx 0.14$ and
$(q,s)=(1,2)$ with square overlap $\approx 0.12$. These terms have periods of
$\approx 0.52\pi/\omega$ and $\approx 0.34\pi/\omega$. For $R_{t}=5a_{\mu}$
$(E_{q_{\rm i}=0},s_{\rm i})\approx(1.077,i\cdot 1.006)$ (solid red line of
the middle panel of Fig. 3) the largest term is the overlaps with
$(q,s)=(0,2)$ with square overlap $\approx 0.58$ and period
$\approx\pi/\omega$. For $R_{t}=10a_{\mu}$ $(E_{q_{\rm i}=0},s_{\rm
i})\approx(1.602,i\cdot 1.006)$ (solid red line of the lower panel of Fig. 3)
the largest term is the overlaps with $(q,s)=(0,2)$ with square overlap
$\approx 0.72$ and period $\approx 1.43\pi/\omega$.

###  III.2 Particle separation

We are not limited to calculating only the Ramsey signal. It is also possible
to calculate the particle separation, $\langle\tilde{R}(t)\rangle$.

The expectation value of $\tilde{R}(t)$ is given

$\displaystyle\langle\tilde{R}(t)\rangle=\bra{\Psi^{\prime}(t)}\tilde{R}\ket{\Psi^{\prime}(t)}$
$\displaystyle=$ $\displaystyle\sum_{j,j^{\prime}}\bra{\Psi_{\rm
i}(0)}\ket{\Psi^{\prime}_{j}}\bra{\Psi^{\prime}_{j^{\prime}}}\ket{\Psi_{\rm
i}(0)}$
$\displaystyle\times\bra{\Psi^{\prime}_{j}}\tilde{R}\ket{\Psi^{\prime}_{j^{\prime}}}e^{-i(E_{j^{\prime}}-E_{j})t/\hbar},$

where $\Psi_{\rm i}$ is the initial pre-quench state with energy $E_{\rm i}$
and $\Psi^{\prime}(t)$ is the post-quench state. $\Psi^{\prime}_{j}$ and
$\Psi^{\prime}_{j^{\prime}}$ are eigenstates of the post-quench system with
eigenenergy $E_{j}$ and $E_{j^{\prime}}$ respectively, with the sums over $j$
and $j^{\prime}$ taken over all post-quench eigenstates.

The COM wavefunction is independent of the interparticle interaction and does
not impact the post-quench dynamics. Due to the hyperangular wavefunction’s
orthogonality in $s$, two sums over $s$ and $s^{\prime}$ collapse into a
single sum over $s$. Hence $\langle\tilde{R}(t)\rangle$ is given

$\displaystyle\langle\tilde{R}(t)\rangle$ $\displaystyle=$
$\displaystyle\sum_{q^{\prime},q}\sum_{s}\bra{F_{q_{\rm i}s_{\rm
i}}\phi_{s_{\rm
i}}}\ket{F_{q^{\prime}s}\phi_{s}}\bra{F_{qs}\phi_{s}}\ket{F_{q_{\rm i}s_{\rm
i}}\phi_{s_{\rm i}}}$ (22)
$\displaystyle\times\bra{F_{q^{\prime}s}\phi_{s}}\tilde{R}\ket{F_{qs}\phi_{s}}e^{-i(E_{qs}-E_{q^{\prime}s})t/\hbar},$

As in Eq. (20) indices with subscript i refer to the initial state and indices
with no subscript refer to the post-quench eigenstates. As with the Ramsey
signal all relevant integrals are presented in the appendix.

Figure 4: $\langle\tilde{R}(t)\rangle$ of a system initially in the ground
state quenched from non-interacting to unitarity. The upper panel corresponds
to $R_{t}=a_{\mu}$, the middle to $R_{t}=5a_{\mu}$ and the lower to
$R_{t}=10a_{\mu}$. Efimov states with $q\geq-1$ are included in the
calculation except $R_{t}=a_{\mu}$ where $E_{q=-1}\approx-566\hbar\omega$ does
not contribute meaningfully. The dot-dashed red line corresponds to $N_{\rm
max}=3$, the dashed green line to $N_{\rm max}=6$, the solid blue line to
$N_{\rm max}=12$ and the dotted black line to $N_{\rm max}=24$. We find that
the summation is convergent.

In Fig. 4 we plot $\langle\tilde{R}(t)\rangle$ for a system initially in the
non-interacting ground state quenched to unitarity with the upper, middle and
lower panels corresponding to $R_{t}=a_{\mu}$, $R_{t}=5a_{\mu}$,
$R_{t}=10a_{\mu}$ respectively. For $R_{t}\geq 5a_{\mu}$ we include states
with $q\geq-1$ in the calculations but for $R_{t}=a_{\mu}$
$E_{q=-1}\approx-566\hbar\omega$ and this state does not meaningfully
contribute so we only include the $q\geq 0$ Efimov states in the calculation.
In Eq. (22) terms with $q=q^{\prime}$ are constants and the $s$ contributions
to the energies cancel out, the angular frequencies depend only on $q$ and
$q^{\prime}$. The universal state terms oscillate with an angular frequency
that is an even integer multiple of $\omega$ because $q-q^{\prime}$ is an
integer but the Efimov state terms oscillate with irrational angular
frequencies because the differences between the Efimov energies, $E_{q}$ and
$E_{q^{\prime}}$, are irrational in general. For each plot we have calculated
$\langle\tilde{R}(t)\rangle$ by summing up to $N_{\rm max}=3,6,12,24$ terms in
each of the three sums in Eq. (22), and we find that the sum is convergent.

For $R_{t}=a_{\mu}$ the largest oscillating terms are
$(q^{\prime},q,s)=(q,q^{\prime},s)=(0,1,i\cdot 1.006\dots)$ with total
coefficient $\approx 0.17$ and $E_{q=1}-E_{q=0}\approx 3.12\hbar\omega$ and
$(q,q^{\prime},s)=(1,2,i\cdot 1.006\dots)$ with total coefficient $\approx
0.09$ and $E_{q=2}-E_{q=1}\approx 2.12\hbar\omega$ . This implies
characteristic periods of $\approx 2\pi/3\omega$ and $\approx\pi/\omega$. For
$R_{t}=5a_{\mu}$ the largest oscillating terms are
$(q,q^{\prime},s)=(0,1,i\cdot 1.006\dots)$ with total coefficient $\approx
0.3$ and $E_{q=1}-E_{q=0}\approx 2.3\hbar\omega$ and
$(q,q^{\prime},s)=(0,1,4.465)$ with total coefficient $\approx 0.054$ and
associated energy difference $2\hbar\omega$. This leads to characteristic
periods of $\approx 0.9\pi/\omega$ and $\pi/\omega$. For $R_{t}=10a_{\mu}$ the
largest oscillating terms are $(q,q^{\prime},s)=(0,1,i\cdot 1.006\dots)$ with
total coefficient $\approx 0.22$ and $E_{q=1}-E_{q=0}\approx 2.3\hbar\omega$
and $(q,q^{\prime},s)=(0,1,4.465)$ with total coefficient $\approx 0.054$ and
associated energy difference $2\hbar\omega$. This leads to characteristic
periods of $\approx 0.9\pi/\omega$ and $\pi/\omega$.

Figure 5: $\langle\tilde{R}(t)\rangle$ of a system following a quench from
unitarity to non-interacting. The upper, middle and lower panels correspond to
$R_{t}=a_{\mu}$, $R_{t}=5a_{\mu}$ and $R_{t}=10a_{\mu}$ and the corresponding
initial states are Efimov states with $q=0$. The dot-dashed red line
corresponds to $N_{\rm max}=3$, the dashed green line to $N_{\rm max}=6$, the
solid blue line to $N_{\rm max}=12$ and the dotted black line to $N_{\rm
max}=24$.

In Fig. 5 we plot $\langle\tilde{R}(t)\rangle$ for the backwards quench where
the system is initially in a variety of Efimov states. Unlike the forwards
quench we find that $\langle\tilde{R}(t)\rangle$ is periodic for the backwards
quench. This is because in the backwards quench the post-quench states are
universal states where the differences between eigenenergies are always even
multiples of $\hbar\omega$, leading to a period of $\pi/\omega$. However
similar to how Ref. Kerin and Martin (2020) found a divergence in
$r=|\vec{r}_{2}-\vec{r}_{1}|$ in the backwards quench we find that
$\langle\tilde{R}(t)\rangle$ also diverges for the backwards quench. In
particular we find that it is logarithmically divergent with the number of
terms in the summation, i.e. $\langle\tilde{R}(t\neq
n\pi/\omega)\rangle\propto\ln(N_{\rm max})$. This divergence is not
exclusively due to the Efimov states as the divergence is present even when
there are no Efimov states Kerin and Martin (2022b).

This divergence is quite unusual, it is not obvious why it occurs nor why it
occurs only for the reverse quench. To investigate further we look at how the
probability distribution of $R$, $P(R,t)$, evolves over time for both
quenches. $P(R,t)$ is given

$\displaystyle
P(R^{\prime},t)=\bra{\Psi^{\prime}(t)}\delta(R^{\prime}-R)\ket{\Psi^{\prime}(t)},$
(23)

where $\Psi^{\prime}(t)$ is the post-quench state.

Figure 6:  The evolution of the probability distribution of the hyperradius,
Eq. (23), for the forwards quench. The upper, middle and lower panels
correspond to $R_{t}/a_{\mu}=1,5$ and $10$ respectively. The horizontal axis
is the hyperradius and the vertical axis is time with dark blue corresponding
to low probability density, and yellow to high density. For all plots the
initial state is $(q_{\rm i},s_{\rm i})=(0,2)$ and each plot is constructed
with $N_{\rm max}=24$. In this case $P(R,t)$ is convergent with $N_{\rm max}$.

In Fig. 6 we plot the evolution of $P(R,t)$ for the forwards quench with
$R_{t}/a_{\mu}=1,5$ and $10$ in the upper, middle and lower panels
respectively. For all values of $R_{t}$ we see a qualitatively similar
evolution, the system oscillates between a broad distribution and a tightly
peaked one with a smaller mean value. This oscillation is only approximately
periodic due to the influence of the irrational Efimov energies. The broad
distribution corresponds to the initial universal state and the tightly peaked
distribution is dominated by the Efimov states, the system oscillates between
these two regimes.

To understand this oscillation it is useful to consider
$\langle\tilde{R}\rangle$ for the initial and post-quench states. For example
for $R_{t}=a_{\mu}$ the states with the largest overlaps with the initial are
$(E_{q=1},s)\approx(2.27\hbar\omega,i\cdot 1.006)$ with a square overlap of
$\approx 0.666$ and $(E_{q=0},s)\approx(-0.85\hbar\omega,i\cdot 1.006)$ with
square overlap $\approx 0.14$, these states have
$\langle\tilde{R}\rangle\approx 1.5$ and $\approx 0.57$ respectively. The
initial state is $(q,s)=(0,2)$ so we have $\langle\tilde{R}(t=0)\rangle\approx
1.66$, hence the position of the initial broad distribution is to the right of
the tightly peaked Efimov distribution. As $R_{t}$ increases
$\langle\tilde{R}\rangle$ of the strongly overlapping Efimov states increases
but they are, on average, still less than $\langle\tilde{R}\rangle$ of the
initial universal state, hence the peak of the narrow distribution moves
rightward with increasing $R_{t}$. Additionally as $R_{t}$ increases the
narrow Efimov distribution broadens because the higher energy Efimov states
are simply broader. Note in the third panel of Fig. 6 one can see small local
peaks in probability near $\tilde{R}=0$. These come from the
$(E_{q=-1},s)\approx(-5.6,i\cdot 1.006)$ state and have a square overlap with
the initial state of $\approx 0.01$ and $\langle\tilde{R}\rangle\approx 0.3$,
the $(E_{q=-1},s)$ state is also accounted for in the $R_{t}=5a_{\mu}$
calculation but the overlap is approximately 50 times smaller.

Figure 7:  The evolution of the probability distribution of the hyperradius,
Eq. (23), for the backward quench plotted at $t=0$ (dot-dashed red line),
$t=0.17\pi/\omega$ (dashed green), $t=0.34\pi/\omega$ (solid blue) and
$t=\pi/\omega$ (dotted black). The initial states in the upper, middle and
lower panels are the $q=0$ Efimov states for $R_{t}/a_{\mu}=1,5$ and $10$
respectively. All calculations are performed with $N_{\rm max}=60$, unlike in
the forwards quench we find that $P(R,t)$ is only convergent for $t=0$.
Figure 8:  The tail of $P(R,t=\pi/2\omega)$ for the reverse quench with the
energy of the initial Efimov state given $E\approx 1.077$ ($q=0$ for
$R_{t}=5a_{\mu}$) for various values of $N_{\rm max}$. The dot-dashed red line
corresponds to $N_{\rm max}=10$, the dashed green line to $N_{\rm max}=20$,
the solid blue line to $N_{\rm max}=30$ and the dotted black line to $N_{\rm
max}=40$.

In Fig, 7 we plot $P(R,t)$ for the backwards quench at
$t=0,0.17\pi/\omega,0.34\pi/\omega$ and $\pi/2\omega$. Unlike the forwards
quench the evolution of $P(R,t)$ here is periodic. The mean of $P(R,t)$
increases with time, reaching a maximum at $\pi/2\omega$ before returning to
its initial shape at $t=\pi/\omega$, this then repeats with period
$\pi/\omega$. Initially $P(R,t)$ is tightly peaked but develops a long tail as
it evolves and in Fig. 8 we present the long tail in detail for various values
of $N_{\rm max}$. The tail behaves approximately like $\tilde{R}^{-2}$ until
ending in an exponential-like “cut-off”. This cut-off occurs at larger $R$ for
larger $N_{\rm max}$ and in the limit of $N_{\rm max}\rightarrow\infty$ the
behaviour of the tail of $P(R,t)$ approaches $\tilde{R}^{-2}$ with no cut-off.
This means that the integral of $P(R,t)$ over $R$ from $R=0$ to
$R\rightarrow\infty$ is finite and properly normalised in the $N_{\rm
max}\rightarrow\infty$ limit, however $RP(R,t)$ has a $\tilde{R}^{-1}$ tail
and so is the integral is not finite, hence the divergence in
$\langle\tilde{R}(t)\rangle$ for the backwards quench.

Physically speaking there are two likely candidates for the source of the
divergence; the zero-range contact interaction and the instantaneous nature of
the quench, in reality atoms interact at some finite range and the quench in
$a_{\rm s}$ occurs over some finite time. These are two non-physical inputs
into this model and may be responsible for the non-physical outputs.

By considering the finite range of the interaction it is possible to estimate
a maximum value of $\langle\tilde{R}(t)\rangle$. The lengthscale of the
interaction provides a justification for a maximum energy and thus a cut-off
in Eq. (22). Specifically, the range of interaction defines a minimum de
Broglie wavelength which defines a maximum energy and thus the cut-off. For
sodium in a 1kHz trap and assuming a van der Waals range of one nanometre we
obtain an energy of $E_{\rm rel}\approx 8.7\times 10^{6}\hbar\omega$ and so we
predict $\langle\tilde{R}(t)\rangle_{\rm max}\approx 21$ for an initial Efimov
energy of $E_{q=0}\approx-0.85$. This is an order of magnitude larger than the
amplitude of oscillations when the system is quenched from the non-interacting
to the strongly interacting regime.

In contrast it is difficult to quantify the effects of a finite duration
quench. In the formalism used here only quenches between the non-interacting
and unitary regimes can be described, meaning a quench to or from the
intermediate regime can’t be elucidated. However in the two-body case a quench
between any two scattering lengths can be considered Kerin and Martin (2020)
so the effects of a finite duration quench can be investigated.

##  IV Conclusion

In this paper we have examined the effects of different Efimov energy spectra
on the time dependent post-quench dynamics of an interacting few-body system.
This was done in the context of three interacting bosons in a spherically
symmetric trap, where the contact interactions were quenched from the non-
interacting regime to the strongly-interacting regime (forwards quench) and
vice-versa (backwards quench). In each case we were able to evaluate the post
quench dynamics of both the Ramsey signal and the expectation value of the
hyperradius.

For the forwards quench we find an irregularly repeating signal for both the
Ramsey signal and $\langle\tilde{R}(t)\rangle$. For the Ramsey signal this is
due to both the Efimov energies and unitary $s$-eigenspectrum being irrational
in general. In the case of the particle separation the contributions from $s$
cancel out and the irregularity is due to the irrationality of the Efimov
energies. In both cases the results are convergent and well defined.

For the backwards quench the magnitude of the Ramsey signal and
$\langle\tilde{R}(t)\rangle$ oscillate with period $\pi/\omega$. This is
because the non-interacting $s$-eigenvalues are even integers and Efimov
states are not present when $a_{\rm s}=0$. The phase of the Ramsey signal is
still irregular due to the influence of the initial irrational Efimov energy.
However we find, analogous to previous results Kerin and Martin (2020), that
the particle separation diverges logarithmically. By enforcing a cut-off on
Eq. (22) motivated by a minimum de Broglie wavelngth derived from the van der
Waals range we expect a maximum $\langle\tilde{R}(t)\rangle\approx 21$. This
estimate of the size of the oscillations is extremely large compared to the
forwards quench case.

##  V Acknowledgements

A.D.K. is supported by an Australian Government Research Training Program
Scholarship and by the University of Melbourne.

With thanks to Victor Colussi for illuminating discussions regarding the
evaluation of the hyperangular integral.

## Appendix

In this work we calculate quench observables of systems where the
wavefunctions of the pre- and post-quench systems are known. To obtain these
observables we need to perform numerous integrals involving these
wavefunctions. In this appendix we present those integrals. Firstly, the
Jacobian in hyperspherical coordinates is given by

$\displaystyle
dV=d\vec{r}_{1}d\vec{r}_{2}d\vec{r}_{3}=\frac{3\sqrt{3}}{32}R^{5}\sin^{2}(2\alpha)dRd\alpha
d\vec{\Omega}_{r}d\vec{\Omega}_{\rho}d\vec{C},$ (24)

and for convenience we define

$\displaystyle\bra{F_{qs}(R)}\ket{F_{q^{\prime}s^{\prime}}(R)}$
$\displaystyle=$
$\displaystyle\int_{0}^{\infty}RF_{qs}(R)^{*}F_{q^{\prime}s^{\prime}}(R)dR,$
(25) $\displaystyle\bra{\phi_{s}(\alpha)}\ket{\phi_{s^{\prime}}(\alpha)}$
$\displaystyle=$
$\displaystyle\int\int\int_{0}^{\pi/2}\phi_{s}(\alpha)^{*}\phi_{s^{\prime}}(\alpha)2\sin^{2}(2\alpha)d\alpha
d\vec{\Omega}_{r}d\vec{\Omega}_{\rho}.$ (26)

To calculate the Ramsey signal, Eq. (20), we need the wavefunction overlaps,
i.e. $\bra{F_{qs}(R)}\ket{F_{q^{\prime}s^{\prime}}(R)}$ and
$\bra{\phi_{s}(\alpha)}\ket{\phi_{s^{\prime}}(\alpha)}$. Whether $s$ is
imaginary or not does not change the functional form of the hyperangular
wavefunction, $\phi_{s}(\alpha)$, unlike the hyperradial wavefunction,
$F_{qs}(R)$. For the hyperangular integral there is only one case, but for the
hyperradial integral there are three; the universal-universal, universal-
Efimov and Efiomv-Efimov. We begin by considering the hyperangular integral.

The presence of the permutation operators makes evaluating the hyperangular
integral directly difficult. To evaluate we transform the permuted terms into
the same Jacobi set as the unpermuted term Nielsen  _et al._ (2001). However
this limits us to the $l=0$ case, if the spherical harmonic term is non-
constant then the coordinate transform is more complicated and the integral
becomes intractable. The hyperangular integral is given Fedorov and Jensen
(1993, 2001); Braaten and Hammer (2006); Thøgersen (2009)

$\displaystyle\bra{\phi_{0s}}\ket{\phi_{0s_{\rm i}}}$ $\displaystyle=$
$\displaystyle
8\pi\int_{0}^{\pi/2}\left(\left(1+\hat{P}_{23}+\hat{P}_{13}\right)\frac{\varphi_{s}(\alpha)}{\sin(2\alpha)}\right)^{*}\left(\left(1+\hat{P}_{23}+\hat{P}_{13}\right)\frac{\varphi_{s_{\rm
i}}(\alpha)}{\sin(2\alpha)}\right)\sin^{2}(2\alpha)d\alpha,$ (27)
$\displaystyle=$ $\displaystyle
24\pi\Bigg{[}\int_{0}^{\pi/2}\varphi_{s}^{*}(\alpha)\varphi_{s_{\rm
i}}(\alpha)d\alpha+\frac{4}{\sqrt{3}}\int_{0}^{\pi/2}\varphi_{s}^{*}(\alpha)\left[\int_{|\pi/3-\alpha|}^{\pi/2-|\pi/6-\alpha|}\varphi_{s_{\rm
i}}(\alpha^{\prime})d\alpha^{\prime}\right]d\alpha\Bigg{]}.$

For $l=0$ we have Werner (2008); Fedorov and Jensen (2001)

$\displaystyle\varphi_{0s}\propto\sin\left(s\left(\frac{\pi}{2}-\alpha\right)\right).$
(28)

Note that evaluating Eq. (27) with Eq. (28) does not give the same result as
Ref. Werner (2008) in general. This is because Ref. Werner (2008) is firstly
concerned with the overlaps with the ground state and secondly combine Eq.
(27) with Eq. (15) and so the results presented here and in the latter
reference agree when $s$ is a unitary eigenvalue and $(q_{\rm i},s_{\rm
i})=(0,2)$. Note that different hyperangular states of the same regime (i.e.
two different unitary values of $s$ or two different non-interacting values of
$s$) are orthogonal, but there is non-zero overlap between unitary and non-
interacting states.

The hyperradial integrals are integrals of products of well understood
functions. The universal-universal integral is given Srivastava  _et al._
(2003)

$\displaystyle\bra{F_{qs}(R)}\ket{F_{q_{\rm i}s_{\rm i}}(R)}=$
$\displaystyle\frac{a_{\mu}^{2}}{2}\binom{q+s}{q}\binom{q_{\rm
i}+\dfrac{s_{\rm i}-s}{2}-1}{q_{\rm i}}\Gamma\left(\frac{s+s_{\rm
i}}{2}+1\right){}_{3}F_{2}\bigg{(}-q,\frac{s+s_{\rm i}}{2}+1,\frac{s-s_{\rm
i}}{2}+1;s+1,\frac{s-s_{\rm i}}{2}-q_{\rm i}+1;1\bigg{)},$ (29)

the Efimov-Efimov Gradshteyn and Ryzhik (2014)

$\displaystyle\bra{F_{qs}(R)}\ket{F_{q_{\rm i}s}(R)}$ $\displaystyle=$
$\displaystyle
a_{\mu}^{2}\real[\frac{\Gamma\left(s+1\right)\Gamma(-s)}{\Gamma\left(\dfrac{1-E_{q_{\rm
i}}/\hbar\omega-s}{2}\right)\Gamma\left(\dfrac{3-E_{q}/\hbar\omega+s}{2}\right)}$
(30) $\displaystyle\times{}_{3}F_{2}\bigg{(}s+1,1,\frac{1-E_{q_{\rm
i}}/\hbar\omega+s}{2};1+s,\frac{3-E_{q}/\hbar\omega+s}{2};1\bigg{)}\Bigg{]},$

and the universal-Efimov Gradshteyn and Ryzhik (2014)

$\displaystyle\bra{F_{qs}(R)}\ket{F_{q_{\rm i}s_{\rm
i}}(R)}=\frac{a_{\mu}^{2}}{4}\frac{(-1)^{q_{\rm i}}}{\Gamma(1+q_{\rm
i})}\Bigg{[}$ $\displaystyle\frac{\Gamma\left(\dfrac{2-s^{*}+s_{\rm
i}}{2}\right)\Gamma\left(\dfrac{2+s^{*}+s_{\rm i}}{2}\right)\Gamma(-s_{\rm
i})}{\Gamma(-q_{\rm i}-s_{\rm i})\Gamma\left(\dfrac{3-E_{q}/\hbar\omega+s_{\rm
i}}{2}\right)}{}_{3}F_{2}\bigg{(}1+\frac{s_{\rm
i}-s^{*}}{2},1+\frac{s^{*}+s_{\rm i}}{2},-q_{\rm i};1+s_{\rm
i},\frac{3-E_{q}/\hbar\omega+s_{\rm i}}{2};1\bigg{)}$
$\displaystyle+\frac{\Gamma\left(\dfrac{2+s^{*}-s_{\rm
i}}{2}\right)\Gamma\left(\dfrac{2-s^{*}-s_{\rm i}}{2}\right)\Gamma(s_{\rm
i})}{\Gamma(-q_{\rm i})\Gamma\left(\dfrac{3-E_{q}/\hbar\omega-s_{\rm
i}}{2}\right)}{}_{3}F_{2}\bigg{(}1+\frac{s^{*}-s_{\rm
i}}{2},1-\frac{s^{*}+s_{\rm i}}{2},-q_{\rm i}-s_{\rm i};1-s_{\rm
i},\frac{3-E_{q}/\hbar\omega-s_{\rm i}}{2};1\bigg{)}\Bigg{]},$ (31)

where we have used the identity

$\displaystyle
L_{n}^{\alpha}(z)=\frac{(-1)^{n}}{n!}e^{z/2}z^{-(n+1)/2}W_{\dfrac{2n+\alpha+1}{2},\dfrac{\alpha}{2}}(z).$
(32)

For the hyperradial integral we have that for $s=s^{\prime}$ the integral is 0
for $q\neq q^{\prime}$.

To calculate the particle separation expectation value, Eq. (22), we again
need to calculate a number of integrals involving the wavefunction. All the
needed integrals except $\bra{F_{q^{\prime}s}}\tilde{R}\ket{F_{qs}}$ are given
above. Previously for the hyperradial integral we had three cases, here we do
not need to consider the universal-Efimov case as $s$ is the same in both bra
and ket due to the orthogonality in $s$ of the hyperangular integral. For the
universal-universal case we have Srivastava  _et al._ (2003)

$\displaystyle\bra{F_{qs}(R)}\tilde{R}\ket{F_{q^{\prime}s}(R)}=\frac{a_{\mu}^{2}}{2}\binom{q+s}{q}\binom{q^{\prime}-\frac{3}{2}}{q^{\prime}}\Gamma\left(s+\frac{3}{2}\right){}_{3}F_{2}\left(-q,s+\frac{3}{2},\frac{3}{2};s+1,\frac{3}{2}-q^{\prime};1\right),$
(33)

and for the Efimov-Efimov case we have Gradshteyn and Ryzhik (2014)

$\displaystyle\bra{F_{qs}(R)}\tilde{R}\ket{F_{q^{\prime}s}(R)}$
$\displaystyle=$ $\displaystyle
a_{\mu}^{2}\real[\frac{\Gamma\left(\dfrac{3}{2}+s\right)\Gamma\left(\dfrac{3}{2}\right)\Gamma(-s)}{\Gamma\left(\dfrac{1-E_{q^{\prime}}/\hbar\omega-s}{2}\right)\Gamma\left(\dfrac{4-E_{q}/\hbar\omega+s}{2}\right)}$
(34)
$\displaystyle\times{}_{3}F_{2}\bigg{(}\frac{3}{2}+s,\frac{3}{2},\frac{1-E_{q^{\prime}}/\hbar\omega+s}{2};1+s,\frac{4-E_{q}/\hbar\omega+s}{2};1\bigg{)}\Bigg{]}.$

## References

  * Efimov (1971) V. Efimov, Soviet Journal of Nuclear Physics 12, 589 (1971).
  * 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, _et al._ , Nature 440, 315 (2006).
  * Kartavtsev and Malykh (2007) O. Kartavtsev and A. Malykh, Journal of Physics B: Atomic, Molecular and Optical Physics 40, 1429 (2007).
  * Kartavtsev and Malykh (2008) O. Kartavtsev and A. Malykh, JETP letters 86, 625 (2008).
  * Endo  _et al._ (2011) S. Endo, P. Naidon, and M. Ueda, Few-Body Systems 51, 207 (2011).
  * Petrov (2012) D. S. Petrov, Many-Body Physics With Ultracold Gases (Les Houches 2010) Lecture Notes of the Les Houches Summer School 94, 109 (2012).
  * Werner (2008) F. Werner, Theses, Université Pierre et Marie Curieâ Paris VI, Paris France (2008).
  * Jonsell  _et al._ (2002) S. Jonsell, H. Heiselberg, and C. Pethick, Physical Review Letters 89, 250401 (2002).
  * Kerin and Martin (2022a) A. Kerin and A. Martin, arXiv preprint arXiv:2204.09205 (2022a).
  * Fletcher  _et al._ (2017) R. J. Fletcher, R. Lopes, J. Man, N. Navon, R. P. Smith, M. W. Zwierlein, and Z. Hadzibabic, Science 355, 377 (2017).
  * Colussi  _et al._ (2018) V. Colussi, J. P. Corson, and J. P. D’Incao, Physical review letters 120, 100401 (2018).
  * Eigen  _et al._ (2018) C. Eigen, J. A. Glidden, R. Lopes, E. A. Cornell, R. P. Smith, and Z. Hadzibabic, Nature 563, 221 (2018).
  * D’Incao  _et al._ (2018) J. P. D’Incao, J. Wang, and V. E. Colussi, Physical Review Letters 121, 023401 (2018).
  * Colussi  _et al._ (2019) V. Colussi, B. van Zwol, J. D’Incao, and S. Kokkelmans, Physical Review A 99, 043604 (2019).
  * Musolino  _et al._ (2019) S. Musolino, V. Colussi, and S. Kokkelmans, Physical Review A 100, 013612 (2019).
  * Musolino  _et al._ (2022) S. Musolino, H. Kurkjian, M. Van Regemortel, M. Wouters, S. Kokkelmans, and V. Colussi, Physical Review Letters 128, 020401 (2022).
  * Serwane  _et al._ (2011) F. Serwane, G. Zürn, T. Lompe, T. Ottenstein, A. Wenz, and S. Jochim, Science 332, 336 (2011).
  * Murmann  _et al._ (2015) S. Murmann, A. Bergschneider, V. M. Klinkhamer, G. Zürn, T. Lompe, and S. Jochim, Physical Review Letters 114, 080402 (2015).
  * Zürn  _et al._ (2013) G. Zürn, A. N. Wenz, S. Murmann, A. Bergschneider, T. Lompe, and S. Jochim, Physical Review Letters 111, 175302 (2013).
  * Zürn  _et al._ (2012) G. Zürn, F. Serwane, T. Lompe, A. N. Wenz, M. G. Ries, J. E. Bohn, and S. Jochim, Physical Review Letters 108, 075303 (2012).
  * Stöferle  _et al._ (2006) T. Stöferle, H. Moritz, K. Günter, M. Köhl, and T. Esslinger, Physical Review Letters 96, 030401 (2006).
  * Cui (2012) X. Cui, Few Body Systems 52, 65 (2012).
  * Blume and Daily (2010) D. Blume and K. M. Daily, Physical review letters 105, 170403 (2010).
  * Kestner and Duan (2007) J. Kestner and L.-M. Duan, Physical Review A 76, 033611 (2007).
  * Liu  _et al._ (2009) X.-J. Liu, H. Hu, and P. D. Drummond, Physical Review Letters 102, 160401 (2009).
  * Liu  _et al._ (2010) X.-J. Liu, H. Hu, and P. D. Drummond, Physical Review A 82, 023619 (2010).
  * Rakshit  _et al._ (2012) D. Rakshit, K. M. Daily, and D. Blume, Physical Review A 85, 033634 (2012).
  * Kaplan and Sun (2011) D. B. Kaplan and S. Sun, Physical Review Letters 107, 030601 (2011).
  * Mulkerin  _et al._ (2012a) B. Mulkerin, C. Bradly, H. Quiney, and A. Martin, Physical Review A 85, 053636 (2012a).
  * Mulkerin  _et al._ (2012b) B. C. Mulkerin, C. J. Bradly, H. M. Quiney, and A. M. Martin, Physical Review A 86, 053631 (2012b).
  * Nascimbène  _et al._ (2010) S. Nascimbène, N. Navon, F. Jiang, K. Chevy, and C. Salomon, Nature 463, 1057 (2010).
  * Ku  _et al._ (2012) M. Ku, A. Sommer, L. Cheuk, and M. Zwierlein, Science 335, 563 (2012).
  * Levinsen  _et al._ (2017) J. Levinsen, P. Massignan, S. Endo, and M. M. Parish, Journal of Physics B: Atomic, Molecular and Optical Physics 50, 072001 (2017).
  * Daily and Blume (2010) K. Daily and D. Blume, Physical Review A 81, 053615 (2010).
  * Bougas  _et al._ (2021) G. Bougas, S. Mistakidis, P. Giannakeas, and P. Schmelcher, New Journal of Physics 23, 093022 (2021).
  * Colussi (2019) V. E. Colussi, Atoms 7, 19 (2019).
  * Enss  _et al._ (2022) T. Enss, N. C. Braatz, and G. Gori, arXiv preprint arXiv:2203.06098 (2022).
  * Bougas  _et al._ (2022) G. Bougas, S. Mistakidis, P. Giannakeas, and P. Schmelcher, arXiv preprint arXiv:2205.02015 (2022).
  * Pecak  _et al._ (2016) D. Pecak, M. Gajda, and T. Sowiński, New Journal of Physics 18, 013030 (2016).
  * Volosniev (2017) A. G. Volosniev, Few-Body Systems 58, 1 (2017).
  * Kehrberger  _et al._ (2018) L. M. A. Kehrberger, V. J. Bolsinger, and P. Schmelcher, Physical Review A 97, 013606 (2018).
  * Sowiński and García-March (2019) T. Sowiński and M. Á. García-March, Reports on Progress in Physics 82, 104401 (2019).
  * Guan  _et al._ (2019) Q. Guan, V. Klinkhamer, R. Klemt, J. H. Becher, A. Bergschneider, P. M. Preiss, S. Jochim, and D. Blume, Physical Review Letters 122, 083401 (2019).
  * Fano (1935) U. Fano, Il Nuovo Cimento 12, 154 (1935).
  * Feshbach (1958) H. Feshbach, Annals of Physics 5, 357 (1958).
  * Tiesinga  _et al._ (1993) E. Tiesinga, B. J. Verhaar, and H. T. C. Stoof, Physical Review A 47, 4114 (1993).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Reviews of Modern Physics 82, 1225 (2010).
  * Cetina  _et al._ (2016) M. Cetina, M. Jag, R. S. Lous, I. Fritsche, J. T. Walraven, R. Grimm, J. Levinsen, M. M. Parish, R. Schmidt, M. Knap, _et al._ , Science 354, 96 (2016).
  * Bethe and Peierls (1935) H. Bethe and R. Peierls, Proceedings of the Royal Society of London. Series A-Mathematical and Physical Sciences 148, 146 (1935).
  * Werner and Castin (2006a) F. Werner and Y. Castin, Physical Review Letters 97, 150401 (2006a).
  * Werner and Castin (2006b) F. Werner and Y. Castin, Physical Review A 74, 053604 (2006b).
  * Kerin and Martin (2020) A. D. Kerin and A. M. Martin, Physical Review A 102, 023311 (2020).
  * Kerin and Martin (2022b) A. Kerin and A. Martin, arXiv preprint arXiv:2207.09091v2 (2022b).
  * Nielsen  _et al._ (2001) E. Nielsen, D. V. Fedorov, A. S. Jensen, and E. Garrido, Physics Reports 347, 373 (2001).
  * Fedorov and Jensen (1993) D. Fedorov and A. Jensen, Physical review letters 71, 4103 (1993).
  * Fedorov and Jensen (2001) D. V. Fedorov and A. Jensen, Journal of Physics A: Mathematical and General 34, 6003 (2001).
  * Braaten and Hammer (2006) E. Braaten and H.-W. Hammer, Physics Reports 428, 259 (2006).
  * Thøgersen (2009) M. Thøgersen, arXiv preprint arXiv:0908.0852 (2009).
  * Srivastava  _et al._ (2003) H. M. Srivastava, H. A. Mavromatis, and R. S. Alassar, Applied Mathematics Letters 16, 1131 (2003).
  * Gradshteyn and Ryzhik (2014) I. S. Gradshteyn and I. M. Ryzhik,  _Table of integrals, series, and products_ (Academic press, 2014).