# Evolution of the unitary Bose gas for broad to narrow Feshbach resonances

D.J.M. Ahmed-Braun<EMAIL_ADDRESS>Eindhoven University of Technology,
P.O. Box 513, 5600 MB Eindhoven, The Netherlands  S. Musolino  Eindhoven
University of Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands  V.
E. Colussi  INO-CNR BEC Center and Dipartimento di Fisica, Università di
Trento, 38123 Povo, Italy  S. J. J. M. F. Kokkelmans  Eindhoven University of
Technology, P.O. Box 513, 5600 MB Eindhoven, The Netherlands

###### Abstract

We study the post-quench dynamics of unitary Bose gases using a two-channel
model, focusing on the effect of variations in the width of the Feshbach
resonance due to density changes. We generally find that increasing the
density leads to a corresponding increase in the production of closed channel
molecules, a decrease in the build up of quantum depletion and a transition
from linear to quadratic early-time growth of the two-body contact as well as
the condensed pair fraction. Motivated by the presence of closed-channel
molecules in the unitary regime, we study the embedded two-body problem
finding a transition from open to closed-channel dominated dimers due to many-
body effects.

###### pacs:

Valid PACS appear here

##  I Introduction

The magnetic tunability of Feshbach resonances Feshbach (1958) in ultracold
gases makes it possible to experimentally control the two-particle interaction
strength as characterized by the s-wave scattering length $a(B)$

$a(B)=a_{\mathrm{bg}}\left(1-\frac{\Delta B}{B-B_{0}}\right),$ (1)

with background scattering length $a_{\mathrm{bg}}$, magnetic resonance width
$\Delta B$ and resonance position $B_{0}$ Chin  _et al._ (2010). By
diabatically quenching the scattering length to unitarity
($\absolutevalue{a(B)}\rightarrow\infty$) Makotyn  _et al._ (2014); Klauss
_et al._ (2017); Eigen  _et al._ (2017, 2018) , it is possible to beat
disastrous per-particle losses in bosonic gases that scale as $n^{2}a^{4}$,
with atomic density $n$. At unitarity the density scales remain finite and
according to the “Universality Hypothesis” all properties of unitary quantum
gases should scale continuously with the Fermi scales
$k_{\mathrm{n}}=(6\pi^{2}n)^{1/3}$,
$E_{\mathrm{n}}=\hbar^{2}k_{\mathrm{n}}^{2}/(2m)$ and
$t_{\mathrm{n}}=\hbar/E_{\mathrm{n}}$, where $m$ is the atomic mass Ho (2004).
This universality and the associated scale invariance relates the properties
of ultracold gases at unitary to other seemingly unrelated strongly correlated
systems, such as the quark-gluon plasma and the inner crust of neutron stars
Schäfer and Teaney (2009); Braaten and Hammer (2003); Castin and Werner
(2012). However, the non-universal scales associated with the Feshbach
resonance in addition to the finite size of Efimov states in a unitary Bose
gases Efimov (1971, 1979); Klauss  _et al._ (2017); D’Incao  _et al._ (2018);
Colussi  _et al._ (2018a) can alter the universal scaling and time-dependence
of system properties.

To describe the interplay between resonance and density scales, the vacuum
classification of a Feshbach resonance must be revisited in the many-body
context Ho  _et al._ (2012). The magnetic resonance width is momentum
dependent and equal to

$\displaystyle\Delta B=\frac{\hbar^{2}k}{\delta\mu mR^{*}},$ (2)

where $\delta\mu$ is the difference in the magnetic moment between two free
atoms and the Feshbach molecule and $R^{*}$ is the resonance width parameter.
In the vacuum classification scheme, the momentum scale $k\sim
1/a_{\mathrm{bg}}$ in estimating $\Delta B$, such that
$\absolutevalue{R^{*}/a_{\mathrm{bg}}}\ll 1$ corresponds to a broad resonance.
In the many-body classification scheme on the other hand, the momentum $k$ is
set by $k_{\mathrm{n}}$ Ho  _et al._ (2012), such that $k_{n}R^{*}\ll 1$
corresponds to a broad resonance. Furthermore, whereas the vacuum
classification is fixed by the specifics of the resonance and the atomic
species, the many-body classification varies with the density of the gas.

In this work, using the many-body classification of the resonance width, we go
beyond the single-channel models of Refs. Musolino  _et al._ (2019); Colussi
_et al._ (2018b); Sykes  _et al._ (2014); Colussi  _et al._ (2018a); D’Incao
_et al._ (2018); Kain and Ling (2014); Corson and Bohn (2015); Ancilotto  _et
al._ (2015); Yin and Radzihovsky (2016); Muñoz de las Heras  _et al._ (2019);
Gao  _et al._ (2020); Colussi  _et al._ (2020); Musolino  _et al._ (2022) and
study the post-quench evolution of the unitary Bose gas using a two-channel
model that explicitly contains the multichannel nature of the Feshbach
resonance. The one and two-body correlation dynamics in the system are
modelled using the cumulant method Musolino  _et al._ (2019); Colussi  _et
al._ (2018b). The correlation dynamics are studied over a range of
$k_{\mathrm{n}}R^{*}$, where we characterize the onset of non-universal
effects in the dynamics of the two-body contact, the atomic, pair and
molecular condensate fractions and the quantum depletion. To describe the pair
formation, the single-channel formalism of Ref. Musolino  _et al._ (2022) is
generalized to the two-channel model. Here, the presence of the molecular
fraction in the two-channel system provides us with an additional probe of the
many-body state at unitarity by using molecular spectroscopy Partridge  _et
al._ (2005). To explain the presence of molecules in the unitary regime, we
study the impact of many-body effects on the embedded few-body problem. This
problem is a generalization of the single-channel formalism of Ref. Colussi
_et al._ (2018b), which considers both the dressing due to the channel
couplings and the impact of the quantum statistics of the medium in the spirit
of the Cooper pair problem Zwerger (2011); Nozières and Schmitt-Rink (1985);
Holland  _et al._ (2001a). We observe that, for all considered values of
$k_{\mathrm{n}}R^{*}$, the embedded dimers become the dominant contribution to
the quantum depletion as the system evolves in the unitary regime.

The paper is outlined as follows. In Sec. II.1 we review the two-channel model
of a Feshbach resonance, which provides a foundation for the many-body model
that follows in Sec. II.2. We then apply the cumulant expansion to include up
to two-body correlations in the dynamics. Within this formalism we derive the
expressions for the pair condensate fraction and the dynamical two-body
contact. Next, we proceed with the formalism of the embedded two body problem.
These concepts are subsequently used to understand the results presented in
Sec. III and to form the conclusions in Sec. IV.

##  II  Two-channel model

In this paper we aim to study the onset of non-universal effects over a range
of resonance widths $k_{\mathrm{n}}R^{*}$. To this extend we develop a two-
channel model that captures the multichannel nature of Feshbach resonances. We
start by reviewing the two-channel model in vacuum in Sec. II.1 and then solve
the many-body version in Sec. II.2 using the method of cumulants. Here we
derive the two-channel generalizations of the pair condensate fraction and the
two-body contact. In Sec. II.3 we include many-body effects in the vacuum two-
channel model in an effort to investigate how the few-body physics is altered
by the medium.

###  II.1  The two-body problem

We employ a two-channel model in order to describe the effect of the resonance
width, classified in terms of $k_{\mathrm{n}}R^{*}$, on the dynamics of the
quenched unitary Bose gas. To that end we begin with a brief review of
Feshbach resonance theory Feshbach (1993) which serves as a basis for the
generalizations to the many-body context in what follows.

Considering an energetically open channel in the subspace $\mathcal{P}$
coupled to an energetically closed channel in the subspace $\mathcal{Q}$ with
a bound state with energy $\nu$ as sketched in Fig. 1, we split the coupled-
channel Schrödinger equation into two components Kokkelmans  _et al._ (2002),
such that

$E\begin{bmatrix}\ket{\Psi_{\mathrm{P}}}\\\
\ket{\Psi_{\mathrm{Q}}}\end{bmatrix}=\begin{bmatrix}\hat{H}_{\mathrm{PP}}&\hat{H}_{\mathrm{PQ}}\\\
\hat{H}_{\mathrm{QP}}&\hat{H}_{\mathrm{QQ}}\end{bmatrix}\begin{bmatrix}\ket{\Psi_{\mathrm{P}}}\\\
\ket{\Psi_{\mathrm{Q}}}\end{bmatrix},$ (3)

where we have employed the open- and closed-channel projection operators
$\hat{P}$ and $\hat{Q}$ that project the total scattering wave function $\Psi$
and the Hamiltonian $\hat{H}=\hat{H}^{0}+\hat{V}$ onto the open- and closed-
channel subspace respectively, such that
$\ket{\Psi_{\mathrm{P}}}\equiv\hat{P}\ket{\Psi}$,$\ket{\Psi_{\mathrm{Q}}}\equiv\hat{Q}\ket{\Psi}$,
$\hat{H}_{\mathrm{PQ}}\equiv\hat{P}\hat{H}\hat{Q}$ etc.

Figure 1: Schematic of the basic two-channel model. A bound state with energy
$\nu$ is in the closed channel subspace $\mathcal{Q}$ (black curve). Atoms
enter in the open channel subspace $\mathcal{P}$ (red curve) with incident
energy $E$ and couple to the closed channel bound state. By varying the
applied magnetic field, the bound state can be tuned to resonance.

Assuming that the energy $\nu$ of the closed channel bound state is close to
the collision energy $E$ of the interacting particles in the subspace
$\mathcal{P}$, we make the single resonance approximation and neglect the
scattering states and other bound states in the closed channel subspace
$\mathcal{Q}$. Under this approximation, we can solve Eq. (3) for the closed-
channel wave function $\ket{\Psi_{\mathrm{Q}}}$, finding

$\ket{\Psi_{\mathrm{Q}}}=\frac{\ket{\phi}\bra{\phi}}{E-\nu}\hat{H}_{\mathrm{QP}}\ket{\Psi_{\mathrm{P}}},$
(4)

where $\ket{\phi}$ is the closed channel bare bound state. Substituting the
previous relation into Eq. (3), we eliminate $\ket{\Psi_{\mathrm{Q}}}$ from
the expression for the open channel scattering wave function
$\ket{\Psi_{\mathrm{P}}}$ such that

$E\ket{\Psi_{\mathrm{P}}}=(\hat{H}^{0}_{\mathrm{PP}}+\hat{V}_{\mathrm{eff}})\ket{\Psi_{\mathrm{P}}},$
(5)

with

$\hat{V}_{\mathrm{eff}}=\hat{V}_{\mathrm{PP}}+\hat{H}_{\mathrm{PQ}}\frac{\ket{\phi}\bra{\phi}}{E-\nu}\hat{H}_{\mathrm{QP}}.$
(6)

Equation (5) can be solved straightforwardly, yielding

$\ket{\Psi_{\mathrm{P}}}=\ket{\psi_{\mathrm{P}}^{+}}+\frac{1}{E^{+}-\hat{H}_{\mathrm{PP}}}\frac{\hat{H}_{\mathrm{PQ}}\ket{\phi}\bra{\phi}\hat{H}_{\mathrm{QP}}}{E-\nu}\ket{\Psi_{\mathrm{P}}},$
(7)

where $\ket{\psi_{\mathrm{P}}^{+}}$ is the eigenstate of the direct open-
channel interaction Hamiltonian $\hat{H}_{\mathrm{PP}}$ and $E^{+}=E+i\delta$,
with $\delta\to 0^{+}$ to avoid singularities.

Throughout this work, we consider separable potentials

$\displaystyle\hat{V}_{\mathrm{PP}}$ $\displaystyle=v\ket{\zeta}\bra{\zeta},$
(8a) $\displaystyle\hat{V}_{\mathrm{PQ}}$
$\displaystyle=\beta\ket{\zeta}\bra{\zeta},$ (8b)

with form factors $\ket{\zeta}$ and open-channel and coupling potential
strengths $v$ and $\beta$ respectively. Analogous to Ref. Kokkelmans  _et al._
(2002), we define the closed channel amplitude
$\Phi_{Q}\equiv\braket{\phi}{\Psi_{\mathrm{Q}}}/\sqrt{2}$ for the system to be
in the bound state. Using Eqs. (4)-(8b), Eq. (3) takes the form

$\displaystyle E\Psi_{\mathrm{P}}(\mathbf{k})$
$\displaystyle=\frac{\hbar^{2}k^{2}}{m}\Psi_{\mathrm{P}}(\mathbf{k})+v\zeta(2\mathbf{k})\sum_{\mathbf{q}}\zeta^{*}(2\mathbf{q})\Psi_{\mathrm{P}}(\mathbf{q})$
$\displaystyle\quad+g\zeta(2\mathbf{k})\Phi_{\mathrm{Q}},$ (9a)
$\displaystyle\vspace{0.1cm}E\Phi_{\mathrm{Q}}$
$\displaystyle=\nu\Phi_{\mathrm{Q}}+\frac{g}{2}\sum_{\mathbf{k}}\Psi_{\mathrm{P}}(\mathbf{k})\zeta^{*}(2\mathbf{k}),$
(9b)

