\section{Introduction}
Since the past few decades there have been significant ongoing efforts directed to understanding the rich properties of vortex matter in type-2 superconductors \cite{blatter:94}.
The glass formation has been intensively studied for vortex systems with pinning disorder \cite{VG1,VG2,blatter:94,VG4}.
One major research theme is investigating the effects of pinning disorder on vortex glass formation, and the corresponding depinning transition and dynamics upon driving. This is important in superconductor applications, as dissipationless current-carrying states require arresting the dynamics of vortices.
In the much less studied multicomponent superconductors,
the vortex interactions
can take complex forms as a result of multiple competing
(attractive and repulsive) interaction length scales \cite{Babaev.Speight:05,Babaev.Carlstrom.ea:10,carlstrom2011type,Silaev.Babaev:11,babaev2017type}.
This has been proposed not only in
superconductors with several superconducting components but also in multilayer systems with different coherence ($\xi_i$) and penetration ($\lambda_i$) lengths \cite{carlstrom2011semi,varney2013hierarchical,Rogelio:Glass}.
In the bulk of intrinsic type-$1.5$ superconductors, intervortex forces are long-range attractive and short-range repulsive \cite{Babaev.Speight:05,Babaev.Carlstrom.ea:10,Silaev.Babaev:11,babaev2017type}. Various multiband materials were proposed to fall into the type-1.5 class \cite{RB1,ray2014muon,PhysRevB.102.144523}.
The double-transition materials are good candidates for type-1.5 behaviour due to
the diverging coherence length at the lower transition  \cite{carlstrom2011length,garaud2018properties}.
For thin films, an additional repulsive interaction arises due to the magnetic stray fields.
Many vortex phases have been predicted for vortices in type-$1.5$ superconducting systems \cite{meng16,meng2014honeycomb,varney2013hierarchical,Rogelio:Glass}. Particularly, the cluster phases are in many ways similar to ones
emerging in a wide range of cluster-forming monodisperse particle systems such as colloidal suspensions and ultracold atoms \cite{MPstripe,cinti14,pc,RR:stripes,RR:phases,RR:depinning,ca:patterns,CCSH,Glaser:SS}.

In type-2 non-cluster-forming vortex systems, glass phases are associated with pinning, while
in contrast, with cluster-forming interactions the situation is different and
quench dynamics can produce a glass also without pinning \cite{Rogelio:Glass}. The effect of pinning is much less understood for type-1.5 superconductors. The especially interesting  question is the effect of pinning on vortex cluster glass formation. In addition, the nonequlibrium depinning transition of monodisperse particles has been intensively studied for particle solids; see \cite{RR:depinning} for a recent review.
However, the depinning of pattern forming systems, in particular the cluster phases \cite{RR:hysteresis} which are experimentally relevant, has not been extensively studied. Indeed, even the nature of the depinning and the main features of the flowing phase have not been explored in detail.

The main purpose of this work is to present a systematic study of cluster glass formation following a thermal quench in presence of pinning disorder, and the corresponding depinning dynamics when a constant external force is applied. We are interested in, e.g., how pinning disorder density and depth affect glass formation, and in fact whether cluster glass itself remains well defined, and for depinning we are interested in the nature of the depinning transition, i.e., elastic depinning or plastic depinning, and the main features of the flowing phase. Here, plastic flow refers to that the flowing clusters change their neighbours, while in elastic flow, they do not.

For traps of a moderate depth, we find that cluster glass remains reasonably well defined, and persists even to the regime of dense pinning.
The depinning transition is found to be plastic, and the flowing phase has a number of interesting features compared with the non-cluster-forming counterpart.
In addition to the cluster plastic flow, the clusters can undergo structural transformations.
Most notably, the clusters can merge or break up due to thermal effects,
and, more importantly, the interplay between the pinning force and the driving force.
The effective polydispersity assists the plastic flow also in that a smaller cluster tends to have a better mobility, and consequently diffuses and escapes from a trap faster. Finally, we have compared the depinning of a cluster glass and a particle glass, finding that the cluster glass exhibits stronger pinning effects. This suggests that cluster vortex glasses might be suitable for certain technological applications where a high critical current is desired.

The paper is organized as follows: the model, observables, and numerical methods are described in Sec.~\ref{mm}.
Numerical results are given in Sec.~\ref{results}.
Finally, a summary and discussion of the main findings is presented in Sec.~\ref{cc}.

\section{Model, observables, and methods}
\label{mm}

\subsection{Model}

We study a system of two-dimensional monodisperse particles interacting via an effective pairwise potential,
and with pinning centers, and an external driving force modelled by the following dimensionless Hamiltonian:
\begin{eqnarray}
H = \sum_{i<j} U_{ij} - \sum_{i\alpha} U_p \exp[-r_{i\alpha}^2/(2\sigma^2)] - F_D\sum_i x_i.
\end{eqnarray}
The indices $i, j$ run from 1 to the number of particles $N$, and
$\alpha$ runs from 1 to the number of pinning centers $M$.
The three terms are for particle-particle interactions, randomly located and quenched Gaussian
pinning centers,
and
a constant driving force along the $x$-direction, respectively.
The particles are placed on an $L \times L$ square with periodic boundary conditions. The density of particles is $n=N/L^2$ and the density of impurities is $m=M/L^2=nM/N$.
$U_p>0$ and $\sigma$ denote the depth and size of the pinning centers, respectively.
$F_D$ denotes the strength of the driving force. We
used $\sigma=0.7$
unless specified otherwise, which is approximately half of the lattice constant of the corresponding cluster crystal
in the absence of disorder and driving.
A point-particle representation is applicable for superconducting vortices where
the interparticle interaction forces are chosen to be consistent with those
obtained in superconductivity models.

The potential, relevant to vortex interactions in
type-1.5 superconductors, is defined as:
\begin{eqnarray}
        U(r) =
        \begin{cases}
        U_0 - \alpha_0(r+0.15)^4, & \text{if } r\leq0.15, \\
        \sum_{l=1}^{3} C_lK_0(\alpha_l r), & \text{if } r>0.15,
        \end{cases}
    \label{potential}
\end{eqnarray}
where $K_0$ is the modified Bessel function of the second kind,
and the parameters are
$U_0=2.6796, \alpha_0=144.3760, C_1=7.2900, \alpha_1=2.3810, C_2=37.2100, \alpha_2=7.1174, C_3=-36.1911, \alpha_3=5.1020$.
The potential is plotted as the red curve in Fig.~\ref{UR}, it is similar to the one used in \cite{Rogelio:Glass} in medium and long ranges, but it features an attractive intermediate part and the unphysical divergence from long-range asymptotic analysis at the origin is regularized with a fourth-order polynomial to match better with the full numerically calculated potential \cite{carlstrom2011semi}.

For comparison, we also use the standard repulsive type-2 vortex
interaction potential from a single modified Bessel function $U(r)=CK_0(\alpha r)$ with parameters $C=1.7511$ and $\alpha=1.4867$,
shown as the blue curve in Fig.~\ref{UR}.
This potential generates single particle crystals with no clustering behavior.

The physical origin of the potential (\ref{potential}) is the following.
In the type-1.5 regime a superconductor has several coherence lengths,
originating from the individual superconducting components. Some of the coherence lengths are
larger and some are smaller than the magnetic field penetration length.
The resulting intervortex forces asymptotically have the form
of a sum of Bessel functions
which give contributions to attractive interactions
at the coherence length scales
and repulsive interaction at the magnetic field
penetration length \cite{Babaev.Speight:05,Babaev.Carlstrom.ea:10,babaev2017type,carlstrom2011type,garaud2018properties}.
In the standard type-1.5 regime the intervortex interaction is short-range
repulsive and long-range attractive.
However, the situations can also occur where the interaction has a local minimum at an intermediate length scale.
This is expected to be realized in
thin films of type-1.5 superconductors where
at the long range the interaction should be repulsive due to stray fields \cite{Rogelio:Glass},
and in layered systems where there are different
magnetic field penetration lengths in different layers \cite{varney2013hierarchical,meng16,Rogelio:Glass}. Multiple repulsive length scales can also arise as an intrinsic feature in certain anisotropic models \cite{winyard2019hierarchies}.

Since we are interested
in not very dense ensembles the only retained attractive Bessel function
corresponds to the interactions of larger cores.
We consider a bilayer system with different penetration lengths that
is approximated by two repulsive Bessel functions. Note that the
potential with such an intermediate minimum also arises
in a thin film of type-1.5 superconductors where instead of the Bessel function with the largest length one includes a power-law repulsion \cite{Rogelio:Glass}.
We emphasize that this does not affect the existence of cluster crystals, which appear at moderately low densities.
Such a regime is the main focus of this work.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{UR.pdf}
\caption{
Pairwise particle-particle, i.e., vortex-vortex interaction potentials.
The red curve is a cluster-forming potential [Eq.~(\ref{potential})] while the blue curve forms single-particle crystals.
They model vortex interactions
in type-1.5 and type-2 superconductors, respectively.
}
\label{UR}
\end{center}
\end{figure}

We briefly compare our potential to other commonly used cluster-forming potentials.
In addition to \cite{Rogelio:Glass}, our potential is also similar to ultrasoft
potentials such as $U(r)=1/(1+r^6)$ \cite{Wang:VC} in the intermediate and long ranges, but differs near the origin.
Therefore the resulting phase diagrams differ at high densities, but all three potentials lead to cluster crystal phases at low densities.
Our potential has a characteristic
medium attractive range from the attractive Bessel term.
This feature is similar to the potential studied in \cite{RR:phases}, which is a sum of a repulsive Coulomb term and an attractive exponential term. The two potentials again differ near the origin and therefore the high-density regimes will be different.

\subsection{Observables and methods}

Both equilibrium and dynamical properties of the system were studied for a large range of particle as well as trap densities
using Monte Carlo (MC) and molecular dynamics (MD) simulations.
The simulations have three different parts.
First, equilibrium sampling was used to map the phase diagram of the clean system with no pinning and no driving, i.e., $U_p = F_D =0$.
The transition temperatures as well as the ordered low-temperature phases
were obtained using parallel tempering MC \cite{ptmc1,ptmc2,Hukushima:PT}.
The calculations were complemented by simulated annealing MC \cite{SA} as a consistency check.
While simulated annealing cannot faithfully maintain thermal equilibrium at low temperatures,
it can nevertheless capture the essential features of the phases, e.g.,
produce typical low energy states \cite{Wang:VC,ca:patterns}.
We used the Metropolis algorithm with a sequential update order of the particles.

