# Dynamical formation of polarons in a Bose-Einstein Condensate: A variational
approach

L. A. Peña Ardila  Institut für Theoretische Physik, Leibniz Universität,
30167 Hannover, Germany

###### Abstract

We investigate the non-equilibrium dynamics of an impurity coupled to a Bose-
Einstein condensate, comparing systematically with recent experimental results
(arXiv:2005.00424). The dynamics of the impurity is tracked down by using a
time-dependent variational coherent ansatz. For weak coupling between the
impurity and the bath, analytical expressions for the time-dependent overlap
(or contrast) are derived, matching quite well with previous findings obtained
within a Master Equation Approach (MEA). For strong coupling instead, the
variational ansatz provides a good quantitative description of the polaron
dynamics, in particular, in signaling the transition from the few to the many-
body correlated regime where polarons are expected to form.

##  I I. Introduction

One of the most challenging problems in current quantum science and technology
has to do with providing a complete description of a many-body system. Besides
the large number of degrees of freedom, interactions and correlations makes
the problem even harder to be tackled. The concept of quasiparticle arises as
an alternative for mapping this complex problem into a more tractable one. In
fact, quasiparticles are crucial for understanding several phenomena ranging
from atomic, nuclear, to high-energy physics. One of the most famous paradigms
of quasiparticles relies on Pekar and Landau’s original proposal of the
polaron concept Landau and Pekar (1948). When an electron travels through a
material, it gets dressed by the low-energy excitations of the lattice
(phonons), forming a polaron with renormalised energy and mass Feynman (1955).
Polarons give insight into transport properties in semiconductors Devreese and
Peters (1984), electrical conduction in polymer chains Mahani et al. (2017),
spin transport in organic materials Gershenson et al. (2006); Watanabe et al.
(2014) and superconductivity Lee et al. (2006). Moreover, they can be used as
a "probe" for highly correlated quantum many-body environments as well, for
instance, 3He impurities immersed into 4He superfluids Baym and Pethick
(1991).

Ultra-cold quantum gases offer a pristine platform where interaction between
atoms can be tuned at will Bloch et al. (2008); Chin et al. (2010), allowing
for the possibility to create polarons by immersing atomic impurities in
quantum gases. Here, we do not rely on a lattice, but the mechanism is
analogous. The impurity gets dressed by the low-energy excitations of the
quantum gas instead. Polarons can be created in a Bose-Einstein condensate
(Bose polarons) Jørgensen et al. (2016); Hu et al. (2016); Yan et al. (2020);
Camargo et al. (2018); Peña Ardila et al. (2019) or in a degenerate Fermi gas
(Fermi polarons) Schirotzek et al. (2009); Ngampruetikorn et al. (2012);
Koschorreck et al. (2012); Kohstall et al. (2012); Cetina et al. (2016a);
Scazza et al. (2017); Sidler et al. (2012). Albeit research on transport
properties using polarons in quantum gases is still premature, these quasi-
particles can be used for probing quantum Christensen et al. (2015) and
thermal effects Mehboudi et al. (2019) as well as topological invariants
Grusdt et al. (2019, 2016) in the host bath. However, the precise nature of
the dopant is fundamental as the many-body environment can be probed using
impurities, rather than polarons Bouton et al. (2020) and hence, it is
relevant to understand the time scales of polaron formation. Contrary to the
condensed matter scenario, the shortest response to collective excitations in
ultracold quantum gases are on the order of microseconds enabling current
state-of-the-art Ramsey techniques to measure these time scales with high
precision. Thus, the non-equilibrium dynamics of polaron in quantum gases has
attracted a lot of attention both in theory and experiments. In the former
case, the dynamics have been investigated using field-theory, renormalization
group, time-dependent variational and perturbative approaches Grusdt et al.
(2018); Volosniev et al. (2015); Shchadilova et al. (2016); Lampo et al.
(2018); Lausch et al. (2018); Nielsen et al. (2019); Liu et al. (2019);
Drescher et al. (2019); Mistakidis et al. (2019a, b), whereas recent
experiments shed light on quasiparticle formation and the non-equilibrium
dynamics of impurities which are quantified by using Ramsey interferometry
protocols Cetina et al. (2016b, a); Skou et al. (2020).

In this work we use a time-dependent variational method Lee et al. (1953);
Shchadilova et al. (2016) to study the real-time evolution of an impurity
embedded in a Bose-Einstein condensate. The non-equilibrium dynamics is
tracked down by measuring the time dependent overlap, namely the probability
to find the polaron at a certain time with respect to its initial non-
interacting state, $\left|\psi(0)\right\rangle$. Formally the overlap is
defined as
$S(t)=\left\langle\psi(0)\right|\exp[-i\hat{H}t]\left|\psi(0)\right\rangle$.
Within the variational approach this quantity is derived for arbitrary times
and coupling strength, which means that one can go beyond the truncated
Bogoliubov-Fröhlich Hamiltonian that is accurate only in the weakly
interacting regime.

The article is organized as follows. In Sec.II we describe the model and we
introduce the full Hamiltonian of the system. By using a low-energy
truncation, we obtain the Bogoliubov-Fröhlich Hamiltoinian which enables a
"system plus reservoir" description in terms of a MEA. In Sec.III we review
this perturbative formalism and introduce the observable of interest, i.e. the
contrast. In Sec.IV we use the time dependent variational approach to compute
the Euler-Lagrange equations of motion (EoM). In Sec.V the EoM can be solved
exactly in the weakly interacting regime and proper benchmark is done with
respect to the perturbative approach. In order to do comparisons with
experiments, additional decoherence effects must be included in the theory, in
Sec.VI we discuss trap dephasing and in Sec.VII decoherence by magnetic-field
fluctuations and losses. In Sec VIII we discuss some results within the
experimental available coupling strengths both in the weak and strongly
interacting regime. Comparisons yield a good agreement between the theory and
the experimental measurements Skou et al. (2020). Finally, conclusions are
drawn in Sec. IX.

##  II II. Model and Hamiltonian

We consider an ultradilute gas of impurities of mass $m_{I}$ immersed in a
Bose-Einstein condensate (BEC) of mass $m_{B}$ and density $n$. The
Hamiltonian of the system reads

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\sum_{\mathbf{p}}\Omega_{\mathbf{p}}\hat{c}_{\mathbf{p}}^{\dagger}\hat{c}_{\mathbf{p}}+\sum_{\mathbf{k}}\epsilon_{\mathbf{k}}\hat{a}_{\mathbf{k}}^{\dagger}\hat{a}_{\mathbf{k}}$
(1) $\displaystyle+$
$\displaystyle\frac{1}{2V}\sum_{\mathbf{k,k^{\prime},q}}V_{B}(\mathbf{q})\hat{a}_{\mathbf{k+q}}^{\dagger}\hat{a}_{\mathbf{k^{\prime}-q}}^{\dagger}\hat{a}_{\mathbf{k^{\prime}}}\hat{a}_{\mathbf{k}}$
$\displaystyle+$
$\displaystyle\frac{1}{V}\sum_{\mathbf{k,k^{\prime},q}}V_{I}(\mathbf{q})\hat{a}_{\mathbf{k+q}}^{\dagger}\hat{c}_{\mathbf{k^{\prime}-q}}^{\dagger}\hat{c}_{\mathbf{k^{\prime}}}\hat{a}_{\mathbf{k}},$