with potential interaction strength $g=\sqrt{2}\beta\braket{\phi}{\zeta}$.
Following Refs. Musolino  _et al._ (2019); Colussi  _et al._ (2018b, 2020);
Musolino  _et al._ (2022), we choose a step-function form factor
$\braket{\mathbf{k}/2}{\zeta}=\zeta(\mathbf{k})=\Theta(\Lambda-\absolutevalue{\mathbf{k}}/2)$,
with relative two-body momentum $\mathbf{k}$ and momentum cut-off $\Lambda$.
The value of the cut-off is calibrated by matching the dimer binding energy of
the two-channel model to the full coupled-channel dimer binding energy tho
111For instance, in the case of the $B_{0}=155.04$ G resonance in
${}^{85}\text{Rb}$ Claussen  _et al._ (2003), we find
$\Lambda=0.6/r_{\mathrm{vdW}}$, with $r_{\mathrm{vdW}}$ the van der Waals
length of ${}^{85}\text{Rb}$ Chin  _et al._ (2010). We obtained the momentum
cut-off $\Lambda$ by comparing the coupled-channels data from tho  to a
polynomial where we keep the zeroth order term in the low-energy expansion of
$a(B)^{-1}$. Other resonances in different atomic species result in different
calibrated values of the cut-off. For example, the calibration of the
$B_{0}=402$ G resonance in 39K yields $\Lambda=0.46/r_{\mathrm{vdW}}$. The
calibrated value depends on the width of the considered resonance and
approaches the value $\Lambda=2\pi/\bar{a}$ Colussi  _et al._ (2018b);
Musolino  _et al._ (2019) with mean-scattering length $\bar{a}=0.955r_{vdw}$
in the single channel limit.. The potential strengths $v$ and $g$ can be
modelled using the following set of renormalization equations Kokkelmans and
Holland (2002); Snyder  _et al._ (2012)

$v=v_{0}\Gamma=\frac{4\pi\hbar^{2}a_{\mathrm{bg}}}{m}\Gamma,$ (10)
$g=g_{0}\Gamma=\frac{(2\pi)^{3/2}\hbar^{2}}{\sqrt{R^{*}}\pi m}\Gamma,$ (11)

with $\Gamma=\left(1-2\Lambda a_{\mathrm{bg}}/\pi\right)^{-1}$ and with
renormalized interaction parameters $v$ and $g$. In order to investigate the
effect of the separable potential interaction on the open channel wave
function $\ket{\Psi_{\mathrm{P}}}$, we multiply Eq. (7) from the left by
$\bra{\mathbf{k}}\hat{V}_{\mathrm{eff}}$, where $\ket{\mathbf{k}}$ is an
unscattered state. Exploiting the relation between the transition operator
$\hat{T}$ and the scattering potential operator
$\hat{T}\ket{\mathbf{k}}=\hat{V}_{\mathrm{eff}}\ket{\Psi_{\mathrm{P}}}$, we
can find the following expression for the coupled channels two-body transition
matrix in the case of a separable potential interaction

$T(E)=v\zeta(2\mathbf{k})\braket{\zeta}{\psi_{\mathrm{P}}^{+}}+\frac{\frac{g^{2}}{2}\absolutevalue{\braket{\zeta}{\psi_{\mathrm{P}}^{+}}}^{2}}{E-\nu-\frac{g^{2}}{2}\braket{\zeta}{\hat{G}_{\mathrm{P}}(E)}{\zeta}},$
(12)

where we have introduced the open channel Green’s operator
$\hat{G}_{\mathrm{P}}(E)=(E-\hat{H}_{0}-v\ket{\zeta}\bra{\zeta})^{-1}$. We can
use Eq. (12) in order to extract important scattering parameters and obtain
the bound state energies from its poles, which we do presently.

####  II.1.1  Vacuum classification of the resonance width

As mentioned in Sec. I, the description of narrower resonances, which in the
vacuum classification corresponds to the limit
$\absolutevalue{R^{*}/a_{\mathrm{bg}}}\gg 1$, requires the inclusion of finite
range scales. These finite range scales affect the universal scaling of the
system, which is clear from the effective range expansion Taylor (2012)

$\displaystyle
k\cot[\delta_{0}(k)]\approx\frac{-1}{a}+\frac{1}{2}k^{2}R_{\mathrm{eff}}+\mathcal{O}(k^{4}),$
(13)

where $\delta_{0}(k)$ is the s-wave phase shift which depends on the
scattering length $a$ and the effective range $R_{\mathrm{eff}}$. The
necessity to consider the energy-dependent correction to the phase shift
follows directly from the inspection of the Breit-Wigner form of the phase
shift for a non-resonant open-channel interaction, which can be expressed as
Ho  _et al._ (2012)

$\text{tan}[\delta_{0}(k)]=-ka_{\mathrm{bg}}-\frac{\Delta E}{E-\nu_{0}},$ (14)

where $\Delta E=\hbar^{2}/ma_{\mathrm{bg}}R^{*}$ is the energy width on the
vacuum length scale $a_{\mathrm{bg}}$ and where $\nu_{0}=\delta\mu(B-B_{0})$
is the energy detuning from the Feshbach resonance. Analogous to Ref.
Kokkelmans  _et al._ (2002), the detuning $\nu_{0}$ can be related to the
closed channel bound state energy $\nu$ as Kokkelmans and Holland (2002);
Snyder  _et al._ (2012)

$\displaystyle\nu=\nu_{0}+\frac{m\Lambda}{4\pi^{2}\hbar^{2}}gg_{0}.$ (15)

Multiplying the numerator and denominator of the second term on the right-hand
side of Eq. (14) by a factor $mR^{*}/\hbar^{2}$, we can recognize that the
term that scales with the square of the momentum and sets the energy scaling
depends on the size of the resonance strength parameter $R^{*}$ 222By applying
Eq. (2), we find that Eq. (14) can be rewritten as
$\text{tan}[\delta_{0}(k)]=-ka_{\mathrm{bg}}-\left(kR^{*}-\frac{B-B_{0}}{\Delta
B}\right)^{-1}$..

In order to find the exact relation that connects the effective range to the
resonance strength parameter $R^{*}$, we use the transition matrix as
presented in Eq. (12) and relate it to the phase shift $\delta_{0}(E)$ as

$k\text{cot}[\delta_{0}(E)]=ik-\frac{4\pi\hbar^{2}}{mT(E)}.$ (16)

In the case of step function separable potential interactions the expansion of
Eq. (16) around $k\rightarrow 0$ and its comparison to Eq. (13) yield the
following expression for the effective range

$R_{\mathrm{eff}}=-2R^{*}+\frac{4}{\pi\Lambda}+\frac{1}{a}\left(4a_{\mathrm{bg}}R^{*}-\frac{2a_{\mathrm{bg}}^{2}R^{*}}{a}\right),$
(17)

which reduces to $R_{\mathrm{eff}}\approx-2R^{*}$ at unitarity and for large
values of the momentum-space cut-off Törmä and Sengstock (2014).

####  II.1.2  Dimer wave function

In the two-channel model, the bound state is dressed by the channel coupling.
The closed-channel amplitude of this dressed bound state is quantified by the
dimer wave function normalization factor, or $Z$-parameter. This parameter can
be introduced consistently with Refs. Duine and Stoof (2004); Cohen-Tannoudji
and Guéry-Odelin (2011) as

$\begin{bmatrix}\ket{\Psi_{\mathrm{P}}}\\\
\ket{\Psi_{\mathrm{Q}}}\end{bmatrix}=\sqrt{Z}\begin{bmatrix}\hat{G}_{\mathrm{P}}(E_{\mathrm{D}})\beta\ket{\zeta}\braket{\zeta}{\phi}\\\
\ket{\phi}\end{bmatrix},$ (18)

such that $\braket{\Psi}{\Psi}=1$ and where $E_{\mathrm{D}}$ is the bound
state energy of the dressed dimer state. In the case of the separable
potential, the $Z$ parameter takes the form

$Z=\left[1+\frac{\frac{g^{2}}{2}\sum_{\mathbf{k}}\frac{\zeta(2\mathbf{k})}{(E_{\mathrm{D}}-\frac{\hbar^{2}k^{2}}{m})^{2}}}{\left(1-v\sum_{\mathbf{k}}\frac{\zeta(2\mathbf{k})}{E_{\mathrm{D}}-\frac{\hbar^{2}k^{2}}{m}}\right)^{2}}\right]^{-1}.$
(19)

The inspection of the previous equation reveals that at unitarity, where
$E_{\mathrm{D}}=0$, $Z=0$ Falco and Stoof (2007); Partridge  _et al._ (2005).
Departing from unitarity, the $Z$-parameter is bounded to a maximum value of
one, which corresponds to the dimer fully in the closed channel subspace. In
Sec. II.3.2, the definition of the $Z$-parameter will be extended to include
many-body effects of the background gas.

###  II.2  The many-body problem

We now proceed to model a uniform gas of identical bosons interacting via the
separable cut-off potential interactions as introduced in Sec. II.1. Using
second quantization, the two-channel many-body Hamiltonian corresponds to

$\displaystyle\hat{H}=$
$\displaystyle\sum_{\mathbf{k}}\frac{\hbar^{2}k^{2}}{2m}\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{k}}+\sum_{\mathbf{k}}\left(\frac{\hbar^{2}k^{2}}{4m}+\nu\right)\hat{b}^{\dagger}_{\mathbf{k}}\hat{b}_{\mathbf{k}}$
(20)
$\displaystyle+\frac{v}{2}\sum_{\mathbf{k},\mathbf{k^{\prime}},\mathbf{q}}\zeta\left(\mathbf{k}-\mathbf{k^{\prime}}+2\mathbf{q}\right)\zeta^{*}\left(\mathbf{k}-\mathbf{k^{\prime}}\right)\hat{a}^{\dagger}_{\mathbf{k}+\mathbf{q}}\hat{a}^{\dagger}_{\mathbf{k^{\prime}}-\mathbf{q}}\hat{a}_{\mathbf{k}^{\prime}}\hat{a}_{\mathbf{k}}$
$\displaystyle+\frac{g}{2}\sum_{\mathbf{k},\mathbf{q}}\left[\zeta^{*}(\mathbf{2k})\hat{b}^{\dagger}_{2\mathbf{q}}\hat{a}_{\mathbf{q}+\mathbf{k}}\hat{a}_{\mathbf{q}-\mathbf{k}}+\text{H.c}\right],$

where we have neglected the interactions among closed-channel molecules
Timmermans  _et al._ (1999); van Abeelen and Verhaar (1999); Holland  _et al._
(2001b); Kokkelmans  _et al._ (2001, 2002). We have introduced the open
channel atomic operators $\hat{a}_{\mathbf{k}}$ and the closed channel
molecular operators $\hat{b}_{\mathbf{k}}$, which can be expressed as Musolino
_et al._ (2019); Altman and Vishwanath (2005)

$\displaystyle\hat{b}^{\dagger}_{\mathbf{q}}=\sum_{\mathbf{k}}\frac{\phi(k)}{\sqrt{2}}\hat{c}^{\dagger}_{-\mathbf{k}+\mathbf{q}/2}\hat{c}^{\dagger}_{\mathbf{k}+\mathbf{q}/2},$
(21)

where $\phi(k)$ is the wave number projection of the closed channel bound
state $\ket{\phi}$ as introduced in Eq. (4) and where
$\hat{c}^{\dagger}_{\mathbf{k}}$ is a closed-channel atomic operator. We note
that we are working in the limit where $\phi(k)$ is localized with respect to
the density scales, such that $\hat{b}^{\dagger}_{\mathbf{q}}$ is a true
bosonic molecular operator, contrary to the composite pair operator introduced
in Ref. Musolino  _et al._ (2022).

We apply the Bogoliubov decomposition to the operators such that
$\braket{a_{0}}=\psi_{a}$ ($\braket{b_{0}}=\psi_{m}$) and
$\braket{a_{\mathbf{k}\neq 0}=0}$ ($\braket{b_{\mathbf{k}\neq 0}=0}$). Here
$\psi_{a}$ ($\psi_{m}$) is the atomic (molecular) wave function. We neglect
excitations of the molecular condensate Timmermans  _et al._ (1999); van
Abeelen and Verhaar (1999); Holland  _et al._ (2001b); Kokkelmans  _et al._
(2001, 2002).

To model correlations in the many-body system, we perform a cumulant expansion
Kira (2015); Colussi  _et al._ (2018b); Musolino  _et al._ (2019). This
expansion separates clusters of correlated atoms and molecules within the
many-body model. We truncate the cumulant expansion at the second order, such
that we obtain the doublet model where only single- and two-particle
correlations are considered. These correlations (or clusters) consist of the
singlets $\psi_{a}$ and $\psi_{m}$, which represent the atomic and molecular
condensate respectively, and the doublets
$\rho_{\mathbf{k}}\equiv\braket{\hat{a}^{\dagger}_{\mathbf{k}}\hat{a}_{\mathbf{k}}}$
and
$\kappa_{\mathbf{k}}\equiv\braket{\hat{a}_{-\mathbf{k}}\hat{a}_{\mathbf{k}}}$,
which are the one-body density and pairing matrices for $k\neq 0$ respectively
Blaizot and Ripka (1986).