For the glass formation simulations with pinning but no driving,
i.e., $U_p \neq 0$ and $F_D=0$, we have used single temperature MC simulations, and
a random update order of the particles
starting from fully disordered configurations.

For the depinning transition, both pinning and driving are included, i.e., $U_p \neq 0$ and $F_D \neq 0$.
Since it is unclear if MC can correctly capture transport dynamics, we
instead used MD for this study.
Details of the MD simulation of the overdamped Langevin dynamics
are discussed in Sec.~\ref{dt}.

The MC simulation methods used here are essentially the same as in \cite{Wang:VC}.
Each MC sweep for a replica is $N$ attempted updates of the particles in either sequential or random order.
The MC trial moves are attempts to shift the particle position randomly within a square box of length $2\sqrt{1/n}$ centered on the particle. We only use single-particle update moves.
For glass formation dynamics this corresponds to the overdamped regime.

The
main observables are the cluster orientational order parameter $\phi_6$, the
heat capacity $C_V=\beta^2 \rm{var}(E)$ where $E$ is the total energy of the system, and the pair correlation function $g(r)$.
The parameter $\phi_6$ quantifies the ordering of the triangular cluster-crystal phase.
The calculation of $\phi_6$ for a given configuration is based on the cluster positions
and we now discuss in detail how to compute it.

We first use a hierarchical clustering technique \cite{DataMining,Rogelio:Cluster} that group particles into clusters deterministically.
Each particle starts as a single cluster, and the two closest clusters are joined together if their center of mass distance is smaller than a chosen cutoff.
We used a cutoff of $0.7 \approx a/2$, where $a$ is the lattice constant of the corresponding cluster crystal in the absence of pinning disorder; see Fig.~\ref{PD}.
The results present in this work are not sensitive to the precise value of the chosen cutoff, e.g., around $0.7\pm 0.1$, as we work in parameter regimes where clusters are reasonably well defined.
When two clusters
are merged, the new center of mass is updated.
The process is repeated until no further grouping is accepted.
Hence the process takes a particle configuration $\{\vec{r}_i, i=1,2,...,N\}$
and outputs the centers of mass $\{\vec{R}_i, i=1,2,...,C\}$ of the $C$ clusters.
The clustering process is followed by a
Voronoi decomposition to identify neighbouring clusters.

The orientational order parameter $\phi_6$ is finally defined as:
\begin{eqnarray}
\phi_6 = \left| \dfrac{1}{C} \sum_{j=1}^{C} \left( \dfrac{1}{N_j} \sum_{\ell=1}^{N_j} e^{i6\theta_{j\ell}} \right) \right|,
\end{eqnarray}
where $\theta_{j\ell}$ is the angle
between the vector $\vec{R}_{\ell} - \vec{R}_j$ and an arbitrary direction, often the $\hat{x}$-axis,
and $N_j$ is the number of neighbours of the $j_{th}$ cluster.
Note that only neighbouring pairs are summed over.
For a perfect triangular lattice, $\phi_6=1$, and with no long-range order and random positions, $\phi_6=0$.
Hence $\phi_6$ is an order parameter.

The pair correlation function is defined as:
\begin{eqnarray}
g(r) = \frac{1}{N} \dfrac{\delta n(r)}{2\pi r \delta r},
\end{eqnarray}
where $\delta n(r)$ is the number of particles in the shell $2\pi r \delta r$, with reference to an arbitrary particle.
Note that the function is normalized as $\int_0^{\infty} g(r) 2\pi r dr =1$.
The lattice constant of a cluster configuration can be extracted from the peaks of $g(r)$.

\section{Numerical results}
\label{results}

\subsection{Cluster glass formation}

In order to study cluster glass formation, it is important to study the underlying pure system without pinning disorder for input, particularly the relevant particle density for clustering and the corresponding transition temperature. The phase diagram of the pure system along with all the ordered phases are depicted in Fig.~\ref{PD}. Here, the transition temperature is estimated from the prominent $C_V$ peak using $N=1000$, which is sufficiently large for a reasonably accurate estimation. The crystalline phases are identified by examining typical low-temperature equilibrium configurations, and they are also verified from ``ground states'' found by simulated annealing. It should be noted that our phase diagram here is only a sketch, and the more subtle features such as the number and the order of these phase transitions are not considered. For the potential in Eq.~(\ref{potential}), we find ordered phases of cluster crystals, stripes, and hole crystals as the particle density is increased up to about $n=18$. This series of ordered phases appears to be rather common in cluster-forming potentials with a repulsive core \cite{RR:phases,Rogelio:Glass}. Note that our potential restores the cluster crystal phase at even higher densities after the hole crystal phase.  This behaviour is different from the potentials with divergent cores near $r=0$; cf. \cite{RR:phases,Rogelio:Glass}. It should be noted that this second cluster phase at high density differs from the first cluster phase at low density with, e.g., a very distinct lattice constant. While pattern formations and the nature of these phase transitions are interesting in their own right, we here focus on the glass formation of the physically relevant low-density cluster phase.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{XY.pdf}
\includegraphics[width=\columnwidth]{PD.pdf}
\caption{
Top panels: Typical equilibrium configurations of the ordered phases of cluster crystals (CC), stripes (SP), hole crystals (HC) and simple crystals in order of increasing density.
Bottom panel: The schematic phase diagram of the potential [Eq.~(\ref{potential})] in absence of pinning disorder estimated from the major peak of the specific heat. The two cluster crystal phases have clustering onset densities at about $n_c=0.65$ and $30.5$, respectively.
}
\label{PD}
\end{center}
\end{figure}

We now present more details of the first triangular cluster phase. The onset density of clustering is at about $n_c=0.65$, below which there is only one particle per unit cell. When clustering occurs, the lattice constant settles to $a \approx 1.34$, which is estimated from the pair correlation function, and is independent of the density. Interestingly, both the onset density and the lattice constant are very similar to that of the ultrasoft potential $U(r)=1/(1+r^6)$ \cite{Wang:VC}.
This is presumably because the two potentials have in some sense similar features in the medium and long ranges. In the following, we work with a typical cluster phase at density $n=2$ where each cluster has approximately $3$ or $4$ particles. The transition temperature at this density is $\beta_c \approx 16.3$. We now turn to the cluster glass formation of this system.

In the absence of pinning disorder, clustering particles form cluster glasses when quenched from a high temperature, e.g., $T=\infty$ to a sufficiently low temperature in a single step \cite{Rogelio:Glass}. The system relaxes following the quench, and it was found that the system can reasonably restore the equilibrium order parameter except when the final quenching temperature $T_f$ is too low. The glass transition temperature is then loosely defined when the resulting order parameter $\phi_6$ departs from that of the equilibrium cluster crystal phase. It was argued in \cite{Rogelio:Glass} that cluster-forming particles are better glass formers compared with non-cluster-forming particles because of the effective polydispersity of the cluster sizes. On the other hand, it is well known that pinning by itself can lead to glass formation in conventional vortex matter \cite{blatter:94}. We next investigate the situation where both mechanisms are present.

For the simulation of thermal quench dynamics in presence of pinning disorder we consider different density and strength of pinning centers, and compare results with the case of no pinning. For each set of parameters we run quench dynamics to various final temperatures $T_f$ from
random configurations corresponding to $\beta=0$, using single-temperature Monte Carlo dynamics.
For all the temperatures that we studied $10^5$ sweeps were enough for a good relaxation, and then
the order parameter $\phi_6$ of the final configuration was measured.
Only one measurement per temperature was made from one disorder realization to prevent correlations, and statistics were collected using averages over disorder realizations. Here, $100$ realizations are simulated for each pinning density and strength.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{Dynamics2.pdf}
\caption{
Effects of density and pinning strength on glass formation after a thermal
quench as a function of the final temperature $T_f$ at density $n=2$ with $N=1000$ particles.
The topmost red curve in each panel is the pure case without disorder which provides an upper bound.
Left panels: Changing the number of pins
$M$ at fixed trapping strengths (a) $U_p=1$ and (b) $U_p=4$.
Right panels: Changing the pinning strength
at fixed number of pins for
(c) $M=20$ and (d) $M=40$.
}
\label{Dynamics}
\end{center}
\end{figure}

Results for dynamical glass formation are summarized in Fig.~\ref{Dynamics}.
In all panels, the red (topmost) curve corresponds to thermal quench of the clean system.
This result is similar to that of the molecular dynamics for other cluster-forming potentials \cite{Rogelio:Glass}.
There are two relevant temperature scales, one is the transition temperature of the pure system where
the system forms global order as temperature is lowered below $T_C$. This matches well with the transition
temperature we find from the
heat capacity.
The other temperature is the crystal-to-glass transition temperature $T_G$ \cite{Rogelio:Glass}, which can be estimated as where $\phi_6$ bends downwards and starts to drop, here, $T_G \approx 0.05$.
This is a characteristic temperature below which dynamics becomes glassy and the system readily runs out of equilibrium.

It is natural that the clean curve provides an upper bound for restoring the global order, since
adding random pinning
and increasing the pining strength generally
enhance glass formation.
However, it is interesting to note that the effect of pinning is far from linear.
Pinning actually has little effect on the order parameter if the disorder is weak or dilute. For example, at $U_p=1$ there is no prominent change of the order parameter up to about $M=160$. Further increasing the number of pins or the pinning strength eventually makes the order depart prominently from the clean limit. In the strong pinning limit, the global order is essentially entirely eliminated.
Another interesting observation that can be drawn from Fig.~\ref{Dynamics} is that addition of pinning gives the strongest suppression of the $\phi_6$ order
parameter at the intermediate temperatures, while at least for relatively dilute pinning the curves merge with the clean curve at lower
temperatures. This suggest that at least for dilute pinning the system falls into a state similar to
that of an interaction-dominated effective polydispersity class pinned by rare pinning centers.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{XY4.pdf}
\includegraphics[width=\columnwidth]{GR.pdf}
\caption{
Top panels: Typical cluster glass configurations corresponding to those of Fig.~\ref{Dynamics} at approximately $T=0.03$. The blue dashed circles are attractive Gaussian centers. Note that a fraction of some large clusters is formed in or near the (overlapping) traps. These configurations suggest that clusters are reasonably well defined in our parameter regimes.
Bottom panel: The clustering feature is further confirmed in the pair correlation function $g(r)$. Particularly, the second largest peak is the inter-cluster peak, or the lattice constant peak of the corresponding cluster crystal. The decay tail ($r \gtrsim 11)$ is a finite-size effect.
}
\label{XY4}
\end{center}
\end{figure}