where the operator $\hat{c}_{\mathbf{p}}$ ($\hat{c}_{\mathbf{p}}^{\dagger}$)
annihilates (creates) an impurity of momentum $\mathbf{p}$ and energy
$\Omega_{\mathbf{p}}=\mathbf{p}^{2}/2m_{I}$, whereas $\hat{a}_{\mathbf{k}}$
($\hat{a}_{\mathbf{k}}^{\dagger}$) annihilates (creates) a boson with momentum
$\hbar\mathbf{k}$ and energy
$\epsilon_{\mathbf{k}}=\hbar^{2}\mathbf{k^{2}/}2m_{B}$. Moreover,
$V_{B}(\mathbf{q})=T_{B}$ and $V_{I}(\mathbf{q})=T_{\nu}$ are the Fourier
transform of the short-ranged boson-boson and impurity-boson potentials
respectively. In addition, we define $T_{\nu}=2\pi\hbar^{2}a/m_{red}$ as the
zero-energy impurity-bath coupling constant depending on the tunable s-wave
scattering length $a$ and $m_{red}=m_{B}^{-1}+m_{I}^{-1}$ the reduced mass. We
perform the Bogoliubov transformation on the bosonic operators, i.e.
$\hat{a}_{\mathbf{k}}=u_{k}\hat{b}_{\mathbf{k}}-v_{k}^{*}\hat{b}_{-\mathbf{k}}^{\dagger}$,
with amplitudes
$u_{k}=1+v_{k}^{2}=\left(\epsilon_{\mathbf{k}}+T_{B}n+\omega_{k}\right)/2\omega_{k}$
and $v_{k}=-T_{B}n/2\omega_{k}$, being the bosonic coupling strength
$T_{B}=4\pi\hbar^{2}a_{B}/m_{B}$ and $a_{B}$ the s-wave boson-boson scattering
length. Moreover,
$\omega_{\mathbf{k}}=\frac{\hbar^{2}k}{2m_{B}\xi}\sqrt{(k\xi)^{2}+2}$ is the
dispersion relation of the bosonic bath written in terms of the healing length
$\xi=1/\sqrt{8\pi na_{B}}$. Thus the transformed hamiltonian in the single
impurity limit (characterized by its position operator $\hat{\mathbf{R}}$ )
yields

$\displaystyle\mathcal{H}$ $\displaystyle=$
$\displaystyle\mathcal{H}_{F}+\frac{T_{\nu}}{V}\sum_{\mathbf{k,q}}V_{\mathbf{k,q}}^{(+)}e^{i\left(\mathbf{k-q}\right)\hat{\cdot\mathbf{R}}}\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{q}}$
(2) $\displaystyle+$
$\displaystyle\frac{T_{\nu}}{V}\sum_{\mathbf{k,q}}V_{\mathbf{k,q}}^{(-)}e^{i\left(\mathbf{k+q}\right)\hat{\cdot\mathbf{R}}}\left(\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{q}}^{\dagger}+\hat{b}_{\mathbf{-k}}\hat{b}_{\mathbf{-q}}\right).$

In order to capture all order of multi-scattering process in the problem,
quadratic terms in the Bogoliubov operators must be also included. Moreover,
$V_{\mathbf{k,q}}^{(\pm)}=\frac{1}{2}\left[W_{\mathbf{k}}W_{\mathbf{q}}\pm\left(W_{\mathbf{k}}W_{\mathbf{q}}\right)^{-1}\right]$
with the coupling matrix
$W_{\mathbf{k}}=\left[\frac{(k\xi)^{2}}{(k\xi)^{2}+2}\right]^{1/4}$ and
$\mathcal{H}_{F}$ reads

$\displaystyle\mathcal{H}_{F}$ $\displaystyle=$
$\displaystyle\frac{\hat{\mathbf{p}}^{2}}{2m_{I}}+\sum_{\mathbf{k}}\hbar\omega_{\mathbf{k}}\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{k}}+nT_{\nu}$
(3) $\displaystyle+$
$\displaystyle\frac{T_{\nu}\sqrt{n}}{\sqrt{V}}\sum_{\mathbf{k}}W_{\mathbf{k}}e^{i\mathbf{k}\hat{\cdot\mathbf{R}}}\left(\hat{b}_{\mathbf{k}}^{\dagger}+\hat{b}_{\mathbf{-k}}\right).$

The truncated Hamiltonian Eq. (3) is known as the Bogoliubov-Fröhlich
Hamiltonian which provides the ground-state properties of a single impurity
weakly interacting with the low-energy excitations of the quantum gas. This
process is represented by the scattering of an impurity with a single-
excitation of the condensate. However, this Hamiltonian is inaccurate for
strong coupling, which is characterized by high energy multi-scattering
process featuring a transition from low-energy excitations to bare bosons and
enabling thus the formation of many-body bound states once the resonance is
crossed Peña Ardila et al. (2019); Shchadilova et al. (2016). Hence, beyond-
Fröhlich terms in Eq. (2) are needed for arbitrary coupling strength in order
to describe the physics properly.

Lee-Low-pines (LLP) transformation.– To simplify the problem further, one can
translate the impurity-bath system into a new set of coordinates by evoking
the celebrated LLP canonical transformation
$\mathcal{H}=\hat{U}^{\dagger}\mathcal{H}\hat{U}$ where
$\hat{U}=\exp\left[i\hat{\mathbf{R}}\cdot\sum_{\mathbf{k}}\mathbf{k}\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{k}}\right]$
Lee et al. (1953); Shchadilova et al. (2016); Chatterjee and mukhopadhyay
(2019). The transformed Hamiltonian depends on the total momentum of the
system,
$U^{\dagger}\mathbf{\hat{P}}U=\mathbf{\hat{P}}-\sum_{\mathbf{k}}\mathbf{k}\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{k}}$.
The Hamiltonian Eq. (2) in the new reference frame reads

$\displaystyle\mathcal{H}$ $\displaystyle=$ $\displaystyle
T_{\nu}n+\frac{1}{2m_{I}}\mathbf{\left({\mathbf{P}}-\sum_{\mathbf{k}}\mathbf{k}\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{k}}\right)}^{2}+\sum_{\mathbf{k}}\omega_{\mathbf{k}}\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{k}}$
(4) $\displaystyle+$ $\displaystyle
T_{\nu}\frac{\sqrt{n}}{\sqrt{V}}\sum_{\mathbf{k}}W_{\mathbf{k}}\left(\hat{b}_{\mathbf{k}}^{\dagger}+\hat{b}_{-\mathbf{k}}\right)+\frac{T_{\nu}}{V}\sum_{\mathbf{k,q}}V_{\mathbf{k,q}}^{(+)}\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{q}}$
$\displaystyle+$
$\displaystyle\frac{T_{\nu}}{V}\sum_{\mathbf{k,q}}V_{\mathbf{k,q}}^{(-)}\left(\hat{b}_{\mathbf{k}}^{\dagger}\hat{b}_{\mathbf{q}}^{\dagger}+\hat{b}_{\mathbf{-k}}\hat{b}_{\mathbf{-q}}\right).$

In the new impurity frame of reference, the total momentum commutes with the
transformed Hamiltonian and is a conserved quantity. Hence the bare impurity
degrees of freedom are fully eliminated and the expectation value of the
impurity momentum can be set to zero, which is consistent with the polaron
ground-state of the system. One more detail that should be taken with care in
the Hamiltonian Eq. (4) is the renormalization of the impurity-boson coupling
strength as the contact interaction potential is modeled by an nonphysical
$\delta$-potential in real space. From the Lippmann-Schwinger equation the
zero-energy coupling constant and the s-wave scattering length are related via