We limit the evolution of the initially pure atomic condensate to a time up to
$t=2t_{\mathrm{n}}$. At longer times $\rho_{\mathbf{k}}$ begins to exceed
unity and the exclusion of strongly-driven higher-order cumulants can no
longer be justified Colussi  _et al._ (2020, 2018b); Musolino  _et al._
(2019).

Applying the cumulant model up to the doublet level and implementing the
Heisenberg equation of motion
$i\hbar\,d\mathcal{\hat{O}}/dt=[\hat{\mathcal{O}},\hat{H}]$, we obtain the
two-channel Hartree-Fock Bogoliubov (HFB) equations of motion Blaizot and
Ripka (1986)

$\displaystyle i\hbar\dot{\psi}_{a}$ $\displaystyle=$ $\displaystyle
v\left(|\zeta(0)|^{2}|\psi_{a}|^{2}+2\sum_{\mathbf{k}\neq
0}|\zeta(\mathbf{k})|^{2}\rho_{\mathbf{k}}\right)\psi_{a}$ (22)
$\displaystyle+$ $\displaystyle v\psi_{a}^{\ast}\sum_{\mathbf{k}\neq
0}\zeta(0)\zeta^{\ast}(2\mathbf{k})\kappa_{\mathbf{k}}+g\zeta(0)\psi_{m}\psi_{a}^{*},$
$\displaystyle i\hbar\dot{\psi}_{m}$ $\displaystyle=$
$\displaystyle\nu\psi_{m}+\frac{g}{2}\psi_{a}^{2}\zeta^{*}(0)+\frac{g}{2}\sum_{\mathbf{k}\neq
0}\zeta^{*}(2\mathbf{k})\kappa_{\mathbf{k}},$ (23)
$\displaystyle\hbar\dot{\rho}_{\mathbf{k}}$ $\displaystyle=$ $\displaystyle
2\mbox{Im}\left[\Delta_{\mathbf{k}}\kappa_{\mathbf{k}}^{\ast}\right],$ (24)
$\displaystyle i\hbar\dot{\kappa}_{\mathbf{k}}$ $\displaystyle=$
$\displaystyle
2h_{\mathbf{k}}\kappa_{\mathbf{k}}+\left(1+2\rho_{\mathbf{k}}\right)\Delta_{\mathbf{k}},$
(25)

where

$h_{\mathbf{k}}=\frac{\hbar^{2}k^{2}}{2m}+2v\left(|\zeta(\mathbf{k})|^{2}|\psi_{a}|^{2}+\sum_{\mathbf{q}\neq
0}|\zeta(\mathbf{k}-\mathbf{q})|^{2}\rho_{\mathbf{q}}\right),$ (26)

and

$\Delta_{\mathbf{k}}=v\zeta(2\mathbf{k})\left(\zeta^{\ast}(0)\psi_{a}^{2}+\sum_{\mathbf{q}\neq
0}\zeta^{\ast}(2\mathbf{q})\kappa_{\mathbf{q}}\right)+g\zeta(2\mathbf{k})\psi_{m},$
(27)

are the Hartree-Fock Hamiltonian and the pairing field, respectively. The HFB
equations of motion conserve the total number of atoms $N=N_{op}+N_{cl}$,
where $N_{op}$ consists of $N_{a}=V\absolutevalue{\psi_{a}}^{2}$ condensate
atoms, with system volume $V=N/n$, and
$N_{exc}=V\sum_{\mathbf{k}}\rho_{\mathbf{k}}$ excitations and where $N_{cl}$
consists of $2N_{m}$ atoms, with $N_{m}$ the number of closed channel
molecules $N_{m}=V\absolutevalue{\psi_{m}}^{2}$.

####  II.2.1  Many-body classification of the resonance width

As previously outlined in Sec. I, the density can be used to derive a set of
Fermi scales that quantify the properties of the unitary Bose gas. Exploiting
the Fermi energy, we can then obtain the many-body analogue of the resonance
width classification in terms of $k_{\mathrm{n}}R^{*}$ Ho  _et al._ (2012).
Similarly to the two-body classification scheme discussed in Sec. II.1.1, the
energy width $\Delta E=\hbar^{2}k_{\mathrm{n}}/mR^{*}$ on the Fermi momentum
scale $k_{\mathrm{n}}$ quantifies how fast the phase-shift varies as a
function of the typical many-body energy. We observe that for
$k_{\mathrm{n}}R^{*}\ll 1$, the phase shift varies slowly over a large range
of energies, such that we can neglect the energy dependent term in Eq. (13).
This corresponds to the case of a broad resonance. Instead, for a rapid
variation, where $k_{\mathrm{n}}R^{*}\gg 1$, we cannot neglect the energy
dependent term in Eq. (13) and we quantify the resonance to be narrow Ho  _et
al._ (2012). In addition to examining the energy width, we can define the
lifetime of the closed-channel molecules $\tau$ relative to the Fermi
timescale $t_{\mathrm{n}}$ as

$\displaystyle\frac{\tau}{t_{\mathrm{n}}}=\frac{E_{\mathrm{n}}}{\Delta
E}=k_{\mathrm{n}}R^{*},$ (28)

where the linear scaling with $k_{\mathrm{n}}R^{*}$ implies longer lifetimes
for narrower resonances in the many-body classification scheme.

As pointed out in Ref. Ho  _et al._ (2012), the vacuum and many-body
classifications of the resonance width are not equivalent for all values of
$a_{\mathrm{bg}}$ and $k_{\mathrm{n}}$. We can still satisfy the vacuum
classification $\absolutevalue{R^{*}/a_{\mathrm{bg}}}\ll 1$ of a broad-
resonance as introduced in Sec. II.1.1 whilst observing the narrow-resonance
behavior quantified by $k_{\mathrm{n}}R^{*}\gg 1$ for sufficiently large gas
densities.

####  II.2.2  Pair condensation

Extending the analysis performed in Ref. Musolino  _et al._ (2022) to the two-
channel model, the closed-channel molecular fraction $\psi_{m}$ and the
presence of the non-zero pairing matrix $\kappa_{\mathbf{k}}$ introduced in
Sec. II.2 signal off-diagonal long-range ordering (ODLRO) and pair
condensation. Interestingly, these pairs could be investigated experimentally
by using a rapid quenching procedure towards the weakly interacting regime
along the lines of Refs. Hodby  _et al._ (2005); Musolino  _et al._ (2019),
which maps the condensed pairs onto true molecules. Isolating the atomic
condensate from the fluctuations in order to omit the ODLRO that arises
trivially due to the presence of the atomic condensate Kira (2014); Yang
(1962), we define the following two-body density matrix as

$\displaystyle\underline{\underline{\rho}}^{(2)}(\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime},\mathbf{r}_{1},\mathbf{r}_{2})=\braket{\delta\hat{\underline{\psi}}^{\dagger}(\mathbf{r}_{1}^{\prime})\delta\hat{\underline{\psi}}^{\dagger}(\mathbf{r}_{2}^{\prime})\delta\hat{\underline{\psi}}(\mathbf{r}_{1})\delta\hat{\underline{\psi}}(\mathbf{r}_{2})},$
(29)

where $\delta\hat{\underline{\psi}}$ is the two-channel vector containing the
open and closed channel fluctuations 333We can generally define
$\delta\hat{\psi}$ as
$\begin{bmatrix}\psi_{P}-\psi_{a},&\psi_{Q}-\psi_{c}\end{bmatrix}^{\text{T}}$,
with $\psi_{c}$ a closed channel atomic condensate. Next, by using the single
resonance approximation and applying Eq. (21) we retain only a small portion
of the closed channel part and rewrite
$\braket{\delta\hat{\psi}_{Q}\delta\hat{\psi}_{Q}(\mathbf{r})}=\sqrt{2}\phi(r)\braket{\hat{b}_{0}}$..
Generalizing Refs. Yang (1962); Musolino  _et al._ (2022) we spectrally
decompose Eq. (29), such that

$\displaystyle\underline{\underline{\rho}}^{(2)}(\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime},\mathbf{r}_{1},\mathbf{r}_{2};t)=$
$\displaystyle\sum_{\nu}N^{(2)}_{\nu}(t)\underline{\varphi}_{\nu}^{(2)}(\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime},t)\underline{\varphi}_{\nu}^{(2)*}(\mathbf{r}_{1},\mathbf{r}_{2},t),$
(30)

where $\underline{\varphi}_{\nu}^{(2)}$ represent the two-component
orthonormal eigenvectors with eigenvalues $N^{(2)}_{\nu}$. In the case of
ODLRO, we expect
$\underline{\underline{\rho}}^{(2)}(\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime},\mathbf{r}_{1},\mathbf{r}_{2};t)$
to be non-zero in the long-range limit where
$\absolutevalue{\sum_{i=1}^{2}\mathbf{r}_{i}-\mathbf{r}_{i}^{\prime}}/2\rightarrow\infty$.
Applying the second order cumulant expansion as presented in Sec. II.2 to the
LR limit of Eq. (29), only the eigenstates with the anomalous contractions
$\braket{\delta\hat{\underline{\psi}}^{\dagger}(\mathbf{r}_{1}^{\prime})\delta\hat{\underline{\psi}}^{\dagger}(\mathbf{r}_{2}^{\prime})}\braket{\delta\hat{\underline{\psi}}(\mathbf{r}_{1})\delta\hat{\underline{\psi}}(\mathbf{r}_{2})}$
remain, since the Hartree and the Fock terms, defined as
$\braket{\delta\hat{\underline{\psi}}^{\dagger}(\mathbf{r}_{1}^{\prime})\delta\hat{\underline{\psi}}(\mathbf{r}_{1})}\braket{\delta\hat{\underline{\psi}}^{\dagger}(\mathbf{r}_{2}^{\prime})\delta\hat{\underline{\psi}}(\mathbf{r}_{2})}$
and
$\braket{\delta\hat{\underline{\psi}}^{\dagger}(\mathbf{r}_{1}^{\prime})\delta\hat{\underline{\psi}}(\mathbf{r}_{2})}\braket{\delta\hat{\underline{\psi}}^{\dagger}(\mathbf{r}_{2}^{\prime})\delta\hat{\underline{\psi}}(\mathbf{r}_{1})}$
respectively, typically vanish for separations that exceed a few Fermi lengths
Leggett (2008). Consequently, we find that

$\displaystyle\underline{\underline{\rho}}^{(2)}(\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime},\mathbf{r}_{1},\mathbf{r}_{2};t)\mathrel{\mathop{=}\limits_{\text{long
range}}}$ $\displaystyle
N^{(2)}_{0}(t)\underline{\varphi}_{0}^{(2)}(\mathbf{r}_{1}^{\prime},\mathbf{r}_{2}^{\prime},t)\underline{\varphi}^{(2)*}_{0}(\mathbf{r}_{1},\mathbf{r}_{2},t),$
(31)

with the associated pair wave function

$\displaystyle\underline{\varphi}_{0}^{(2)}(\mathbf{r},t)=\frac{1}{\sqrt{N_{0}^{(2)}(t)}}\begin{bmatrix}\sum_{\mathbf{k}}e^{i\mathbf{k}\cdot\mathbf{r}}\kappa_{\mathbf{k}}\\\
\sqrt{2}\phi(\mathbf{r})\psi_{m}\end{bmatrix},$ (32)

and the macroscopic eigenvalue

$\displaystyle
N_{0}^{(2)}(t)=\sum_{\mathbf{k}}\absolutevalue{\kappa_{\mathbf{k}}}^{2}+2\absolutevalue{\psi_{m}}^{2}.$
(33)

Contrary to the atomic condensate fraction that can be derived from the one-
body density matrix, the eigenvalue $N_{0}^{(2)}$ cannot be directly related
to the pair condensate fraction. This is a consequence of the background gas
of excitations that Bose-enhance interactions and thereby violate the bosonic
commutation relations as discussed in Ref. Musolino  _et al._ (2022). As a
result, we need to renormalize the eigenvalue $N_{0}^{(2)}$ in order to avoid
over counting the pair condensate fraction. Following the renormalization
procedure as presented in App. B, we find that the pair condensate fraction
$n_{0}^{(2)}$ can be expressed as

$\displaystyle
n_{0}^{(2)}(t)=\frac{N_{0}^{(2)}(t)}{1+\left[2/N_{0}^{(2)}(t)\right]\sum_{\mathbf{k}\neq
0}\absolutevalue{\kappa_{\mathbf{k}}^{2}}\rho_{\mathbf{k}}}.$ (34)

We note that the renormalization factor in the denominator of Eq. (34) is
identical to the renormalization factor presented in Ref. Musolino  _et al._
(2022), where a similar analysis is performed for the single channel Bose gas.
This is due to the presence of the background gas in the open channel. The
evolution of the pair fraction for various values of $k_{n}R^{*}$ will be
discussed in Sec. III.1.