It is important to confirm that the quenched glasses are genuine cluster glasses in presence of substrate disorder. To this end, we have checked the final configurations. Several typical cases are illustrated in the top panels of Fig.~\ref{XY4}. The blue circles are Gaussian traps, it is clear from observation that clusters are still reasonably well defined, despite that the attractive traps distort the cluster order as well as the distribution of the cluster sizes. For example, it is more likely to find larger clusters in or near the (overlapping) traps. Interestingly, these clusters do not necessarily all sit at the centers of the traps.  In the dense case $M=640$, there are a number of large clusters of sizes $5, 6$ and even $7$, leaving at the same time some small clusters in the trap interstitial voids.

Cluster glass manifests itself when the corresponding cluster crystal lattice constant peak is relevant in the pair correlation function. This is shown in the bottom panel of Fig.~\ref{XY4} for various $M$ values. It is seen that the second largest peak, which is the cluster crystal lattice constant peak, is present for all the cases studied. The one near $r=0$ is from the neighbouring particles within clusters, note that a small shoulder peak is developed for the dense cases $M=320$ and $640$, as some clusters are larger compared with the typical size of the free case. By contrast, the lattice constant of the particle triangular crystal $3^{-\frac{1}{4}}=0.7598$ at the same density is not (at least yet) relevant.

It is important to note that we should not take the weak and strong traps here too literally, as in our cases cluster glasses are \textit{all} reasonably well defined. It is clearly possible to further increase the trap depth, and new phases, e.g., particle glass may appear. It is expected that fully particle glasses should form in certain trap settings, depending on the trap size, depth, and density. Our main conclusion of this section is that interaction-driven cluster glasses are stable with respect to inclusion of quenched disorder.

\subsection{Cluster glass depinning dynamics}
\label{dt}