$\displaystyle
T_{\nu}^{-1}=\frac{m_{red}}{2\pi\hbar^{2}a}-\frac{1}{V}\sum_{\mathbf{k}}^{\Lambda\rightarrow\infty}\frac{2m_{red}}{\hbar^{2}\mathbf{k}^{2}},$
(5)

##  III III. Master Equation approach

At $t=0$ the total density matrix of the whole system can be factorized as the
product of the density matrix of a small subsystem (impurities) and the host
reservoir (bath). In the spirit of the Born approximation, when the
interaction system-bath is turned on, correlations play an important role
giving place to decoherence of the polaron and hence the density matrix of the
system displays deviations on the order of the system-reservoir coupling
strength Carmichael (1999). This situation can be described accurately within
a MEA, provided that interactions (impurity-bath) are weak enough Nielsen et
al. (2019). In fact, Eq. (3) follows a system-reservoir kind of Hamiltonian.
In order to characterize the coherence of the system we study the time-
dependent overlap (also known as Ramsey overlap), $S(t)=\left\langle
0\left|\hat{c}(t)\rho_{I}(t)\right|0\right\rangle$, where $\rho_{I}(t)$ is the
reduced density matrix of the impurity obtained explicitly with the MEA. In
reference Nielsen et al. (2019) a detailed derivation of $\rho_{I}(t)$ is
discussed. Moreover, $\hat{c}(t)$ depicts the impurity operators and
$\left|0\right\rangle$ is the vacuum of phonons. Thus within the Fröhlich
model, the decoherence for an uniform system yields Nielsen et al. (2019)

$\displaystyle
S^{h}(t)=\exp\left[-\frac{i}{\hbar}E_{pol}^{F}t+i\frac{1}{8}\frac{a_{B}}{a_{0}}F(t/t_{n})\left(\frac{a}{a_{B}}\right)^{2}\right],$
(6)

with $t_{n}=m/8\pi\hbar na_{B}$ and the function
$F(t)=t+i\left(it+3\right)\left(\frac{1+i}{2}\right)\sqrt{\pi
t}\exp\left(it/2\right)\mathrm{erfc}\left[\left(\frac{1+i}{2}\right)\sqrt{t}\right]$.
The upper index $``h"$ indicates the contrast computed for a homogeneous
system and

$E_{pol}^{F}=\frac{8\pi}{(6\pi^{2})^{2/3}}(na_{B}^{3})^{1/3}\left[\frac{a}{a_{B}}+\frac{a_{B}}{a_{0}}\left(\frac{a}{a_{B}}\right)^{2}\right]\frac{\hbar^{2}k_{n}^{2}}{2m}$
(7)

is the polaron energy in the Fröhlich regime. Here
$k_{n}=\left(6\pi^{2}n\right)^{1/3}$ and
$\frac{a_{B}}{a_{0}}=\frac{32}{3\sqrt{\pi}}\sqrt{na_{B}^{3}}$. In addition,
for short-time, namely $t\ll\hbar/E_{\mathrm{MF}}$ (being
$E_{\mathrm{MF}}=T_{\nu}n$ the mean-field polaron energy) the contrast yields
a particular form

$\displaystyle S^{0}(t)\sim
1-(i+1)\sqrt{\frac{t}{t_{\Omega}}}-\frac{i}{\hbar}E_{\mathrm{MF}}t,$ (8)

with the characteristic time scale $t_{\Omega}=\frac{m_{B}}{32\pi\hbar
n^{2}a^{4}}$ which depends on the shape of the impurity-bath potential. Both
theoretical approaches addressed in this work assumes a quenched impurity.
Universality of the very short-time dynamics and the corresponding time scales
are discussed in Skou et al. (2020) and they are not within the scope of the
present work.

##  IV IV. VARIATIONAL ANSATZ

In this section we employ a time-dependent variational ansatz to study polaron
dynamics as previously used in references Shchadilova et al. (2016); Lee et
al. (1953); Drescher et al. (2019). On one hand we benchmark the results
obtained within a variational ansatz against a MEA. On the other hand the
variational formalism enables us to access to the strongly interacting regime
ruled by the full Hamiltonian Eq. (4). The ansatz is based on a coherent state
product

$\displaystyle\left|\psi_{c}(t)\right\rangle=\exp[-i\phi(t)]\hat{O}(t)\left|0\right\rangle,$
(9)

where
$\hat{O}(t)=\exp\left(\frac{1}{\sqrt{V}}\sum_{\mathbf{k}}\beta_{\mathbf{k}}(t)\hat{b}_{\mathbf{k}}^{\dagger}-h.c\right)$
and $\left|0\right\rangle=\left|0_{\mathbf{k}},\mathbf{p}\right\rangle$ is the
total state at $t=0$ before the impurity-boson interaction is quenched, i.e.
an impurity with momentum $\mathbf{p}$ and the vacuum of phonons. Here
$\gamma_{\mathbf{k}}(t)=\left\\{\beta_{\mathbf{k}}(t),\phi(t)\right\\}$ plays
the role of time-dependent variational parameters. The displacement operator
$\hat{O}$ obeys
$\hat{O}^{\dagger}(t)\hat{b}_{\mathbf{k}}\hat{O}(t)=\hat{b}_{\mathbf{k}}+\beta_{\mathbf{k}}$.
These properties are heavily used to compute the equations of motion and
expectation values. Note that, the ansatz Eq. (9) is exact in the limit of an
immobile impurity and non-interacting bosons. Following Shchadilova et al.
(2016), the Euler Lagrange equations,
$\frac{d}{dt}\frac{\partial\mathcal{L}}{\partial\dot{\gamma}}-\frac{\partial\mathcal{L}}{\partial\gamma}=0$,
with the Lagrangian density,
$\mathcal{L}=\left\langle\psi_{c}(t)\left|i\hbar\partial_{t}-\mathcal{H}\right|\psi_{c}(t)\right\rangle$
and for the manifold of variational parameters $\gamma_{\mathbf{k}}(t)$ reads

$\displaystyle i\hbar\dot{\beta}_{\mathbf{k}}$ $\displaystyle=$ $\displaystyle
T_{\nu}W_{\mathbf{k}}\sqrt{\frac{n}{V}}+\Gamma(\mathbf{k})\beta_{\mathbf{k}}+\frac{T_{\nu}}{2V}\left[W_{\mathbf{k}}\sum_{\mathbf{k}^{\prime}}W_{\mathbf{k^{\prime}}}\left(\beta_{\mathbf{k}^{\prime}}+\beta_{\mathbf{k}^{\prime}}^{*}\right)\right.$
(10) $\displaystyle+$
$\displaystyle\left.\frac{1}{W_{\mathbf{k}}}\sum_{\mathbf{k}^{\prime}}\frac{1}{W_{\mathbf{k^{\prime}}}}\left(\beta_{\mathbf{k}^{\prime}}-\beta_{\mathbf{k}^{\prime}}^{*}\right)\right]$

and