####  II.2.3  The dynamical two-body contact

When two bosons in a quantum gas separated by a distance
$r_{12}=\absolutevalue{\mathbf{r}_{1}-\mathbf{r}_{2}}$ approach each other,
such that $\Lambda^{-1}\ll
r_{12}\ll\\{n^{-1/3},\absolutevalue{a},\lambda_{\mathrm{dB}},\text{etc.}\\}$,
the many-body wave function
$\Psi_{\mathrm{MB}}(\mathbf{r}_{1}\sigma_{1},\mathbf{r}_{2}\sigma_{2},...,\mathbf{r}_{N}\sigma_{N})$
can be factorized as Werner  _et al._ (2009)

$\displaystyle\Psi_{\mathrm{MB}}(\mathbf{r}_{1}\sigma_{1},..,\mathbf{r}_{N}\sigma_{N})\approx\phi_{2B}(r_{12})\mathcal{A}(\mathbf{R}_{12}\sigma_{P},..,\mathbf{r}_{N}\sigma_{N}),$
(35)

with center of mass coordinate
$\mathbf{R}_{12}=\frac{m_{1}\mathbf{r}_{1}+m_{2}\mathbf{r}_{2}}{m_{1}+m_{2}}$
and channel index $\sigma_{i}$. The above factorization is a starting point in
the derivation of a set of universal relations related to the two-body, or Tan
contact Tan (2008a, b, c); Werner and Castin (2012a); Corson and Bohn (2015)
$\mathcal{C}_{2}$, defined as 444We have normalized our momentum distribution
$\sum_{\mathbf{k}}\rho_{\mathbf{k}}$ using the particle density $n$ instead of
the particle number $N$. As a result, we are analysing the contact density
$\mathcal{C}_{2}$.

$\mathcal{C}_{2}\mathrel{\mathop{\equiv}\limits_{\scriptstyle\small\Lambda\gg
k\gg\\{k_{\mathrm{n}},a^{-1},\lambda^{-1}_{\mathrm{dB}},\text{etc.}\\}}}k^{4}\rho_{\mathbf{k}}.$
(36)

Since the two-body contact effectively measures the probability for pairs of
atoms to be close together, it is a valuable parameter in the analysis of the
unitary Bose gas. In our doublet model the asymptotic scaling behavior
described by Eq. (36) emerges for momenta $k$ in the universal regime compared
to the system specific length scales $a$, $k_{\mathrm{n}}$ and
$R_{\mathrm{eff}}$. The set of HFB equations presented in Sec. II.2 then allow
us study the time evolution of $\rho_{\mathbf{k}}$ and hence to compute the
two-body contact. Alternatively, by integration over the regularized part of
the wave function
$\mathcal{A}(\mathbf{R}_{12}\sigma_{P},..,\mathbf{r}_{N}\sigma_{N})$
introduced in Eq. (35), it is possible to express the contact in terms of the
open-open channel component of the two-body correlation function
$g_{\mathrm{PP}}^{(2)}(\mathbf{r},\mathbf{r}^{\prime})=\braket{\hat{\psi}_{\mathrm{P}}^{\dagger}(\mathbf{r})\hat{\psi}_{\mathrm{P}}^{\dagger}(\mathbf{r}^{\prime})\hat{\psi}_{\mathrm{P}}(\mathbf{r})\hat{\psi}_{\mathrm{P}}(\mathbf{r}^{\prime})}$,
such that

$\displaystyle g_{\mathrm{PP}}^{(2)}(\mathbf{r},\mathbf{r}^{\prime})$
$\displaystyle\mathrel{\mathop{\approx}\limits_{\Lambda^{-1}\ll
r_{12}\ll\\{n^{-1/3},\absolutevalue{a},\lambda_{\mathrm{dB}},\text{etc.}\\}}}\absolutevalue{\phi_{2B}(r\rightarrow
0)}^{2}\frac{\mathcal{C}_{2}}{(4\pi)^{2}}.$ (37)

Relating the two-body wave function to the open-channel wave function as
$\phi_{2B}(r)=-a^{-1}\lim_{E\rightarrow 0}\Psi_{P}(r)$ Werner and Castin
(2012b), it is possible to compute the well-defined zero range
($\Lambda\rightarrow\infty)$ and zero energy limit of the two-body wave
function 555we would have to consider energy corrections in order to obtain
the sub leading order contact. Using
$\hat{T}\ket{\mathbf{k}}=\hat{V}_{\mathrm{eff}}\ket{\Psi_{P}}$ with the
effective potential as defined in Eq. (6) we obtain

$\displaystyle\left(v-\frac{g^{2}}{2\nu}\right)\phi_{2B}(r\rightarrow
0)\underset{\stackrel{{\scriptstyle E\rightarrow
0}}{{\Lambda\rightarrow\infty}}}{\approx}-\frac{4\pi\hbar^{2}}{m},$ (38)

such that the two-body contact can be expressed as

$\displaystyle\mathcal{C}_{2}=\frac{m^{2}}{\hbar^{4}}\braket{\left(v-\frac{g^{2}}{2\nu}\right)^{2}\hat{\psi}_{P}^{\dagger}\hat{\psi}_{P}^{\dagger}\hat{\psi}_{P}\hat{\psi}_{P}(r\rightarrow
0)}.$ (39)

As derived in App. C, Eq. (39) is the bosonic version of the generalized Tan
relation introduced in Ref. Braaten  _et al._ (2008). Contrary to its single
channel analogue, which is retrieved for $g=0$ Braaten and Platter (2008) and
only holds in the broad resonance limit, the expression for the generalized
contact is valid for all considered values of $k_{\mathrm{n}}R^{*}$. In
addition, as outlined in App. C, the zero-range limit of Eq. (39) can be
directly related to the set of cumulants as presented in Sec. II.2. This
facilitates the analysis of the contact as extracted from the tail of the
momentum distribution as well as computed using the generalized Tan relation.
The dynamics of this quantity are the subject of Sec. III.2.

###  II.3  Embedded two-body interactions

Having discussed the two-body as well as the many-body problem in terms of a
two channels model, we now proceed to study the resonance-width dependent
physics of two-body interactions embedded in a many-body environment. This
allows us to probe the effect of the medium on two-body interactions and
relate these effects to the observation of non-zero values of the Z-parameter
in the unitary Fermi gas Koetsier  _et al._ (2009); Partridge  _et al._
(2005).

####  II.3.1  The embedded two-body transition matrix

Following Ref. Colussi  _et al._ (2018b), we start our analysis of the
embedded two-body interactions by extending the two-body transition matrix in
vacuum to its embedded analogue. To this extend we decompose
$\kappa_{\mathbf{k}}$ and $\psi_{m}$ in terms of a complete basis set with
open- and closed channel wave functions
$\Psi^{R}_{\mathrm{P},\mu}(\mathbf{k})$ and $\Phi_{\mathrm{Q},\mu}$
respectively, such that we find

$\begin{bmatrix}\kappa_{\mathbf{k}}\\\
\psi_{m}\end{bmatrix}=\sum_{\mu}c_{\mu}(t)\begin{bmatrix}\Psi^{R}_{\mathrm{P},\mu}(\mathbf{k})\\\
\Phi_{\mathrm{Q},\mu}\end{bmatrix},\vspace{-0.15cm}$ (40)

with $c_{\mu}(t)=c_{\mu}e^{-iE_{\mu}t/\hbar}$ and
$\Phi_{\mathrm{Q},\mu}=\braket{\phi}{\Psi_{\mathrm{Q},\mu}}/\sqrt{2}$,
analogous to the definition of the amplitude $\Phi_{\mathrm{Q}}$ in Sec. II.1.
In this model we treat the density effects as quasi-stationary Kira (2015);
Kira and Koch (2011); Colussi  _et al._ (2018b). The subscript $R$ introduced
in the open channel wave function $\Psi_{\mathrm{P},\mu}^{R}(\mathbf{k})$
indicates the usage of right eigenvectors. The origin of the asymmetry will be
discussed shortly. Using Eq. (21), we recognize that $\psi_{m}$ can be
interpreted as the closed channel analogue of the pairing matrix. The
coefficients $c_{\mu}$ then tell us how much of the total pairing matrix,
consisting of open- and closed channel contributions, is contained in a dimer-
basis state $\mu$.

Applying the quasi-stationary approximation and neglecting the source terms
related to the presence of the atomic condensate, we use Eq. (40) to obtain
the following set of two channel eigenvalue equations

$\displaystyle E_{\mu}\Psi^{R}_{\mathrm{P},\mu}(\mathbf{k})$
$\displaystyle=2h_{\mathbf{k}}\Psi^{R}_{\mathrm{P},\mu}(\mathbf{k})+(1+2\rho_{\mathbf{k}})\zeta(2\mathbf{k})\left(v\sum_{\mathbf{q}}\Psi^{R}_{\mathrm{P},\mu}(\mathbf{q})\zeta^{*}(2\mathbf{q})+g\Phi_{\mathrm{Q},\mu}\right),$
(41a) $\displaystyle E_{\mu}\Phi_{\mathrm{Q},\mu}$
$\displaystyle=\nu\Phi_{\mathrm{Q},\mu}+\frac{g}{2}\sum_{\mathbf{k}}\Psi^{R}_{\mathrm{P},\mu}(\mathbf{k})\zeta^{*}(2\mathbf{k}).$
(41b)

Equations (41a) and (41b) should be compared to the two-body eigenvalue
equations as presented in Eqs. (9) and (9b). Here the effective interaction
potential operator $\hat{V}_{\mathrm{eff}}$ that has been introduced in Eq.
(6) is replaced by the interaction potential operator
$\hat{\mathcal{V}}_{\mathrm{eff}}$, where
$\mathcal{\hat{V}}_{\mathrm{eff}}=\hat{B}\hat{V}_{\mathrm{eff}}$, with
$\hat{B}$ the Bose-enhancement operator
$\bra{\mathbf{k},\mathbf{k^{\prime}}}\hat{B}=(1+\rho_{\mathbf{k}}+\rho_{\mathbf{k^{\prime}}})\bra{\mathbf{k},\mathbf{k^{\prime}}}$
Colussi  _et al._ (2018b). The Bose-enhancement of open channel excitations
causes the asymmetry of the open channel eigenvalue equation.

Additionally, the kinetic energy term in Eq. (9) is represented by the
Hartree-Fock term $h(\mathbf{k})\approx\hbar^{2}k^{2}/2m+2vn_{op}$ in Eq.
(41a), meaning that the energy in the open channel is mean-field shifted.
Therefore, analogously to Ref. Colussi  _et al._ (2018b), we quantify the
binding energy and the detuning relative to the mean-field energy shifted
threshold, such that $\mathcal{E}_{\mu}\equiv E_{\mu}-4vn_{op}$ and
$\mathpzc{v}\equiv\nu-4vn_{op}$ respectively.

By using the effective potential interaction operator
$\mathcal{\hat{V}}_{\mathrm{eff}}$, we can straightforwardly introduce the
embedded transition operator as
$\mathcal{\hat{T}}\ket{\mathbf{k}}=\mathcal{\hat{V}}_{\mathrm{eff}}\ket{\Psi_{P,\mu}^{R}}$
and, for the separable potential introduction, obtain the embedded transition
matrix

$\mathcal{T}=\braket{\mathbf{k}}{\hat{B}v}{\zeta}\braket{\zeta}{\psi_{\mathrm{P},\mu}^{R,+}}+\frac{\frac{g^{2}}{2}\braket{\psi_{\mathrm{P},\mu}^{R,-}}{\hat{B}}{\zeta}\braket{\zeta}{\psi_{\mathrm{P},\mu}^{R,+}}}{\mathcal{E}_{\mu}-\mathpzc{v}-\frac{g^{2}}{2}\braket{\zeta}{\mathcal{\hat{G}}_{\mathrm{P}}(\mathcal{E}_{\mu})\hat{B}}{\zeta}},$
(42)

where $\mathcal{\hat{G}}_{\mathrm{P}}$ is the embedded open channel Green’s
operator
$\mathcal{\hat{G}}_{\mathrm{P}}(\mathcal{E}_{\mathrm{\mu}})=(\mathcal{E}_{\mathrm{\mu}}-\hat{B}v\ket{\zeta}\bra{\zeta})^{-1}$.
Equation (42) is related to the many-body transition operator $\hat{T}^{MB}$
Stoof  _et al._ (2009) according to $\mathcal{\hat{T}}=\hat{B}\hat{T}^{MB}$.

One of the dimer-basis states $\mu=D$ corresponds to the embedded dimer. The
energy of this dimer $\mathcal{E}_{\mathrm{D}}$ can be extracted from the pole
of Eq. (42). Through the analysis of the zero energy limit of the embedded
transition matrix, we can extract an embedded analogue of the scattering
length $\mathpzc{a}$ and the effective range $\mathpzc{R_{\mathrm{eff}}}$,
such that 666Equation (43) is valid for interactions with constant form
factors.