In this section, we use molecular dynamics to study the
depinning dynamics of the cluster glasses we obtained. We use the overdamped Langevin dynamics following the stochastic ordinary differential equations:
\begin{eqnarray}
\gamma\dot{\vec{r}}_i &=& -\nabla U_i + \sqrt{2T\gamma} \vec{\eta}_i(t) + F_D \hat{x},
\end{eqnarray}
where
without loss of generality we set the friction coefficient $\gamma=1$.
The Gaussian stochastic force $\eta$ has zero mean and are uncorrelated in time, $\langle \eta_{i\alpha}(t) \eta_{i\gamma}(t') \rangle = \delta(t-t') \delta_{\alpha\gamma}$, where $\alpha, \gamma = x, y$ denote components in different spatial directions. $U_i$ is the potential energy of particle $i$ interacting with other particles and the traps, and $F_D$ is applied along the $\hat{x}$-axis.
For the thermal noise in each direction, it is straightforward to show from the correlation function that $\int_t^{t+dt} \sqrt{2T} \eta_{i\alpha}(t') dt' = \sqrt{2Tdt}\zeta$, where $\zeta$ follows the standard Gaussian distribution $n(0,1)$.
In the large driving limit, we expect a drift velocity $v = \langle v_x \rangle = F_D$ as in the case of single particles \cite{RR:depinning}.
The force on a particle for the symmetric potential can be evaluated exactly using the relations $F(r) = -\frac{dV}{dr} \hat{r}$ and $d K_0(r)/dr = - K_1(r)$ of the modified Bessel functions of the second kind. The $K_1$ function is similar in shape to $K_0$ but it is further shifted away from $r=0$. The integration is done using the fourth-order Runge-Kutta method. We study the cluster glass depinning for a typical case of $U_p=1$ at $T=T_C/2$, where various features of the cluster flow can be studied. In certain cases, we have also performed
simulations at $T=0$ to suppress thermal fluctuations for a more clear visualization of the dynamics.

In the pure
case with no pinning, the clusters slide together along the direction of the driving force, and the drift velocity is given by the linear relation $v=F_D$ for any $F_D$; see the (topmost) red line of Fig.~\ref{Drive}. At finite temperatures, the particles will fluctuate while sliding, but the underlying cluster glass order is unchanged. When the temperature is quenched to $T=0$, the internal dynamics of the clusters will be rapidly frozen and the cluster glass essentially slides rigidly afterwards. In either case, the cluster glass remains intact.
These observations are expected since
without pinning the particles motion is given by a
superposition of the internal dynamics when driving is absent and a constant translational motion, i.e., in the absence of pinning the state can be simply identified as a {\it sliding cluster glass}.
A typical movie of such dynamics for $F_D=0.2$ at both finite and zero temperatures is presented in \cite{comment:mvmd1}.

In presence of pinning, a dynamical depinning transition should occur as the driving force is increased. We have conducted a large-scale MD simulation to measure the disorder-averaged drift velocity at various driving forces. For each run, we integrate up to a sufficiently long time $t=1000$, and the drift velocity is measured using the latter half of the simulation, i.e., the initial transit dynamics is ignored and we measure the asymptotic steady state drift velocity. The results are shown in Fig.~\ref{Drive}. For moderate and dense pinnings, there is a prominent pinned phase. The critical depinning force is naturally very small but nonzero when the pinning is dilute. Detailed examination near $F_D=0$ reveals that there are pinned phases as well for dilute pinnings.
For a depinning transition, a critical exponent $\beta$ can be defined as $v \sim (F-F_C)^\beta$ when $F \geq F_C$ in the vicinity of the depinning critical force $F_C$ \cite{Fisher:beta}. This is analogous to the order parameter scaling at a second order phase transition. As the curves are concave down near the transition, it suggests that $\beta>1$, i.e., the depinnings are likely plastic \cite{RR:depinning}. In plastic flows, particles or
clusters change their neighbours during the flow, contrary to the elastic flow where neighbouring particles
and clusters are preserved. A close examination
of the dynamics shows that the cluster flows are indeed plastic, with a number of features distinct from the non-cluster-forming particle systems, e.g., clusters undergo merging and breaking dynamics.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{Drive2.pdf}
\caption{
Drift velocity of the cluster glass subject to an external driving force at $T=T_C/2$ for, from left to right, the pure system, $M=10, 40, 160, 320$, and $640$, respectively.
The inset panel is a zoom of the main panel near the origin.
For the pure case, there is a linear force-velocity relation.
For dense pinning there is a prominent depinning transition and the drift velocity converges slowly to the free case.
For
dilute pinning the critical force is much smaller and the drift velocity converges rapidly to the pure case. The concave down of the depinning curves near the transition suggests plastic depinnings, which are subsequently confirmed in the details of the flowing dynamics; see the text for details.
}
\label{Drive}
\end{center}
\end{figure}

We discuss the dense pinning first, where the flowing phase is clearly plastic. The flowing dynamics is quite complex, and here we focus on the main features and typical dynamical processes. Considering the complexity of the flow, a movie showing both the pinned and flowing phases are shown in \cite{comment:mvmd2},
where the driving forces are $F_D=0.0508$ and $0.2$, respectively, for a disorder realization of $M=640$. There are a number of interesting observations for the flowing phase. First, cluster chaotic flow occurs instead of the particle chaotic flow in non-cluster-forming ensembles. Here, the clusters have a strong tendency to flow through or near the traps, forming rivers and changing neighbours. Interestingly, the pinning centers are like stepping stones for the clusters. In addition, the clusters are no longer as robust as in the free case. There are frequent clusters merging and breaking dynamics. For example, unusually large clusters may form as a result of merging and there are also a number of single-particle clusters as a result of
particle emission from clusters. As a result, the clusters are much more heterogeneous than the pure case due to the interplay of pinning and driving. It is important to emphasize that this increased heterogeneity is \textit{not} a result of thermal noise as we are working in the low-temperature glass phase, despite that it also contributes. This is clearly demonstrated as the above processes also occur at $T=0$ as shown in the movie in \cite{comment:mvmd3}. A typical cluster merging process and also a cluster splitting process are depicted in Fig.~\ref{pf}. It is interesting that trap centers play a ``catalytic'' role in these dynamics. Indeed, these dynamics are much suppressed in the dilute pinning regime.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{Plasticflow.jpg}
\caption{
Typical clusters merging (top panels) and splitting (bottom panels) dynamics induced by the trap centers ($T=0$), the pertinent clusters are in blue and pink while other clusters are in red, and the green circles are pinning centers. In this dense pinning regime, such processes are frequent, rendering the flow highly plastic.
}
\label{pf}
\end{center}
\end{figure}

In the opposite dilute pinning regime, the flow is much more ordered. After close inspections of $M=10$, we find that the flow is still plastic, but it is very weak.
Upon depinning the glass mostly slides over the pinning centers. The plastic flow and cluster merging and breaking are much suppressed. Nevertheless, the effective polydispersity of the clusters can still undergo trap-induced neighbour-changing dynamics because of their different mobilities. Similar to a cluster liquid \cite{Wang:VC}, small clusters tend to have better mobilities, and moreover larger clusters are better pinned by impurity traps.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{Plasticweak.jpg}
\caption{
Example of a plastic neighbour-changing process for a dilute pinning case with $M=10$. The clusters can change their neighbours due to their polydispersity and consequently their mobilities. Here, the smaller black cluster depins faster and has a better mobility. As a result, when the blue cluster depins, the black cluster is no longer a neighbour of the blue cluster, inducing weak plastic flows.
}
\label{pw}
\end{center}
\end{figure}

Figure~\ref{pw} illustrates a neighbour-changing plastic flow process in the dilute pinning case. Here, the two relevant clusters are highlighted in blue and black. The blue cluster has four particles while the black cluster has three particles. The black cluster depins at around $t=450$ and afterwards the blue cluster arrives at the pinning center and is temporarily pinned. This pinning is stronger than that of the black one as there are more particles, and also importantly a smaller cluster has a higher mobility and is easier to diffuse around. As a result, although the blue cluster eventually depins also at around $t=1000$, the black cluster has already diffused away and is thus no longer a neighbour of the blue cluster. This process makes the flow weakly plastic. We emphasize, however, that such event does not occur frequently in dilute pinning; compare the time scale with Fig.~\ref{pf}.

To summarize this section,
all the depinning flows that we studied are plastic cluster flows.  Clusters tend to flow along pinning
channels when pinning is dense. The effective polydispersity assists plastic flows in two aspects. First, it generates new forms of plastic flow processes such as clusters merging and breaking. In addition, their different mobilities can also induce neighbour-changing plastic flows when flowing through the pinning centers. These are quite distinct features from the depinning of ordinary non-cluster-forming glasses.

\subsection{Dynamical reordering, critical exponent, and comparison with particle glass}

In this section we first demonstrate that cluster glasses upon depinning can undergo dynamical reordering when the driving force is increased. Next, we measure the critical exponent $\beta$ of the depinning transition. Finally, the cluster glass pinning strength is compared with that of ordinary particle glass, finding that cluster glass can have much stronger pinning effects than particle glass. This is potentially useful in superconductor applications.

It is well known that the particle glass and the stripe pattern-forming systems can exhibit dynamical reordering as the driving force is increased \cite{RR:stripes,RR:depinning}. For example, stripes may break into pieces in presence of pinning disorder, but upon depinning stripes may reform when the driving force gets larger. Here, we investigate this phenomenon for the cluster glass. To study the dynamical reordering, we look at an order statistic $P_6$. This is the fraction of clusters having $6$ cluster neighbours in the asymptotic steady state as a function of the driving force. The result of $M=160$ is shown in Fig.~\ref{P6}.
It is interesting that $P_6$ drops very sharply at the onset of the plastic flow, the critical force here is in good agreement with the depinning transition of Fig.~\ref{Drive}. This indicates that the system in the flowing phase is most chaotic when is just depins, and pinning and driving compete most strongly near the depinning transition. As the driving force increases, it is remarkable that $P_6$ is gradually restored and the effects of pinning are increasingly washed out; note that the force scale is larger than that of Fig.~\ref{Drive}. When $F_D \gtrsim 0.8$, $P_6$ is almost the same as that of the pinned phase. However, it seems that the reordering is by no means perfect, it does not achieve a perfect order for the largest force we studied, which is approximately $50$ times larger than the depinning force. In addition, the reordering occurs gradually and appears to be a crossover rather than a genuine phase transition.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{P6.pdf}
\caption{
Cluster glass plastic flow exhibits dynamical reordering when the driving force is increased. Here the fraction $P_6$ of clusters with $6$ neighbouring clusters is shown as a function of the driving force for $M=160$. The fraction decreases significantly at the onset of the plastic flow around $F_D=0.02$, but it then gradually restores when $F_D$ is increased as the pinning effects are increasingly washed out. Note that this dynamical reordering appears to be a crossover rather than a genuine phase transition.
}
\label{P6}
\end{center}
\end{figure}

We now study the critical exponent $\beta$ of the depinning transitions. The critical exponent $\beta$ is defined as $v \sim (F-F_C)^\beta$ when $F \geq F_C$ in the vicinity of the critical force $F_C$. When the pinning is dilute, the range of $F_D$ where a power-law behaviour applies is limited; see Fig.~\ref{Drive}. By contrast, the depinning in the dense pinning regime has a much wider range of forces pertinent to scaling, and we therefore compute $\beta$ using the cases $M=320$ and $640$. We use a cubic function fit in the flowing phase but keep away from the transition to first locate the critical force $F_C$ and then make a power-law fit. In both cases, we find good power-law behaviours and the scaling plots are shown in Fig.~\ref{beta}. The two exponents are close but slightly different, $\beta=1.049(5)$ and $1.351(13)$, respectively. Adding two intermediate cases $M=400$ and $500$, we find that there appears to be a slow drift of the exponent $\beta$.

This exponent is unfortunately in general very poorly understood, and it is yet unclear whether there is a universality class for plastic depinnings. Our case appears to be even more complicated as our starting configurations are by design quenched glasses rather than slowly annealed states. Hence, it is not clear whether different pinning densities, corresponding to different degrees of glassiness, necessarily have the same critical exponent. More pinning centers introduce more defects into the quenched cluster glass, and therefore the depinning could be more plastic for a denser pinning, yielding a slightly larger exponent.
On the other hand, our computed exponents also fall within the broad range of exponents
found for plastic depinnings in many different systems \cite{RR:depinning}.
Therefore, this could be a finite-size effect or finite-temperature effect and there is a universality class, or alternatively quenched glasses may be different from annealed glasses. This intriguing question requires further systematic studies in its own right, and we shall not discuss it further here.

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{beta.pdf}
\caption{
Power-law scaling of the velocity $v$ versus $F_D-F_C$ near the transition for $M=320, 400, 500$, and $640$. While the power-law behaviour appears to be reasonably good in all cases, the exponent or the slope becomes slightly larger with the increasing number of pinning centers, yielding $\beta=1.049(5), 1.121(5), 1.233(8)$ and $1.351(13)$, respectively. See the text for more discussions.
}
\label{beta}
\end{center}
\end{figure}

\begin{figure}[htb]
\begin{center}
\includegraphics[width=\columnwidth]{ST.jpg}
\caption{
Velocity-force curve for particle glass (left) and cluster glass (right) for $M=320$ and $U_P=1$. The trap size is $\sigma=0.3799$, which is half the lattice constant of the particle crystal. Both glass formation and molecular dynamics are performed at their respective $T_C/2$. Note that the depinning force of the cluster glass is much larger than of the particle glass, suggesting that vortex clusters can have much better pinning properties than ordinary single vortices.
}
\label{ST}
\end{center}
\end{figure}

Finally, we
compare the depinning forces of cluster glasses and particle glasses. To this end, we have simulated a group of particles at $N=1000$ and $n=2$ but interacting with a single Bessel function potential as shown in Fig.~\ref{UR}. The numerical setup is essentially the same for particle glasses, except that the temperature scales differ. Here, the transition temperature of the pure system is much smaller $T_C\approx0.01077$. Similarly, we prepare glasses and do molecular dynamics in presence of pinning at its own $T_C/2$.

We consider $M=320$ pins and
$N=1000$ particles. We have similarly checked that the quenched states are particle glasses with $\phi_6 \lesssim 0.1$.
There are two natural choices of the trap size $\sigma$ for a fair comparison, i.e., we can use either half the lattice constant of the cluster crystal or the particle crystal. For the same trap depth $U_p=1$, they yield wider or narrower traps, or weaker or stronger pinnings, respectively. Fortunately, they yield qualitatively the same conclusion. When wider traps are used, the particle glass has a very small depinning force that is not easily measurable, i.e., the traps are too wide to pin the particle glass effectively.
Hence we use smaller traps for \textit{both} systems, and the results are shown in Fig.~\ref{ST}.
It is found that the cluster glass again pins much better than the particle glass, the critical forces are approximately $F_C=0.01$
for particle glass and $0.3$ for cluster glass. There are many parameters including the interaction potentials that one can tune to explore this more systematically.
However, considering that the depinning force scales are very different in the cases studied here, it is quite suggestive that cluster-forming type-1.5 superconductors can have much better pinning properties than system of ordinary vortices.

\section{Conclusions}
\label{cc}

We studied an ensemble of two-dimensional monodisperse particles interacting via a cluster-forming potential, relevant to vortices in type-1.5 superconductors in the presence of demagnetizing fields and layered type-1.5 systems.
Vortices in such systems were previously demonstrated to form glasses in the absence of pinning, purely due to
peculiarities of intervortex forces.

We examined cluster glass formation following a thermal quench in presence of pinning centers of moderate depth, and the corresponding depinning transition. We find that weak and dilute pinnings have little effects on the interaction-dominated glass formation: the disordered
configuration adapts to the presence of dilute pinning sites, preserving the intervortex-interaction-driven effective polydispersity.

However, moderate, strong and dense pinnings can substantially assist glass formation, especially at intermediate temperatures, and indeed can entirely remove  order. For the pinning strength we have studied, the depinning transition is plastic and the cluster plastic flow exhibits features distinct from particle plastic flow. This includes cluster chaotic flow, clusters merging and breaking, and their different mobilities. Therefore, the effective polydispersity of clusters assists both glass formation and plastic flows. The cluster plastic flow is also found to exhibit dynamical reordering when the driving force is large.

Finally, a conclusion of practical importance is that
cluster glasses can have significantly better pinning properties than systems
of non-cluster-forming particles.
Since vortex pinning is important for transport properties of superconductors, it suggests that these
properties of  type-1.5 superconductors can be used in applications. Cluster-forming potentials can be engineered for example by layering systems \cite{meng16} or by using superconductors where there is an additional phase transition associated with the time-reversal symmetry breakdown \cite{carlstrom2011length,garaud2018properties}.

\acknowledgments
We thank Cynthia Reichhardt for helpful exchanges on depinning transitions, and Alexander Zyuzin for helpful discussions.
W.W. and E.B. acknowledge support from the G\"{o}ran Gustafsson Foundation for Research in Natural Sciences and Medicine. E.B. acknowledges support from Swedish Research Council Grants No. 642-2013-7837, 2016-06122, 2018-03659, and the Olle Engkvists Stiftelse. W.W. also acknowledges support from the Fundamental Research Funds for the Central Universities, China. The computations were performed on resources provided by the Swedish National Infrastructure for Computing (SNIC) at the National Supercomputer Center in Link\"{o}ping, Sweden, and the High Performance Computing Center North (HPC2N) partially funded by the Swedish Research Council through grant agreement No. 2018-05973


\section{Introduction}
\label{sec:I}

Three-body collisions are commonly associated with atom loss and heating in ultracold atomic gases.
However, also interesting but subtle three-body phenomena such as the Efimov effect \cite{Efimov:1970} appear in the recombination rate and modify its behavior.
The now fifty years old prediction of the Efimov effect was experimentally discovered only in 2006 \cite{Kraemer:2006} in an ultracold gas of cesium atoms.
Key to this breakthrough observation is the
tunability of the scattering length $a$ that parametrizes the two-body interaction strength at ultralow temperatures \cite{Chin:2010}.
The Efimov regime is determined by $|a|/r_\mathrm{vdW} \gg 1$, where $r_\mathrm{vdW}$ characterizes the range of the van der Waals attraction between the atoms.
This regime is accessible close to a Feshbach resonance where $a$ goes through a pole.

In the Efimov regime the three-body system shows universal behavior that does
not depend on the details of the two-body interaction \cite{Efimov:1970, braaten:2006, Greene:2017, Naidon:2017, D_Incao:2018}.
Remarkably an infinite sequence of loosely bound three-body states emerges
for resonant two-body interactions.
On resonance the binding energies $E_n$ of these trimers follow the universal scaling relations $E_{n+1} / E_n = \mathrm{e}^{-2 \pi / s_0}$ with $s_0 \approx 1.00624$ for identical bosons \cite{braaten:2006, Greene:2017, Naidon:2017, D_Incao:2018}.
Those scaling relations also transfer to related quantities like the scattering lengths $a_-^{(n)}$ at which the $n$-th Efimov trimer state hits the three-body threshold and causes an Efimov resonance in the three-body recombination rate.
In the universal regime the position of Efimov features is determined by a single three-body parameter, which is often determined experimentally from the position of the lowest Efimov resonance $a_-^{(0)}$ in ultracold atomic systems.
Following the pioneering work in 2006 \cite{Kraemer:2006} the $a_-^{(0)}$ has been measured over a wide range of species \cite{Kraemer:2006,pollack:2009,Gross:2009,Gross:2010,zaccanti:2009,Wild:2012,Ferlaino:2011,Berninger:2012}.
Surprisingly many of the
three-body parameters were found to have roughly the value $a_-^{(0)} / r_\mathrm{vdW} \approx - 9$
\cite{Kraemer:2006,Gross:2009,Gross:2010,Wild:2012,Ferlaino:2011,Berninger:2012}.
The discovery of the origin of this van der Waals universal behavior was a major theoretical success in recent years \cite{Wang:2012,Naidon:2014a,Naidon:2014l}.
However, there are still some experimental results \cite{Roy:2013,Chapurin:2019} which pose an exception to universality.
To describe those results complex numerical models taking also the atomic spin structure into account
are necessary \cite{Chapurin:2019}.

Advanced numerical models are also needed outside the universal Efimov regime, when $|a|/r_\mathrm{vdW} \lesssim 1$.
In this regime three-body recombination persists to constitute a major loss mechanism in an ultracold Bose gas,
but the universal expressions fixed by a three-body parameter are no longer valid.
In addition to the total recombination rate experiments can now also reveal partial recombination rates by identifying the recombination products and are thus ranging in the realm of ultracold chemistry \cite{Haerter:2013, Wolf:2017}.
Even in the regime where $a / r_\mathrm{vdW} \approx 0 $ elastic three-body effects have recently been proposed to determine the phase diagram of a Bose-Einstein condensate \cite{Zwerger:2019,Mestrom:2019,Mestrom:2020}.

The above mentioned examples substantiate that the fast experimental progress creates a demand to advance state-of-the-art theoretical models to calculate three-body effects ranging from $a/r_\mathrm{vdW} \approx 0 $ to $|a|/r_\mathrm{vdW} \gg 1$.
Therefore new numerical approaches are needed, which allow to calculate the three-body problem in an efficient way.
The Alt-Grassberger-Sandhas (AGS) equations \citep{Alt:1967} are one way to formulate the three-body problem.
To solve the AGS equations numerically in momentum representation, many partial wave components of the two-body $t$-operator need to be calculated for a large number of energy points, which puts constraints on calculation time and accuracy.
Fortunately, there are many coordinate space methods  \cite{Feit:1982,Monovasilis:2007,Tal:1984,Light:1985,Fattal:1996,Willner:2004,Karman:2014} to solve the two-body problem efficiently, whose capacity has been demonstrated in
two-body scattering calculations.
A for our purpose advantagous category of methods is based on
a Discrete Variable Representation (DVR) \cite{Light:1985} or a mapped DVR \cite{Fattal:1996,Willner:2004}.
These approaches lead to an approximate finite dimensional matrix representation of the Hamiltonian, which can be directly used to calculate the $t$-operator via
\begin{equation}
t(z^{2b})=V+\sum_{i} V \frac{|\psi_{i}\rangle\langle \psi_{i}|}{z^{2b}-E^{2b}_{i}} V, \label{tintro}
\end{equation}
where $|\psi_{i}\rangle$ and $E^{2b}_{i}$ are eigenvectors and eigenvalues of the Hamiltonian matrix, respectively and $V$ denotes the pairwise interaction potential.
It should be noted that the diagonalization of the Hamiltonian matrix only needs to be done once, after which the $t$-operator can be calculated at any two-body energy $z^{2b}$ using Eq.~(\ref{tintro}).

In numerical
practice we consider a
finite relative distance between the particles
for this two-body problem.
The (mapped) DVR brings the free-space two-body system into a finite distance region with specific boundary conditions.
A hard wall boundary condition is frequently chosen for some numerical benefits, especially in the mapped case \cite{Willner:2004}.
Even though the finite distance region with hard wall boundary condition has only minor impact on traditional calculations of bound state energies and wave functions, its influence on the $t$-operator which includes off-shell scattering
properties remains to be determined and is subject of study in this paper.
A study of the effect of restricting to a finite distance region on the two-body off-shell $t$-matrix and its consequences on three-body quantities is critical to clarify whether most aformentioned numerical methods \cite{Feit:1982,Tal:1984,Light:1985,Fattal:1996,Willner:2004} can facilitate the three-body calculation.

This paper is organized as follows.
In section \ref{sec:IIA} we review the AGS equations related to three-body recombination of identical bosons interacting via pairwise interaction potentials.
In section \ref{sec:IIB} we introduce an analytic model including the finite distance region approximation with a hard wall boundary condition for the two-body $t$-operator based on a square-well interaction to analyze the validity and convergence properties of the approximation.
In section \ref{sec:IIC} we then review the mapped DVR method that we apply to a Lennard-Jones van der Waals potential with realistic long range interaction properties.
In section \ref{sec:III} we present our results for the three-body recombination rate for the square-well and Lennard-Jones potential.
We compare the results for different sizes of the finite distance region over a wide range of scattering lengths and analyze the influence of the finite distance approximation on the position of Efimov features.
Finally we analyze the convergence properties of our approach with respect to a momentum cut-off in the AGS equations for Lennard-Jones potentials that support almost 4 and 6 $s$-wave bound states.

\section{Theory}
\label{sec:II}

\subsection{Three-body recombination}
\label{sec:IIA}

We consider a system of three identical bosonic alkali metal atoms.
The interaction in the system is described by pairwise interaction potentials $V_\alpha$, where the index $\alpha=(ij)$ indicates that the interaction takes place between particles $i$ and $j$.
We calculate the three-body recombination rate $K_3$ at zero kinetic energy.
This quantity is relevant since it is directly related to the loss rate in a sample of ultracold atoms \cite{
Braaten:2008}.

To calculate $K_3$ we start from the AGS equation for three identical particles that define the transition operator $U_{\alpha 0}(z)$ related to three-body recombination into a dimer state of particles $\alpha = (ij)$ and a free particle $k$ \cite{Alt:1967,Mestrom:2019b}
\begin{align} \label{eq:AGSeqrecom}
U_{\alpha 0} (z) & = G_0^{-1}(z) \left[1 + P_+ + P_-\right]/3 \nonumber \\
& \phantom{=} + \left[ P_+ + P_- \right] \mathcal{T}_\alpha (z) G_0(z) U_{\alpha 0} (z) \, .
\end{align}
The operators $U_{\alpha 0}$, $G_0$ and $\mathcal{T}_\alpha$ depend on the complex energy $z$. $G_0$ denotes the Green's operator of the free three-body system and is defined as
\begin{equation}
G_0 (z) = (z - H_0)^{-1} \, ,
\end{equation}
whereas $\mathcal{T}_\alpha $ is related to the two-body  $t$-operator and given by
\begin{equation}
\mathcal{T}_\alpha (z) = (1 - V_\alpha G_0 (z))^{-1} V_\alpha \, .
\end{equation}
In the following we will omit the explicit dependence on $z$ for notational compactness unless it is needed.
The operators $P_+$ and $P_-$ are the cyclic and anticylclic permutation operators, respectively.

We restrict ourselves to the case of zero total angular momentum in the system, which is suitable for the low collision energy limit in an ultracold system.
For recombination from a free incoming state $\Psi_\mathrm{in}$ of energy $E$ into a $\alpha$-dimer state labeled by $d$ with wave function $\varphi_d$ and energy $E_d$ plus a free atom of absolute momentum $q_d$ relative to the dimer center-of-mass, one needs to evaluate the transition operator element
\begin{align}
& U_{\mathrm{rc}}\left( \{ \Psi_\mathrm{in} \} \to \{  q_d, \varphi_d \} \right) \nonumber \\
& \equiv {}_\alpha \langle q_d, \varphi_d | U_{\alpha 0} (z) | \Psi_\mathrm{in} \rangle
\end{align}
on the energy shell. This leads to $E = 3 q_d^2 / 4 m + E_d $ with $m$ the mass of an atom and $z = E + i0$, which means that we take the limit in $z$ from the upper half of the complex energy plane.
In this case the relative angular momentum  between the atom and dimer in the final state is determined by the angular momentum of the dimer and the requirement that the total angular momentum needs to be zero.

We consider the limit of zero kinetic energy in the incoming state as is common for systems of ultracold atoms.
The recombination rate at zero energy is then given by \cite{Moerdijk:1995,Lee:2007,Smirne:2007}
\begin{align}
K_3 (0)
&= \frac{24 \pi m}{ \hbar} (2 \pi \hbar)^6 \nonumber \\
& \phantom{=} \sum_{d} q_d |U_{\mathrm{rc}}\left( \{ \Psi_\mathrm{in} \} \to \{  q_d, \varphi_d \} \right) |^2 \, .
\end{align}
Note that we follow the conventions of \cite{Braaten:2008} in defining $K_3$ which deviates from \cite{Moerdijk:1995,Lee:2007,Smirne:2007} by a factor of $2$.

The on-shell transition operator elements can be rewritten as
\begin{align}\label{eq:transopelms}
& {}_\alpha \langle q_d, \varphi_d | U_{\alpha 0} | \Psi_\mathrm{in} \rangle \nonumber \\
& = {}_\alpha \langle q_d, \varphi_d | V_\alpha G_0 \left( P_+ + P_- \right) \mathcal{T}_\alpha  G_0  U_{\alpha 0} | \Psi_\mathrm{in} \rangle \, ,
\end{align}
since the inhomogeneous term in Eq.~(\ref{eq:AGSeqrecom}) evaluates to zero in the on-shell limit and ${}_\alpha \langle q_d, \varphi_d | = {}_\alpha \langle q_d, \varphi_d | V_\alpha G_0$.
In our numerical treatment it is advantageous to consider the operator
\begin{equation}
A_\alpha = 3 G_0 \left( P_+ + P_- \right) \mathcal{T}_\alpha  G_0  U_{\alpha 0} \, .
\end{equation}
With Eq.~(\ref{eq:AGSeqrecom}) we obtain
\begin{equation}
A_\alpha = G_0 \left( P_+ + P_- \right)  \mathcal{T}_\alpha \left[ 1 + P_+ + P_- + A_\alpha \right]\, .
\end{equation}

We expand $\mathcal{T}_\alpha = \int dq q^{2}\sum_i \tau (q,i) | q, i \rangle_{\alpha \alpha} \langle q, i |$ and use the incoming state $\Psi_\mathrm{in}$, such that we arrive at the linear system
\begin{align}\label{eq:3bodyint}
& {}_\alpha \langle q', i | A_\alpha | \Psi_\mathrm{in} \rangle  \\
& = \int dq q^{2}\sum_j {}_\alpha\langle q',i | G_0 (P_+ + P_-) | q , j \rangle_\alpha \tau (q,j) \nonumber \\
& \phantom{=} \left[ {}_\alpha \langle q, j| (1 + P_+ + P_-) | \Psi_\mathrm{in} \rangle +  {}_\alpha \langle q, j | A_\alpha | \Psi_\mathrm{in} \rangle \right] \, . \nonumber
\end{align}
The on-shell transition operator element in Eq.~(\ref{eq:transopelms}) is then directly related to the single components  ${}_\alpha \langle q_d , i | A_\alpha | \Psi_\mathrm{in} \rangle$, since the expansion base ${}_\alpha \langle q, i  |$ naturally includes terms ${}_\alpha \langle q_d, \varphi_d |V_\alpha$.
More details on the linear system can be found in appendix \ref{B}.

\subsection{Off-Shell Scattering}
\label{sec:IIB}
\begin{figure*}[htb!]
	  \centering
	\begin{minipage}[c]{1 \columnwidth}
	  \includegraphics[width = 1 \textwidth]{figure1}
	\end{minipage}
	\centering
	\begin{minipage}[c]{1 \columnwidth}
	\includegraphics[width = 1 \textwidth]{figure2}
	\end{minipage}
\caption{We compare $t_\ell^{r_\mathrm{b}}(z,0,0)$ for
$r_\mathrm{b} / r_\mathrm{sw}$ equal to $\infty$ (solid line), $5000$ (dashed line), $2000$ (dash-dotted line) and $1000$ (dotted line) for negative (left) and positive (right) scattering lengths.}
\label{fig:1}
\end{figure*}

The operator $\mathcal{T}_{\alpha}$ in the AGS equation contains all information about the interaction between the atoms.
It is directly related to the two-body transition operator $t$ by
\begin{align}
\mathcal{T}_\alpha(z) &= \int d \mathbf{q} d \mathbf{k} d \mathbf{k}' \langle \mathbf{k}' |t(z - 3 q^2 / 4 m ) | \mathbf{k} \rangle  \\
& \phantom{=} \qquad \times |  \mathbf{k}', \mathbf{q}  \rangle_\alpha {}_\alpha  \langle  \mathbf{k}, \mathbf{q} | \, , \nonumber
\end{align}
with $\mathbf{k}$ or $\mathbf{k}'$ the relative momentum between atoms $i$ and $j$ and $\mathbf{q}$ the momentum of atom $k$ relative to the center-of-mass of the pair $\alpha=(ij)$.
We intend to compute the $t$-matrix directly using a mapped DVR
approach described in Section \ref{sec:IIC} in combination with Eq.~(\ref{tintro}).
However, in DVR practice a finite region $[0,r_\mathrm{b}]$ with a hard wall boundary condition is introduced in the relative separation of the two atoms.
Therefore we shall analyze how well the full  $t$-operator of a two particle system can be approximated by the $t$-operator of a system with a hard wall boundary condition at a finite particle separation of $r_\mathrm{b}$.

This can be done in a clear and easy way
by analyzing a square-well interaction potential, since $t(z^{2b})$ can be worked out analytically and main features like locality and the finite range of the atomic interaction are maintained.
We define the square-well interaction potential by
\begin{equation}
V_\mathrm{sw} =
\begin{cases}
V_0, & r < r_\mathrm{sw}\\
0 , & r \ge r_\mathrm{sw}
\end{cases} \, ,
\end{equation}
with $-V_0$ the depth and $r_\mathrm{sw}$ the range of the potential.
In case of a spherically symmetric potential $t(z^{2b})$ can be split into its partial wave components as
\begin{align}
 t(z^{2b})
& = \int d\mathbf{k}' d\mathbf{k} \,| \mathbf{k}' \rangle \langle \mathbf{k} |\\
& \phantom{=} \times \left[\sum_{\ell, m_\ell} \mathcal{Y}_{\ell, m_\ell}(\hat{\mathbf{k}}') t_\ell (z^{2b},k',k) \mathcal{Y}^*_{\ell, m_\ell}(\hat{\mathbf{k}}) \right]  \nonumber \, ,
\end{align}
with $\mathcal{Y}_{\ell, m_\ell}(\hat{\mathbf{k}})$ a spherical harmonic function in direction $\hat{\mathbf{k}}$.
To indicate the size of the finite distance region $r_\mathrm{b}$ we switch to the notation $ t_\ell^{r_\mathrm{b}} (z^{2b},k',k)$ with $ t_\ell^{\infty} (z^{2b},k',k)$ corresponding to the free-space case.
Changing to units where $\hbar$, the particle mass $m$ and the square-well radius $r_\mathrm{sw}$ are equal to one we get
\begin{widetext}
\begin{align} \label{eq:topsymmetric}
 t_\ell^{r_\mathrm{b}} (z^{2b}, k', k) &= \frac{2}{\pi} \frac{V_0}{k' k} \left[ \frac{-(k'^2 + k^2)(q_z^2 + q_0^2) + 2 k'^2 k^2 - 2 q_z^2 q_0^2}{2 ( k^2 - k'^2 )( k'^2 - q_0^2) (k^2 - q_0^2)} \mathrm{W}[\mathrm{S}_\ell (k r),\mathrm{S}_\ell (k' r)] \right. \\
 & \phantom{=} \left. - \frac{V_0}{2 (k'^2 - q_0^2)(k^2 - q_0^2) } \left( \frac{\mathrm{W}[\mathrm{S}_\ell(k r),  \phi_\ell^{r_\mathrm{b}}( q_z, r)] \mathrm{W}[\mathrm{S}_\ell(q_0 r), \mathrm{S}_\ell(k' r) ]}{\mathrm{W}[\mathrm{S}_\ell(q_0 r), \phi_\ell^{r_\mathrm{b}}( q_z, r))]} + k \leftrightarrow k' \right) \right]_{r = r_\mathrm{sw} = 1} \, , \nonumber
\end{align}
\end{widetext}
where the symbol $k \leftrightarrow k'$ represents the expression it is in brackets with, but with $k$ and $k'$ interchanged. We also used $q_z \equiv \sqrt{z^{2b}}$ and $q_0 \equiv \sqrt{z^{2b} - V_0} $, while $\mathrm{W}[\cdot, \cdot]$ denotes the Wronskian and $\mathrm{S}_\ell (\xi) \equiv \xi \mathrm{j}_\ell (\xi)  $ is a Riccati-Bessel function that we define via the spherical Bessel function of the first kind $\mathrm{j}_\ell$.
$\phi_\ell^{r_\mathrm{b}}( q_z, r)$ is the outer solution to the Hamiltonian differential equation $(q_z^2 - H_\ell)\phi_\ell^{r_\mathrm{b}} = 0$ with boundary condition $\phi_\ell^{\infty}( q_z, r)$ exponentially decaying as $r \rightarrow \infty$ or $\phi_\ell^{r_\mathrm{b}}( q_z, r_\mathrm{b}) = 0$ in case of finite $r_\mathrm{b}$.
$H_\ell$ is the partial-wave two-body Hamiltonian in the relative separation $r$ between the particles.
A derivation can be found in appendix \ref{A}.

We analyze the quality of the finite $r_\mathrm{b}$ approximation by considering the difference
\begin{align} \label{eq:topdifference}
& t_\ell^\infty(z^{2b}, k', k) - t_\ell^{r_\mathrm{b}}(z^{2b}, k', k) \nonumber \\
&= h_\ell(z^{2b}, k') h_\ell(z^{2b}, k) f_\ell(z^{2b}, r_\mathrm{b}) \, ,
\end{align}
which separates in the momenta $k, k'$ and the size of the finite distance region $r_\mathrm{b}$ (see appendix \ref{A} for more details).
We also define the functions
\begin{equation}
h_\ell (z^{2b}, k) = \left. \sqrt{\frac{2}{\pi}} \frac{V_0 \mathrm{W}[\mathrm{S}_\ell (q_0 r), \mathrm{S}_\ell (k r)]}{k (k^2 - q_0^2)} \right|_{r = 1} \, ,
\end{equation}
which are well behaved in $k$ for fixed $z^{2b}$
and
\begin{align}
&f_\ell (z^{2b}, r_\mathrm{b})  \\
&= \left. \frac{ - \mathrm{W}[\phi_\ell^{r_\mathrm{b}}( q_z, r), \phi_\ell^\infty( q_z, r)]}{\mathrm{W}[\mathrm{S}_\ell( q_0 r), \phi_\ell^\infty( q_z, r)]\mathrm{W}[\mathrm{S}_\ell( q_0 r), \phi_\ell^{r_\mathrm{b}}( q_z, r)]} \right|_{r = 1} \, . \nonumber
\end{align}
For $\mathrm{Im}(z^{2b}) \geq 0$ we find that the Wronskian in the numerator
\begin{equation}
\left. \mathrm{W}[\phi_\ell^{r_\mathrm{b}}( q_z, r), \phi_\ell^\infty( q_z, r)] \right|_{r = 1} \propto \mathrm{k}_\ell(r_\mathrm{b} \sqrt{-z^{2b}})
\end{equation}
is proportional to a modified spherical Bessel function $\mathrm{k}_\ell$ that behaves like $\sim \mathrm{e}^{-r_\mathrm{b} \sqrt{-z^{2b}}} /( r_\mathrm{b} \sqrt{-z^{2b}}) $ in the limit $r_\mathrm{b} \sqrt{-z^{2b}} \rightarrow \infty $ and thus guaranties
\begin{equation}
t_\ell^{r_\mathrm{b}} (z^{2b}) \xrightarrow[r_\mathrm{b} \rightarrow \infty]{} t_\ell^\infty (z^{2b})
\,
\end{equation}
for all $z^{2b} \notin \mathbb{R}^+_0$ with $\mathrm{Im}(z^{2b}) \geq 0$.
We can conclude that the convergence in $r_\mathrm{b}$ will be slowest for $z^{2b}$ close to $\mathbb{R}^+_0$.
We focus on the regime close to $z^{2b} = 0$ since for zero energy three-body scattering $z^{2b} \le 0$.
In Fig.~\ref{fig:1} we compare $t_0^{r_\mathrm{b}}(z^{2b},0,0)$ for different sizes of the finite distance region and scattering lengths between the first and second potential resonance.
We find similar good convergence properties for positive and negative scattering lengths, while the deviation between the finite distance approximation and the free-space case increases  almost quadratically with $|a|$.

In Fig.~\ref{fig:2} we show the relative deviation between the free-space and the finite distance region $t$-operator between the first and second potential resonance to quantify the quality of the approximation also for higher values of the angular momentum quantum number $\ell$.
We find that the relative deviation decreases with increasing $\ell$.
Also the expected behavior $\sim \mathrm{e}^{-r_\mathrm{b} \sqrt{-z^{2b}}} /( r_\mathrm{b} \sqrt{-z^{2b}}) $ in $q_z = \sqrt{-z^{2b}}$ can be observed.

\begin{figure}[htb]
	  \centering{}%
\includegraphics[width = 1
	  \columnwidth]{figure3}
\caption{Relative deviation of $t_\ell^\infty$ and $t_\ell^{r_\mathrm{b}}$ for $r_\mathrm{b} / r_\mathrm{sw} = 1000$ at $a / r_\mathrm{sw} = -1$ (full line), $a/ r_\mathrm{sw} = -10$ (dashed line) and $a / r_\mathrm{sw} = -100$ (dotted line) and $\ell = 0,1,2,3,4$ (black to light gray).}
\label{fig:2}
\end{figure}

\subsection{Mapped DVR}
\label{sec:IIC}

To compute the $t$-operator we use a mapped grid method in combination with a sine function basis as discussed in \cite{Willner:2004}.
We consider central interaction potentials with a $C_6/r^6$ van der Waals long range behavior. We focus on the Lennard-Jones van der Waals potential
\begin{equation}
V_\mathrm{LJ} = - \frac{C_6}{r^6}\left( 1 - \frac{\lambda^6}{r^6} \right) \, .
\end{equation}
The range of the van der Waals attraction is then defined as $r_\mathrm{vdW} = (m \, C_6 / \hbar^2)^{1/4} / 2 $ \cite{Chin:2010}.
For energies close to threshold the wave function will oscillate with a very big wavelength at large separation, where the potential almost vanishes and with short wavelength close to the potential minimum.
This wide range of wavelengths can increase the numerical cost when using standard grid representation methods.
Therefore we want to transform to a new relative coordinate $x$ in which the wave function oscillates with a more regular frequency.
In the following we drop the index in the interaction potential $V_\mathrm{LJ}$, since the procedure applies generally also to potentials of similar shape.
A coordinate transformation with the desired properties is given by \cite{Willner:2004}
\begin{equation}
x(r) = \frac{\sqrt{m}}{p_\mathrm{max}} \int_{r_\mathrm{in} = 0}^r \mathrm{d}r' \sqrt{E_\mathrm{max} - V^\mathrm{env}(r')} \, ,
\end{equation}
where $p_\mathrm{max}$ and $E_\mathrm{max}$ are specifying the transformation and with the enveloping potential $V^\mathrm{env}$ whose value is at all separation lower or equal than the one of the potential $ V $ and defined as
\begin{equation}
V^\mathrm{env}(r) = \mathrm{min}_{r' \geq r}\left(V (r') \right) \, .
\end{equation}
To see that the transformation has the desired properties, we look at the semi-classical approximation of the $s$-wave Schr\"odinger equation for the phase function $\varphi$ after transforming coordinates
\begin{align}
\partial_x \varphi & = \pm \frac{\sqrt{m [E - V(r(x))]}}{\hbar (\partial_r x(r))} \\
& =  \pm \frac{p_\mathrm{max} \sqrt{ E - V(r(x))}}{ \sqrt{ E_\mathrm{max} - V^\mathrm{env}(r(x))}} \nonumber
\end{align}
which is indeed approximately constant for $E \sim E_\mathrm{max}$ and $V^\mathrm{env}(r(x)) \sim V(r(x))$.
We absorb the change in volume element $\partial_x r(x) = J(x)$ into the wave function $\langle x | \psi \rangle = \bar{\psi} (x) = \sqrt{J(x)} \psi(r(x))$ and choose an equally spaced grid $x_i = i x(r_\mathrm{b})/N  $ in $x$ with $N+1$ grid points in the region $[0,x(r_\mathrm{b})]$ together with a sine or particle in a box base
\begin{equation}
s_k(x) = \sqrt{\frac{2}{N}}\mathrm{sin}(k \pi x / x(r_\mathrm{b}))
 \quad (k = 1, ... N-1)
 \end{equation}
on the grid.
This choice for the bases introduces a hard wall boundary condition in the system.
We define the transformation to the base of $x$-grid points
\begin{equation}
S_{jk} = s_k (x_j) = \langle x_j | s_k \rangle \, ,
\end{equation}
which leads us to the representation of the $x$-grid points $x_j$ in terms of the sine base
\begin{equation}
\tilde{s}_j(x) = \sum_k (S^\dag)_{jk}  s_k(x) \, .
\end{equation}
The radial part of the kinetic energy operator $T = -\hbar^2 \partial_r^2 / m$ can then be obtained in the $\tilde{s}_j$-base
\begin{align}
\langle \tilde{s}_i | T |\tilde{s}_j \rangle & = - \frac{\hbar^2}{m}\frac{\pi^2}{x(r_\mathrm{b})^2}\sum_{k = 0}^N J(x_i)^{-1/2} (D^\dag)_{ik} \nonumber \\
& \phantom{=} \qquad \times J(x_k)^{-1} D_{kj} J(x_j)^{-1/2}
\end{align}
with
\begin{equation}
D_{ij} = \frac{x(r_\mathrm{b})}{\pi} (\partial_x \tilde{s}_j)(x_i) \, .
\end{equation}

As a consequence we can obtain the radial Hamiltonian $H_\ell$ of the two-atom system as a finite dimensional matrix on the $x$-grid using the $\tilde{s}_j$-base functions
\begin{equation}
\langle \tilde{s}_i | H_\ell |\tilde{s}_j \rangle = \langle \tilde{s}_i | T |\tilde{s}_j \rangle + \delta_{ij} V_\ell (x_i) \, .
\end{equation}
Here we use the potential including the angular momentum barrier $ V_\ell (r) = V(r) + \hbar^2 \ell(\ell +1)/ m r^2$ and note that we introduce a cut-off in $V_\ell$ at large positive energy to prevent the potential to diverge at zero distance in numerical practice.
We can find the eigenvalues $E_\ell^i$ and eigenstates $| E_\ell^i \rangle$ of the resulting matrix and the two-body $t$-operator can be obtained
in momentum representation by
\begin{align}
t_\ell(z^{2b},k',k) &= \sum_i \langle \ell,k' | x_i \rangle V(x_i) \langle x_i | \ell, k \rangle \nonumber \\
& \phantom{=} + \sum_{ijn} \langle \ell, k' | x_i \rangle \frac{\langle x_i | V | E_\ell^n \rangle \langle E_\ell^n| V | x_j \rangle}{z - E_\ell^n} \nonumber \\
& \phantom{=} \qquad \times \langle x_j | \ell, k \rangle \, ,
\end{align}
with $\langle x_i | \ell,  k \rangle = [2 J(x_i) x(r_\mathrm{b})/\pi N]^{1/2} r(x_i) \mathrm{j}_\ell [k r(x_i)] $.

Diagonalizing $t_\ell(z - 3 q_\alpha^2/4m,k',k)$ for a discrete set of momenta $k$, $k'$
we can find the expansion
$\mathcal{T}_\alpha = \int dq q^{2}\sum_i \tau (q,i)| q, i \rangle_{\alpha \alpha} \langle q, i |$ that we need in the three-body calculation similar to \cite{Mestrom:2019b}.
In the following section we present three-body results for $V_\mathrm{LJ}$
that we obtained using the mapped DVR approach presented here.

\section{Three-body results}
\label{sec:III}
We performed three-body recombination calculations using the analytic expressions for $t$ we obtained for the square-well interaction
and the values for $t$ that we obtained numerically for the
Lennard-Jones potential with the mapped DVR method.
In Fig.~\ref{fig:3} we compare the three-body recombination rate for different values of $r_\mathrm{b}$ over a wide range of scattering lengths $a$ between the first and second potential resonance.
The results for the free-space case $r_\mathrm{b} = \infty $ have been performed using the Weinberg expansion for the $t$-operator \cite{Mestrom:2019b} and have been partly published earlier \cite{Mestrom:2019,Mestrom:2020}.
We find that the hard wall boundary condition has a minor influence on the three-body recombination rate when $|a|/r_\mathrm{b} \ll 1$.
For the square-well interaction with $r_\mathrm{b} = 2000 \, r_\mathrm{sw}$ in the regime $|a|/r_\mathrm{sw}< 1.5$ we find a small relative deviation in $K_3$ of $\sim 1 \%$.
For the Lennard-Jones potential with $r_\mathrm{b} = 2000 \, r_\mathrm{vdW}$ in the regime $|a|/r_\mathrm{vdW} < 1.5$ we find a still small but larger relative deviation of $\sim 4 \%$ (see Fig.~\ref{fig:5} in appendix \ref{app:C}).
We attribute this increase in relative deviation to the numerical error in the mapped DVR approach.
However, when $|a| \sim r_\mathrm{b}$ the deviation becomes more significant.
We note that the relative deviation of the analytic $t_0$ is also increasing significantly in this regime (see Fig.~\ref{fig:2}).

To quantify the accuracy of the finite distance approximation in the Efimov regime we investigate the ground and first excited Efimov resonance peaks and the first and second excited Efimov recombination minima for various values of $r_\mathrm{b}$.
We obtain the resonance peak positions $a_-^{(n)}$ and widths $\eta^{(n)}$ as well as the recombination minima $a_+^{(n)}$ by fitting the universal expressions \cite{Efimov:1979,Esry:1999,Bedaque:2000,Braaten:2001,Braaten:2002,Braaten:2004,braaten:2006,D_Incao:2018} for the recombination rate to our calculation close to the peak or minimum position.
We note that for fitting recombination minimum positions we introduce an overall scaling factor as an additional fit parameter to improve the quality of the fit.
The results can be found in Table \ref{tab:1} and \ref{tab:2} for different values of $r_\mathrm{b}$.
We find good agreement for both the square-well and Lennard-Jones potential in $a_-^{(0)}$ and $a_+^{(1)}$ with $r_\mathrm{b} = 2000 \, r_\mathrm{sw}$ or $r_\mathrm{b} = 2000 \, r_\mathrm{vdW}$, respectively, with a deviation of $\lesssim 1 \%$ in all cases.

\begin{table}[h!]
\centering
\caption{Width and position of the ground and first excited Efimov resonances along with first and second excited recombination minima for a square-well interaction between the first and second potential resonance with different sizes of the finite distance region.
Our results for $a_\pm^{(n)}$ are converged in at least three significant figures.
The results for $r_\mathrm{b} / r_\mathrm{sw} = \infty $ have been obtained using the Weinberg expansion for the $t$-operator and have been partly published earlier \cite{Mestrom:2019}. \label{tab:1}}
\begin{tabularx}{\columnwidth}{>{\raggedright \arraybackslash}X|>{\centering \arraybackslash}X >{\centering \arraybackslash}X| >{\centering \arraybackslash}X >{\centering \arraybackslash}X|>{\centering \arraybackslash}X|>{\centering \arraybackslash}X}
\hline
\hline
$\frac{r_\mathrm{b}}{r_\mathrm{sw}}$ & $\eta^{(0)}$ & $\frac{a_-^{(0)}}{r_\mathrm{sw}}$ & $\eta^{(1)}$ & $\frac{a_-^{(1)}}{r_\mathrm{sw}}$ & $\frac{a_+^{(1)}}{r_\mathrm{sw}}$ & $\frac{a_+^{(2)}}{r_\mathrm{sw}}$ \\
\hline
 $1000$ & $0.059$ & $-17.12$ & $0.062$ & $-279.4$ & $-$ & $323$ \\
 $2000$ & $0.060$ & $-17.27$ & $0.049$ & $-314.6$ & $9.58$ & $284$ \\
 $4000$ & $0.061$ & $-17.34$ & $-$ & $-$ & $-$ & $-$ \\
 $5000$ & $0.061$ & $-17.36$ & $0.058$ & $-346.6$ & $9.55$ & $263$ \\
 $\infty$  & $0.061$ & $-17.42$ & $0.067$ & $-371.9$ & $9.530$ & $249.7$ \\
\hline
\hline
\end{tabularx}
\end{table}

\begin{table}[h!]
\centering
\caption{Width and position of the ground and first excited Efimov resonances along with first and second excited recombination minima for a Lennard-Jones interaction between the first and second potential resonance with different sizes of the finite distance region. Our results for $a_\pm^{(n)}$ are converged in at least two significant figures.
The results for $r_\mathrm{b} / r_\mathrm{vdW} = \infty $ have been obtained using the Weinberg expansion for the $t$-operator and have been partly published earlier \cite{Mestrom:2020}. \label{tab:2}}
\begin{tabularx}{\columnwidth}{>{\raggedright \arraybackslash}X|>{\centering \arraybackslash}X >{\centering \arraybackslash}X| >{\centering \arraybackslash}X >{\centering \arraybackslash}X|>{\centering \arraybackslash}X|>{\centering \arraybackslash}X}
\hline
\hline
$\frac{r_\mathrm{b} }{r_\mathrm{vdW}} $ & $\eta^{(0)}$ & $\frac{a_-^{(0)}}{r_\mathrm{vdW}}$ & $\eta^{(1)}$ & $\frac{a_-^{(1)}}{r_\mathrm{vdW}}$ & $\frac{a_+^{(1)}}{r_\mathrm{vdW}}$ & $\frac{a_+^{(2)}}{r_\mathrm{vdW}}$ \\
\hline
 $500$ & $0.022$ & $-9.46$ & $-$ & $-$& $-$ & $-$ \\
 $1000$ & $0.022$ & $-9.55$ & $0.030$ & $-139$ & $27.9$ & $954$ \\
 $2000$ & $0.021$ & $-9.59$ & $0.032$ & $-150$ & $27.5$ & $968$ \\
 $\infty$  & $0.020$ & $-9.67$ & $0.038$ & $-163$ & $27.2$ & $722$ \\
 $\infty$ \cite{Mestrom:2017}  & $-$ & $-9.74$  & $-$ & $-164$  & $27.2$  & $-$ \\
\hline
\hline
\end{tabularx}
\end{table}

\begin{figure*}[htb!]
	  \centering
	      \includegraphics[width=\textwidth]{figure4}

\caption{\label{fig:3} Three-body recombination rates for  the square-well (top) and Lennard-Jones potential (bottom). We compare $K_3(0)$ for sizes of the finite distane region $r_\mathrm{b} / r_{\mathrm{sw}} $ or $r_\mathrm{b} / r_\mathrm{vdW} $ of $1000$ (green circles), $2000$ (blue squares), $5000$ (red diamonds),  $\infty$ (black line). The data for $r_\mathrm{b} = \infty$ are calculated using the Weinberg expansion method (see \cite{Mestrom:2019b} for more details) and have partly been published earlier \cite{Mestrom:2019,Mestrom:2020}.
}
\end{figure*}

\begin{figure*}[htb]
	  \centering
	\begin{minipage}[c]{1 \columnwidth}
	  \includegraphics[width = 1 \textwidth]{figure5}
	\end{minipage}
	\begin{minipage}[c]{1 \columnwidth}
	\includegraphics[width = 1 \textwidth]{figure6}
	\end{minipage}
\caption{ We show the dependence of the lowest Efimov resonance position $a_-^{(0)}$ and width $\eta^{(0)}$ on the integration range of the atom-dimer momentum $[0,q_\mathrm{max}]$ for a Lennard-Jones van der Waals potential with almost four $s$-wave bound states (black dots) and almost 6 $s$-wave bound states (red dots). The vertical lines indicate the positions of two-body bound states for almost four (gray) and 6 (red)
$s$-wave bound states when the system is tuned close to the lowest Efimov resonance. }
\label{fig:4}
\end{figure*}

The mapped DVR approach allows us to obtain $t$ even for deep potentials without any complications.
However, in the equations for three-body recombination, Eq.~(\ref{eq:3bodyint}), we have to limit to a finite integration range $[0, q_\mathrm{max}]$ in $q$ and to a finite expansion in relative atom-dimer partial waves $\ell$.
In the following we consider the lowest Efimov resonance of $V_\mathrm{LJ}$ close to the 4th and 6th potential resonance.
We include up to $\ell_\mathrm{max} = 20$ to guarantee the convergence in $\ell$
and analyze the convergence in $q_\mathrm{max}$.
The results are shown in Fig.~\ref{fig:4}.
We find that the convergence in $q_\mathrm{max}$ of $a_-^{(0)}$ is similar in both cases.
For a value of $a_-^{(0)}$ accurate up to $2\%$ it is sufficient to choose $q_\mathrm{max} \approx  20 \, \hbar / r_\mathrm{vdW}$.
For $q_\mathrm{max} \approx  40 \, \hbar / r_\mathrm{vdW}$ we find $a_-^{(0)}$ to be accurate up to $0.2\%$.
The resonance width $\eta^{(0)}$ is however converging more slowly in $q_\mathrm{max}$.
To obtain converged results we find that $q_\mathrm{max}$ needs to be chosen such that $3 q_\mathrm{max}^2 / 4 m$ is larger than the lowest two-body binding energy.
We note that when $3 q_\mathrm{max}^2 / 4 m$ is lower than some dimer binding energies recombination into those dimer states will be neglected.
The values of $a_-^{(0)}$ for the 4th and 6th potential resonance of a Lennard-Jones potential have been calculated earlier by an adiabatic hyperspherical coordinate approach in \cite{Wang:2012} to be $a_-^{(0)} /r_\mathrm{vdW} \approx -10.7$ and $a_-^{(0)}/r_\mathrm{vdW} \approx -10.4 $ respectively.
Comparing to our results we find a relative deviation of about  $\sim 10\%$.
We checked the convergence of our results in momentum grids, partial waves, separable expansion terms of the $t$-operator and the number of two-body $x$-grid points, but were not able to explain this deviation.
We hope to be able to resolve this discrepancy in future work.

\section{Conclusion and Outlook}
\label{sec:IV}

We extended a well established mapped DVR method to calculate off-shell scattering properties of a two-atom system that can be used directly as input for three-body scattering calculations in an AGS momentum space approach.
Using an analytic example we demonstrated that results from a method restricting to a finite distance region with a hard wall boundary condition such as the mapped DVR can be used to approximate the free space result on the two-body level.
We performed three-body recombination calculations using the mapped DVR and compared them with results calculated with a standard Weinberg expansion of the $t$-operator.
The corresponding results are in good agreement when the magnitude of the scattering length is much smaller than the size of the finite distance region in the mapped DVR scheme.
Using the mapped DVR method we are able to perform three-body calculations even for deep interaction potentials that are more difficult to access with standard momentum space treatments.
We find that applying a cut-off in the relative atom-dimer momentum for the three-body recombination equations can lead to accurate results within a few percent when determining the three-body parameter $a_-^{(0)}$.
However, this approximation is much less accurate for the Efimov resonance width $\eta^{(0)}$, since effects of recombination into dimer states beyond the integration range are neglected.

Our results for the atom-dimer cut-off momentum dependence suggest that the method can easily be generalized to deep realistic potentials without significant numerical complications.
Since mapped DVR methods have been applied to multichannel systems as well, the presented numerical procedure allows for a straightforward generalization to a multichannel version.
So far multichannel methods including the full spin structure of the atomic system just exist for effective non-local potentials with a few separable terms \cite{Li:2019}.
For local van der Waals potentials on the other hand the spin structure is usually just approximated by effective models to be able to perform the calculations.
The method we present here allows for calculating three-body scattering amplitudes involving local potentials including the full spin-structure of the atomic system.
In addition we find that the regime $|a|/r_\mathrm{vdW} \lesssim 1$ is particularly well represented in the mapped DVR approach.
This is promising since in this regime partial recombination rates have been determined experimentally \cite{Haerter:2013, Wolf:2017}.
Also the prospects for studying three-body elastic collisions around $|a| / r_\mathrm{vdW} \approx 0$ are good with the method presented here. The corresponding elastic cross sections affect
the phase diagram of a Bose-Einstein condensate, as has been demonstrated recently \cite{Zwerger:2019, Mestrom:2020}.

\section*{Acknowledgements}
We thank Denise Ahmed-Braun, Victor Colussi, Gijs Groeneveld, and Silvia Musolino for discussions.
This research is financially supported by the
Netherlands Organisation for Scientific Research (NWO)
under Grant No. 680-47-623 and by the Foundation for
Fundamental Research on Matter (FOM).

\clearpage