$\displaystyle\hbar\partial_{t}\phi=T_{\nu}n+\frac{T_{\nu}\sqrt{n}}{2\sqrt{V}}\sum_{\mathbf{k}}W_{\mathbf{k}}\left(\beta_{\mathbf{k}}+\beta_{\mathbf{k}}^{*}\right)-\frac{\mathbf{P}_{B}^{2}}{2m_{I}}.$
(11)

Here
$\Gamma(\mathbf{k})=\omega_{\mathbf{k}}+\Omega_{\mathbf{k}}+\hbar\mathbf{k}\cdot\mathbf{P}_{B}/m_{I}$
and
$\mathbf{P}_{B}=\sum_{\mathbf{k}}\mathbf{k}\beta_{\mathbf{k}}^{*}\beta_{\mathbf{k}}$.
As mentioned previously and without loss of generality, we restrict to the
case of slow motion impurities, i.e. $\mathbf{P}=0$. Moreover, note that
$\mathbf{P_{B}}$ is always zero for all times due to spherical symmetry.
Solving the coupled differential equations will allow us to compute quantities
of interest such as the contrast

$S^{0}(t)=\left\langle\psi_{c}(0)|\psi_{c}(t)\right\rangle,$ (12)

where $\left|\psi_{c}(t)\right\rangle$ is given by Eq. (9). Carrying out the
calculation and after some algebra one ends up with a simple expression

$\displaystyle
S^{0}(t)=\exp\left[-i\phi(t)-\frac{1}{2}\sum_{\mathbf{k}}\left|\beta_{\mathbf{k}}(t)\right|^{2}\right].$
(13)

In the stationary case, i.e. $\dot{\gamma}_{\mathbf{k}}(t)=0$, Eq. (10) can be
split into a real and an imaginary part and hence one obtains

$\displaystyle\left(\begin{array}[]{c}\mathrm{Re}[\beta_{\mathbf{k}}]\\\
\mathrm{Im}[\beta_{\mathbf{k}}]\end{array}\right)=\left(\begin{array}[]{c}-\frac{W_{\mathbf{k}}\sqrt{n}}{\Omega_{\mathbf{k}}}\left(T_{\nu}^{-1}+\sum_{\mathbf{k^{\prime}}}\frac{\left|W_{\mathbf{k^{\prime}}}\right|^{2}}{\Omega_{\mathbf{k^{\prime}}}}\right)^{-1}\\\
0\end{array}\right).$ (18)

Thus the polaron energy in the stationary case can be computed by
$E_{pol}=\left\langle\psi_{st}\left|\mathcal{\widehat{H}}\right|\psi_{st}\right\rangle$
yielding

$\displaystyle
E_{pol}=\frac{8\pi(na_{B}^{3})^{1/3}}{\left(6\pi^{2}\right)^{2/3}}\left[\frac{a_{B}}{a}-\frac{a_{B}}{a_{0}}\right]^{-1}\frac{\hbar^{2}k_{n}^{2}}{2m_{B}}.$
(19)

Alternatively, one can check that for long times
$\left|S(t\rightarrow\infty)\right|\rightarrow
Z=\exp\left(-\frac{1}{2V}\sum_{\mathbf{k}}\left|\beta_{\mathbf{k}}\right|^{2}\right)$,
which coincides with the amplitude of the contrast for the homogeneous case
and in particular at weak coupling, one recovers the quasiparticle residue
predicted by perturbative approaches Christensen et al. (2015). As
aforementioned the expression for the energy Eq. (19) is in principle valid
for arbitrary coupling strength, in fact one recovers the perturbative result
Eq. (7) for $a/\xi\ll 1$ . In order to check the quality of the ansatz against
another theories and the experiment, in Fig. 1 we compare the polaron energy
using recent experimental data and quantum Monte-Carlo methods Peña Ardila et
al. (2019) against the variational result in Eq. (19). The agreement is quite
good within statistical error, even for large values of $1/k_{n}a$.
Nonetheless, near to the unitary regime differences might be significant. In
fact, the variational ansatz relies on the Bogoliubov approximation which is
well known to be less accurate close to the resonance, however for time scales
far enough from polaron stationary state the variational method still provides
a quantitatively good description of the impurity decoherence.

Figure 1: (a) polaron energy at equilibrium as a function of the coupling
strength $1/k_{n}a$ using quantum-Monte Carlo (blue square) and experimental
data (orange circle) from Peña Ardila et al. (2019), second-order perturbation
theory (red curve) and time-dependent variational ansatz (green curve). (b)
example of instantaneous polaron energy
$E_{Pol}(t)=\left\langle\psi_{c}(t)\left|\mathcal{H}\right|\psi_{c}(t)\right\rangle$
as a function of time for a coupling strength $1/k_{n}a=-5$. The polaron is
expected to reach equilibrium with an energy $E_{Pol}$ (see Eq. (19)).

##  V V. Decoherence-weak coupling

The coupled system of Eqs. (10) and (11) are computed using the full
Hamiltonian in Eq. (4) and should be solved numerically. However, for weak-
coupling one can use the truncated Fröhlich Hamiltonian and the EoM are
reduced to

$\displaystyle i\hbar\partial_{t}\beta_{\mathbf{k}}$ $\displaystyle=$
$\displaystyle\frac{T_{\nu}}{\sqrt{V}}\sqrt{n}W_{\mathbf{k}}+\Gamma(\mathbf{k})\beta_{\mathbf{k}}$
$\displaystyle\hbar\partial_{t}\phi$ $\displaystyle=$ $\displaystyle
T_{\nu}n+\frac{T_{\nu}\sqrt{n}}{2\sqrt{V}}\sum_{\mathbf{k}}W_{\mathbf{k}}\left(\beta_{\mathbf{k}}+\beta_{\mathbf{k}}^{*}\right).$
(20)

The previous equations can be solved exactly. Notice that, in the mean-field
case the solution for the contrast is trivial since $\beta_{\mathbf{k}}(t)=0$
and $\phi(t)=\frac{T_{\nu}n}{\hbar}t$, hence the contrast reads

$S^{h}(t)=\exp\left[-iT_{\nu}nt/\hbar\right].$ (21)

In order to recover beyond mean-field terms, the strategy is to solve the
first equation in Eq. (20) for $\beta_{\mathbf{k}}$. It yields
$\beta_{\mathbf{k}}(t)=\frac{T_{\nu}}{\sqrt{V}}\sqrt{n}\frac{W_{k}}{\Gamma(\mathbf{k})}\left[\exp\left(-i\Gamma(\mathbf{k})t/\hbar\right)-1\right]$
for the initial condition $\beta_{\mathbf{k}}(0)=0$ and subsequently this
solution is plugged into the $\phi-$equation, obtaining thus,
$\phi=\frac{1}{\hbar}\int_{0}^{t}ds\left[T_{\nu}n+\frac{T_{\nu}\sqrt{n}}{2\sqrt{V}}\sum_{\mathbf{k}}W_{\mathbf{k}}\left(\beta_{\mathbf{k}}(s)+\beta_{\mathbf{k}}^{*}(s)\right)\right]$.
In addition, one should take into account the proper renormalization of the
impurity-boson coupling strength, $T_{\nu}$, by means of Eq. (5). After some
algebra a closed form for both $\beta_{\mathbf{k}}(t)$ and $\phi(t)$ is found.
Plugging the previous solutions into the definition of the contrast Eq. (13),
one finally gets
$S^{h}(t)=\exp\left[-iE_{pol}^{F}t+i\frac{1}{8}\frac{a_{B}}{a_{0}}\left(\frac{a}{a_{B}}\right)^{2}F(t/t_{n})\right]$,
which coincides with the very same expression found by the MEA in Eq. (6). The
variational approach captures surprisingly the quantitative behaviour of the
polaron dynamics out of equilibrium. This conclusion yields interesting
consequences. In reference Nielsen et al. (2019) the master equation for the
impurity problem was derived upon two strong assumptions. First, the Born
approximation, which was genuinely justified on the fact that the impurity
weakly perturbs the host medium. Second, for this system a highly degree of
Markovianity is expected Lampo et al. (2018). Within the MEA at zero
temperature, this latter approximation is valid due to the sufficiently large
difference in time scales. Thus, the variational approach indirectly justifies
the use of the Markovian-approximation in this regime.