$\mathcal{T}_{E\rightarrow
0}\approx\frac{4\pi\hbar^{2}\mathpzc{a}}{m}\left(1-ik\mathpzc{a}+\frac{\mathpzc{a}\mathpzc{R_{\mathrm{eff}}}-2\mathpzc{a}^{2}}{2}k^{2}+\mathcal{O}(k^{3})\right).$
(43)

The evolution of these quantities as a function of the resonance width (or
equivalently $k_{\mathrm{n}}R^{*}$) will be the investigated in more detail in
Sec. III.3.

####  II.3.2  The embedded dimer wave function normalization factor

As mentioned in Sec. II.1.2, the Z-parameter is zero at unitarity for the two-
body model. However, a finite value of the Z-parameter was previously
predicted and observed in the unitary Fermi gas Koetsier  _et al._ (2009);
Partridge  _et al._ (2005). This motivates us to analyze the embedded version
of the Z-parameter. This parameter quantifies the division of the embedded
dimer with energy $\mathcal{E}_{\mathrm{D}}$ amongst the open- and closed
channel subspaces.Analogous to Eq. (19), we define the $\mathcal{Z}$-parameter
as

$\begin{bmatrix}\ket{\Psi^{R}_{\mathrm{P,D}}}\\\
\ket{\Psi_{\mathrm{Q,D}}}\end{bmatrix}=\sqrt{\mathcal{Z}}\begin{bmatrix}\mathcal{\hat{G}}_{\mathrm{P}}(E_{\mathrm{D}})\hat{B}\beta\ket{\zeta}\braket{\zeta}{\phi}\\\
\ket{\phi}\end{bmatrix},$ (44)

such that $\braket{\Psi_{\mathrm{D}}}{\Psi_{\mathrm{D}}}=1$ and we find that

$\mathcal{Z}=\left[1+\frac{\frac{g^{2}}{2}\sum_{\mathbf{k}}\zeta(2\mathbf{k})\frac{(1+2\rho_{\mathbf{k}})}{(\mathcal{E}_{\mathrm{D}}-\frac{\hbar^{2}k^{2}}{m})^{2}}}{\left(1-v\sum_{\mathbf{k}}\zeta(2\mathbf{k})\frac{(1+2\rho_{\mathbf{k}})}{\mathcal{E}_{\mathrm{D}}-\frac{\hbar^{2}k^{2}}{m}}\right)^{2}}\right]^{-1}.$
(45)

As the $\mathcal{Z}$-parameter can be linked to the change in $\delta\mu$
Jochim  _et al._ (2003); Falco and Stoof (2007) as well as to the observed
atom-loss in molecular probe experiments as investigated in Ref. Partridge
_et al._ (2005). The relation between the $\mathcal{Z}$-parameter and
experimental observables makes it a valuable quantity to connect the theory of
embedded dimers to experiments. As such, we will analyze its resonance width
dependent evolution in Sec. III.4.

####  II.3.3  The relative dimer state occupation

In addition to analyzing how the dimer is distributed amongst the two
subspaces, we also wish to quantify the relative importance of this dimer
state to the total pairing field. Therefore, we introduce the relative dimer
state occupation $\mathcal{F}_{\mathrm{D}}$ Kira (2015), defined as

$\mathcal{F}_{\mathrm{D}}=\frac{\absolutevalue{c_{D}}^{2}}{\sum_{\mu}\absolutevalue{c_{\mu}}^{2}},$
(46)

in terms of the coefficients $c_{\mu}$ of the dimer basis in Eq. (3). In order
to compute the numerator of Eq. (46), we multiply Eq. (40) from the left side
by
$\begin{bmatrix}[\Psi_{\mathrm{P,D}}^{L}(\mathbf{k})]^{*}&2[\Phi_{\mathrm{Q,D}}]^{*}/\sqrt{V}\end{bmatrix}$,
sum over all $\mathbf{k}$ and compute the square of the absolute value, such
that we find

$\absolutevalue{c_{\mathrm{D}}}^{2}=\absolutevalue{\sum_{\mathbf{k}}\kappa_{\mathbf{k}}\frac{[\Psi_{\mathrm{P,D}}^{R}(\mathbf{k})]^{*}}{1+2\rho_{\mathbf{k}}}+\sqrt{2\mathcal{Z}}\psi_{m}}^{2},$
(47)

where we have used the normalization condition
$\sum_{\mathbf{k}}[\Psi_{P,\lambda}^{L}(\mathbf{k})]^{*}\Psi_{\mathrm{P},\mu}^{R}(\mathbf{k})+2[\Phi_{Q,\lambda}]^{*}\Phi_{\mathrm{Q},\mu}=\delta_{\lambda,\mu}$,
the relation between the right- and left-eigenvectors of the open channel
subspace
$\Psi_{\mathrm{P},\mu}^{R}(\mathbf{k})=\left(1+2\rho_{\mathbf{k}}\right)\Psi_{\mathrm{P},\mu}^{L}(\mathbf{k})$
and where we have applied Eq. (44).

Applying a similar strategy, we next compute the denominator of Eq. (46) by
multiplying Eq. (40) from the left side by
$\begin{bmatrix}[c_{\lambda}\Psi_{P,\lambda}^{L}(\mathbf{k})]^{*}&2[c_{\lambda}\Psi_{Q,\lambda}]^{*}\end{bmatrix}$,
summing over all values of $\mathbf{k}$ and computing the square of the
absolute value. Once more we exploit the normalization condition and find

$\absolutevalue{c_{\mu}}^{2}=\sum_{\mathbf{k}}\frac{\left[c_{\mu}\Psi_{\mathrm{P},\mu}(\mathbf{k})\right]^{*}}{1+2\rho_{\mathbf{k}}}\kappa_{\mathbf{k}}+2\left[c_{\mu}\Psi_{\mathrm{Q},\mu}\right]^{*}\psi_{m}.$
(48)

Summing the previous expression over all basis-states $\mu$ and using Eq. (40)
in order to rewrite the dimer-basis states in terms of cumulants, we obtain
the following expression for the denominator of $\mathcal{F}_{\mathrm{D}}$

$\sum_{\mu}\absolutevalue{c_{\mu}}^{2}=\sum_{\mathbf{k}}\frac{\absolutevalue{\kappa_{\mathbf{k}}}^{2}}{1+2\rho_{\mathbf{k}}}+2\absolutevalue{\psi_{m}}^{2}.$
(49)

The expression for $\mathcal{F}_{\mathrm{D}}$ as derived here is normalized at
every time step and can be applied to two channel systems. Consequently, we
can analyze how the relative dimer state occupation evolves as a function of
the time for various values of the resonance width, presenting our results in
Sec. III.5.

##  III  Results

Having outlined our model, we now analyze the results for the dynamics of the
quenched unitary Bose gas over a range of $k_{\mathrm{n}}R^{*}$. In this
study, we follow the procedure as outlined by Refs. Musolino  _et al._ (2019);
Colussi  _et al._ (2018b), starting with a non-interacting pure atomic
condensate and perform an effectively sudden quench to unitarity 777The
analysis is not limited to a fixed atomic species or resonance and can be
applied generally for varying values of $k_{\mathrm{n}}R^{*}$. In our
simulations, we consider a condensate of ${}^{85}\text{Rb}$ atoms at an
experimentally relevant density of $n=4\times 10^{12}\,\text{cm}^{-3}$, such
that $t_{\mathrm{n}}=69.84\,\mu s$. We model the Feshbach resonance located at
$B_{0}=155.04$ G with magnetic field width $\Delta B=10.7$ G and background
scattering length $a_{\mathrm{bg}}=-443\,a_{0}$ Claussen  _et al._ (2003) and
quench the magnetic field to resonance in a time span of
$\tau_{in}=0.072\,t_{\mathrm{n}}$, such that the quench is effectively sudden.
We mark the completion of the quench as $t=0$ and start our analysis with this
broad-resonance system, where $k_{\mathrm{n}}R^{*}\approx 10^{-3}$, in order
to verify the consistency of the coupled-channels model with the single-
channel model as presented in Ref. Colussi  _et al._ (2018b) and then
gradually increase the value of $k_{\mathrm{n}}R^{*}$.. Whereas we vary the
value of $k_{\mathrm{n}}R^{*}$, we keep the value of $a_{\mathrm{bg}}$ fixed
in order to satisfy the diluteness criterion $na^{3}_{\mathrm{bg}}\ll 1$. This
means that our results can be universally extended to different atomic species
in the dilute regime with the same values of $k_{\mathrm{n}}R^{*}$. In the
following sections, we express all our results in terms of Fermi-units unless
mentioned otherwise and indicate the dimensionless resonance widths in terms
of the many-body classification $k_{\mathrm{n}}R^{*}$.

###  III.1  Population fractions

Figure 2 illustrates the (normalized) population dynamics as the time spent at
unitarity progresses.

(a)

(b)

(c)

(d)

Figure 2: Population fractions at unitarity: (a) atomic condensate fraction
$n_{a}$, (b) excited state fraction $n_{exc}$, (c) molecular condensate
fraction $n_{m}$ and (d) pair condensate fraction $n_{0}^{(2)}$ as a function
of the time $t/t_{\mathrm{n}}$ for various values of the dimensionless
resonance width parameter $k_{\mathrm{n}}R^{*}$. The inset in (c) represents
the molecular condensate fraction as a function of the time, with the time-
axis rescaled with the embedded transition time
$t_{*}=\sqrt{t_{\mathrm{n}}\tau}$. The legend in the inset indicates the
resonance width in multiples of $k_{\mathrm{n}}R^{*}$. The dashed black line
in (d) represents the single channel data as presented in Ref. Musolino  _et
al._ (2022). Our results are consistent with this calculation.

In particular Fig. 2(a) shows how the atomic condensate is quantum depleted.
The atoms that leave the condensate can either form excitations (Fig. 2(b)) or
closed channel molecules (Fig. 2(c)). Since the formed closed channel
molecules have a shorter lifetime for broader resonances
($k_{\mathrm{n}}R^{*}\ll 1$), more excitations are formed at early times with
respect to resonances with larger $k_{\mathrm{n}}R^{*}$. Additionally, the
presence of background excitations Bose-enhances the production of extra
excitations. Consequently, the excited state fraction as presented in Fig.
2(b) grows more rapidly for small $k_{\mathrm{n}}R^{*}$ at early times.

Aiming to analyze the early-time dynamics of the molecular condensate on a
more quantitative level, we consider the narrow-resonance
($k_{\mathrm{n}}R^{*}\gg 1$) and early-time limit ($t\ll t_{\mathrm{n}}$) of
Eq. (23). Since all atoms start out in the atomic condensate and Fig. 2(b)
reveals that the excited state fraction remains limited at early times for
resonances with large $k_{\mathrm{n}}R^{*}$, we approximate the atomic wave
function as $\psi_{a}\approx\sqrt{n}$ and neglect the terms scaling with
$\psi_{m}$ and $\kappa_{\mathbf{k}}$ in Eq. (23). Under these approximations,
we can integrate Eq. (23) with respect to the time and obtain

$2\absolutevalue{\psi_{m}}^{2}\propto(t/t_{\ast})^{2},$ (50)

where $t_{*}=\sqrt{t_{\mathrm{n}}\tau}$ is the geometric mean of two relevant
time scales. Here, $\tau=2mR^{*}/k_{\mathrm{n}}\hbar$ is the lifetime
associated to a molecule on the Fermi scale as discussed in Sec. II.1.
Physically $t_{*}$ represents the density-averaged time for open channel atoms
to transition to closed channel molecules. Hence, we refer to $t_{*}$ as the
mean transition time. The scaling of the early time dynamics of resonances in
terms of the mean transition time $t_{*}$ is supported by the inset of Fig.
2(c), where the rescaling of the time-axis with $t_{*}$ results in the
collapse of the molecular condensate fraction curves for sufficiently large
values of $k_{\mathrm{n}}R^{*}$.

Following the stage of rapid initial growth, Fig. 2(c) shows the saturation of
the molecular condensate fraction at later times. The time it takes to reach
the saturated value increases as a function of $k_{\mathrm{n}}R^{*}$. We
interpret the delay in the observed saturation time to be a result of the
increased mean transition time of closed channel molecules for narrower
resonances.

Once the molecular condensate fraction has saturated, the depletion of the
atomic condensate fraction effectively solely results in the formation of new
excitations. This process can be observed in Fig. 2(b), where we recognize
that, at later times, the excited state fraction grows more rapidly for narrow
resonances with respect to broader resonances. Furthermore, Fig. 2(a) shows
how the difference between the atomic condensate fraction for broad and narrow
resonances starts to decrease at later times.

The dynamical evolution of the excited state fraction and the molecular
condensate is echoed by the pair condensate fraction as presented in Fig.
2(d). We recognize that, for all considered resonance widths, a considerable
fraction of the gas contributes to the pair condensate within the simulated
time frame. In addition, this fraction is relatively insensitive to
$k_{\mathrm{n}}R^{*}$ compared to the other fractions in Fig. 2 and reflects
the early time scaling law change observed in the molecular condensate
fraction.

###  III.2  Two-body contact