##  VI VI. Trap dephasing

Figure 2: (a) Contrast amplitude for $1/k_{n}a=-5.0$. The experimental points
from Skou et al. (2020) are represented by orange circles. The blue crosses
are the second-order perturbative result Eq. (23) which overlaps fully with
the green curve depicting the variational ansatz. The purple dashed line
depicts for the short-time behaviour characterised by the time scale
$\sqrt{t_{n}}$ (defined in the mean text). On top of the trap dephasing
decoherence, the variational approach also includes decoherence by losses
represented by the red dashed line and finally magnetic-field fluctuations
depicted by the dashed black line (MFF) (b) zoom-in of Fig. 2(a) in order to
distinguish better all theories for short and intermediate times. The MFF
fitting experimental parameter here is $\Delta=1800$Hz.

If we consider the simplest case at the mean-field level, the amplitude of the
contrast Eq. (21) is one for all times. In order to do proper comparison with
the experiments, the contrast must be averaged over a non-homogeneous density
profile, i.e. dephasing by the trap. In the mean-field case and using Eq. (21)
one finds

$\displaystyle S(t)=\frac{1}{N_{B}}\int
d^{3}\mathbf{r}n(\mathbf{r})\exp\left[-iT_{\nu}n(\mathbf{r})\frac{t}{\hbar}\right],$
(22)

being $n(\mathbf{r})=n_{\mathrm{TF}}(\mathbf{r})$ the condensate density (see
Appendix A), whereas the first beyond-mean field correction reads

$\displaystyle S(t)=\frac{1}{N_{B}}\int
d^{3}\mathbf{r}n(\mathbf{r})S(\mathbf{r},t),$ (23)

with

$S(\mathbf{r},t)=\exp\left[-\frac{i}{\hbar}E_{pol}^{F}\left[n(\mathbf{r})\right]t+i\frac{1}{8}\frac{a_{B}}{a_{0}}F(t/t_{n(\mathbf{r})})\left(\frac{a}{a_{B}}\right)^{2}\right],$
(24)

containing the inhomogeneity of the condensate in terms of the local time
scale $t_{n(\mathbf{r})}=\frac{m_{B}}{8\pi n(\mathbf{r})a_{B}\hbar}$.

##  VII VII. Dephasing by losses and magnetic field fluctuations (MFF)

Inelastic losses arising from three-body recombination process can strongly
influence the contrast and give rise to decoherence by losses. The
experimental impurity loss rate as a function of the coupling strength has
been measured in Skou et al. (2020) using the function
$\frac{\hbar\Gamma_{\mathrm{loss}}}{E_{n}}=A\exp(\frac{B}{k_{n}a})+C$, with
fitting parameters $A=0.1788$, $B=1.0130$ and $C=0.0061$. In addition,
$\hbar/E_{n}=2m/\left(\hbar^{3/2}6\pi^{2}n\right)^{2/3}$. On the other hand,
the Ramsey interferometry protocol is sensible to magnetic field variations
for each experimental measurement and hence MFF is another important source of
decoherence Skou et al. (2020). Therefore, the total contrast needs to include
also these additional inherent effects from the experiment. Thus, the total
contrast is modified as

$\displaystyle S(t)=\frac{1}{N_{B}}S^{\mathrm{Loss}}(t)S^{\mathrm{MFF}}(t)\int
d^{3}\mathbf{r}n(\mathbf{r})S(\mathbf{r},t),$ (25)

with

$\displaystyle S^{\mathrm{Loss}}(t)=\exp\left[-\Gamma_{\mathrm{Loss}}t\right]$
(26)

and

$\displaystyle
S^{\mathrm{MFF}}(t)=\frac{\int_{-\infty}^{\infty}d\phi\exp\left[-i\phi-\frac{\phi^{2}}{8\pi^{2}\varDelta^{2}t^{2}}\right]}{\sqrt{8\pi^{3}\varDelta^{2}t^{2}}},$
(27)

being the decoherence associated to losses and MFF respectively. Instead,
$S(\mathbf{r},t)$ is the contrast including only the trap dephasing. From the
experiments, one observe that dephasing due to the trap and by losses are the
most important contributions on top of the dephasing by collisions. Additional
decoherence effects, such as the one due to magnetic field fluctuations (MFF)
are quite relevant in the weakly interacting regime and becomes less important
for strong coupling, albeit never negligible.

##  VIII VIII. Results

Figure 3: (a-b) Amplitude and phase of $S(t)$ for $1/(k_{n}a)=-1$. Orange
points represent the experimental measurements in Skou et al. (2020). The
green line depicts for the variational calculation including only trap
dephasing, whereas the dashed line includes the losses dephasing as well.
(c-d) Amplitude and phase of $S(t)$ at the unitary regime. In addition, the
cyan line indicates the fitting proposed in Eq. (29) and the black dashed line
depicts for the short-time universal behavior. The MFF fitting parameter here
is $\Delta=4600$ Hz.

Weak-coupling.–The time evolution of the contrast is plotted in Fig. 2. The
circles with error bars depict the recent experimental data Skou et al. (2020)
for a fixed coupling strength $1/k_{n}a=-5.0$, featuring a typical value in
the weakly interacting regime. The computed contrast using the perturbative
MEA is shown by the blue crosses (see also Eq. (23)), whereas the variational
ansatz calculation including only trap dephasing is represented by the green
solid curve. As expected both curves overlap according to our discussion in
Sec.V. In addition, the dashed purple line shows the power-law short-time
behavior (see also Eq. (8)) valid for $t_{d}<t\ll t_{n}=\hbar/E_{MF}\sim
50\mu$s for this particular coupling strength. Note that, the short time
prediction in Fig. 2 represents the initial high energy two-body dynamics
Nielsen et al. (2019) which is independent of the host bath interaction, i.e.
the impurity does not experience the condensate presence. Therefore
decoherence is originated only by impurity-boson interactions. This regime is
preceded by the polaron relaxation and formation stage, where the impurity
interact with the condensate excitations and many-body effects start to take
place ($t\approx t_{n}$). Instead, the formation time is expected for times
$t\gg t_{B}\sim 500\mu$s as revealed by the instantaneous polaron energy
defined as
$E_{Pol}(t)=\left\langle\psi_{c}(t)\left|\mathcal{H}\right|\psi_{c}(t)\right\rangle=\hbar\dot{\phi}$
and plotted in Fig. 1(b). A discussion on the degeneracy time scale
$t_{d}=ma^{2}/\hbar$ can be found in Braaten et al. (2010); Skou et al.
(2020). This time scale is set by the coupling strength and for this
particular case is on the order of $t_{d}=0.08\mu$s.

For very short-times $t_{d}<t<t_{n}$, the mean-field prediction already
provides a good qualitatively agreement with respect to the experimental
contrast, however the effect of quantum fluctuations (originated from
impurity-bath collisions) tends to shift the mean-field contrast amplitude for
larger times, following thus the experimental downward trend as shown in
Figures 2(a)-2(b), but still far from a fully quantitatively agreement with
the experimental data. For intermediate and larger times $t_{n}>100\mu$s the
contrast decays faster in the experiment with respect to the theoretical
prediction which only includes trap dephasing. In fact, additional decoherence
effects which are inherent to the experiment are absent in the theory and need
to be included, for instance, decoherence by losses and by MFF Skou et al.
(2020). We have included the losses dephasing as well as the MFF decoherence
as illustrated in Sec.VII. A quantitatively agreement with respect to the
experimental data is reached as shown by the black dashed line in Figures 2(a)
and 2(b). The effect of including losses dephasing is more evident in the
strongly interacting regime, in contrast to the dephasing due to MFF which is
less relevant.

Strong-coupling.–We investigate the time-dependent overlap in the strongly
interacting regime, i.e. $-1<\left(k_{n}a\right)^{-1}<0$. Here we rely solely
on numerical solutions of equations Eq. (10) and Eq. (11) with the
experimental parameters associated to each $1/k_{n}a$ Skou et al. (2020), in
particular for two characteristic coupling strengths, $1/k_{n}a=-1$ and at
unitarity. In the former case, a very good agreement between the variational
results including losses decoherence and MFF and the experimental data is
achieved as shown in Figures 3(a) and 3(b). In particular, for this regime,
MFF decoherence is small in comparison with the other experimental decoherence
effects, however it provides a quantitatively agreement with the experimental
contrast, whereas at unitarity is practically negligible.

Figure 4: Initial departure of the instantaneous polaron energy computed with
the variational approach (red dashed line) against the experimental
measurements (orange points) Skou et al. (2020) at unitarity. The blue region
is the QMC prediction for the polaron stationary energy Peña Ardila et al.
(2019).

In Figures 3(c) and 3(d) we show the results at the unitary regime. We observe
that both the initial relaxation obtained within the variational ansatz (green
line) and the short-time universal behavior Braaten et al. (2010), valid for
$t\ll t_{d}$ (black dashed line) agree with each other and it follows that

$\displaystyle S_{\mathrm{2body}}(t)\approx
1-(1-i)\frac{16}{9\pi^{3/2}}\left(\frac{t}{t_{n}}\right)^{3/2}.$ (28)

The time scale $t_{n}^{3/2}$ is linked with the universal two-body scattering
physics, where the impurity interacts with highly energetic bosons of the
condensate. Note that the same universal scaling is recovered for impurities
immersed in a Fermi sea Parish and Levinsen (2016) with a characteristic time
scale set by the Fermi time ($\hbar/E_{F}$). In order to bridge the two-body
with the many-body physics where the polaron starts its formation, we attempt
to find a functional form of the contrast for arbitrary times. Taking
advantage of the functional form of Eq. (28), we propose the following
expression for arbitrary time

$\displaystyle S(t)=\exp\left[S_{\mathrm{2body}}(t)-1\right],$ (29)

which recovers exactly Eq. (28) for short-times and spans all over one-body,
two-body and all high $n-$body correlations for larger times, thus signaling
the transition from few-body to the many-body correlated regime (around $t\sim
5\mu$s). The agreement between Eq. (29) depicted by the cyan curve in Fig.
3(c) and the full numerical calculations obtained with the variational ansatz
(green curve) is quite remarkable. Interestingly, in Fig. 4 one observes that
the polaron is expected to enter the equilibrium regime for times $t\gg 5\mu$s
with a polaron energy near to $E_{Pol}\approx-E_{n}$ at equilibrium. In
addition, the proposed function is heavily based on the initial form of the
ansatz Eq. (9) which encodes all correlations of the system. In fact, by
expanding out the wave-function Eq. (9) one has

$\displaystyle\left|\psi(t)\right\rangle$ $\displaystyle=$
$\displaystyle\exp\left(-i\phi(t)\sum_{\mathbf{k}}\left[\beta_{\mathbf{k}}(t)\hat{b}_{\mathbf{k}}^{\dagger}-\beta_{\mathbf{k}}^{*}(t)\hat{b}_{\mathbf{k}}\right]\right)\left|0\right\rangle$
$\displaystyle=$
$\displaystyle\sum_{n=0}^{\infty}\frac{\left(-i\phi(t)\left[\sum_{\mathbf{k}}\beta_{\mathbf{k}}(t)\hat{b}_{\mathbf{k}}^{\dagger}-\beta_{\mathbf{k}}^{*}(t)\hat{b}_{\mathbf{k}}\right]\right)^{n}}{n!}\left|0\right\rangle,$

and taking into account only the first terms of the series, one obtains

$\left|\psi(t)\right\rangle\simeq\left(1-i\phi(t)\sum_{\mathbf{k}}\beta_{\mathbf{k}}(t)\hat{b}_{\mathbf{k}}^{\dagger}\right)\left|0\right\rangle,$
(30)

which resembles the Chevy ansatz Li and Das Sarma (2014) where only a single
excitation on top of the unperturbed condensate is included. Eq. (30) is
accurate in the weakly interacting regime. Hence, perturbation theory and the
variational approach again consistently agree in this regime.

Finally we would like to stress that in the neighborhood of the unitarity
regime and for very long time (essentially the polaron state at equilibrium),
the variational approach might provide quantitatively results as in this
regime the Bogoliubov approximation breaks down. However, for our current work
the time scales are way far from equilibrium at unitarity and hence good
qualitatively results in the short-time dynamics holds.

##  IX IX. Conclusions

In this work we have employed a time-dependent variational approach to
investigate the out-of-equilibrium dynamics of impurities interacting with a
Bose-Einstein condensate. The main results are two-fold. First, the method is
benchmarked against time-dependent perturbative approaches such as MEA, valid
in the weakly interacting regime. Second, the method is compared with state-
of-the-art experimental data. Besides impurity-bath decoherence, agreement
between theory and experiment relies on the inclusion of decoherence effects
in the theoretical approaches due to trap inhomogeneity, losses and
fluctuations of magnetic field inherently incorporated in the experiment.
Close to the unitarity regime, we also bridged in a consistent and
quantitatively way the few-body physics characterized by the universal two-
body scattering process to the many-body process where polarons starts its
formation. Note that, apart from the initial experimental preparation of the
system, the contrast might be altered due to the comparable Efimov states
energy that are also relevant for polaron physics in bosonic systems Levinsen
et al. (2015). Also, how the polaron formation happens once the resonance is
crossed is an open question since a new time scale related with the two-body
bound state physics should play an important role. Real-time formation of
quasiparticles emerging from impurity-bath decoherence might be important not
only for gases, but also to other exotic many-body environments such as ultra-
dilute quantum liquids or recently investigated long ranged polarons such as
Rydberg, dipolar and ionic polarons. Bisset et al. (2020); Ardila and Pohl
(2018); Kain and Ling (2014); Camargo et al. (2018); Astrakharchik et al.
(2020).

## Acknowledgments