By applying Eqs. (36) and (71) to the doublet model we can compute the two-
body contact, which encodes the probability of finding clustered particles at
short distances, for different values of $k_{\mathrm{n}}R^{*}$.

Figure 3: Two-body contact rescaled with the density $n^{-4/3}$ fitted from
the tail of the momentum distribution (full lines) and computed using the
generalized Tan relation (points) as a function of the time for various values
of the dimensionless resonance width $k_{\mathrm{n}}R^{*}$. The black line
indicates the linear fit as presented in Eq. (51). In the inset, the half-way
time $\tau_{c}$ to reach the maximum value of the contact is plotted as a
function of $k_{\mathrm{n}}R^{*}$.

Figure 3 reveals that the contact obtained using the generalized Tan relation
agrees with the contact obtained using the tail of the momentum distribution
over the whole range of considered resonance widths. Furthermore, Fig. 3 shows
that the contact for a broad resonance interaction initially grows linearly
and follows the analytic expression derived in Ref. Corson and Bohn (2015)
according to

$n^{-4/3}\mathcal{C}_{2}(t)=\frac{128\pi}{(6\pi^{2})^{2/3}}\frac{t}{t_{\mathrm{n}}}.$
(51)

In addition, consistent with the derivation presented in App. A, the contact
is observed to be proportional to the molecular condensate fraction in the
broad resonance limit Werner  _et al._ (2009).

Considering resonances with large $k_{\mathrm{n}}R^{*}$ on the other hand, we
recognize that the initial growth of the contact is gradually becoming less
rapid as $k_{\mathrm{n}}R^{*}$ increases. In order to quantify the time it
takes for the contact to evolve, we have computed the half-way time
$\tau_{c}$. This time scale is defined as the time it takes to reach half the
maximum value of the contact. The inset of Fig. 3 reveals how $\tau_{c}$
increases as a function of the resonance width. Similarly to the initial
growth of the molecular lifetime, we observe the initial growth of the contact
to evolve from linear to quadratic. In addition, the saturation time increases
for larger values of $k_{\mathrm{n}}R^{*}$. We relate this to the increase in
the mean transition time $t^{*}$ of the closed-channel molecules. The increase
in $t^{*}$ slows down the formation of excitations from closed-channel
molecules. Equation 36 then reveals how a delay in the saturation of the large
momentum modes of the excited state fraction delays the saturation of the
contact. As the large momentum modes of the excited state fraction saturate
more rapidly than the smaller momentum modes Makotyn  _et al._ (2014); Sykes
_et al._ (2014), we expect the contact in Fig. 3 to saturate before the
excited state fraction in Fig. 2(b).

###  III.3  Dimer

In this section we focus on how the evolution of the embedded dimer energy
$\mathcal{E}_{\mathrm{D}}$ as investigated in the single channel limit in
Refs. Musolino  _et al._ (2019); Colussi  _et al._ (2018b) is affected by the
gradual increase of $k_{\mathrm{n}}R^{*}$. This allows us to investigate how
the size of the dimer and the universal scaling of the dimer energy with the
scattering length is affected by the medium as well as the resonance width.
The results of this analysis are presented in Figs. 4 and 5.

Figure 4: Dimer binding energy rescaled with the Fermi energy as a function of
the time spent at unitarity for three different values of the dimensionless
resonance width $k_{\mathrm{n}}R^{*}$. The dimer shrinks progressively with
increasing time at unitarity, as schematically illustrated by the cartoon. The
dashed orange line is consistent with the single channel result (orange
circles) as presented in Ref. Colussi  _et al._ (2018b).

We observe that, contrary to vacuum two-body interactions, the dimer energy
decreases towards more deeply bound values as the time spent at unitarity
progresses. We relate this to the reduction of the magnitude of the embedded
scattering length as presented in Fig. 5(a), which indicates a decrease in the
size of the embedded dimers as illustrated by the sketch in Fig. 4.

The sketch indicates that the localization is more rapid for increasing
$k_{\mathrm{n}}R^{*}$, consistent with the production of more excitations for
broader resonances, as discussed in Sec. III.1. Since these atoms Bose-enhance
open-channel interactions, the change in the effectively experienced open
channel potential $\hat{\mathcal{V}}_{\mathrm{eff}}$ is more drastic for broad
resonances, such that the system is effectively pushed away from resonance
more rapidly, corresponding to a more swift localization of the dimer.

(a)

(b) Figure 5: The effective scattering length $\mathpzc{a}$ (5a) and the
effective range $\mathpzc{R_{\mathrm{eff}}}$ (5b) as a function of the time
spent at unitarity for various values of the dimensionless resonance width
$k_{\mathrm{n}}R^{*}$. The inset in (a) presents the inverse of the scattering
lengths $a$ and $a(k)$ as a function of the dimer wavenumber rescaled with the
Fermi wavenumber. Dashed black lines have been added to (b) in order to
indicate the vacuum limit of the effective range as computed using Eq. (17) .

Besides commenting on the size of the embedded dimer, we notice that Fig. 5(a)
shows that, before the completion of the quench, the scattering length has
already assumed finite values. This is a result of the, albeit limited,
depletion of the initial state during the quench Colussi  _et al._ (2018b).

Similarly to two-body physics, we expect that the universal relation between
the embedded scattering length and the dimer energy
$\absolutevalue{\mathcal{E}_{\mathrm{D}}}\approx\hbar^{2}/(m\mathpzc{a}^{2})$
only holds for sufficiently large values of the scattering length.
Furthermore, the correction to the universal quadratic relation becomes more
important at larger values of the scattering length for increasingly narrow
resonances. Therefore, in addition to analyzing the embedded scattering
length, we present the evolution of the effective range
$\mathpzc{R_{\mathrm{eff}}}$ as a function of the time spent at unitarity in
Fig. 5(b). Initially, the embedded effective range correctly reduces to the
expected vacuum result as extracted from Eq. (17). However, as the time spent
at unitarity progresses, the value of the embedded effective range starts to
increase.

Especially interesting is the observation that, at later times, the effective
range assumes non-negligible values even for broad resonances and cannot be
ignored. This effect is visualized in the inset of Fig. 5, where the inverse
scattering length $\mathpzc{a}^{-1}$ and the function
$\mathpzc{a}(k)^{-1}=\mathpzc{a}^{-1}+\frac{1}{2}\mathpzc{R_{\mathrm{eff}}}k^{2}$
are plotted versus the rescaled dimer wavenumber
$k_{\mathrm{D}}/k_{\mathrm{n}}=\sqrt{\mathcal{E}_{\mathrm{D}}/2E_{\mathrm{n}}}$
and the curves are observed to differ for larger values of
$k_{\mathrm{D}}/k_{\mathrm{n}}$.

###  III.4  Z-parameter

We now aim to analyze how the embedded dimer is distributed amongst the
closed- and open-channel subspaces through the computation of the dynamical
$\mathcal{Z}$-parameter. Consistent with Eq. (19) for a two-body system, our
analysis of the $\mathcal{Z}$-parameter as presented in Fig. 6 reveals that
initially $\mathcal{Z}=0$ at unitarity.

Figure 6: The embedded Z-parameter (Eq. (45)) as a function of the time spent
at unitarity for various values of the dimensionless resonance width
$k_{\mathrm{n}}R^{*}$.

We interpret this as being a result of the dimer extending over an infinite
length scale at unitarity, such that the overlap of this infinite quantity
with the localized closed-channel molecule vanishes. However, in agreement
with the picture of the localization of the dimer as elaborated in Sec. III.3
and as predicted in Ref. Koetsier  _et al._ (2009), the
$\mathcal{Z}$-parameter assumes non-zero values as the time spent at unitarity
increases. Figure 6 shows that the finite value of the $\mathcal{Z}$ increases
as a function of $k_{\mathrm{n}}R^{*}$, consistent with the dimer becoming
progressively less open channel dominated for narrow resonances.

###  III.5  Relative dimer state occupation

In the past two subsections, we have focused on the characteristics of the
embedded dimer state. However, we are yet to quantify the relative importance
of this dimer state with respect to the total pairing matrix. As outlined in
Sec. II.3.3, this importance is gauged by the quantity
$\mathcal{F}_{\mathrm{D}}$ as presented in Eq. (46). The evolution of this
parameter for various values of the resonance width is presented Fig. 7.

Figure 7: The relative dimer state occupation (Eq. (46)) for various values of
the dimensionless resonance width $k_{\mathrm{n}}R^{*}$ as a function of the
time spent at unitarity.

Consistent with the results found in Ref. Kira (2015), the embedded dimer
state dominates the quantum depletion in the case of relatively broad
resonances, displaying evolution of an initially fast growth followed by a
slower relaxation towards a finite value. Furthermore, Fig. 7 reveals that
this complicated dynamics is affected by the resonance width. In the simulated
time-frame, the relatively broad-resonance case presents an upper-bound to the
value of $\mathcal{F}_{\mathrm{D}}$. This means that the embedded dimer state
is comparatively less relevant for pair correlations in narrow resonance
systems.

In order to gain an intuitive picture of this observation, we refer back to
Fig. 4, where we observe that dimers formed using narrow resonance
interactions are less-deeply bound compared to dimers formed using broad
resonance interactions. Consequently, narrow resonance interactions result in
comparatively large dimers, which don’t localize to sizes comparable to the
Fermi scale during the considered time frame. As the pair correlation physics
occurs on the level of the Fermi scales, the relatively extended dimers extend
over a larger range than the typical two-body pairing physics. Hence, the
observation that the dimer state is less relevant for narrower resonances,
corresponding to lower values of the $\mathcal{F}_{\mathrm{D}}$.

##  IV Conclusion

In this work, we employ a two-channel model with finite range pairwise
interactions in order to study the effect of the resonance width on the
dynamics of quenched unitary Bose gases. Using a two-channel model we vary the
value of $k_{\mathrm{n}}R^{*}$, effectively changing the width of the
resonance. This allows us to analyse the onset of non-universal effects. For
increasing $k_{\mathrm{n}}R^{*}$ the lifetime of the closed channel molecules
is increased relative to the Fermi scale, and the pairwise production of
excitations in the open channel is decreased. Consequently the early-time
growth of the dynamical two body contact becomes more gradual, transitioning
from linear- to quadratic behaviors. In the limit $k_{\mathrm{n}}R^{*}\gg 1$,
the early-time growth of the molecular condensate fraction also shifts from
linear to quadratic and is set by the mean molecular lifetime
$t_{*}=\sqrt{t_{\mathrm{n}}\tau}$. Here the system scales as combinations of
resonance parameters and the density. Within the embedded few-body problem, we
study the localization of dressed embedded dimers that are bound purely by
many-body effects and are dressed by multichannel couplings. This analysis
reveals their dominant role in the quantum depletion and how they become
increasingly closed channel dominated as $k_{\mathrm{n}}R^{*}$ is increased.

Our analysis of the closed-channel contribution to the dressed embedded dimers
and the pair condensate fraction opens up additional experimental probes in
the unitary regime Partridge  _et al._ (2005); Musolino  _et al._ (2019);
Hodby  _et al._ (2005). Furthermore, including third-order cumulants in the
two-channel model along the lines of Refs. Colussi  _et al._ (2018b, 2020);
Musolino  _et al._ (2022) represents an interesting prospect for future
research, allowing for the investigation of Efimov physics, the triple
condensate fraction and atom-molecule interactions in the unitary Bose gas for
various values of the resonance width.

##  V Acknowledgments

The authors thank Thomas Secker for supplying a full coupled-channels
calculation data set in order to compute the momentum-space cut-off and
acknowledge Paul Mestrom and Jing-Lun Li for useful discussions. D.J.M.A.B,
S.M. and S.J.J.M.F.K acknowledge financial support by the Netherlands
Organisation for Scientific Research (NWO) under Grant No. 680-47-623. V.E.C.
acknowledges financial support from Provincia Autonoma di Trento and the
Italian MIUR under the PRIN2017 projectCEnTraL.

## References

  * Feshbach (1958) H. Feshbach, Annals of Physics 5, 357 (1958).
  * Chin  _et al._ (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010).
  * Makotyn  _et al._ (2014) P. Makotyn, C. E. Klauss, D. L. Goldberger, E. Cornell, and D. S. Jin, Nat. Phys. 10, 116 (2014).
  * Klauss  _et al._ (2017) C. E. Klauss, X. Xie, C. Lopez-Abadia, J. P. D’Incao, Z. Hadzibabic, D. S. Jin, and E. A. Cornell, Phys. Rev. Lett. 119, 143401 (2017).
  * Eigen  _et al._ (2017) C. Eigen, J. A. P. Glidden, R. Lopes, N. Navon, Z. Hadzibabic, and R. P. Smith, Phys. Rev. Lett. 119, 250404 (2017).
  * Eigen  _et al._ (2018) C. Eigen, J. A. P. Glidden, R. Lopes, E. A. Cornell, R. P. Smith, and Z. Hadzibabic, Nature 556, 221 (2018).
  * Ho (2004) T.-L. Ho, Phys. Rev. Lett. 92, 090402 (2004).
  * Schäfer and Teaney (2009) T. Schäfer and D. Teaney, Rep. Prog. Phys. 72, 126001 (2009).
  * Braaten and Hammer (2003) E. Braaten and H.-W. Hammer, Phys. Rev. Lett. 91, 102002 (2003).
  * Castin and Werner (2012) Y. Castin and F. Werner,  _The BCS-BEC Crossover and the Unitary Fermi Gas_ , edited by W. Zwerger (Springer Berlin Heidelberg, Berlin, Heidelberg, 2012) pp. 127–191.
  * Efimov (1971) V. Efimov, Sov. J. Nucl.Phys. 12 (1971).
  * Efimov (1979) V. Efimov, Sov. J. Nucl.Phys. 29 (1979).
  * D’Incao  _et al._ (2018) J. P. D’Incao, J. Wang, and V. E. Colussi, Phys. Rev. Lett. 121, 023401 (2018).
  * Colussi  _et al._ (2018a) V. E. Colussi, J. P. Corson, and J. P. D’Incao, Phys. Rev. Lett. 120, 100401 (2018a).
  * Ho  _et al._ (2012) T.-L. Ho, X. Cui, and W. Li, Phys. Rev. Lett. 108, 250401 (2012).
  * Musolino  _et al._ (2019) S. Musolino, V. E. Colussi, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 100, 013612 (2019).
  * Colussi  _et al._ (2018b) V. E. Colussi, S. Musolino, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 98, 051601(R) (2018b).
  * Sykes  _et al._ (2014) A. G. Sykes, J. P. Corson, J. P. D’Incao, A. P. Koller, C. H. Greene, A. M. Rey, K. R. A. Hazzard, and J. L. Bohn, Phys. Rev. A 89, 021601(R) (2014).
  * Kain and Ling (2014) B. Kain and H. Y. Ling, Phys. Rev. A 90, 063626 (2014).
  * Corson and Bohn (2015) J. P. Corson and J. L. Bohn, Phys. Rev. A 91, 013616 (2015).
  * Ancilotto  _et al._ (2015) F. Ancilotto, M. Rossi, L. Salasnich, and F. Toigo, Few-Body Systems 56, 801 (2015).
  * Yin and Radzihovsky (2016) X. Yin and L. Radzihovsky, Phys. Rev. A 93, 033653 (2016).
  * Muñoz de las Heras  _et al._ (2019) A. Muñoz de las Heras, M. M. Parish, and F. M. Marchetti, Phys. Rev. A 99, 023623 (2019).
  * Gao  _et al._ (2020) C. Gao, M. Sun, P. Zhang, and H. Zhai, Phys. Rev. Lett. 124, 040403 (2020).
  * Colussi  _et al._ (2020) V. E. Colussi, H. Kurkjian, M. Van Regemortel, S. Musolino, J. van de Kraats, M. Wouters, and S. J. J. M. F. Kokkelmans, Phys. Rev. A 102, 063314 (2020).
  * Musolino  _et al._ (2022) S. Musolino, H. Kurkjian, M. Van Regemortel, M. Wouters, S. J. J. M. F. Kokkelmans, and V. E. Colussi, “Bose-einstein condensation of efimovian triples in the unitary bose gas,”  (2022).
  * Partridge  _et al._ (2005) G. B. Partridge, K. E. Strecker, R. I. Kamar, M. W. Jack, and R. G. Hulet, Phys. Rev. Lett. 95, 020404 (2005).
  * Zwerger (2011) W. Zwerger, ed.,  _The BCS-BEC Crossover and the Unitary Fermi Gas_, Lecture Notes in Physics Vol. 836 (Springer, 2011).
  * Nozières and Schmitt-Rink (1985) P. Nozières and S. Schmitt-Rink, Journal of Low Temperature Physics 59, 195 (1985).
  * Holland  _et al._ (2001a) M. Holland, S. J. J. M. F. Kokkelmans, M. L. Chiofalo, and R. Walser, Phys. Rev. Lett. 87, 120406 (2001a).
  * Feshbach (1993) H. Feshbach,  _Theoretical Nuclear Physics, Nuclear Reactions_ , Theoretical Nuclear Physics (Wiley, 1993).
  * Kokkelmans  _et al._ (2002) S. J. J. M. F. Kokkelmans, J. N. Milstein, M. L. Chiofalo, R. Walser, and M. J. Holland, Phys. Rev. A 65, 053617 (2002).
  * (33) T. Secker (private communication).
  * Note (1) For instance, in the case of the $B_{0}=155.04$ G resonance in ${}^{85}\text{Rb}$ Claussen  _et al._ (2003), we find $\Lambda=0.6/r_{\mathrm{vdW}}$, with $r_{\mathrm{vdW}}$ the van der Waals length of ${}^{85}\text{Rb}$ Chin  _et al._ (2010). We obtained the momentum cut-off $\Lambda$ by comparing the coupled-channels data from tho  to a polynomial where we keep the zeroth order term in the low-energy expansion of $a(B)^{-1}$. Other resonances in different atomic species result in different calibrated values of the cut-off. For example, the calibration of the $B_{0}=402$ G resonance in 39K yields $\Lambda=0.46/r_{\mathrm{vdW}}$. The calibrated value depends on the width of the considered resonance and approaches the value $\Lambda=2\pi/\bar{a}$ Colussi  _et al._ (2018b); Musolino  _et al._ (2019) with mean-scattering length $\bar{a}=0.955r_{vdw}$ in the single channel limit.
  * Kokkelmans and Holland (2002) S. J. J. M. F. Kokkelmans and M. J. Holland, Phys. Rev. Lett. 89, 180401 (2002).
  * Snyder  _et al._ (2012) V. D. Snyder, S. J. J. M. F. Kokkelmans, and L. D. Carr, Phys. Rev. A 85, 033616 (2012).
  * Taylor (2012) J. Taylor,  _Scattering Theory: The Quantum Theory of Nonrelativistic Collisions_, Dover Books on Engineering (Dover Publications, 2012).
  * Note (2) By applying Eq.+.2777em(2), we find that Eq.+.2777em(14) can be rewritten as $\text{tan}[\delta_{0}(k)]=-ka_{\mathrm{bg}}-\left(kR^{*}-\frac{B-B_{0}}{\Delta B}\right)^{-1}$.
  * Törmä and Sengstock (2014) P. Törmä and K. Sengstock, eds.,  _Quantum gas experiments-exploring many-body states_ (Imperial College Press, London, 2014) Chap. 4.
  * Duine and Stoof (2004) R. Duine and H. Stoof, Physics Reports 396, 115 (2004).
  * Cohen-Tannoudji and Guéry-Odelin (2011) C. Cohen-Tannoudji and D. Guéry-Odelin,  _Advances in atomic physics: an overview_ (World Scientific Publishing, 2011).
  * Falco and Stoof (2007) G. M. Falco and H. T. C. Stoof, Phys. Rev. A 75, 023612 (2007).
  * Timmermans  _et al._ (1999) E. Timmermans, P. Tommasini, R. Côté, M. Hussein, and A. Kerman, Phys. Rev. Lett. 83, 2691 (1999).
  * van Abeelen and Verhaar (1999) F. A. van Abeelen and B. J. Verhaar, Phys. Rev. Lett. 83, 1550 (1999).
  * Holland  _et al._ (2001b) M. Holland, J. Park, and R. Walser, Phys. Rev. Lett. 86, 1915 (2001b).
  * Kokkelmans  _et al._ (2001) S. J. J. M. F. Kokkelmans, H. M. J. Vissers, and B. J. Verhaar, Phys. Rev. A 63, 031601 (2001).
  * Altman and Vishwanath (2005) E. Altman and A. Vishwanath, Phys. Rev. Lett. 95, 110404 (2005).
  * Kira (2015) M. Kira, Annals of Physics 356, 185 (2015).
  * Blaizot and Ripka (1986) J.-P. Blaizot and G. Ripka,  _Quantum theory of finite systems_ , Vol. 3 (MIT press Cambridge, 1986).
  * Hodby  _et al._ (2005) E. Hodby, S. T. Thompson, C. A. Regal, M. Greiner, A. C. Wilson, D. S. Jin, E. A. Cornell, and C. E. Wieman, Phys. Rev. Lett. 94, 120402 (2005).
  * Kira (2014) M. Kira, Annals of Physics 351, 200 (2014).
  * Yang (1962) C. N. Yang, Rev. Mod. Phys. 34, 694 (1962).
  * Note (3) We can generally define $\delta\hat{\psi}$ as $\begin{bmatrix}\psi_{P}-\psi_{a},&\psi_{Q}-\psi_{c}\end{bmatrix}^{\text{T}}$, with $\psi_{c}$ a closed channel atomic condensate. Next, by using the single resonance approximation and applying Eq.+.2777em(21) we retain only a small portion of the closed channel part and rewrite $\mathinner{\delimiter 69632778{\delta\hat{\psi}_{Q}\delta\hat{\psi}_{Q}(\mathbf{r})}\delimiter 86414091}=\sqrt{2}\phi(r)\mathinner{\delimiter 69632778{\hat{b}_{0}}\delimiter 86414091}$.
  * Leggett (2008) A. J. Leggett,  _Quantum Liquids: Bose condensation and Cooper pairing in condensed-matter systems_, Vol. 9780198526438 (Oxford University Press, United Kingdom, 2008).
  * Werner  _et al._ (2009) F. Werner, L. Tarruell, and Y. Castin, The European Physical Journal B 68, 401 (2009).
  * Tan (2008a) S. Tan, Annals of Physics 323, 2952 (2008a).
  * Tan (2008b) S. Tan, Annals of Physics 323, 2971 (2008b).
  * Tan (2008c) S. Tan, Annals of Physics 323, 2987 (2008c).
  * Werner and Castin (2012a) F. Werner and Y. Castin, Phys. Rev. A 86, 053633 (2012a).
  * Note (4) We have normalized our momentum distribution $\sum@\slimits@_{\mathbf{k}}\rho_{\mathbf{k}}$ using the particle density $n$ instead of the particle number $N$. As a result, we are analysing the contact density $\mathcal{C}_{2}$.
  * Werner and Castin (2012b) F. Werner and Y. Castin, Phys. Rev. A 86, 013626 (2012b).
  * Note (5) We would have to consider energy corrections in order to obtain the sub leading order contact.
  * Braaten  _et al._ (2008) E. Braaten, D. Kang, and L. Platter, Phys. Rev. A 78, 053606 (2008).
  * Braaten and Platter (2008) E. Braaten and L. Platter, Phys. Rev. Lett. 100, 205301 (2008).
  * Koetsier  _et al._ (2009) A. Koetsier, P. Massignan, R. A. Duine, and H. T. C. Stoof, Phys. Rev. A 79, 063609 (2009).
  * Kira and Koch (2011) M. Kira and S. W. Koch,  _Semiconductor Quantum Optics_ (Cambridge University Press, 2011).
  * Stoof  _et al._ (2009) H. T. Stoof, K. B. Gubbels, and D. B. Dickerscheid,  _Ultracold Quantum Fields_, Theoretical and Mathematical Physics (Springer, Berlin, Germany, 2009).
  * Note (6) Equation (43) is valid for interactions with constant form factors.
  * Jochim  _et al._ (2003) S. Jochim, M. Bartenstein, A. Altmeyer, G. Hendl, C. Chin, J. H. Denschlag, and R. Grimm, Phys. Rev. Lett. 91, 240402 (2003).
  * Note (7) The analysis is not limited to a fixed atomic species or resonance and can be applied generally for varying values of $k_{\mathrm{n}}R^{*}$. In our simulations, we consider a condensate of ${}^{85}\text{Rb}$ atoms at an experimentally relevant density of $n=4\times 10^{12}\tmspace+{.1667em}\text{cm}^{-3}$, such that $t_{\mathrm{n}}=69.84\tmspace+{.1667em}\mu s$. We model the Feshbach resonance located at $B_{0}=155.04$ G with magnetic field width $\Delta B=10.7$ G and background scattering length $a_{\mathrm{bg}}=-443\tmspace+{.1667em}a_{0}$ Claussen  _et al._ (2003) and quench the magnetic field to resonance in a time span of $\tau_{in}=0.072\tmspace+{.1667em}t_{\mathrm{n}}$, such that the quench is effectively sudden. We mark the completion of the quench as $t=0$ and start our analysis with this broad-resonance system, where $k_{\mathrm{n}}R^{*}\approx 10^{-3}$, in order to verify the consistency of the coupled-channels model with the single-channel model as presented in Ref.+.2777emColussi  _et al._ (2018b) and then gradually increase the value of $k_{\mathrm{n}}R^{*}$.
  * Claussen  _et al._ (2003) N. R. Claussen, S. J. J. M. F. Kokkelmans, S. T. Thompson, E. A. Donley, E. Hodby, and C. E. Wieman, Phys. Rev. A 67, 060701 (2003).

##  Appendix A  Broad resonance limit

In the vacuum broad resonance limit, the two-body interactions can be
described by a single channel model, since the number of closed channel
molecules remains small. In the many-body environment, the broad resonance
limit of Eqs. (22)-(25) can be obtained by eliminating the molecular wave
function $\psi_{m}$ from this set of equations, meaning that
$\partial_{t}\psi_{m}\approx 0$, such that