We gratefully acknowledge to Kristian K. Nielsen for discussions at earlier
stage of this work. We also thank to Magnus G. Skou, Nils Byg Jørgensen and
Jan Artl for insightful discussions and providing the experimental data. The
authors also acknowledge to Luis Santos, Fabian Grusdt, Artem Volosniev and
Giacomo Bighin for reading the manuscript and for the critical feedback.
Special Acknowledgments to Roberta Giusteri for critical comments on the
manuscript. This research was funded by the DFG Excellence Cluster
QuantumFrontiers.

##  X Appendix A: contrast for an inhomogeneous trap

The density of the condensate can be written in a very simple form as the
interaction energy scales are larger with respect to the kinetic energy one,
hence the density of the condensate within the Thomas-Fermi (TF) approximation
can be written as

$\displaystyle
n_{\mathrm{TF}}(\mathbf{r})=\left\\{\begin{array}[]{c}\frac{\mu_{\mathrm{TF}}-V(\mathbf{r})}{T_{B}},\quad\mu_{\mathrm{TF}}\geq
V(\mathbf{r})\\\ \\\ 0,\quad\mathrm{Otherwise},\end{array}\right.$

with $\mu_{\mathrm{TF}}$ the chemical potential in TF regime. In experiments
the trapping potential is modeled by
$V(\mathbf{r})=\frac{1}{2}m_{B}\sum_{i=1}^{3}\omega_{i}^{2}x_{i}^{2}=\frac{1}{2}m_{B}\overline{\omega}^{2}\sum_{i=1}^{3}\overline{x_{i}}^{2}$
where the geometric frequency is defined as
$\overline{\omega}=(\omega_{x}\omega_{y}\omega_{z})^{1/3}$ and
$\overline{x}_{i}=x_{i}\omega_{i}/\overline{\omega}$. Upon the condition
$\mu_{\mathrm{TF}}=V(\mathbf{r})$, one finds the size of the condensate
confined in this potential, namely the Thomas-Fermi radius
$\overline{R}=\sqrt{\frac{2\mu_{\mathrm{TF}}}{m_{B}\overline{\omega}^{2}}}$
and the chemical potential depending on the number of particles of the
condensate. Thus, from the normalization condition $N_{B}=\int
d^{3}\mathbf{r}n_{\mathrm{TF}}(\mathbf{r})$ one obtains
$\mu_{\mathrm{TF}}=\frac{\hbar\overline{\omega}}{2}\left[15N_{B}\frac{a_{B}}{a_{\mathrm{ho}}}\right]^{2/5}$.
The contrast for an inhomogeneous distribution can be written as

$\displaystyle S(t)$ $\displaystyle=$ $\displaystyle\frac{1}{N_{B}}\int
d^{3}\mathbf{r}n_{\mathrm{TF}}(\mathbf{r})S(\mathbf{r},t)$ (32)
$\displaystyle=$ $\displaystyle\frac{1}{N_{B}}\int
d^{3}\mathbf{r}\left[\frac{\mu_{\mathrm{TF}}-V(\mathbf{r})}{T_{B}}\right]S(\mathbf{r},t).$

The contrast $S(\mathbf{r},t)$ depends on time and density in each point of
space. In other words, the spatial dependence of the contrast is implicitly
encoded in the density, thus we substitute the homogeneous density
$n\rightarrow n_{\mathrm{TF}}(\mathbf{r})$ using the local density
approximation.

## References

  * Landau and Pekar (1948) L. Landau and S. Pekar, J. Exp. Theor. Phys 18, 419 (1948).
  * Feynman (1955) R. P. Feynman, Phys. Rev. 97, 660 (1955).
  * Devreese and Peters (1984) J. Devreese and F. Peters, _Polarons and Excitons in Polar Semiconductors and Ionic Crystals_ (Plenum Press, New York, 1984).
  * Mahani et al. (2017) M. R. Mahani, A. Mirsakiyeva, and A. Delin, J. Phys. Chem. C 121, 10317 (2017), URL http://dx.doi.org/10.1021/acs.jpcc.7b02368.
  * Gershenson et al. (2006) M. E. Gershenson, V. Podzorov, and A. F. Morpurgo, Rev. Mod. Phys. 78, 973 (2006), URL https://link.aps.org/doi/10.1103/RevModPhys.78.973.
  * Watanabe et al. (2014) S. Watanabe, K. Ando, K. Kang, S. Mooser, Y. Vaynzof, H. Kurebayashi, E. Saitoh, and H. Sirringhaus, Nat. Phys. 10, 308 (2014).
  * Lee et al. (2006) P. A. Lee, N. Nagaosa, and X.-G. Wen, Rev. Mod. Phys. 78, 17 (2006), URL https://link.aps.org/doi/10.1103/RevModPhys.78.17.
  * Baym and Pethick (1991) G. Baym and C. Pethick, _Landau Fermi-Liquid Theory: Concepts and Applications_ (Wiley-VCH-New York., 1991).
  * Bloch et al. (2008) I. Bloch, J. Dalibard, and W. Zwerger, Rev. Mod. Phys. 80, 885 (2008), URL https://link.aps.org/doi/10.1103/RevModPhys.80.885.
  * Chin et al. (2010) C. Chin, R. Grimm, P. Julienne, and E. Tiesinga, Rev. Mod. Phys. 82, 1225 (2010), URL https://link.aps.org/doi/10.1103/RevModPhys.82.1225.
  * Jørgensen et al. (2016) N. B. Jørgensen, L. Wacker, K. T. Skalmstang, M. M. Parish, J. Levinsen, R. S. Christensen, G. M. Bruun, and J. J. Arlt, Phys. Rev. Lett. 117, 055302 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.055302.
  * Hu et al. (2016) M.-G. Hu, M. J. Van de Graaff, D. Kedar, J. P. Corson, E. A. Cornell, and D. S. Jin, Phys. Rev. Lett. 117, 055301 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.055301.
  * Yan et al. (2020) Z. Z. Yan, Y. Ni, C. Robens, and M. W. Zwierlein, Science 368, 190 (2020), ISSN 0036-8075, eprint https://science.sciencemag.org/content/368/6487/190.full.pdf, URL https://science.sciencemag.org/content/368/6487/190.
  * Camargo et al. (2018) F. Camargo, R. Schmidt, J. D. Whalen, R. Ding, G. Woehl, S. Yoshida, J. Burgdörfer, F. B. Dunning, H. R. Sadeghpour, E. Demler, et al., Phys. Rev. Lett. 120, 083401 (2018), URL https://link.aps.org/doi/10.1103/PhysRevLett.120.083401.
  * Peña Ardila et al. (2019) L. A. Peña Ardila, N. B. Jørgensen, T. Pohl, S. Giorgini, G. M. Bruun, and J. J. Arlt, Phys. Rev. A 99, 063607 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.99.063607.
  * Schirotzek et al. (2009) A. Schirotzek, C.-H. Wu, A. Sommer, and M. W. Zwierlein, Phys. Rev. Lett. 102, 230402 (2009), URL https://link.aps.org/doi/10.1103/PhysRevLett.102.230402.
  * Ngampruetikorn et al. (2012) V. Ngampruetikorn, J. Levinsen, and M. M. Parish, EPL (Europhysics Letters) 98, 30005 (2012), URL https://doi.org/10.1209%2F0295-5075%2F98%2F30005.
  * Koschorreck et al. (2012) M. Koschorreck, D. Pertot, E. Vogt, B. Fröhlich, M. Feld, and M. Köhl, Nature 485, 619 (2012).
  * Kohstall et al. (2012) C. Kohstall, M. Zaccanti, M. Jag, A. Trenkwalder, P. Massignan, G. M. Bruun, F. Schreck, and R. Grimm, Nature 485, 615 (2012).
  * Cetina et al. (2016a) M. Cetina, M. Jag, R. S. Lous, I. Fritsche, J. T. Walraven, R. Grimm, J. Levinsen, M. M. Parish, R. Schmidt, M. Knap, et al., Science 354, 96 (2016a).
  * Scazza et al. (2017) F. Scazza, G. Valtolina, P. Massignan, A. Recati, A. Amico, A. Burchianti, C. Fort, M. Inguscio, M. Zaccanti, and G. Roati, Phys. Rev. Lett. 118, 083602 (2017), URL https://link.aps.org/doi/10.1103/PhysRevLett.118.083602.
  * Sidler et al. (2012) M. Sidler, P. Back, O. Cotlet, A. Srivastava, T. Fink, M. Kroner, E. Demler, and A. Imamoglu, Nature 485, 615 (2012).
  * Christensen et al. (2015) R. S. Christensen, J. Levinsen, and G. M. Bruun, Phys. Rev. Lett. 115, 160401 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.160401.
  * Mehboudi et al. (2019) M. Mehboudi, A. Lampo, C. Charalambous, L. A. Correa, M. A. García-March, and M. Lewenstein, Phys. Rev. Lett. 122, 030403 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.030403.
  * Grusdt et al. (2019) F. Grusdt, N. Y. Yao, and E. A. Demler, Phys. Rev. B 100, 075126 (2019), URL https://link.aps.org/doi/10.1103/PhysRevB.100.075126.
  * Grusdt et al. (2016) F. Grusdt, N. Yao, and D. Abanin, Nat Commun 7, 043614 (2016).
  * Bouton et al. (2020) Q. Bouton, J. Nettersheim, D. Adam, F. Schmidt, D. Mayer, T. Lausch, E. Tiemann, and A. Widera, Phys. Rev. X 10, 011018 (2020), URL https://link.aps.org/doi/10.1103/PhysRevX.10.011018.
  * Grusdt et al. (2018) F. Grusdt, K. Seetharam, Y. Shchadilova, and E. Demler, Phys. Rev. A 97, 033612 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.97.033612.
  * Volosniev et al. (2015) A. G. Volosniev, H.-W. Hammer, and N. T. Zinner, Phys. Rev. A 92, 023623 (2015), URL https://link.aps.org/doi/10.1103/PhysRevA.92.023623.
  * Shchadilova et al. (2016) Y. E. Shchadilova, R. Schmidt, F. Grusdt, and E. Demler, Phys. Rev. Lett. 117, 113002 (2016), URL https://link.aps.org/doi/10.1103/PhysRevLett.117.113002.
  * Lampo et al. (2018) A. Lampo, C. Charalambous, M. A. García-March, and M. Lewenstein, Phys. Rev. A 98, 063630 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.98.063630.
  * Lausch et al. (2018) T. Lausch, A. Widera, and M. Fleischhauer, Phys. Rev. A 97, 023621 (2018), URL https://link.aps.org/doi/10.1103/PhysRevA.97.023621.
  * Nielsen et al. (2019) K. K. Nielsen, L. A. P. Ardila, G. M. Bruun, and T. Pohl, New Journal of Physics (2019), URL http://iopscience.iop.org/10.1088/1367-2630/ab0a81.
  * Liu et al. (2019) W. E. Liu, J. Levinsen, and M. M. Parish, Phys. Rev. Lett. 122, 205301 (2019), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.205301.
  * Drescher et al. (2019) M. Drescher, M. Salmhofer, and T. Enss, Phys. Rev. A 99, 023601 (2019), URL https://link.aps.org/doi/10.1103/PhysRevA.99.023601.
  * Mistakidis et al. (2019a) S. I. Mistakidis, G. C. Katsimiga, G. M. Koutentakis, T. Busch, and P. Schmelcher, Phys. Rev. Lett. 122, 183001 (2019a), URL https://link.aps.org/doi/10.1103/PhysRevLett.122.183001.
  * Mistakidis et al. (2019b) S. I. Mistakidis, F. Grusdt, G. M. Koutentakis, and P. Schmelcher, New Journal of Physics 21, 103026 (2019b), URL https://doi.org/10.1088/1367-2630/ab4738.
  * Cetina et al. (2016b) M. Cetina, M. Jag, R. S. Lous, I. Fritsche, J. T. M. Walraven, R. Grimm, J. Levinsen, M. M. Parish, R. Schmidt, M. Knap, et al., Science 354, 96 (2016b), ISSN 0036-8075, eprint https://science.sciencemag.org/content/354/6308/96.full.pdf, URL https://science.sciencemag.org/content/354/6308/96.
  * Skou et al. (2020) M. G. Skou, T. G. Skov, N. B. Jørgensen, K. K. Nielsen, A. Camacho-Guardian, T. Pohl, G. M. Bruun, and J. J. Arlt, arXiv e-prints arXiv:2005.00424 (2020), eprint 2005.00424.
  * Lee et al. (1953) T. D. Lee, F. E. Low, and D. Pines, Phys. Rev. 90, 297 (1953), URL https://link.aps.org/doi/10.1103/PhysRev.90.297.
  * Chatterjee and mukhopadhyay (2019) A. Chatterjee and S. mukhopadhyay, _Polarons and Bipolarons_ (Taylor and Francis Group LLC, 2019).
  * Carmichael (1999) H. J. Carmichael, _Statistical Methods in Quantum Optics 1_ (Springer-Verlag Berlin Heidelberg, 1999).
  * Braaten et al. (2010) E. Braaten, D. Kang, and L. Platter, Phys. Rev. Lett. 104, 223004 (2010), URL https://link.aps.org/doi/10.1103/PhysRevLett.104.223004.
  * Parish and Levinsen (2016) M. M. Parish and J. Levinsen, Phys. Rev. B 94, 184303 (2016), URL https://link.aps.org/doi/10.1103/PhysRevB.94.184303.
  * Li and Das Sarma (2014) W. Li and S. Das Sarma, Phys. Rev. A 90, 013618 (2014), URL https://link.aps.org/doi/10.1103/PhysRevA.90.013618.
  * Levinsen et al. (2015) J. Levinsen, M. M. Parish, and G. M. Bruun, Phys. Rev. Lett. 115, 125302 (2015), URL https://link.aps.org/doi/10.1103/PhysRevLett.115.125302.
  * Bisset et al. (2020) R. N. Bisset, L. A. Peña Ardila, and L. Santos, arXiv e-prints arXiv:2007.00404 (2020), eprint 2007.00404.
  * Ardila and Pohl (2018) L. A. P. Ardila and T. Pohl, Journal of Physics B: Atomic, Molecular and Optical Physics 52, 015004 (2018), URL https://doi.org/10.1088%2F1361-6455%2Faaf35e.
  * Kain and Ling (2014) B. Kain and H. Y. Ling, Phys. Rev. A 89, 023612 (2014), URL https://link.aps.org/doi/10.1103/PhysRevA.89.023612.
  * Astrakharchik et al. (2020) G. E. Astrakharchik, L. A. Peña Ardila, R. Schmidt, K. Jachymski, and A. Negretti, arXiv e-prints arXiv:2005.12033 (2020), eprint 2005.12033.