$\displaystyle\psi_{m}=-\frac{g}{2\nu}\psi_{a}^{2}\zeta^{*}(0)-\frac{g}{2\nu}\sum_{\mathbf{k}\neq
0}\zeta^{*}(2\mathbf{k})\kappa_{\mathbf{k}}.$ (52)

The reduced set of equations then corresponds to the single channel HFB
equations as presented in Refs. Colussi  _et al._ (2018b); Musolino  _et al._
(2019), with

$v_{\mathrm{1ch}}=v-\frac{\absolutevalue{g}^{2}}{2\nu},$ (53)

where $v_{\mathrm{1ch}}$ is the renormalized single channel potential
interaction strength Kokkelmans  _et al._ (2002); Kokkelmans and Holland
(2002). Additionally, in a single-channel model, the two-body transition
matrix can be expressed as

$T_{\mathrm{1ch}}=v_{\mathrm{1ch}}\zeta(2\mathbf{k})\braket{\zeta}{\psi_{\mathrm{P}}^{+}},$
(54)

Apart from the difference in the potential strength, the previous expression
corresponds to the uncoupled open channel part of the transition matrix as
presented in Eq. (12). In the broad resonance limit, Eq. (54) must match the
expression for the two channel transition matrix as given in Eq. (12), such
that

$v_{\mathrm{1ch}}\zeta(2\mathbf{k})\braket{\zeta}{\psi_{P,1ch}^{+}}=v\zeta(2\mathbf{k})\braket{\zeta}{\psi_{P}^{+}}+\frac{\frac{g^{2}}{2}\absolutevalue{\braket{\zeta}{\psi_{P}^{+}}}^{2}}{E-\nu-\frac{g^{2}}{2}\braket{\zeta}{\hat{G}_{P}(E)}{\zeta}},$
(55)

where $\braket{\zeta}{\psi_{P}^{+}}$ can be expressed as

$\braket{\zeta}{\psi_{P}^{+}}=\zeta(2\mathbf{k})^{*}\left(1+\frac{v\xi(E)}{1+4\pi
v\xi(E)}\right),$ (56)

with

$\xi(E)=\frac{1}{(2\pi)^{3}}\int\absolutevalue{\zeta(2\mathbf{k}^{\prime})}^{2}\frac{k^{\prime
2}}{E-\hbar^{2}k^{\prime 2}/m}dk^{\prime}.$ (57)

Furthermore, we can write the term $\braket{\zeta}{\hat{G}_{P}(E)}{\zeta}$
that appears in the denominator of Eq. (55) as

$\braket{\zeta}{\hat{G}_{P}(E)}{\zeta}=4\pi\xi(E)\left(1+\frac{v\xi(E)}{1+4\pi
v\xi(E)}\right).$ (58)

Substituting the expression for the potential strength $v$ as presented in Eq.
(10) into Eq. (55), defining $v_{\mathrm{1ch}}$ analogously to $v$ and
considering the low-energy limit, we can express the single channel scattering
length $a_{\mathrm{1ch}}$ as

$a_{\mathrm{1ch}}=a_{\mathrm{bg}}-\frac{m}{4\pi\hbar^{2}}\frac{g_{0}^{2}}{2\nu_{0}}.$
(59)

Comparing Eq. (59) to the definition of $a_{\mathrm{eff}}$ in Ref. Snyder  _et
al._ (2012), we recognize that the correction factor $g_{c}$ used in that work
has been replaced by a factor $2$ in our definition of $a_{\mathrm{1ch}}$.
There, the correction factor $g_{c}$ was set to a value of 1.816 in order to
match the binding energy of the contact potential model as closely as
possible. A similar parametrization and correction factor was found in Ref.
Kokkelmans and Holland (2002). We recognize that, using the separable
potential model, our analytically obtained factor of $2$ is in relatively
close agreement to the calibrated correction factors used in these works.

Extending the broad-resonance limit analysis to the two-body contact, we
realize that this limit allows for the computation of the contact using the
adiabatic sweep theorem Tan (2008b); Werner  _et al._ (2009)

$\displaystyle\mathcal{C}_{2}=\frac{4\pi m}{\hbar^{2}}\frac{dE}{d(-1/a)}.$
(60)

Applying the dispersive relation between the scattering length and the
magnetic field as presented in Eq. (1), the previous expression can be
rewritten into the form

$\displaystyle\mathcal{C}_{2}=\frac{4\pi m}{\delta\mu
R^{*}(1-a_{\mathrm{bg}}/a)^{2}}\frac{dE}{dB}.$ (61)

Since only the bare bound-state energy of the closed-channel molecule depends
on the magnetic field, the application of the Helmann-Feynmann theorem Werner
_et al._ (2009); Braaten  _et al._ (2008) reveals that the two-body contact
can be expressed as Werner  _et al._ (2009); Castin and Werner (2012)

$\displaystyle\mathcal{C}_{2}=\frac{4\pi\absolutevalue{\psi_{m}}^{2}}{R^{*}}\left(1-\frac{a_{\mathrm{bg}}}{a}\right)^{-2}.$
(62)

The previous expression implies that the two-body contact scales linearly with
the molecular condensate fraction in the broad-resonance limit.

##  Appendix B  Renormalization of the pair condensate

As outlined in Sec. II.2.2, we cannot directly relate the eigenvalue
$N_{0}^{(2)}$ as presented in Eq. (33) to the number of condensed pairs due to
the Bose enhancement that alters the bosonic commutation relation and results
in an over counting of the number of pairs. As presented in Fig. 8, the
unrenormalized number of pairs exceeds the number of particles available for
pairing $N-N_{a}$, such that $N_{0}^{(2)}+N_{a}>N$. The overcounting of bosons
at early times is larger for broader resonances due to the increased value of
the open-channel pairing matrix $\kappa_{\mathbf{\mathrm{k}}}$.

Figure 8: The unrenormalized pair fraction $N_{0}^{(2)}$ added to the atomic
condensate fraction as a function of the time $t/t_{\mathrm{n}}$ spent at
unitarity.

In order to retrieve the number of condensed pairs and the correct bosonic
commutation relations, the eigenvalue $N_{0}^{(2)}$ has to be renormalized. To
this extend we define the following momentum space composite operator that
annihilates a condensed pair

$\displaystyle\hat{d}_{0}=\frac{1}{\sqrt{2}}\sum_{k\neq
0}\begin{bmatrix}\hat{a}_{k}\hat{a}_{-k}\varphi_{P,0}(k)\\\
\sqrt{2}\phi(k)\hat{b}_{0}\varphi_{Q,0}(k)\end{bmatrix},$ (63)

such that $\braket{d_{0}^{\dagger}d_{0}}=N_{0}^{(2)}/2$. By calculating the
expectation value of the commutator $[\hat{d}_{0},\hat{d}^{\dagger}_{0}]$, we
can compute the overcounting factor required for the renormalization. Using
Eq. (63) we find that

$\displaystyle\left[\hat{d}_{0},\hat{d}^{\dagger}_{0}\right]$
$\displaystyle=\frac{1}{2}\sum_{k,k^{\prime}\neq
0}\left(\varphi_{P,0}(k)\varphi_{P,0}(k^{\prime})[\hat{a}_{k}\hat{a}_{-k},\hat{a}^{\dagger}_{k^{\prime}}\hat{a}^{\dagger}_{-k^{\prime}}]\right.$
$\displaystyle\left.+2\varphi_{Q,0}(k)\varphi_{Q,0}(k^{\prime})\phi(k)\phi(k^{\prime})[\hat{b}_{0},\hat{b}^{\dagger}_{0}]\right),$
(64)

where we have used that the open and closed channel operators commute. Since
both operators satsify the bosonic commutation relation and using the
orthonormality of the bare closed channel wave function, the expectation value
of Eq. (B) can be expressed as

$\displaystyle\left\langle\left[\hat{d}_{0},\hat{d}^{\dagger}_{0}\right]\right\rangle=\sum_{k\neq
0}\frac{\absolutevalue{\braket{\hat{a}_{k}\hat{a}_{-k}}}^{2}}{N_{0}^{(2)}}(1+2\hat{a}_{k}^{\dagger}\hat{a}_{k})+2\frac{\absolutevalue{\braket{\hat{b}_{0}}}^{2}}{N_{0}^{(2)}},$
(65)

where we have used the momentum space analogue of the pair wave function
presented in Eq. (32) in order to express the eigenvector components in terms
of the expectation value of the cumulants introduced in Sec. II.2. Recognizing
the factor $N_{0}^{(2)}$ as defined in Eq. (33), we can rewrite Eq. (66) as

$\displaystyle\left\langle\left[\hat{d}_{0},\hat{d}^{\dagger}_{0}\right]\right\rangle=1+\frac{2}{N_{0}^{(2)}}\sum_{k\neq
0}\absolutevalue{\braket{\hat{a}_{k}\hat{a}_{-k}}}^{2}\hat{a}_{k}^{\dagger}\hat{a}_{k}.$
(66)

Since we neglect the presence of background molecules in the closed channel,
Eq. (66) is identical to its single channel analogue presented in Ref.
Musolino  _et al._ (2022), despite the presence of the closed channel
molecular fraction in the definition of $N_{0}^{(2)}$ presented in Eq. (33).
In order to retreive the desired bosonic commutation for the pair condensate,
Eq. (66) inspires us to define the following renormalized analogue of the pair
operator presented in Eq. (63)

$\displaystyle\hat{D}_{0}=\frac{\hat{d}_{0}}{\sqrt{1+\left[2/N_{0}^{(2)}\right]\sum_{k\neq
0}\absolutevalue{\braket{a_{k}a_{-k}}}^{2}a_{k}^{\dagger}a_{k}}}$ (67)

Using this renormalized operator, we can now compute the number of condensed
pairs using $\braket{D_{0}^{\dagger}D_{0}}$ for various values of the
resonance width (which in the many-body system can be parametrized in terms of
$k_{n}R^{*}$).

##  Appendix C  Relating the contact to cumulants

In order to relate Eq. (39) to the set of cumulants introduced in Sec. II.2
such that we can compute its value as a function of the time spent at
unitarity, we consider the zero energy limit of Eq. (23), finding that

$\displaystyle\hat{\phi}_{\mathrm{Q}}=-\frac{g^{2}}{2\nu}\hat{\psi}_{\mathrm{P}}\hat{\psi}_{\mathrm{P}},$
(68)

where we have used
$\braket{\hat{\psi}_{P}}=\zeta(0)\psi_{a}^{2}+\sum_{\mathbf{k}}\zeta(2\mathbf{k})\kappa_{\mathbf{k}}$
and $\braket{\hat{\phi}_{\mathrm{Q}}}=\psi_{m}$. Analogous to Ref. Braaten
_et al._ (2008) we can then define the following compound operator
$\hat{\Phi}$

$\displaystyle\hat{\Phi}=v\hat{\psi}_{\mathrm{P}}\hat{\psi}_{\mathrm{P}}+g\hat{\phi}_{\mathrm{Q}},$
(69)

such that we can obtain the following expression for the generalized Tan
contact

$\displaystyle\mathcal{C}_{2}=\braket{\hat{\Phi}^{\dagger}\hat{\Phi}}.$ (70)

In terms of our set of cumulants the previous expression shows that the two-
channel contact can be computed as

$\begin{split}\mathcal{C}_{2}(t)=&\frac{m^{2}v^{2}}{\hbar^{4}}\left(\absolutevalue{\psi_{a}}^{4}\absolutevalue{\zeta(0)}^{2}+4\absolutevalue{\psi_{a}}^{2}\sum_{\mathbf{k}}\absolutevalue{\zeta(\mathbf{k})}^{2}\rho_{\mathbf{k}}+2\sum_{\mathbf{k},\mathbf{k}^{\prime}}\absolutevalue{\zeta(\mathbf{k}-\mathbf{k}^{\prime})}^{2}\rho_{\mathbf{k}}\rho_{\mathbf{k}^{\prime}}+\sum_{\mathbf{k}}\left[\zeta(2\mathbf{k})\zeta^{*}(0)\kappa_{\mathbf{k}}^{*}\psi_{a}^{2}+\text{h.c.}\right]\right.\\\
&\left.+\absolutevalue{\sum_{\mathbf{k}}\zeta^{*}(\mathbf{k})\kappa_{\mathbf{k}}}^{2}+\frac{1}{\sqrt{2\pi
a_{\mathrm{bg}}^{2}R^{*}}}\left[(\psi^{\dagger}_{a})^{2}\psi_{m}\zeta(0)+\sum_{\mathbf{k}}\zeta(2\mathbf{k})\kappa_{\mathbf{k}}^{*}\psi_{m}+\text{h.c.}\right]+\frac{1}{2\pi
a_{\mathrm{bg}}^{2}R^{*}}\absolutevalue{\psi_{m}}^{2}\right).\end{split}$ (71)

Equation (71) is used to obtain the results as presented in Sec. III.2.