# Dynamics of a trapped ion in a quantum gas: effects of particle statistics

Lorenzo Oghittu, Melf Johannsen  Zentrum für Optische Quantentechnologien,
Universität Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany  Rene
Gerritsma  Van der Waals Zeeman Institute, Institute of Physics, University of
Amsterdam, Science Park 904, 1098 XH Amsterdam, The Netherlands  Antonio
Negretti  Zentrum für Optische Quantentechnologien, Universität Hamburg,
Luruper Chaussee 149, 22761 Hamburg, Germany

###### Abstract

We study the quantum dynamics of an ion confined in a radiofrequency trap in
interaction with either a Bose or spin-polarized Fermi gas. To this end, we
derive quantum optical master equations in the limit of weak coupling and the
Lamb-Dicke approximations. For the bosonic bath, we also include the so-called
“Lamb-shift” correction to the ion trap due to the coupling to the quantum gas
as well as the extended Fröhlich interaction within the Bogolyubov
approximation that have been not considered in previous studies. We calculate
the ion kinetic energy for various atom-ion scattering lengths as well as gas
temperatures by considering the intrinsic micromotion and we analyse the
damping of the ion motion in the gas as a function of the gas temperature. We
find that the ion’s dynamics depends on the quantum statistics of the gas and
that a fermionic bath enables to attain lower ionic energies.

##  I Introduction

The interest in studying both experimentally and theoretically quantum
mixtures of ions and ultracold gases is increasing rapidly. Indeed, such a
compound system offers various fascinating perspectives both on fundamental
quantum few- and many-body physics research and on technological applications
that the two systems separately cannot afford. For instance, the exploration
of novel polaronic states [1, 2, 3] and quantum simulation of the electron-
phonon coupling [4, 5, 6]. For an extensive overview on atom-ion physics
research, we refer to the review [7, 8, 9]. Experimentally, a considerable
effort has been undertaken in the last few years in order to cool the atom-ion
compound system down to the quantum regime, namely when only s-wave atom-ion
collisions take place. To this end, three experimental approaches have been
pursed so far: ionization of a highly excited Rydberg atom in a Bose-Einstein
condensate [10, 11]; an ion confined in a radiofrequency (rf) trap interacting
with an optically trapped atomic gas [12, 13]; sympathetic cooling of ions and
atoms in optical dipole traps [14, 15]. In the former approach, the ion is not
trapped after ionisation and a controlled momentum kick via external electric
fields is imparted in order to investigate charge transport in a bosonic
medium [16, 17]. The second relies on the well-established laser cooling and
manipulation techniques of trapped ions, which are confined by means of a
combination of time-dependent and time-independent electric fields. The
exquisite control of the ion motion enables to prepare various non-classical
states [18, 19] and, in principle, to infer on environment properties by
reading out the ion quantum state. The time-dependent fields, however, can
seriously jeopardise the attainment of the ultracold atom-ion collisional
regime, a notorious issue named micromotion. Finally, the third approach is
somehow a hybrid of the former two, as it reduces drastically the impact of
ion micromotion in Paul traps and, at least in principle, it enables to
control the ion motion optically. At the same time, however, since the optical
trap is not as deep as the rf-trap, unavoidable stray fields can reduce the
ion lifetime in the trap [20, 21]. Moreover, an intense laser light can induce
undesired chemical reactions between the ion and the particles of the gas.

Here, we focus our attention on the second approach and investigate the
cooling dynamics of a trapped ion immersed in either a bosonic or fermionic
environment. Albeit the s-wave regime in hybrid traps has been attained with
only fermionic atoms so far [12, 22], several ongoing experiments involve
bosonic ensembles [23, 24, 25, 26, 27, 28, 14]. Up until now, however, a few
theoretical studies have been undertaken in order to assess the impact of ion
micromotion on the atom-ion quantum dynamics: a quantum mechanical calculation
in one dimension (1D) [29, 30], a semiclassical investigation of confinement-
induced resonances in quasi-1D [31, 32], and a 3D master equation analysis
[33]. Apart from the latter, all others concerned a single trapped atom and
ion, and therefore only the emerging two-body physics has been investigated.
In Ref. [33], however, the Fröhlich model was assumed, while the Lamb-shift
and, specifically for the bosonic case, the contribution of the non-condensed
fraction were not considered. Here, by Lamb-shift we refer to the
renormalisation of the ion trap parameters due to the coupling to the quantum
gas. In condensed-matter and for a free impurity such a shift is named polaron
shift, whereas in this work we use the quantum optics terminology. The
aforementioned studies confirmed that the impact of the ion micromotion can be
mitigated by choosing a small atom-ion mass ratio. With the present study, we
first aim at developing a formalism for mobile quantum impurities based on an
open quantum system approach, which does not rely on the (linear) Fröhlich
impurity-bath interaction model and rotating-wave approximation. Moreover, we
want to understand the role of the gas quantum statistics on the ion cooling
dynamics and whether atom-ion pairs different than Li/Yb+ can reach the s-wave
limit. As far as the first objective is concerned, we find that the
corrections causing the Lamb-shift yield additional dissipative contributions
in the master equation, as a consequence of the non-applicability of the
rotating-wave approximation. Furthermore, the quantum statistics of the
ultracold gas affects significantly the ion dynamics. While for temperatures
larger than the Fermi temperature $\mathcal{T}_{F}$ and the critical
temperature of condensation $\mathcal{T}^{0}_{\mathrm{c}}$ the ion dynamics
reproduces essentially the one corresponding to the interaction with a buffer
gas, at low temperature the fermionic and bosonic nature of the gas is
observable in a distinct gas temperature dependence of the ion energy. Both
for the fermionic and bosonic gas we observe a minimum in the ion energy for a
temperature $\mathcal{T}_{\min}$ close to $\mathcal{T}^{0}_{\mathrm{c}}$.
Below $\mathcal{T}_{\min}$ the ion energy increases again, but only marginally
for fermions. On the other hand, for a bosonic gas the ion energy dependence
on the gas temperature $\mathcal{T}$ exhibits a more rich structure. Indeed,
after the enhancement of the ion energy for
$\mathcal{T}^{0}_{\mathrm{c}}<\mathcal{T}<\mathcal{T}_{\min}$, below
$\mathcal{T}^{0}_{\mathrm{c}}$ it presents a maximum due to the interplay
between the condensate contribution and that of the normal part of the gas.
Moreover, the damping rate of the ion motion exhibits the same dependence on
the gas temperature as the condensate fraction for
$\mathcal{T}<\mathcal{T}^{0}_{\mathrm{c}}$. In addition, we elucidate the role
of the long-range character of the atom-ion polarisation potential on the ion
quantum dynamics by comparing our findings with those of the pseudopotential.
Contrary to neutral impurities, the non-equilibrium dynamics of the ion in the
quantum gas is non-universal, that is, it cannot be uniquely characterised by
the impurity-gas scattering length like for neutral impurities. Since there is
no clear separation of length scales in the many-body problem at typical gas
densities, the tail of the atom-ion interaction plays a crucial role in the
ion dynamics and the effective range corrections cannot be neglected.

The paper is organised as follows: In Sec. II we introduce the atom-ion
potential, while in Sec. III the system plus bath description is outlined,
which is kept on purpose quite general. A master equation including the
contribution of the non-condensed fraction is obtained in a way that is valid
for both neutral and charged impurities. In Sec. IV we focus on the trapped
ion case and make use of the well-known Lamb-Dicke approximation, which
enables us to further simplify the description of the ion dynamics. We
continue with Sec. V by providing the equations of motion of observables of
interest, while in Sec. VI we present our results. In Sec. VII we draw our
conclusions and provide an outlook for future work.

##  II Atom-ion interaction potential

Figure 1:  (Color online). Schematic view of the open quantum system with the
main length scales. The width $l$ of the Paul trap potential (represented by
the saddle) corresponds to the size of the ion in the ground state (large blue
wave) and is (much) smaller than the two other lengths: the average distance
$\bar{d}$ among the gas particles (small red waves), which is defined by
atomic density $n$; the characteristic length $R^{\star}$ of the atom-ion
potential (black thick line).

In this section we discuss briefly the interaction between an atom and an ion
and how we model it for the master equation calculation.

Polarisation potential. The interaction between an atom and an ion in free
space is described asymptothically by ($r\equiv|\mathbf{r}|$)

$\displaystyle V_{\mathrm{ai}}(\mathbf{r})=-\frac{C_{4}}{r^{4}}$ (1)

with $C_{4}=\frac{\alpha e^{2}}{2}\frac{1}{4\pi\epsilon_{0}}$ (in SI units)
111We note that for lithium additive and nonadditive interaction coefficients
for the three-body Li-Li-Li+ system have been recently computed [67]., where
$\alpha$ is the static polarisability of the atom, $e$ is the elementary
electronic charge, and $\epsilon_{0}$ the vacuum permittivity. Here, $r$
denotes the separation between the atom and the ion. The potential is
characterised by the length $R^{\star}=(2\mu C_{4}/\hbar^{2})^{1/2}$ and
energy $E^{\star}=\hbar^{2}/[2\mu(R^{\star})^{2}]$ scales, with $\mu=mM/(m+M)$
the reduced mass, $m$ the atom mass, and $M$ the ion mass.

Length scales. In our setting there are several relevant length scales (see
Fig. 1). Firstly, the aforementioned $R^{\star}$, which is typically in the
range of hundreds of nanometers and gives, as a rule of thumb, the order of
magnitude of the 3D zero-energy s-wave atom-ion scattering length (see also
Refs. [35, 36]). For instance, for the atom-ion pair 7Li/174Yb+ we have
$R^{\star}\simeq 75.15\,$nm, for 23Na/174Yb+ $R^{\star}\simeq 129.85\,$nm, and
for 87Rb/174Yb+ $R^{\star}\simeq 307.23\,$nm. Secondly, the atom-atom
(background) scattering length $a_{\mathrm{aa}}^{s}$ which is typically on the
order of a few nanometers. Thirdly, for a condensate, we have the healing
length, which is defined as $\xi=(8\pi na_{\mathrm{aa}}^{s})^{-1/2}$ with $n$
being the gas density. For a typical gas density $n=10^{14}/\mathrm{cm}^{3}$
we have, e.g., $\xi\simeq 273.61\,$nm for 87Rb, and $\xi\simeq 380.38\,$nm for
23Na. Instead, for a Fermi gas, the inverse of the Fermi wave vector
$\lambda_{F}=2\pi/k_{F}=(3\pi^{2}n)^{-1/3}$ gives another length scale. For
$n=10^{14}/\mathrm{cm}^{3}$, we have $\lambda_{F}\simeq 437.56$ nm. Fourthly,
the mean path length $\bar{d}$ that at the aforementioned typical gas density
is about $215.44\,$nm. Hence, all lengths are comparable and therefore there
is no separation of length scales as in the neutral counterpart. As a
consequence, non-universal behaviour in the ion statical and dynamical
properties is expected (see, for instance, Refs. [37, 38, 39] for a static ion
analysis). Thus, either very low atomic densities are considered, where a
universal behaviour is expected, or else the long-range tail of the atom-ion
interaction strongly matters.
Finally, the last relevant length for a trapped ion system is the ion trap
length $l$, which corresponds to the ion ground state size. This length is
about a few tens of nms for 174Yb+, as we discuss in Sec. IV. While $l$ is
rather small compared to $R^{\star}$ for heavy atoms, for lithium it is
roughly half of the corresponding $R^{\star}$. This means that scattering of
an atom and an ion should not be analysed in free space, as we do in the next
paragraph, as the confinement affects the atom-ion collision as for neutrals
in waveguides. Here, however, we neglect effects like confinement-induced
resonances as a dedicated study of such a phenomenology is required.

Quantum regime condition. It is important to note that the energy $E^{\star}$
sets the onset of s-wave atom-ion collisions, namely for energies smaller than
$E^{\star}$ the quantum regime can be attained 222We note that this is a
necessary condition, but not yet a sufficient one, as trap-shaped induced
resonances can still occur.. Indeed, $E^{\star}$ corresponds to the height of
the centrifugal barrier for the $\ell=1$ partial wave from threshold (see,
e.g., Fig. 1 of Ref. [9]). Assuming that the kinetic energy of the atom is
negligible, since it is ultracold compared to the ion, the collisional energy
in the relative atom-ion coordinate frame is given by [41, 42]

$\displaystyle
E_{\mathrm{coll}}=k_{\mathrm{B}}\mathcal{T}_{\mathrm{coll}}\simeq\frac{\mu}{M}E_{\mathrm{kin}}$
(2)

with $E_{\mathrm{kin}}=k_{\mathrm{B}}\mathcal{T}_{\mathrm{kin}}$ the ion’s
average kinetic energy in the laboratory frame of reference and
$k_{\mathrm{B}}$ being the Boltzmann constant 333The relation (2) is obtained
easily as follows. Since what matters is the relative atom-ion kinetic energy,
the collision energy is thus given by
$E_{\mathrm{coll}}=\mu\upsilon_{\mathrm{rel}}^{2}/2=\mu(\bm{\upsilon}_{\mathrm{a}}-\bm{\upsilon}_{\mathrm{i}})^{2}/2=\mu/2(\upsilon_{\mathrm{i}}^{2}+\upsilon_{\mathrm{i}}^{2}-2\bm{\upsilon}_{\mathrm{a}}\cdot\bm{\upsilon}_{\mathrm{i}})$.
Here, $\upsilon_{\mathrm{i}}\equiv|\bm{\upsilon}_{\mathrm{i}}|$
($\upsilon_{\mathrm{a}}\equiv|\bm{\upsilon}_{\mathrm{a}}|$) denotes the
magnitude of the ion (atom) velocity. Now, if
$\upsilon_{\mathrm{i}}\gg\upsilon_{\mathrm{a}}$, that is, the atoms are very
slow (i.e. ultracold) compared to the ion micromotion, then
$E_{\mathrm{coll}}\simeq\mu\upsilon_{\mathrm{i}}^{2}/2=\mu E_{\mathrm{kin}}/M$
with $E_{\mathrm{kin}}=M\upsilon_{\mathrm{i}}^{2}/2$.. Hence, in order to
enter the quantum regime of s-wave collisions, the ion’s kinetic energy must
be smaller than $E^{\star}$, that is,

$\displaystyle E_{\mathrm{kin}}\ll
E_{\mathrm{s}}=\frac{M}{\mu}E^{\star}=\left(1+\frac{M}{m}\right)E^{\star}.$
(3)

In case of a light atom and a heavy ion we have $\mu\simeq m$ so that
$M/\mu\gg 1$, and thus a significant gain in the upper limit for s-wave
collisions is obtained. For example, for 7Li/174Yb+ we have
$E_{\mathrm{s}}/k_{\mathrm{B}}\simeq 164.26\,\mu$K and for 23Na/174Yb+ we find
$E_{\mathrm{s}}/k_{\mathrm{B}}\simeq 6.07\,\mu$K. This shows that there is a
rather broad range of temperatures and that these are at least an order of
magnitude smaller than those of ultracold neutral collisions (on the mK
scale).

Regularized potential. Given the fact that the potential (1) is singular and
that later in the master equation we need to compute the Fourier transform of
the atom-ion potential, we introduce the regularization [33]

$\displaystyle
V_{\mathrm{ai}}^{r}(\mathbf{r})=-C_{4}\frac{r^{2}-c^{2}}{r^{2}+c^{2}}\frac{1}{(b^{2}+r^{2})^{2}}.$
(4)

Here, $b$ and $c$ are tuneable parameters that have units of a length and
control the energy spectrum of the potential as well as the atom-ion
scattering length. The Fourier transform of (4) is linked to the scattering
amplitude in the first-order Born approximation, which is defined as

$\displaystyle
f(q)=-\frac{\mu}{2\pi\hbar^{2}}\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{r}\,e^{i\mathbf{q}\cdot\mathbf{r}}V_{\mathrm{ai}}^{r}(\mathbf{r}).$
(5)

By using spherical coordinates and by integrating out the angular part, we
obtain

$\displaystyle f(q)$
$\displaystyle=-\frac{2\mu}{q\hbar^{2}}\int_{\mathbb{R}^{+}}\mathrm{d}r\,r\sin(qr)V_{\mathrm{ai}}^{r}(r)$
$\displaystyle=\frac{c^{2}\pi(R^{\star})^{2}}{(b^{2}-c^{2})^{2}q}\left\\{e^{-bq}\left[1+\frac{(b^{4}-c^{4})q}{4bc^{2}}\right]-e^{-cq}\right\\},$
(6)

where we used the fact that $(R^{\star})^{2}=2\mu C_{4}/\hbar^{2}$. The
determination of $b$ and $c$ is discussed in appendix A.

##  III Impurity master equation

In this section we provide an open system description of an impurity in a
quantum gas of either bosons or fermions by following the approach of Ref.
[44], where the impurity is described in first quantisation, whereas the
quantum bath in the second one. We focus mainly on the bosonic case, for which
we apply Bogolyubov theory, but we consider also the quadratic terms of the
bosonic field operators, which result in an extended Fröhlich interaction
Hamiltonian. The inclusion of this interaction has been proven to be crucial
in the description of the many-body response function of Rydberg [45] and Bose
polarons in vicinity of a Feshbach resonance [46]. The fermionic case is
considered only for a normal gas, i.e. not superfluid BCS theory, and it is
obtained as a special case of the master equation for a bosonic bath for gas
temperatures above the critical temperature of condensation. We note that in
the literature a master equation treatment of an impurity in a degenerate Bose
gas has already been undertaken [47, 48, 33, 49], but (i) only the (linear)
Fröhlich interaction has been considered and (ii) the Lamb-shift have been not
taken into account. Moreover and specifically for the ionic impurity, the
fermionic bath has been not investigated in Ref. [33].

###  III.1 System plus bath Hamiltonian

The total Hamiltonian of the system, the atomic impurity in interaction with a
bosonic bath, is given by $\hat{H}=\hat{H}_{S}+\hat{H}_{B}+\hat{H}_{BS}$,
where $\hat{H}_{S}$ is the impurity time-dependent Hamiltonian 444For the
specific case of the trapped ion, $\hat{H}_{S}\equiv\hat{H}_{I}^{trap}$ with
$\hat{H}_{I}^{trap}$ given in Appendix C.2,

$\displaystyle\hat{H}_{B}=\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{r}_{b}\,\hat{\Psi}_{b}^{\dagger}(\mathbf{r}_{b})\left[\frac{\hat{\mathbf{p}}_{b}^{2}}{2m}+\frac{g}{2}\hat{\Psi}_{b}^{\dagger}(\mathbf{r}_{b})\hat{\Psi}_{b}(\mathbf{r}_{b})\right]\hat{\Psi}_{b}(\mathbf{r}_{b}),$
(7)

and

$\displaystyle\hat{H}_{BS}=\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{r}_{b}\,\hat{\Psi}_{b}^{\dagger}(\mathbf{r}_{b})V_{\mathrm{ib}}(\mathbf{r}_{b}-\hat{\mathbf{r}})\hat{\Psi}_{b}(\mathbf{r}_{b}).$
(8)

Here, $V_{\mathrm{ib}}$ denotes the interaction between the impurity and a
particle of the bath. Besides, we assume that the bath is confined in a box of
length $L$ 555$L$ is assumed to be much larger than other length scale
involved in the system description. and that the interaction between the
bosons is given by a contact potential with coupling constant
$g=4\pi\hbar^{2}a_{\mathrm{bb}}^{s}/m$ and $a_{\mathrm{bb}}^{s}$ the 3D s-wave
atom-atom scattering length, while for the fermions we assume a spin-polarised
gas with no intra-particle interaction. Let us note that at this stage the
models describing the bosonic and fermionic baths are different, i.e. the
former is interacting while the latter is not. Nonetheless, and specifically
for the ion in a Paul trap, it turns out that the interaction among bosons is
much smaller than other energies involved, so that we can safely neglect it
(see Sec. IV for details). Hence, the differences we observe in the ionic
dynamics in the two baths are owed to their quantum statistics and not to the
interaction. For the sake of completeness, however, we keep here the
derivation of the master equation as general as possible.

The bosonic quantum field is expanded as

$\displaystyle\hat{\Psi}_{b}(\mathbf{r}_{b})=\sqrt{n_{0}}+\delta\hat{\Psi}_{b}(\mathbf{r}_{b})$
(9)

where $n_{0}=N_{0}/L^{3}$ is the density of the condensate, i.e. the zero
momentum component, while $N_{0}$ is the atom number. Fluctuations around the
condensate mode are described in terms of Bogolyubov modes

$\displaystyle\delta\hat{\Psi}_{b}(\mathbf{r}_{b})=L^{-3/2}\sum_{\mathbf{q}}u_{\mathbf{q}}\hat{b}_{\mathbf{q}}e^{i\mathbf{q}\cdot\mathbf{r}_{b}}+v_{\mathbf{q}}\hat{b}_{\mathbf{q}}^{\dagger}e^{-i\mathbf{q}\cdot\mathbf{r}_{b}}$
(10)

where
$[\hat{b}_{\mathbf{q}},\hat{b}_{\mathbf{q}^{\prime}}^{\dagger}]=\delta_{\mathbf{q},\mathbf{q}^{\prime}}$.
Using this expression in Eq. (7), we arrive at

$\displaystyle\hat{H}_{B}=E_{0}+\sum_{\mathbf{q}}\hbar\omega_{\mathbf{q}}\hat{b}_{\mathbf{q}}^{\dagger}\hat{b}_{\mathbf{q}}.$
(11)

Here, $E_{0}=gN_{0}^{2}/(2L^{3})$ is the ground state energy of the condensate
and $\mu_{\mathrm{G}}=\partial_{N_{0}}E_{0}=gn_{0}$ is the chemical potential
at zero temperature. The corresponding dispersion relation is given by [52]

$\displaystyle\varepsilon(\mathbf{q})=\hbar\omega_{\mathbf{q}}=\sqrt{\left(\frac{\hbar^{2}q^{2}}{2m}\right)^{2}+(\hbar
c_{s}q)^{2}}$ (12)

with $c_{s}=(gn_{0}/m)^{1/2}$ being the speed of sound, and the Bogolyubov
amplitudes $u_{\mathbf{q}}$ and $v_{\mathbf{q}}$ are given in Ref. [52]. Given
this, the atomic density operator is

$\displaystyle\hat{\Psi}_{b}^{\dagger}(\mathbf{r}_{b})\hat{\Psi}_{b}(\mathbf{r}_{b})=n_{0}+\Delta\hat{n}(\mathbf{r}_{b}).$
(13)

The first term provides a constant energy term in Eq. (8) for a homogeneous
gas, as we consider here, and it can be neglected, since it shifts merely the
energy minimum. The second term is given by

$\displaystyle\Delta\hat{n}(\mathbf{r}_{b})$
$\displaystyle=\hat{\Psi}_{b}^{\dagger}(\mathbf{r}_{b})\hat{\Psi}_{b}(\mathbf{r}_{b})-n_{0}$
$\displaystyle=\sqrt{n_{0}}[\delta\hat{\Psi}_{b}(\mathbf{r}_{b})+\delta\hat{\Psi}_{b}^{\dagger}(\mathbf{r}_{b})]+\delta\hat{\Psi}_{b}^{\dagger}(\mathbf{r}_{b})\delta\hat{\Psi}_{b}(\mathbf{r}_{b})$
$\displaystyle=\delta\hat{n}(\mathbf{r}_{b})+\delta^{2}\hat{n}(\mathbf{r}_{b})$
(14)

with
$\delta\hat{n}(\mathbf{r}_{b})=\sqrt{n_{0}}[\delta\hat{\Psi}_{b}(\mathbf{r}_{b})+\delta\hat{\Psi}_{b}^{\dagger}(\mathbf{r}_{b})]$.
Hence, we have

$\displaystyle\delta\hat{n}(\mathbf{r}_{b})=\sqrt{\frac{n_{0}}{L^{3}}}\sum_{\mathbf{q}}(u_{\mathbf{q}}+v^{*}_{\mathbf{q}})\hat{b}_{\mathbf{q}}e^{i\mathbf{q}\cdot\mathbf{r}_{b}}+(u^{*}_{\mathbf{q}}+v_{\mathbf{q}})\hat{b}_{\mathbf{q}}^{\dagger}e^{-i\mathbf{q}\cdot\mathbf{r}_{b}},$
(15)

and

$\displaystyle\delta^{2}\hat{n}(\mathbf{r}_{b})=L^{-3}\sum_{\mathbf{q},\mathbf{q}^{\prime}}\left[u^{*}_{\mathbf{q}}u_{\mathbf{q}^{\prime}}e^{-i(\mathbf{q}-\mathbf{q}^{\prime})\cdot\mathbf{r}_{b}}\hat{b}_{\mathbf{q}}^{\dagger}\hat{b}_{\mathbf{q}^{\prime}}\right.$
$\displaystyle\left.+u^{*}_{\mathbf{q}}v_{\mathbf{q}^{\prime}}e^{-i(\mathbf{q}+\mathbf{q}^{\prime})\cdot\mathbf{r}_{b}}\hat{b}_{\mathbf{q}}^{\dagger}\hat{b}_{\mathbf{q}^{\prime}}^{\dagger}+v^{*}_{\mathbf{q}}u_{\mathbf{q}^{\prime}}e^{i(\mathbf{q}+\mathbf{q}^{\prime})\cdot\mathbf{r}_{b}}\hat{b}_{\mathbf{q}}\hat{b}_{\mathbf{q}^{\prime}}\right.$
$\displaystyle\left.+v^{*}_{\mathbf{q}}v_{\mathbf{q}^{\prime}}e^{i(\mathbf{q}-\mathbf{q}^{\prime})\cdot\mathbf{r}_{b}}\hat{b}_{\mathbf{q}}\hat{b}_{\mathbf{q}^{\prime}}^{\dagger}\right].$
(16)

Given this, the system-bath Hamiltonian becomes

$\displaystyle\hat{H}_{BS}=\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{r}_{b}\,V_{\mathrm{ib}}(\mathbf{r}_{b}-\hat{\mathbf{r}})\Delta\hat{n}(\mathbf{r}_{b})$
$\displaystyle=\hbar\sum_{\mathbf{q}}\left(\hat{S}_{\mathbf{q}}\hat{\Gamma}_{\mathbf{q}}+\hat{S}^{{\dagger}}_{\mathbf{q}}\hat{\Gamma}_{\mathbf{q}}^{\dagger}\right)+\hbar\sum_{\mathbf{q},\mathbf{q}^{\prime}}\left(\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(u,u^{\prime})}\hat{\Gamma}_{\mathbf{q}}^{\dagger}\hat{\Gamma}_{\mathbf{q}^{\prime}}\right.$
$\displaystyle\left.+\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(u,v^{\prime})}\hat{\Gamma}_{\mathbf{q}}^{\dagger}\hat{\Gamma}_{\mathbf{q}^{\prime}}^{\dagger}+\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(v,u^{\prime})}\hat{\Gamma}_{\mathbf{q}}\hat{\Gamma}_{\mathbf{q}^{\prime}}+\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(v,v^{\prime})}\hat{\Gamma}_{\mathbf{q}}\hat{\Gamma}_{\mathbf{q}^{\prime}}^{\dagger}\right)$
$\displaystyle=\hat{H}_{BS}^{(1)}+\hat{H}_{BS}^{(2)},$ (17)

where we used the notation of Ref. [44]:
$\hat{\Gamma}_{\mathbf{q}}\equiv\hat{b}_{\mathbf{q}}$, and

$\displaystyle\hat{S}_{\mathbf{q}}$
$\displaystyle=\frac{\sqrt{nL^{3}}}{\hbar}(u_{\mathbf{q}}+v^{*}_{\mathbf{q}})e^{i\mathbf{q}\cdot\hat{\mathbf{r}}}c_{\mathbf{q}},$
$\displaystyle\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(u,u^{\prime})}$
$\displaystyle=\frac{u^{*}_{\mathbf{q}}u_{\mathbf{q}^{\prime}}}{\hbar}e^{i(\mathbf{q}^{\prime}-\mathbf{q})\cdot\hat{\mathbf{r}}}c_{\mathbf{q}^{\prime}-\mathbf{q}},$
$\displaystyle\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(u,v^{\prime})}$
$\displaystyle=\frac{u^{*}_{\mathbf{q}}v_{\mathbf{q}^{\prime}}}{\hbar}e^{-i(\mathbf{q}^{\prime}+\mathbf{q})\cdot\hat{\mathbf{r}}}c^{*}_{\mathbf{q}^{\prime}+\mathbf{q}},$
$\displaystyle\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(v,u^{\prime})}$
$\displaystyle=\frac{v^{*}_{\mathbf{q}}u_{\mathbf{q}^{\prime}}}{\hbar}e^{i(\mathbf{q}^{\prime}+\mathbf{q})\cdot\hat{\mathbf{r}}}c_{\mathbf{q}^{\prime}+\mathbf{q}},$
$\displaystyle\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(v,v^{\prime})}$
$\displaystyle=\frac{v^{*}_{\mathbf{q}}v_{\mathbf{q}^{\prime}}}{\hbar}e^{-i(\mathbf{q}^{\prime}-\mathbf{q})\cdot\hat{\mathbf{r}}}c^{*}_{\mathbf{q}^{\prime}-\mathbf{q}},$
$\displaystyle c_{\mathbf{q}}$
$\displaystyle=\frac{1}{L^{3}}\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{y}\,e^{i\mathbf{q}\cdot\mathbf{y}}V_{\mathrm{ib}}(\mathbf{y}).$
(18)

In addition, $\hat{H}_{BS}^{(1)}$ is the first sum over $\mathbf{q}$ in Eq.
(III.1) and it denotes the so-called Fröhlich model Hamiltonian in the context
of the electron-phonon coupling in solid-state [53], while
$\hat{H}_{BS}^{(2)}$ refers to the double sum over $\mathbf{q}$,
$\mathbf{q}^{\prime}$, which describes the so-called extended Fröhlich model
[54, 45]. Specifically for the ionic impurity, the coefficient
$c_{\mathbf{q}}$ is linked to the scattering amplitude (5) as

$\displaystyle c_{\mathbf{q}}=-\frac{2\pi\hbar^{2}}{\mu L^{3}}f(q).$ (19)

Apart from the Bogolyubov approximation, the expression (III.1) is exact for a
bosonic bath. For a normal Fermi gas, the interaction Hamiltonian reduces to

$\displaystyle\hat{H}_{BS}=\hbar\sum_{\mathbf{q},\mathbf{q}^{\prime}}\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}\hat{\Gamma}_{\mathbf{q}}^{\dagger}\hat{\Gamma}_{\mathbf{q}^{\prime}}$
(20)

with
$\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}=e^{i(\mathbf{q}^{\prime}-\mathbf{q})\cdot\hat{\mathbf{r}}}c_{\mathbf{q}^{\prime}-\mathbf{q}}/\hbar$
and $\hat{\Gamma}_{\mathbf{q}}^{\dagger}\equiv\hat{c}_{\mathbf{q}}^{\dagger}$
($\hat{\Gamma}_{\mathbf{q}}\equiv\hat{c}_{\mathbf{q}}$) being the creation
(annihilation) operator of a free fermion of momentum ${\mathbf{q}}$ with
$\\{\hat{c}_{\mathbf{q}},\hat{c}_{\mathbf{q}^{\prime}}^{\dagger}\\}=\delta_{\mathbf{q},\mathbf{q}^{\prime}}$.

If the impurity-bath interaction is described by the pseudopotential, as for
neutral impurities, we have 666Here we neglect the differential operator
$\partial_{r}(r\cdot)_{|_{r=0}}$, since its application to a plane wave has
unity effect.

$\displaystyle V_{\mathrm{ib}}(\mathbf{r})=g_{\mathrm{ib}}\delta(\mathbf{r}),$
(21)

where $g_{\mathrm{ib}}=2\pi\hbar^{2}a_{\mathrm{ib}}^{s}/\mu$ with
$a_{\mathrm{ib}}^{s}$ being the 3D $s$-wave impurity-bath scattering length.
Thus, Eq. (19) becomes

$\displaystyle c_{\mathbf{q}}=\frac{2\pi\hbar^{2}}{\mu
L^{3}}a_{\mathrm{ib}}^{s}.$ (22)

Specifically for the atom-ion case, we see that $f(q)\rightarrow-
a_{\mathrm{ib}}^{s}$. Thus, if we would replace the atom-ion interaction (4)
by Eq. (21), in all subsequent equations one has simply to replace the
scattering amplitude in the first-Born approximation with the scattering
length $a_{\mathrm{ib}}^{s}$.

###  III.2 Validity requirements and chemical potential

Let us comment on the validity of the Bogolyubov approximation, which implies
that both the quantum and thermal depletion must be small [56]. As we shall
discuss in Sec. IV, we are mostly interested in the high temperature regime,
that is, $k_{\mathrm{B}}\mathcal{T}\gg gn_{\mathrm{t}}$, which means that the
intra-particle interactions are essentially negligible. Here, $n_{\mathrm{t}}$
is the total gas density that in the zero temperature limit corresponds to the
condensate density $n_{0}$. In this regime, the Bogolyubov approximation is
valid only if
$n_{\mathrm{t}}\lambda_{\mathrm{dB}}(\mathcal{T})^{3}\gg\zeta(3/2)$. Here,
$\lambda_{\mathrm{dB}}(\mathcal{T})=[2\pi\hbar^{2}/(mk_{\mathrm{B}}\mathcal{T})]^{1/2}$
denotes the thermal de Broglie wavelength, $\zeta(x)$ is the zeta-Riemann
function. The last inequality can be also rewritten as
$\mathcal{T}\ll\mathcal{T}_{\mathrm{c}}^{0}=\frac{2\pi\hbar^{2}n_{\mathrm{t}}^{2/3}}{mk_{\mathrm{B}}[\zeta(3/2)]^{2/3}}$,
where $\mathcal{T}_{\mathrm{c}}^{0}$ is the critical temperature of
condensation of a non-interacting and untrapped Bose gas. Thus, in order to
perform the Bogolyubov approximation all conditions have to be fulfilled
simultaneously, namely the non-condensed fraction
$n_{\mathrm{t}}(\mathcal{T}/\mathcal{T}_{\mathrm{c}}^{0})^{3/2}\ll 1$. Here,
$n_{\mathrm{t}}=n_{0}+n_{\mathrm{n}}$ is the total density of the gas with
$n_{\mathrm{n}}$ being the normal (i.e. non-condensed) component, from which
we retrieve the condensate density as

$\displaystyle
n_{0}=n_{\mathrm{t}}-n_{\mathrm{n}}=n_{\mathrm{t}}\left[1-(\mathcal{T}/\mathcal{T}_{\mathrm{c}}^{0})^{3/2}\right].$
(23)

In table 1 we provide some values of the critical temperature of condensation
at typical quantum gas densities.

Boson $n_{\mathrm{t}}=10^{12}\,\mathrm{cm}^{-3}$ $10^{13}\,\mathrm{cm}^{-3}$
$10^{14}\,\mathrm{cm}^{-3}$ 7Li 229 1063 4934 23Na 70 324 1506 87Rb 18 86 398
Fermion $n_{\mathrm{t}}=10^{12}\,\mathrm{cm}^{-3}$ $10^{13}\,\mathrm{cm}^{-3}$
$10^{14}\,\mathrm{cm}^{-3}$ 6Li 613 2844 13198 40K 92 428 1987

Table 1:  Top: Critical temperature for condensation
$\mathcal{T}_{\mathrm{c}}^{0}$ of a non-interacting gas for three bosonic
species and densities. Bottom: Fermi temperature $\mathcal{T}_{F}$ for two
fermionic species and three densities. Temperature values are given in units
of nK.

Afterwards, it will be important to compute the chemical potential for a non-
interacting and homogenous Bose and Fermi gas at temperature $\mathcal{T}$. To
this end, let us remind that the chemical potential for the bosons reads [57]

$\displaystyle\mu_{\mathrm{G}}=\left\\{\begin{array}[]{l}\text{0, if
$\mathcal{T}\leq\mathcal{T}_{\mathrm{c}}^{0}$},\\\ \\\ \text{root of
$n_{\mathrm{t}}\lambda_{\mathrm{dB}}^{3}=g_{3/2}(z)$, if
$\mathcal{T}>\mathcal{T}_{\mathrm{c}}^{0}$}.\end{array}\right.$ (27)

Here, $z=e^{\frac{\mu_{\mathrm{G}}}{k_{\mathrm{B}}\mathcal{T}}}$ is the so-
called fugacity and $g_{3/2}(z)=\sum_{l=1}^{\infty}z^{l}l^{-3/2}$. We note
that the chemical potential for $\mathcal{T}>\mathcal{T}_{\mathrm{c}}^{0}$ is
negative or else $g_{3/2}(z)$ does not converge.

For the fermions, the chemical potential is obtained by solving numerically
the equation of state [57]

$\displaystyle n_{\mathrm{t}}\lambda_{\mathrm{dB}}^{3}=f_{3/2}(z),$ (28)

where $f_{3/2}(z)=\sum_{l=1}^{\infty}(-1)^{l+1}z^{l}l^{-3/2}$. At
$\mathcal{T}=0$, the chemical potential corresponds to the Fermi energy
$E_{F}$:

$\displaystyle\mu_{\mathrm{G}}\equiv
E_{F}=\frac{\hbar^{2}}{2m}\left(6\pi^{2}n_{\mathrm{t}}\right)^{2/3}.$ (29)

We note that for $\mathcal{T}>\mathcal{T}_{F}$, $\mu_{\mathrm{G}}$ is
negative, similarly to the bosonic case, where
$\mathcal{T}_{F}=E_{F}/k_{\mathrm{B}}$ is the Fermi temperature. In the high
temperature limit
$\mathcal{T}\gg\mathcal{T}_{\mathrm{c}}^{0},\,\mathcal{T}_{F}$, the chemical
potential of both the bosons and the fermions are well described by that of
the Boltzmann gas

$\displaystyle\mu_{\mathrm{G}}=k_{\mathrm{B}}\mathcal{T}\ln(n_{\mathrm{t}}\lambda_{\mathrm{dB}}^{3}).$
(30)

###  III.3 Markovian master equation

In this section we describe the relevant steps of the derivation of the master
equation for the bosons, while for the fermionic bath we simply provide the
final result, since the derivation is analogous.

We start from the full system-bath density matrix $\hat{\chi}(t)$, which obeys
the von Neumann equation

$\displaystyle\frac{d}{dt}\hat{\chi}(t)=-\frac{i}{\hbar}[\hat{H},\hat{\chi}].$
(31)

Writing the density operator in the interaction picture as

$\displaystyle\tilde{\chi}(t)=\hat{U}^{\dagger}(0,t)e^{i\hat{H}_{B}t/\hbar}\hat{\chi}(t)e^{-i\hat{H}_{B}t/\hbar}\hat{U}(0,t),$
(32)

where

$\displaystyle\hat{U}(t_{1},t_{2})=\mathsf{T}\exp\left[-\frac{i}{\hbar}\int_{t_{1}}^{t_{2}}\mathrm{d}t\,\hat{H}_{S}(t)\right]$
(33)

with $\mathsf{T}$ the time-ordered evolution operator, we have

$\displaystyle\frac{d}{dt}\tilde{\chi}(t)=-\frac{i}{\hbar}[\tilde{H}_{BS}(t),\tilde{\chi}(t)].$
(34)

Here, $\tilde{H}_{BS}$ is the interaction Hamiltonian in the interaction
picture, which is defined similarly to Eq. (32). The formal solution of Eq.
(34) is

$\displaystyle\tilde{\chi}(t)=\tilde{\chi}(0)-\frac{i}{\hbar}\int_{0}^{t}\mathrm{d}t^{\prime}\,[\tilde{H}_{BS}(t^{\prime}),\tilde{\chi}(t^{\prime})]$
(35)

and substituting it into the commutator (34) we obtain

$\displaystyle\frac{d}{dt}\tilde{\chi}(t)$
$\displaystyle=-\frac{i}{\hbar}[\tilde{H}_{BS}(t),\tilde{\chi}(0)]$
$\displaystyle-\frac{1}{\hbar^{2}}\int_{0}^{t}\mathrm{d}t^{\prime}\,[\tilde{H}_{BS}(t),[\tilde{H}_{BS}(t^{\prime}),\tilde{\chi}(t^{\prime})]].$
(36)

Thus, we assume that initially, at $t=0$, the system and the bath are
uncorrelated, namely
$\hat{\chi}(0)=\tilde{\chi}(0)=\hat{\rho}(0)\otimes\hat{B}_{0}$, where
$\hat{B}_{0}$ is the initial bath density matrix. This is a reasonable
assumption if the impurity and the bath are initially well-separated such that
no interaction occurs. By tracing over the bath degrees of freedom in Eq.
(III.3) we arrive to the equation 777Here,
$\mathrm{Tr}_{R}[\tilde{\chi}(t)]=\tilde{\rho}(t)$ and we eliminated the term
Tr${}_{B}\\{[\tilde{H}_{BS}(t),\tilde{\chi}(0)]\\}$ by assuming
Tr${}_{B}\\{\tilde{H}_{BS}(t)\hat{B}_{0}\\}=0$ (see Ref. [44] for details).

$\displaystyle\frac{d}{dt}\tilde{\rho}(t)$
$\displaystyle=-\frac{1}{\hbar^{2}}\int_{0}^{t}\mathrm{d}t^{\prime}\,\mathrm{Tr}_{B}\\{[\tilde{H}_{BS}(t),[\tilde{H}_{BS}(t^{\prime}),\tilde{\chi}(t^{\prime})]]\\}.$
(37)

The next step consists in performing the so-called Born approximation, namely
we assume that the impurity-bath coupling is weak and that the bath is so
large that
$\tilde{\chi}(t^{\prime})\simeq\tilde{\rho}(t^{\prime})\otimes\hat{B}_{0}$
$\forall t^{\prime}$ holds. Thus, Eq. (37) becomes

$\displaystyle\frac{d}{dt}\tilde{\rho}(t)$
$\displaystyle=-\frac{1}{\hbar^{2}}\int_{0}^{t}\mathrm{d}t^{\prime}\,\mathrm{Tr}_{B}\\{[\tilde{H}_{BS}(t),[\tilde{H}_{BS}(t^{\prime}),\tilde{\rho}(t^{\prime})\otimes\hat{B}_{0}]]\\}.$
(38)

In order to further simplify this equation, we make the Markov approximation,
namely we replace $\tilde{\rho}(t^{\prime})$ by $\tilde{\rho}(t)$ in order to
obtain a time-local master equation

$\displaystyle\frac{d}{dt}\tilde{\rho}(t)$
$\displaystyle=-\frac{1}{\hbar^{2}}\int_{0}^{t}\mathrm{d}t^{\prime}\,\mathrm{Tr}_{B}\\{[\tilde{H}_{BS}(t),[\tilde{H}_{BS}(t^{\prime}),\tilde{\rho}(t)\otimes\hat{B}_{0}]]\\}.$
(39)

This equation is known in the literature as the Redfield equation [59]. The
Hamiltonian $\tilde{H}_{BS}(t)$ keeps the original structure of the
Schrödinger picture, but with time-dependent system and bath operators:

$\displaystyle\hat{S}_{\mathbf{q}}(t)$
$\displaystyle=\hat{U}^{\dagger}(0,t)\hat{S}_{\mathbf{q}}\hat{U}(0,t)$
$\displaystyle\hat{\Gamma}_{\mathbf{q}}(t)$
$\displaystyle=e^{i\hat{H}_{B}t/\hbar}\hat{\Gamma}_{\mathbf{q}}e^{-i\hat{H}_{B}t/\hbar}=e^{-\frac{i}{\hbar}\varepsilon({\mathbf{q}})t}\hat{\Gamma}_{\mathbf{q}}$
$\displaystyle\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(u,u^{\prime})}(t)$
$\displaystyle=\hat{U}^{\dagger}(0,t)\hat{S}_{\mathbf{q},\mathbf{q}^{\prime}}^{(u,u^{\prime})}\hat{U}(0,t).$
(40)

Next, we need to perform the partial trace over the bath degrees of freedom,
namely we need to assess

$\displaystyle\mathrm{Tr}_{B}\\{[\tilde{H}_{BS}(t),[\tilde{H}_{BS}(t^{\prime}),\tilde{\rho}(t)\otimes\hat{B}_{0}]]\\}=$
$\displaystyle\mathrm{Tr}_{B}\\{[\tilde{H}_{BS}^{(1)}(t),[\tilde{H}_{BS}^{(1)}(t^{\prime}),\tilde{\rho}(t)\otimes\hat{B}_{0}]]\\}$
$\displaystyle+\mathrm{Tr}_{B}\\{[\tilde{H}_{BS}^{(1)}(t),[\tilde{H}_{BS}^{(2)}(t^{\prime}),\tilde{\rho}(t)\otimes\hat{B}_{0}]]\\}$
$\displaystyle+\mathrm{Tr}_{B}\\{[\tilde{H}_{BS}^{(2)}(t),[\tilde{H}_{BS}^{(1)}(t^{\prime}),\tilde{\rho}(t)\otimes\hat{B}_{0}]]\\}$
$\displaystyle+\mathrm{Tr}_{B}\\{[\tilde{H}_{BS}^{(2)}(t),[\tilde{H}_{BS}^{(2)}(t^{\prime}),\tilde{\rho}(t)\otimes\hat{B}_{0}]]\\}.$
(41)

Thus, we consider the bath’s thermal density matrix

$\displaystyle\\!\\!\\!\hat{B}_{0}=\frac{e^{-\beta(\hat{H}_{B}-\mu_{\mathrm{G}}\hat{N})}}{\mathcal{Z}},\qquad\mathcal{Z}=\mathrm{Tr}_{B}\\{e^{-\beta(\hat{H}_{B}-\mu_{\mathrm{G}}\hat{N})}\\},$
(42)

where $\hat{N}$ is the bath number operator. The mixed terms in the third and
fourth line of Eq. (III.3) are zero, since they contain an odd number of bath
operators. Thus, only the terms of the second (i.e. with
$\tilde{H}_{BS}^{(1)}(t)$ only) and the last line (i.e. with
$\tilde{H}_{BS}^{(2)}(t)$ only) of Eq. (III.3) remain.

First, we consider the thermal average of the double commutator involving
$\tilde{H}_{BS}^{(1)}(t)$, which includes the averages of two bath operators
only, for example,
$\langle\tilde{\Gamma}_{\mathbf{q}}(t)\tilde{\Gamma}_{\mathbf{q}^{\prime}}(t^{\prime})\rangle_{B_{0}}$=
Tr${}_{B}\\{\hat{B}_{0}\tilde{\Gamma}_{\mathbf{q}}(t)\tilde{\Gamma}_{\mathbf{q}^{\prime}}(t^{\prime})\\}$.
On the other hand,
$\langle\tilde{\Gamma}_{\mathbf{q}}(t)\tilde{\Gamma}_{\mathbf{q}^{\prime}}(t^{\prime})\rangle_{B_{0}}=\langle\tilde{\Gamma}^{\dagger}_{\mathbf{q}}(t)\tilde{\Gamma}^{\dagger}_{\mathbf{q}^{\prime}}(t^{\prime})\rangle_{B_{0}}=0$,
while

$\displaystyle\langle\tilde{\Gamma}_{\mathbf{q}}(t)\tilde{\Gamma}^{\dagger}_{\mathbf{q}^{\prime}}(t^{\prime})\rangle_{B_{0}}=e^{-\frac{i}{\hbar}\varepsilon({\mathbf{q}})(t-t^{\prime})}(n_{\mathbf{q}}+1)\delta_{\mathbf{q},\mathbf{q}^{\prime}}.$
(43)

Here,
$n_{\mathbf{q}}=\langle\hat{b}^{\dagger}_{\mathbf{q}}\hat{b}_{\mathbf{q}}\rangle=[e^{\beta_{\mathcal{T}}(\varepsilon({\mathbf{q}})-\mu_{\mathrm{G}})}-1]^{-1}$
is the Bose-Einstein occupation number and
$\beta_{\mathcal{T}}=1/(k_{\mathrm{B}}\mathcal{T})$.

The average of the double commutator with only $\tilde{H}_{RS}^{(2)}(t)$ has
terms that vanish when the number of raising and lowering bath operators is
not the same, while the non-zero contributions are given in appendix B.
Putting all together into Eq. (39) and performing the change of variable
$\tau=t-t^{\prime}$, and finally transforming back to the Schrödinger picture,
we arrive at the following final master equation

$\displaystyle\frac{d}{dt}\hat{\rho}(t)=-\frac{i}{\hbar}[\hat{H}_{S},\hat{\rho}]-\sum_{\mathbf{q}}\int_{0}^{t}\mathrm{d}\tau\,\Omega^{2}_{\mathbf{q}}\left\\{(n_{\mathbf{q}}+1)[\hat{Z}_{\mathbf{q}},\hat{W}_{\mathbf{q}}(t,\tau)\hat{\rho}(t)]+n_{\mathbf{q}}[\hat{\rho}(t)\hat{W}_{\mathbf{q}}(t,\tau),\hat{Z}_{\mathbf{q}}]+\text{H.
c.}\right\\}$
$\displaystyle-\sum_{\mathbf{q},\mathbf{q}^{\prime}}\int_{0}^{t}\mathrm{d}\tau\,\left\\{n_{\mathbf{q}}(n_{\mathbf{q}^{\prime}}+1)[\hat{Z}_{\mathbf{q}^{\prime}-\mathbf{q}},\hat{W}_{\mathbf{q}^{\prime}-\mathbf{q}}(t,\tau)\hat{\rho}(t)]\Omega_{u,u^{\prime}}^{u,u^{\prime}}(\mathbf{q}^{\prime}-\mathbf{q})+2n_{\mathbf{q}}(n_{\mathbf{q}^{\prime}}+1)[\hat{Z}_{\mathbf{q}^{\prime}-\mathbf{q}},\hat{W}_{\mathbf{q}^{\prime}-\mathbf{q}}(t,\tau)\hat{\rho}(t)]\Omega_{u,u^{\prime}}^{v,v^{\prime}}(\mathbf{q}^{\prime}-\mathbf{q})\right.$
$\displaystyle\left((1+n_{\mathbf{q}}+n_{\mathbf{q}^{\prime}}+n_{\mathbf{q}}n_{\mathbf{q}^{\prime}})[\hat{Z}_{\mathbf{q}^{\prime}+\mathbf{q}},\hat{W}_{\mathbf{q}^{\prime}+\mathbf{q}}(t,\tau)\hat{\rho}(t)]+n_{\mathbf{q}}n_{\mathbf{q}^{\prime}}[\hat{Z}_{\mathbf{q}^{\prime}+\mathbf{q}},\hat{\rho}(t)\hat{W}_{\mathbf{q}^{\prime}+\mathbf{q}}(t,\tau)]\right)\left(\Omega_{v,u^{\prime}}^{u,v^{\prime}}(\mathbf{q}^{\prime}+\mathbf{q})+\tilde{\Omega}_{v,u^{\prime}}^{u,v^{\prime}}(\mathbf{q}^{\prime}+\mathbf{q})\right)$
$\displaystyle\left.n_{\mathbf{q}^{\prime}}(n_{\mathbf{q}}+1)[\hat{Z}_{\mathbf{q}-\mathbf{q}^{\prime}},\hat{W}_{\mathbf{q}-\mathbf{q}^{\prime}}(t,\tau)\hat{\rho}(t)]\Omega_{v,v^{\prime}}^{v,v^{\prime}}(\mathbf{q}-\mathbf{q}^{\prime})+\text{H.
c.}\right\\}.$ (44)

Here, we have defined the operators

$\displaystyle\hat{Z}_{\mathbf{q}}=e^{i\mathbf{q}\cdot\hat{\mathbf{r}}}\qquad\hat{W}_{\mathbf{q}}(t,\tau)=e^{-\frac{i}{\hbar}\varepsilon(\mathbf{q})\tau}e^{-i\mathbf{q}\cdot\hat{\mathbf{r}}(t,\tau)},$
$\displaystyle\hat{W}_{\mathbf{q}^{\prime}-\mathbf{q}}(t,\tau)=e^{-\frac{i}{\hbar}(\varepsilon(\mathbf{q}^{\prime})-\varepsilon(\mathbf{q}))\tau}e^{-i(\mathbf{q}^{\prime}-\mathbf{q})\cdot\hat{\mathbf{r}}(t,\tau)},$
$\displaystyle\hat{W}_{\mathbf{q}^{\prime}+\mathbf{q}}(t,\tau)=e^{-\frac{i}{\hbar}(\varepsilon(\mathbf{q}^{\prime})+\varepsilon(\mathbf{q}))\tau}e^{-i(\mathbf{q}^{\prime}+\mathbf{q})\cdot\hat{\mathbf{r}}(t,\tau)},$
(45)

and the coefficients

$\displaystyle\Omega^{2}_{\mathbf{q}}$
$\displaystyle=\frac{|u_{\mathbf{q}}+v^{*}_{\mathbf{q}}|^{2}}{\hbar^{2}}|c_{\mathbf{q}}|^{2}n_{0}L^{3},$
$\displaystyle\Omega_{u,u^{\prime}}^{u,u^{\prime}}(\mathbf{q}^{\prime}-\mathbf{q})$
$\displaystyle=\frac{|u_{\mathbf{q}}|^{2}|u_{\mathbf{q}^{\prime}}|^{2}}{\hbar^{2}}|c_{\mathbf{q}^{\prime}-\mathbf{q}}|^{2},$
$\displaystyle\Omega_{v,v^{\prime}}^{v,v^{\prime}}(\mathbf{q}^{\prime}-\mathbf{q})$
$\displaystyle=\frac{|v_{\mathbf{q}}|^{2}|v_{\mathbf{q}^{\prime}}|^{2}}{\hbar^{2}}|c_{\mathbf{q}^{\prime}-\mathbf{q}}|^{2},$
$\displaystyle\Omega_{u,u^{\prime}}^{v,v^{\prime}}(\mathbf{q}^{\prime}-\mathbf{q})$
$\displaystyle=\frac{u_{\mathbf{q}}^{*}v_{\mathbf{q}}^{*}u_{\mathbf{q}^{\prime}}v_{\mathbf{q}^{\prime}}}{\hbar^{2}}|c_{\mathbf{q}^{\prime}-\mathbf{q}}|^{2},$
$\displaystyle\Omega_{v,u^{\prime}}^{u,v^{\prime}}(\mathbf{q}^{\prime}+\mathbf{q})$
$\displaystyle=\frac{u_{\mathbf{q}}^{*}v_{\mathbf{q}}^{*}u_{\mathbf{q}^{\prime}}v_{\mathbf{q}^{\prime}}}{\hbar^{2}}|c_{\mathbf{q}^{\prime}+\mathbf{q}}|^{2},$
$\displaystyle\tilde{\Omega}_{v,u^{\prime}}^{u,v^{\prime}}(\mathbf{q}^{\prime}+\mathbf{q})$
$\displaystyle=\frac{|v_{\mathbf{q}}|^{2}|u_{\mathbf{q}^{\prime}}|^{2}}{\hbar^{2}}|c_{\mathbf{q}^{\prime}+\mathbf{q}}|^{2},$
(46)

whereas

$\displaystyle\hat{\mathbf{r}}(t,\tau)=\hat{U}(0,t)\hat{U}^{\dagger}(0,t-\tau)\,\hat{\mathbf{r}}\,\hat{U}(0,t-\tau)\hat{U}^{\dagger}(0,t).$
(47)

This relation describes the impurity position evolution in absence of the gas.
Equation (III.3) is not yet a Markovian master equation, even though the time
development of the system density matrix relies only on the density matrix at
time $t$. Indeed, the impurity density matrix in Eq. (III.3) still depends on
the specific choice for the system preparation at $t=0$ via the impurity’s
trajectory in $W_{\mathbf{q}-\mathbf{q}^{\prime}}(t,\tau)$ or, in a more
precise mathematical language, it is not yet a dynamical semigroup [59].
Hence, to render (III.3) a Markovian master equation, we let the upper limit
of the integral go to infinity, which is permissible if the integrand
disappears sufficiently fast for
$\tau\gg\tau_{R}=\hbar/(k_{\mathrm{B}}\mathcal{T})$. This is justifiable if
the time scale of the system, $\tau_{S}$, over which the system density matrix
$\hat{\rho}(t)$ varies appreciably, is much larger than the time scale of the
bath, $\tau_{R}$. In other words, we require the bath correlation functions to
decay much faster than $\tau_{S}$ [59]. Hence, the Markov approximation is
justified when the bath correlation functions, e.g.,
$\langle\tilde{\Gamma}_{\mathbf{q}}(t)\tilde{\Gamma}_{\mathbf{k}}(t^{\prime})\rangle_{B_{0}}$=Tr${}_{B}\\{\hat{B}_{0}\tilde{\Gamma}_{\mathbf{q}}(t)\tilde{\Gamma}_{\mathbf{k}}(t^{\prime})\\}$,
are proportional to $\delta(t-t^{\prime})$. Instead, the Born approximation is
fulfilled if the dissipative damping rate is smaller than the relevant
system’s transition frequencies. We shall come back to this point later in the
paper. Note that with the upper limit of the integral going to infinity, the
first line of Eq. (III.3) is equal to Eq. (23) of Ref. [33]. However, while in
the latter the time integral in the definition of
$W_{\mathbf{k},\mathbf{k}^{\prime}}(t)$ is included, we prefer here to write
it explicitly.

Finally, we underline that up until now Eq. (III.3) is valid for any impurity
in a condensate (not only for an ion), provided that the Fourier transform
(19) can be computed. Indeed, the solution to Eq. (47) depends on the impurity
dynamics only and the Hamiltonian $\hat{H}_{S}(t)$ can also represent the free
evolution of a not trapped ion in a BEC or an impurity atom in an optical
lattice. For a normal Fermi gas, however, the master equation (III.3) reduces
to the double sum only, namely the sum over $\mathbf{q}$ in the first line
disappears. Moreover, only one term of the double sum contributes, as we have
single-particle energy states and not Bogolyubov modes. In practice, one sets
in Eq. (III.3) the Bogolyubov amplitudes $u\equiv 1$ and $v\equiv 0$ and
replaces $n_{\mathbf{q}^{\prime}}+1$ by $1-n_{\mathbf{q}^{\prime}}$, because
of the anticommutation relations of the fermionic field operators. Here,
$n_{\mathbf{q}}=[e^{\beta_{\mathcal{T}}(\varepsilon({\mathbf{q}})-\mu_{\mathrm{G}})}+1]^{-1}$
is the Fermi-Dirac occupation number with $\mu_{\mathrm{G}}$ being the
chemical potential obtained from solving Eq. (28).

##  IV Trapped ion master equation

The Markovian master equation (III.3) can be further simplified for an ion in
a radiofrequency trap, because of the separation of energy and length scales
between the atomic ensemble and the trapped ion system. Indeed, we are going
to make two further approximations:

  * (a)

The particle-like approximation;

  * (b)

The Lamb-Dicke approximation.

The former concerns the bosonic energy dispersion (12). Because of the large
energy difference between the bosonic bath and the ion system, only particle-
like excitations couple to the ion motion. This implies that the Bogolyubov
dispersion relation (12) is essentially quadratic in the wave vector $q$. For
a linear Paul trap, for instance, we have (see Appendix C.1 for details):
$a_{x,y}=-0.001$, $a_{z}=0.002$, $q_{x}=0.2$, $q_{y}=-0.2$, $q_{z}=0$, and
$\Omega_{rf}=2\pi\times$2 MHz, as obtained from the trap design of Ref. [60].
These parameters yield the following reference trap frequencies for an
ytterbium ion (see appendix C.2): $\nu_{x}\simeq 2\pi\times$169 kHz,
$\nu_{y}\simeq 2\pi\times$112 kHz, $\nu_{z}\simeq 2\pi\times$45 kHz (the $z$
axis is the longitudinal direction, where a system with two or more ions would
form a linear crystal). For the sake of convenience, we rescale the dispersion
relation (12) as

$\displaystyle\frac{\omega_{\mathbf{q}}}{\nu_{\xi}}=\sqrt{\left(\frac{\bar{q}^{2}}{2}\right)^{2}+(\bar{c}_{s}\bar{q})^{2}}\qquad\xi=x,y,z$
(48)

where $\bar{q}=\ell_{\xi}q$ with $\ell_{\xi}=\sqrt{\hbar/(m\nu_{\xi})}$, and
$\bar{c}_{s}=4\pi(a_{\mathrm{aa}}^{s}/\ell_{\xi})n\ell_{\xi}^{3}$. For sodium
atoms with a density $n_{0}=10^{14}$cm-3 we have: $\bar{c}_{s}\simeq 0.009$
($\xi\equiv x$), $\bar{c}_{s}\simeq 0.014$ ($\xi\equiv y$), and
$\bar{c}_{s}\simeq 0.034$ ($\xi\equiv z$). The bosons speed of sound
($c_{s}=\bar{c}_{s}\ell_{\xi}\nu_{\xi}$) is therefore quite small compared to
the ion motion in the secular trap such that only phonons of quite low
$\bar{q}\sim\bar{c}_{s}$ (i.e., large wavelength) yield an appreciable
difference in the dispersion relation compared to the free particle energy
$\bar{q}^{2}/2$. On the other hand, only phonons in the condensate of
comparable energy as $\hbar\nu_{\xi}$ will couple to the ion motion, so that
we can safely assume a particle-like dispersion relation,
$\varepsilon(\mathbf{q})=\hbar^{2}q^{2}/(2m)$, and set $u_{\mathbf{q}}\simeq
1$, $v_{\mathbf{q}}\simeq 0$, namely the bosonic bath can be treated as a non-
interacting Bose gas. This corresponds to an atom velocity of
$\sqrt{2\hbar\nu_{\xi}/m}$ in the $\xi$-th direction. For example,
$v_{x}\simeq 0.077\,$m/s or, in rescaled units, $\bar{v}_{x}\simeq 1.414$,
which is much larger than $\bar{c}_{s}$, and therefore the atom is moving at
supersonic velocities. To such an atomic velocity it is associated the
wavelength $\lambda_{x}\simeq 226.63\,$nm. As a consequence of the energy
separation, several terms of the quadratic corrections of the atom-ion
interaction in Eq. (III.3) can be safely discarded.

Figure 2: (Color online). Schematic representation of the Lamb-Dicke
approximation: the size $l$ of the ion trap (corresponding to the size of the
ion in its ground state) is much smaller than the de Broglie wavelength
$\lambda_{\mathrm{dB}}$ of the atoms in the gas.

On the other hand, the approximation (b) implies that the typical wavelength
of the bosons or of the fermions, i.e. the thermal de Broglie wavelength
$\lambda_{\mathrm{dB}}(\mathcal{T})$, is much larger than the $\xi$-th width
of the ion ground state, $l_{\xi}=\sqrt{\hbar/(M\nu_{\xi})}$ (see Fig. 2). Let
us still consider the example of bosonic sodium atoms at a temperature of
$\mathcal{T}=200$ nK. Thus, we have $\lambda_{\mathrm{dB}}(\mathcal{T})\simeq
814.18\,$nm, while for a trapped ytterbium ion the width of the ground state
in the secular trap is $l_{z}\simeq 36.05\,$nm. This yields a ratio
$l_{z}/\lambda_{\mathrm{dB}}(\mathcal{T})\simeq 0.044$ (similarly for the
other directions). Even if we consider the previously estimated supersonic
atom velocity, we get $l_{\xi}/\lambda_{\xi}\simeq 0.045$ for lithium atoms,
$l_{\xi}/\lambda_{\xi}\simeq 0.082$ for sodium atoms, and
$l_{\xi}/\lambda_{\xi}\simeq 0.159$ for rubidium atoms. Hence, the ion
spreading is quite localised compared to that of the bath’s particle, and
therefore the Lamb-Dicke approximation holds very well in the regime we are
interested in.

###  IV.1 Simplified master equation

Under the applicability conditions of the Lamb-Dicke approximation, we can
legitimately expand the exponential functions appearing in Eq. (III.3) up to
second-order in $\mathbf{q}\cdot\hat{\mathbf{r}}$. For instance, the
commutator reduces to

$\displaystyle[\hat{Z}_{\mathbf{q}},\hat{W}_{\mathbf{q}}(t,\tau)\hat{\rho}(t)]\simeq
e^{-\frac{i}{\hbar}\varepsilon(\mathbf{q})\tau}\left\\{i[\mathbf{q}\cdot\hat{\mathbf{r}},\hat{\rho}(t)]\right.$
$\displaystyle\left.+[\mathbf{q}\cdot\hat{\mathbf{r}},\mathbf{q}\cdot\hat{\mathbf{r}}(t,\tau)\hat{\rho}(t)]-\frac{1}{2}[(\mathbf{q}\cdot\hat{\mathbf{r}})^{2},\hat{\rho}(t)]\right\\}.$
(49)

Here, the operators
$\hat{W}_{\mathbf{q}}(t,\tau)\equiv\hat{W}_{\mathbf{q}}(t,t-\tau)$ and
$\hat{\mathbf{r}}(t,\tau)\equiv\hat{\mathbf{r}}(t-\tau)$ [see also Eq. (47)].
Given this, the three directions are decoupled from each other, because odd
powers of the wave vectors $\mathbf{q}$ vanish, as a consequence of the
symmetric summation in the master equation, that is, since the bath is
homogeneously confined. Therefore, Eq. (III.3) can be rewritten as follows

$\displaystyle\frac{d}{dt}\hat{\rho}(t)=-\frac{i}{\hbar}[\hat{H}_{S},\hat{\rho}]-\sum_{\mathbf{q},\xi}\int_{0}^{\infty}\mathrm{d}\tau\,\Omega^{2}_{\mathbf{q}}q_{\xi}^{2}\left\\{i\sin\left(\frac{\varepsilon(\mathbf{q})\tau}{\hbar}\right)[\hat{r}_{\xi}^{2},\hat{\rho}(t)]+e^{-i\frac{\varepsilon(\mathbf{q})\tau}{\hbar}}[\hat{r}_{\xi},\hat{r}_{\xi}(t,\tau)\hat{\rho}(t)]-e^{i\frac{\varepsilon(\mathbf{q})\tau}{\hbar}}[\hat{r}_{\xi},\hat{\rho}(t)\hat{r}_{\xi}(t,\tau)]\right.$
$\displaystyle\left.+2\,n_{\mathbf{q}}\cos\left(\frac{\varepsilon(\mathbf{q})\tau}{\hbar}\right)\left([\hat{r}_{\xi},\hat{r}_{\xi}(t,\tau)\hat{\rho}(t)]-[\hat{r}_{\xi},\hat{\rho}(t)\hat{r}_{\xi}(t,\tau)]\right)\right\\}-\sum_{\mathbf{q},\mathbf{q}^{\prime}}\sum_{\xi}\int_{0}^{\infty}\mathrm{d}\tau\,n_{\mathbf{q}}(n_{\mathbf{q}^{\prime}}+1)(q_{\xi}^{\prime}-q_{\xi})^{2}\Omega_{u,u^{\prime}}^{u,u^{\prime}}(\mathbf{q}^{\prime}-\mathbf{q})$
$\displaystyle\times\left\\{i\sin\left(\frac{\varepsilon(\mathbf{q}^{\prime})-\varepsilon(\mathbf{q})}{\hbar}\tau\right)[\hat{r}_{\xi}^{2},\hat{\rho}(t)]+e^{-i\frac{[\varepsilon(\mathbf{q}^{\prime})-\varepsilon(\mathbf{q})]\tau}{\hbar}}[\hat{r}_{\xi},\hat{r}_{\xi}(t,\tau)\hat{\rho}(t)]-e^{i\frac{[\varepsilon(\mathbf{q}^{\prime})-\varepsilon(\mathbf{q})]\tau}{\hbar}}[\hat{r}_{\xi},\hat{\rho}(t)\hat{r}_{\xi}(t,\tau)]\right\\}$
(50)

with $\xi=x,y,z$. The first sum over $\mathbf{q}$ in the first two lines of
Eq. (IV.1) refers to the Frölich model, while the double sum in the last two
lines refers to the extended Frölich model (III.1), that is,
$\hat{H}_{RS}^{(2)}$. Because of the particle-like approximation, however,
only the term $\Omega_{u,u^{\prime}}^{u,u^{\prime}}$ remains, while for the
fermions no Frölich interaction appears.

In order to make further progress, we need to explicitly use the solution of
the ion dynamics in the Paul trap in the absence of the gas. The full solution
$\hat{r}_{\xi}(t)$ is provided with various details in Ref. [33], which is
given by

$\displaystyle\hat{r}_{\xi}(t,\tau)$
$\displaystyle=\sum_{s,s^{\prime}}C_{s}^{\xi}C_{s^{\prime}}^{\xi}\left[\hat{r}_{\xi}\left(\frac{\beta_{\xi}}{2}+s^{\prime}\right)\frac{\Omega_{rf}}{\nu_{\xi}}\cos(\mathcal{I}^{\xi}_{s,s^{\prime}}(t,\tau))\right.$
$\displaystyle\left.-\frac{\hat{p}_{\xi}}{\nu_{\xi}M}\sin(\mathcal{I}^{\xi}_{s,s^{\prime}}(t,\tau))\right]$
(51)

for $\xi=x,y,z$, $[\hat{\mathbf{r}},\hat{\mathbf{p}}]=i\hbar$
($\hat{\mathbf{p}}$ is the ion momentum operator), and

$\displaystyle\mathcal{I}^{\xi}_{s,s^{\prime}}(t,\tau)=\Omega_{rf}\left[\left(\frac{\beta_{\xi}}{2}+s\right)\tau-(s-s^{\prime})t\right].$
(52)

The coefficients $C_{s}^{\xi}$, the trap parameters $\beta_{\xi}$, and the
frequencies $\nu_{\xi}$ are introduced in the appendix C.1.

Since we consider a gas confined in a box of volume $L^{3}$, the wave vector
$\mathbf{q}$ assumes quantised values: $2\pi s_{\xi}/L$ with
$s_{\xi}\in\mathbb{Z}$ (periodic boundary conditions) and $\xi=x,y,z$. Hence,
in the limit $L\rightarrow+\infty$ the allowed values of $\mathbf{q}$ in
momentum space become closely spaced, and since their density is
d$\mathbf{s}=L^{3}/(2\pi)^{3}$d$\mathbf{q}$, we make the replacement

$\displaystyle\sum_{\mathbf{q}}\rightarrow\frac{L^{3}}{(2\pi)^{3}}\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{q}.$
(53)

Such a continuum limit approximation is reasonable for a large bath.
Furthermore, we use the identity [44]

$\displaystyle\int_{0}^{\infty}\mathrm{d}\tau\,e^{-i(\omega-\omega_{0})\tau}=\pi\delta(\omega-\omega_{0})-i\mathcal{P}\left(\frac{1}{\omega-\omega_{0}}\right),$
(54)

where $\mathcal{P}$ denotes the Cauchy principal value (CPV), whose action on
a test function $\varphi(\omega)$ is:

$\displaystyle\mathcal{P}\left(\frac{1}{\omega}\right)(\varphi)=\lim_{\epsilon\rightarrow
0^{+}}\int_{\mathbb{R}/[-\epsilon,\epsilon]}\mathrm{d}\omega\,\frac{\varphi(\omega)}{\omega}.$
(55)

Now, we apply these results and we focus first our attention on the Frölich
contribution to the master equation, namely the incoherent term in Eq. (IV.1)
involving the summation over $\mathbf{q}$. Besides, at the moment, we neglect
the contribution due to the CPV and look at the $\delta$-contributions only.
Thus, $\omega\equiv\varepsilon(\mathbf{q})/\hbar$, while $\omega_{0}\equiv 0$
or $\omega_{0}\equiv\Omega_{rf}(\beta_{\xi}/2+s)$ in Eq. (54). When
$\omega_{0}\equiv 0$, however, the contribution of the term coming from the
sine function in the first line of Eq. (IV.1) vanishes, as it can be verified
by performing the integration (53) in spherical coordinates. On the other
hand, when $\omega_{0}\equiv\Omega_{rf}(\beta_{\xi}/2+s)$ and after having
moved to spherical coordinates, we first write the Dirac’s delta as

$\displaystyle\delta(\omega-\omega_{0})=\frac{m}{\hbar q}\delta(q-q_{s,\xi})$
(56)

with $q=|\mathbf{q}|$, and 888We note that the definition of $q_{s,\xi}$ is a
mere definition originated by the fact that
$q_{s,\xi}^{2}=2m\Omega_{rf}(\beta_{\xi}/2+s)/\hbar$ has been introduced when
manipulating the exponential and trigonometric functions in the master
equation (IV.1). Importantly, $q_{s,\xi}^{2}$ can assume negative or positive
values, upon the sign of $\beta_{\xi}/2+s$.

$\displaystyle
q_{s,\xi}=\sqrt{\frac{2m\Omega_{rf}}{\hbar}\left|\frac{\beta_{\xi}}{2}+s\right|}.$
(57)

Hence, the integration in momentum space yields

$\displaystyle\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{q}\,\frac{f^{2}(q)q^{2}_{\xi}}{q}\delta(q-q_{s,\xi})=\frac{4}{3}\pi
q_{s,\xi}^{3}f^{2}(q_{s,\xi})$ (58)

with $f$ being the scattering amplitude (II) evaluated in $q=q_{s,\xi}$.

As far as the contribution of the extended Fröhlich model is concerned, namely
the terms due to the double summation over $\mathbf{q}$ and
$\mathbf{q}^{\prime}$ in Eq. (IV.1), we proceed in a very similar manner with
the exception that we now have to assess a double integration in momentum
space. In the appendix D.1 we provide details of this calculation. On the
other hand, in order to assess the contribution due the Lamb-shift we have to
compute the integral (55), whose details are outlined in the appendix D.2. We
note, however, that we performed such a Lamb-shift calculation only for the
linear terms of Eq. (IV.1), since the contribution of the extended Fröhlich
model is much smaller, and therefore it can be neglected.

Under the above outlined approximations, we arrive at

$\displaystyle\frac{d}{dt}\hat{\rho}(t)=-\frac{i}{\hbar}[\hat{H}_{S}+\delta\hat{H}_{S},\hat{\rho}]-\Gamma\sum_{\xi=x,y,z}\left\\{\Lambda_{\xi}[\hat{r}_{\xi},\hat{\rho}(t)\hat{p}_{\xi}]\right.$
$\displaystyle-\Lambda_{\xi}^{*}[\hat{r}_{\xi},\hat{p}_{\xi}\hat{\rho}(t)]-\Phi_{\xi}[\hat{r}_{\xi},\hat{\rho}(t)\hat{r}_{\xi}]\left.+\Phi_{\xi}^{*}[\hat{r}_{\xi},\hat{r}_{\xi}\hat{\rho}(t)]\right\\}.$
(59)

Here, $\Gamma=\frac{2}{3}\frac{mn_{0}\pi\hbar}{\mu^{2}}$ and
$\delta\hat{H}_{S}$ is the correction to the free ion Hamiltonian due to the
interaction to the quantum gas (i.e., the Lamb-shift), which is given by

$\displaystyle\delta\hat{H}_{S}=(1-\phi)\frac{M}{2}\sum_{\xi=,x,y,z}\delta
W_{\xi}(t)\hat{r}_{\xi}^{2}$ (60)

with $\phi=0$ for the bosons and $\phi=1$ for the fermions,

$\displaystyle\delta W_{\xi}(t)=\frac{\Omega_{rf}^{2}}{4}\left[\delta
a_{\xi}-2\delta q_{\xi}\cos(\Omega_{rf}t)-2\delta
q_{\xi}^{\prime}\,g_{\xi}(t)\right],$ (61)

and

$\displaystyle g_{\xi}(t)=\\!\sum_{s,s^{\prime}\notin
S_{\mathrm{i}}}\\!F_{s,s^{\prime}}^{\xi}\cos[(s-s^{\prime})\Omega_{rf}t]\left[{\mathcal{J}}_{+}^{\prime}(q_{s,\xi})-{\mathcal{J}}_{-}^{\prime}(q_{s,\xi})\right],$
$\displaystyle\delta
a_{\xi}=-Q{\mathcal{J}}_{-}^{\prime}(0)+Q\sum_{s}\left\\{F^{\xi}_{s,s}\left[{\mathcal{J}}_{-}^{\prime}(q_{s,\xi})-{\mathcal{J}}_{+}^{\prime}(q_{s,\xi})\right]\right\\},$
$\displaystyle\delta
q_{\xi}=\frac{Q}{2}\sum_{|s-s^{\prime}|=1}\left\\{F^{\xi}_{s,s^{\prime}}\left[{\mathcal{J}}_{+}^{\prime}(q_{s,\xi})-{\mathcal{J}}_{-}^{\prime}(q_{s,\xi})\right]\right\\},$
(62)

where ${\mathcal{J}}_{\pm}^{\prime}(q_{s,\xi})$ are defined in Eq. (154) and
Eq. (155), $\delta q_{\xi}^{\prime}=Q/2$, and

$\displaystyle S_{\mathrm{i}}$
$\displaystyle=\\{(s,s^{\prime}):|s-s^{\prime}|=0\,\,\text{or}\,\,1\\},$
$\displaystyle F_{s,s^{\prime}}^{\xi}$
$\displaystyle=C_{s}^{\xi}C_{s^{\prime}}^{\xi}\left(\frac{\beta_{\xi}}{2}+s^{\prime}\right)\frac{\Omega_{rf}}{2\nu_{\xi}},\qquad
S_{s,s^{\prime}}^{\xi}=i\frac{C_{s}^{\xi}C_{s^{\prime}}^{\xi}}{M\nu_{\xi}2},$
$\displaystyle Q$
$\displaystyle=\frac{32}{3}\frac{m}{M}\frac{\hbar^{2}n_{0}}{\mu^{2}\Omega_{rf}^{2}}.$
(63)

We see that the coupling to the quantum gas renormalises the geometric Paul
trap parameters as $a_{\xi}\mapsto a_{\xi}+\delta a_{\xi}$ and $q_{\xi}\mapsto
q_{\xi}+\delta q_{\xi}$ [see also Eq. (118)] and it yields additional time-
dependent driving terms [i.e., $g_{\xi}(t)$]. Moreover, we introduced the
functions
$\Phi_{\xi}(t)=\Phi_{\xi}^{\delta}(t)+(1-\phi)\Phi_{\xi}^{\mathcal{P}}(t)$,
$\Lambda_{\xi}(t)=\Lambda_{\xi}^{\delta}(t)+(1-\phi)\Lambda_{\xi}^{\mathcal{P}}(t)$,
where

$\begin{split}\Phi_{\xi}^{\delta}(t)=&\sum_{s,s^{\prime}}F_{s,s^{\prime}}^{\xi}\Bigg{\\{}|q_{s,\xi}|^{3}f(q_{s,\xi})^{2}(1-\phi)\\\\[4.30554pt]
\times&\bigg{[}\cos\big{[}(s-s^{\prime})\Omega_{rf}t\big{]}(1+2n_{q_{s,\xi}})\\\\[4.30554pt]
&+i\sin\big{[}(s-s^{\prime})\Omega_{rf}t\big{]}\mathrm{sgn}\bigg{(}\frac{\beta_{\xi}}{2}+s\bigg{)}\bigg{]}\\\\[4.30554pt]
&+\eta^{-}_{s,\xi}e^{-i(s-s^{\prime})\Omega_{rf}t}+\eta^{+}_{s,\xi}e^{i(s-s^{\prime})\Omega_{rf}t}\Bigg{\\}},\end{split}$
(64)
$\begin{split}\Lambda_{\xi}^{\delta}(t)=&\sum_{s,s^{\prime}}S_{s,s^{\prime}}^{\xi}\Bigg{\\{}|q_{s,\xi}|^{3}f(q_{s,\xi})^{2}(1-\phi)\\\\[4.30554pt]
\times&\bigg{[}i\sin\big{[}(s-s^{\prime})\Omega_{rf}t\big{]}(1+2n_{q_{s,\xi}})\\\\[4.30554pt]
&+\cos\big{[}(s-s^{\prime})\Omega_{rf}t\big{]}\mathrm{sgn}\bigg{(}\frac{\beta_{\xi}}{2}+s\bigg{)}\bigg{]}\\\\[4.30554pt]
&-\eta^{+}_{s,\xi}e^{-i(s-s^{\prime})\Omega_{rf}t}+\eta^{-}_{s,\xi}e^{i(s-s^{\prime})\Omega_{rf}t}\Bigg{\\}}\end{split}$
(65)

with

$\eta^{\pm}_{s,\xi}=\frac{m}{16\pi^{2}\hbar^{2}\beta_{\mathcal{T}}n_{0}}\Big{(}\mathcal{F}_{s,\xi}^{(1),\pm}+(-1)^{\phi}\mathcal{F}_{s,\xi}^{(2)}\Big{)}.$
(66)

Here, $\mathcal{F}_{s,\xi}^{(1),\pm}$ and $\mathcal{F}_{s,\xi}^{(2)}$ are
defined in Eq. (143) and Eq. (150), respectively, and

$\displaystyle\Phi_{\xi}^{\mathcal{P}}(t)=\frac{2}{\pi}\sum_{s,s^{\prime}}F_{s,s^{\prime}}^{\xi}\sin[(s-s^{\prime})\Omega_{rf}t]\left[{\mathcal{J}}_{+}^{\prime}(q_{s,\xi})+{\mathcal{J}}_{-}^{\prime}(q_{s,\xi})\right.$
$\displaystyle+\left.2\left({\mathcal{J}}_{+}^{n_{q}}(q_{s,\xi})+{\mathcal{J}}_{-}^{n_{q}}(q_{s,\xi})\right)\right]\left[2\Theta(\mathrm{sgn}(\beta_{\xi}/2+s))-1\right],$
$\displaystyle\Lambda_{\xi}^{\mathcal{P}}(t)=i\frac{2}{\pi}\sum_{s,s^{\prime}}S_{s,s^{\prime}}^{\xi}\left\\{\left[1-2\Theta(\mathrm{sgn}(\beta_{\xi}/2+s))\right]\left[{\mathcal{J}}_{+}^{\prime}(q_{s,\xi})\right.\right.$
$\displaystyle\left.+{\mathcal{J}}_{-}^{\prime}(q_{s,\xi})+2\left({\mathcal{J}}_{+}^{n_{q}}(q_{s,\xi})+{\mathcal{J}}_{-}^{n_{q}}(q_{s,\xi})\right)\right]\cos[(s-s^{\prime})\Omega_{rf}t]$
$\displaystyle\left.-i\sin\left[(s-s^{\prime})\Omega_{rf}t\right]\left[{\mathcal{J}}_{-}^{\prime}(q_{s,\xi})-{\mathcal{J}}_{+}^{\prime}(q_{s,\xi})\right]\right\\}.$
(67)

The labels $\delta$ and $\mathcal{P}$ indicate the origin of the contribution,
namely with $\delta$ from the Dirac’s delta in Eq. (54), whereas with
$\mathcal{P}$ from the Cauchy principal value, i.e. the Lamb-shift. On the
other hand, the terms involving the functions $\mathcal{F}_{s,\xi}^{(1,2)}$
stem from the extended Fröhlich model. Besides, the integrals
${\mathcal{J}}_{\pm}^{n_{q}}(q_{s,\xi})$ are defined in Eq. (D.2) and they
stem from the first term in the second line of Eq. (IV.1), i.e. the one
involving the cosine function. Interestingly, this term does not contribute to
the renormalisation of the Paul trap parameters (IV.1), as the sum of each
single exponential function, when the cosine is written as a linear
combination of exponentials, vanishes for the part that concerns the unitary
dynamics. Furthermore, because of the Lamb-Dicke approximation we made, we
note that the three directions are decoupled. Thus, Eq. (IV.1) can be split
into three components, each representing a direction of the ion motion. In
other words, we have a master equation per each direction.

A few remarks, however, are in order. First, we note that even the
contribution stemming from the CPV yields incoherent terms, as indicated by
the functions $\Phi_{\xi}^{\mathcal{P}}(t)$ and
$\Lambda_{\xi}^{\mathcal{P}}(t)$. Second, the contribution stemming from the
Dirac’s delta yields coherent (i.e., unitary) dynamics as well, since the
functions $\Phi_{\xi}(t)$ and $\Lambda_{\xi}(t)$ are complex. This is because
we did not perform the rotating-wave approximation (RWA), that is, we did not
neglect the non-secular terms with $\omega\neq\Omega_{rf}(\beta_{\xi}/2+s)$,
which is typically applied in quantum optics. Given the fact that the most
relevant Paul trap coefficients $C_{s}^{\xi}$ for a linear geometry are those
with $s=0,\pm 1$, the contributions in $\Phi_{\xi}^{\mathcal{P}}(t)$ and
$\Lambda_{\xi}^{\mathcal{P}}(t)$ with
$\sin\left[(s-s^{\prime})\Omega_{rf}t\right]$ are in general small, and
therefore one could in principle apply the RWA also in this context.
Nonetheless, we find that the term in $\Lambda_{\xi}^{\mathcal{P}}(t)$
involving $\cos\left[(s-s^{\prime})\Omega_{rf}t\right]$ has a non-negligible
contribution. This means that an application of the RWA would underestimate
the overall dissipative dynamics. A similar reasoning applies for the coherent
dynamics stemming from the imaginary contributions of the
$\Phi_{\xi}^{\delta}(t)$ and $\Lambda_{\xi}^{\delta}(t)$ functions. This is
also the reason why we cannot transform Eq. (IV.1) in Lindbland form (see also
appendix F), as an essential requisite is the RWA [59]. We note that up until
now, these effects have been not taken into account in investigations in the
context of impurity physics for settings similar to ours [47, 33, 62, 48].

##  V Ion energy and first order moments

The energy of the ion at time $t$ is given by the expectation value of the ion
Hamiltonian (119):

$\displaystyle\langle\hat{H}_{I}^{trap}(t)\rangle=\langle\hat{H}_{I}^{kin}(t)\rangle+\frac{M}{2}{W}_{\xi}^{\prime}(t)\langle\hat{r}_{\xi}^{2}\rangle,$
(68)

where ${W}_{\xi}^{\prime}(t)=W_{\xi}(t)+\delta W_{\xi}(t)$, with $W_{\xi}(t)$
given by

$\displaystyle
W_{\xi}(t)=\frac{\Omega_{rf}^{2}}{4}\left[a_{\xi}-2q_{\xi}\cos(\Omega_{rf}t)\right],$
(69)

$\delta W_{\xi}(t)$ defined in Eq. (61) and

$\langle\hat{H}_{I}^{\mathrm{kin}}(t)\rangle=\sum_{\xi=x,y,z}\frac{\langle\hat{p}_{\xi}^{2}\rangle}{2M}$
(70)

In order to calculate it, we have to determine the expectation values
$\langle\hat{p}_{\xi}^{2}(t)\rangle=$
Tr$\\{\hat{p}_{\xi}^{2}\hat{\rho}_{\xi}(t)\\}$ and
$\langle\hat{r}_{\xi}^{2}(t)\rangle=$
Tr$\\{\hat{r}_{\xi}^{2}\hat{\rho}_{\xi}(t)\\}$ $\forall\,\xi=x,y,z$ with
$\hat{\rho}_{\xi}(t)$ being the ion density matrix corresponding to the
$\xi$-th direction, whose equation of motion is obtained by considering only
the pertinent direction in the sum appearing in Eq. (IV.1). Instead of solving
the full master equation, however, it is computationally less expensive to
solve the corresponding differential equations for the expectation values of
the square of the position and momentum, which are coupled to the covariance
$\langle\hat{c}_{\xi}\rangle=\langle(\hat{r}_{\xi}\hat{p}_{\xi}+\hat{p}_{\xi}\hat{r}_{\xi})(t)\rangle$.
Therefore, by using the definition of the expectation value of an observable,
$\langle\hat{\mathcal{O}}\rangle=$Tr$\\{\hat{\mathcal{O}}\hat{\rho}(t)\\}$,
and the master equation (IV.1), one arrives at the set of coupled differential
equations:

$\displaystyle\frac{d}{dt}\langle\hat{r}_{\xi}^{2}\rangle=\frac{\langle\hat{c}_{\xi}\rangle}{M},$
$\displaystyle\frac{d}{dt}\langle\hat{p}_{\xi}^{2}\rangle\\!=\\!\left\\{2\hbar\Gamma\Im[\Phi_{\xi}(t)]-M{W}_{\xi}^{\prime}(t)\right\\}\langle\hat{c}_{\xi}\rangle-4\hbar\Gamma\Im[\Lambda_{\xi}(t)]\langle\hat{p}_{\xi}^{2}\rangle$
$\displaystyle+2\hbar^{2}\Gamma\Re[\Phi_{\xi}(t)],$
$\displaystyle\frac{d}{dt}\langle\hat{c}_{\xi}\rangle=2\left\\{2\hbar\Gamma\Im[\Phi_{\xi}(t)]-M{W}_{\xi}^{\prime}(t)\right\\}\langle\hat{r}_{\xi}^{2}\rangle+2\frac{\langle\hat{p}_{\xi}^{2}\rangle}{M}$
$\displaystyle-2\hbar\Gamma\Im[\Lambda_{\xi}(t)]\langle\hat{c}_{\xi}\rangle+2\hbar^{2}\Gamma\Re[\Lambda_{\xi}(t)].$
(71)

In the limit for which the Lamb-shift and the extended Fröhlich model are not
considered, the equations of motion (58) of Ref. [33] are retrieved. We note
that the set of equations (V) holds for both the bosonic and fermionic bath,
but with different $\Phi_{\xi}$ and $\Lambda_{\xi}$ functions.

The radiofrequency fields set the smallest time scale in the open quantum
system. It is therefore useful to evaluate the time averaged energy over a rf-
period, namely

$\displaystyle\langle\langle\hat{H}_{I}^{\mathrm{kin}}(t)\rangle\rangle_{T_{rf}}=\frac{1}{T_{rf}}\int_{t}^{t+T_{rf}}\mathrm{d}t^{\prime}\langle\hat{H}_{I}^{\mathrm{kin}}(t^{\prime})\rangle$
(72)

where $T_{rf}=2\pi/\Omega_{rf}$. Here, the notation
$\langle\langle\dots\rangle\rangle_{T_{rf}}$ denotes the time average over
$T_{rf}$. In this way we average out the fast oscillations due to the rf-field
and the final ion energy at thermal equilibrium with the atomic gas can be
assessed.

Finally, we provide equations of motion of the first order moments of the ion
position and momentum operators. Exactly with the same procedure that we
outlined previously for Eq. (V), the coupled differential equations for the
first order moments read:

$\displaystyle\frac{d}{dt}\langle\hat{r}_{\xi}\rangle=\frac{\langle\hat{p}_{\xi}\rangle}{M},$
(73)
$\displaystyle\frac{d}{dt}\langle\hat{p}_{\xi}\rangle\\!=\\!\left\\{2\hbar\Gamma\Im[\Phi_{\xi}(t)]-M{W}_{\xi}^{\prime}(t)\right\\}\\!\langle\hat{r}_{\xi}\rangle\\!-\\!2\hbar\Gamma\Im[\Lambda_{\xi}(t)]\langle\hat{p}_{\xi}\rangle.$

The first equation of motion of $\langle\hat{r}_{\xi}\rangle$ is simply the
definition of the ion momentum in the $\xi$-direction, while the second one
provides the average force acting on the ion. The latter is, on the one side,
due to the ion trap, i.e. the term proportional to ${W}_{\xi}^{\prime}(t)$,
and, on the other side, to the atom-ion interaction, namely the term
proportional to $\Im[\Phi_{\xi}(t)]$. Moreover, the term proportional to
$\Im[\Lambda_{\xi}(t)]$ relies on the ion momentum, that is, it corresponds to
a velocity-dependent force, which results in a damped ion motion because of
the presence of the gas, unless the initial conditions for
$\langle\hat{r}_{\xi}\rangle$ and $\langle\hat{p}_{\xi}\rangle$ vanish. It is
interesting to note that the form of the equations of motion (73) resembles
that of $\langle\hat{r}_{\xi}^{2}\rangle$ and
$\langle\hat{p}_{\xi}^{2}\rangle$, where $\langle\hat{c}_{\xi}\rangle$ is
replaced by $\langle\hat{p}_{\xi}\rangle$ or $\langle\hat{r}_{\xi}\rangle$ and
$\langle\hat{p}_{\xi}^{2}\rangle$ by $\langle\hat{p}_{\xi}\rangle$, but Eq.
(73) does not have source terms.

##  VI Numerical Results

Before we present and discuss our numerical findings, we note that hereafter
they are based on the parameters $b$ and $c$ computed as explained in appendix
A. Moreover, the bare Paul trap parameters are chosen as: $a_{x,y}=-0.001$,
$a_{z}=0.002$, $q_{x}=-q_{y}=0.2$, $q_{z}=0$, and $\Omega_{rf}=2\pi\times$2
MHz. Such parameters correspond to a linear trap geometry, whose frequencies
are: $\nu_{x}\simeq 2\pi\times$169 kHz, $\nu_{y}\simeq 2\pi\times$112 kHz,
$\nu_{z}\simeq 2\pi\times$45 kHz (a 174Yb ion is always assumed). Thus, we
have $\beta_{x,y}\simeq 0.1389$ and $\beta_{z}\simeq 0.0447$. Details on the
choice of the initial density matrix are provided in the appendix C.2.

###  VI.1 Renormalised Paul trap parameters

To begin with, we analyse the impact of the coupling to the quantum gas on the
ion dynamics by considering the renormalised trap $a$\- and $q$-parameter. The
absolute amount of change of those parameters from their bare values, i.e.
without the gas, provides us a rule of thumb to assess how strong can be the
coupling, that is, how large can be the atomic density such that the master
equation can yield a faithful description of the ion dynamics. The
modification of those parameters relies on the particular atom-ion species via
the mass ratio as well as on the condensate density.

Figure 3:  (Color Online). Corrections to the Paul trap parameters along the
$x$ direction for sodium and lithium atoms for two densities as a function of
the atom-ion scattering length. Top panel: relative correction to the
$q$-parameter. Bottom panel: relative correction to the $a$-parameter.

In Fig. 3 we show the corrections to the Paul trap parameters $\delta a_{\xi}$
and $\delta q_{\xi}$ relatively to their bare values as a function of the
atom-ion scattering length for $n_{0}=10^{12}/$cm3 and $n_{0}=10^{13}/$cm3 for
sodium and lithium, respectively, along the transverse direction $\xi\equiv
x$. Note that for lower densities the values of $\delta a_{\xi}$ and $\delta
q_{\xi}$ are reduced by $n_{0}/(10^{12}/$cm3) for sodium and by
$n_{0}/(10^{13}/$cm3) for lithium due to their definition. As it can be seen,
the $q$-parameter, namely that of the driving rf-field, is very weakly
affected by the coupling to the gas (top panel). The $a$-parameter, instead,
assumes larger values, especially for positive scattering lengths (bottom
panel). Furthermore, we see that the heavier the atom, the larger is the
impact on the trap (for equal densities), as expected (in the figure the
result for lithium has to be divided by ten to be compared with that of
sodium). These results show that, while for lithium densities up to
$n_{0}=10^{13}/$cm3 can be considered (ideally the ratio should be smaller
than unity), it is better not to exceed $n_{0}=10^{12}/$cm3 for sodium atoms
because of the strong modification to the $a$-parameter.

###  VI.2 Ion in a lithium gas

Figure 4:  (Color online). Ion temperature obtained from the averaged energy
(72) for $b\simeq 0.0780\,R^{\star}$, $c\simeq 0.2239\,R^{\star}$, which
correspond to $a_{\mathrm{ia}}^{s}\simeq R^{\star}$. The square and circle
symbols refer to the bosons, while the diamonds to the fermions. The black
vertical line indicates the position of the critical temperature of
condensation $\mathcal{T}_{\mathrm{c}}^{0}$, while the Fermi temperature is
not indicated. Panel (a): $n_{\mathrm{t}}=10^{12}\,\text{cm}^{-3}$,
$T=50\,\mathrm{ms}$, $\mathcal{T}_{\mathrm{F}}=0.61\,\mu\mathrm{K}$. Panel
(b): $n_{\mathrm{t}}=10^{13}\,\text{cm}^{-3}$, $T=6\,\mathrm{ms}$,
$\mathcal{T}_{\mathrm{F}}=2.48\,\mu\mathrm{K}$.

In Fig. 4 we show the averaged ion kinetic energy expressed as a temperature,
$\mathcal{T}_{\mathrm{kin}}=2\langle\langle\hat{H}_{I}^{\mathrm{kin}}(T)\rangle\rangle_{T_{rf}}/3k_{\mathrm{B}}$,
at the final time $T$, namely when the ion has thermalised with the quantum
gas, where the averaged energy is given by Eq. (72). The definition of the
kinetic temperature $\mathcal{T}_{\mathrm{kin}}$ includes the secular motion
and the micromotion of the ion and the factor $2/3$ is due to the
equipartition theorem [42]. We note that the value of the thermalisation time
$T$ relies on various system parameters, especially on the atomic density. For
a fixed scattering length, decreasing the density by an order of magnitude
implies an increase of the thermalisation time by roughly the same amount. As
a consequence and for $\mathcal{T}<\mathcal{T}_{\mathrm{c}}^{0}$, $T$ strongly
depends on the gas temperature, since the latter determines the density of the
condensed fraction, i.e. $n_{0}$. A first estimation of the thermalisation
time for each plot is found by looking at $\mathcal{T}_{\mathrm{kin}}$ as a
function of time for a single value of the temperature (see e.g. the bottom
panel of Fig. 7). Then, the values of $\mathcal{T}_{\mathrm{kin}}$ at all
temperatures are computed at the estimated time and at several larger times up
to the one at which the difference between the curves becomes negligible.
General remarks. Let us first focus on the bosonic case below
$\mathcal{T}_{c}^{0}$. As we can see from the insets of Fig. 4, the behavior
of the ion kinetic energy is the result of the interplay among the different
contributions. The empty blue circles correspond to the Fröhlich model, i.e.
they show the contribution of the condensed part of the gas only. The final
temperature of the ion is basically independent of the gas temperature in this
approximation. This result is consistent with the fact that the density barely
affects the final energy of the ion. Indeed, in the Fröhlich model the
dependence of the equations on the gas temperature is almost entirely carried
by the condensate density $n_{0}$, as the temperature-dependent factor
$n_{q_{s,\xi}}$ in the definitions of $\Phi(t)$ and $\Lambda(t)$ is always
much lower than unity. This behavior has to be traced back to the nature of
the condensate, in which all the particles occupy the same single particle
state. For this reason, the cooling effect of the condensate does not change
when its temperature changes, as the latter only affects the fraction of
condensed particles. A similar reasoning can be applied to the Lamb-shift,
whose contribution can be observed in the purple empty squares of Fig. 4. The
additional cooling effect is stronger at temperatures much lower than
$\mathcal{T}_{\mathrm{c}}^{0}$. This phenomenon is related to the condensate
density $n_{0}$ increasing when the temperature is decreased, which implies a
stronger coupling to the gas, and it is in agreement with what discussed in
section VI.1 about the dependence of the $\delta a_{\xi}$ and $\delta q_{\xi}$
parameters on $n_{0}$. On the other hand, the contribution of the extended
Fröhlich model shows the opposite trend. As it can be seen from the full blue
circles of the aforementioned insets, when $\mathcal{T}$ approaches the
critical temperature the contribution of the interaction with the normal part
of the gas bends the ion temperature downward. The extended Fröhlich model,
which is the only one contributing at
$\mathcal{T}>\mathcal{T}_{\mathrm{c}}^{0}$ in the bosonic case and at every
gas temperature $\mathcal{T}$ for the fermions, is responsible for a minimum
in the final ion temperature, i.e. kinetic energy. In order to understand its
nature, we studied the temperature dependence of some characteristic
quantities involved in the equations such as $\eta_{s,\xi}^{\pm}$ [see Eq.
(66)]. For simplicity, let us now focus on the fermionic case, where no other
contributions have to be considered. In this case, all the temperature
dependence relies on $\eta_{s,\xi}^{\pm}$, which monotonically increases with
$\mathcal{T}$ for all $s$ and $\xi=x,y,z$ (not shown). Moreover, by looking at
the kinetic energy along the three directions at $t=T$ (also not shown), i.e.
when thermalisation is achieved, we found that the temperature dependence
along $\xi=x$ and $\xi=y$ presents the same minimum of
$\mathcal{T}_{\mathrm{kin}}$ while along $\xi=z$ it is monotonic. We can
therefore attribute the emergence of the minimum to the presence of the trap
and, in particular, to the interference between terms with different $s$ due
to the radiofrequency-induced micromotion. Finally, in Fig. 5 we observe that
when the regularized polarization potential is substituted with the
pseudopotential, the depth of the minimum is strongly decreased for
$n_{\mathrm{t}}=10^{12}\mathrm{cm}^{-3}$ (orange data), and even barely
visible at $n_{\mathrm{t}}=10^{13}\mathrm{cm}^{-3}$ (light blue data). We thus
conclude that the long-range character of the atom-ion potential renders the
occurrence of the minimum in the kinetic energy more pronounced.
Density dependence. While the overall behavior is not substantially affected
by the value of the density, there are some differences that it is worth to
remark. First, for large densities, the ion temperature in the case of 6Li
(green diamonds in Fig. 4) is slightly lower at very low $\mathcal{T}$. This
difference, though, is not substantial and is definitely negligible compared
to the scale of $s$-wave energy threshold. Another difference concerns the
contribution of the Lamb-shift (purple squares), which is enhanced at large
densities, thus confirming what we discussed in section VI.1 and in the
previous paragraph. Both of these differences, though, are only visible at
very low temperatures. At high temperatures, nor the density and neither the
statistics of the gas influence the result in a sensible way, a part from the
time required for thermalization that, as anticipated, increases linearly with
the decrease of the density.

Figure 5:  (Color online). Ion temperature obtained from the averaged energy
(72) with the pseudopotential and $a_{\mathrm{ia}}^{s}=R^{\star}$. The dot-
dashed light blue lines correspond to
$n_{\mathrm{t}}=10^{13}\,\text{cm}^{-3}$, while the solid orange ones to
$n_{\mathrm{t}}=10^{12}\,\text{cm}^{-3}$. The simulation time is
$T=10\,\mathrm{ms}$ and $T=100\,\mathrm{ms}$, respectively.

Damping rates. We investigated the temperature dependence of the damping
rates, $\gamma_{x}$, of $\langle\hat{r}_{x}\rangle$. In Fig. 6 they are shown
for 6Li and 7Li with a density $n_{\mathrm{t}}=10^{13}\mathrm{cm}^{-3}$.
Before commenting on the result, let us briefly explain the procedure we
followed in order to calculate the values of $\gamma_{x}$. We started from a
non zero initial condition for $\langle\hat{r}_{x}\rangle$. Its value is not
particularly important, because the damping rates do not rely on it anyway. By
solving the system in Eq. (73), we obtain $\langle\hat{r}_{x}(t)\rangle$,
whose behavior is a damped oscillation with zero average value. We calculated
the curve enveloping the oscillation and we fitted it with the exponential
function $\alpha_{x}\,e^{-\gamma_{x}t}$ (the same procedure was applied to all
other directions with similar findings). Interestingly, the temperature
dependence of the coefficients $\gamma_{\xi}$ does not reproduce that of
$\mathcal{T}_{\mathrm{kin}}$. In the case of 6Li (fermion), the curve is
monotonic and so it is for 7Li (boson) above $\mathcal{T}_{\mathrm{c}}^{0}$.
In the bosonic case, when the gas temperature is reduced below
$\mathcal{T}_{\mathrm{c}}^{0}$, the damping rates increase with the density of
the condensate. The grey circles in the inset of Fig. 6 show that in the
Fröhlich model there is a one-to-one correspondence between the damping rates
and the condensate density. The extended Fröhlich model (brown squares)
enhances the damping rates when the contribution of the normal part of the gas
becomes stronger. This relation between the condensate density and the values
of $\gamma_{\xi}$ strongly underlines the difference between bosonic and
fermionic baths at low temperatures. Moreover, it could be exploited in
experiments, where the condensate fraction may be extracted from the
measurement of the damping rates.

Figure 6:  (Color online). Damping rates of $\langle\hat{r}_{\xi}\rangle$ as a
function of the gas temperature with a total density
$n_{\mathrm{t}}=10^{13}\mathrm{cm}^{-3}$. The dotted line in the inset
represents a fit with the condensate density dependence on the ratio
$\mathcal{T}/\mathcal{T}_{\mathrm{c}}^{0}$.

Scattering length dependence. Figure 7 shows the dependence of the ion energy
on the atom-ion scattering length at $\mathcal{T}=100\,\mathrm{nK}$ in the
case of 6Li. As it can be seen in the top panel, for some values of the
scattering length convergence is already achieved at $T=6\,\mathrm{ms}$. When
the value of $a_{\mathrm{ia}}^{s}$ approaches roughly $-0.5R^{*}$, the
thermalization time strongly increases, as shown in the bottom panel (red
dotted line). Such (numerical) observation suggests an instability that could
be related to the occurrence of a resonance, as experimentally observed
recently [15]. Indeed, although the master equation does not entail
information about two-particle bound states, the behavior of the final energy
could still give some hint about the microscopic dynamics, the latter emerging
through the parameters of the atom-ion potential. As a reference, the grey
dashed line shows the values of the converged ion temperature in the case of
the pseudopotential approximation, where the only parameter involved is the
scattering length. The curve is symmetric because $f(k)\propto-
a_{\mathrm{ia}}$ and the latter enters only via $|f(k)|^{2}$. Hence, the long-
ranged character of the polarisation potential together with the fact that
there is no separation of length scales in the impurity system crucially
affects the ion dynamical properties.
Let us conclude by mentioning that recently the polaronic properties of a free
ion in a condensate have been investigated [2]. There, quite different
polaronic states have been identified on the basis that the atom-ion
polarisation potential supports either one bound state or none. Without a two-
body bound state, a polaron resembling that of a neutral impurity, albeit with
a larger effective mass, has been found. Here, we have also investigated the
scenario for which the potential does not support any bound state. In this
case, however, we found for a broad range of negative atom-ion scattering
lengths that the ion does not thermalise, that is, its dynamics is very
unstable in the Paul trap. Thereby, even though in the framework of the master
equation we cannot make quantitative predictions on the dynamical formation of
many-body bound states such as its size, the presence of two-body bound states
and the inclusion of the extended Fröhlich model in the master equation
description is of paramount importance for stabilising the ion dynamics in the
rf-trap.

Figure 7:  (Color online). Top panel: Ion temperature vs. the atom-ion
scattering length for 6Li for $n_{\mathrm{t}}=10^{13}\mathrm{cm}^{-3}$ and a
gas temperature $\mathcal{T}=0.1\mu\mathrm{K}$. Bottom panel: Ion temperature
vs. time for two different scattering lengths.

###  VI.3 Ion in a sodium gas

We have also investigated the ion dynamics in a heavier bosonic ensemble. In
Fig. 8 we illustrate the result of this analysis. With the linear Paul trap
that we have assumed so far, the ion energy as a function of the gas
temperature is shown by the purple squares. As it can be seen, the ion energy
is always above the s-wave threshold and therefore no ultracold atom-ion
collisions can be expected. Nonetheless, by reducing the rf-frequency by an
order of magnitude, that is, $\Omega_{rf}=2\pi\times$ 200 kHz, and by reducing
the $q$-parameter by one half, i.e. $q_{x}=-q_{y}=0.1$, which result in the
trap frequencies $\nu_{x}\simeq 2\pi\times$ 7 kHz, $\nu_{y}\simeq 2\pi\times$
6 kHz, $\nu_{z}\simeq 2\pi\times$ 4 kHz and $\beta_{x,y}\simeq 0.0634$,
$\beta_{z}\simeq 0.0447$, we find that at low gas temperatures (i.e. $\lesssim
400\,$nK) the $s$-wave limit can be beaten, as shown by the grey squares with
dotted line in Fig. 8. The resulting trap frequencies indicate an almost
isotropic ion trap. It has to be noticed, however, that with such a shallow
Paul trap it will be experimentally challenging to suppress excess micromotion
and to keep under control the ion heating due to the background electric
noise.

Figure 8:  (Color online). Ion temperature obtained from the averaged energy
(72) for $b\simeq 0.0780\,R^{\star}$, $c\simeq 0.2239\,R^{\star}$, which
correspond to $a_{\mathrm{ia}}^{s}\simeq R^{\star}$, and a total density
$n_{\mathrm{t}}=10^{12}\,\text{cm}^{-3}$. The grey squares with dotted line
correspond to the shallow trap (see text for parameters). The black vertical
lines indicate the position of the critical temperature of condensation
$\mathcal{T}_{\mathrm{c}}^{0}$, while the dashed horizontal one to the
$s$-wave threshold.

##  VII Summary of the main results and conclusions

We have investigated the quantum non-equilibrium dynamics of an ion in an rf-
trap superimposed to a quantum gas of either bosons or fermions. To this end,
in Sec. IV, we developed a quantum master equation by including the
contribution of the Lamb-shift and the extended Fröhlich model, that is, the
non-condensate fraction. The final master equation can be found in Eq. (IV.1),
where the definitions of the corresponding functions discern the case of the
fermionic and bosonic bath. The equations of motion for the second and first
moments are calculated in Sec. V directly from the master equation and are
given in Eq. (V) and Eq. (73), respectively. These systems of equations were
numerically solved resulting in the findings exposed in Sec. VI. We also note
that the master equation (III.3) is also an important result of our study, as
it can be the starting point for other investigations, such as a free ion in a
Bose-Einstein condensate. We found significant differences in the ion dynamics
between the bosonic and fermionic bath at low temperatures
($\mathcal{T}\lesssim\mathcal{T}_{\mathrm{c}}^{0}$), where the quantum nature
of the gas emerges. As shown in Fig. 4, in this regime, a fermionic
environment seems to ensure better cooling of the ion compared to the bosonic
one, while at higher temperatures the difference gets less and less pronounced
and it vanishes when $\mathcal{T}\gg\mathcal{T}_{\mathrm{c}}^{0}$. Similarly,
the nature of the gas affects the damping of the ion below
$\mathcal{T}_{\mathrm{c}}^{0}$ (see Fig. 6): the temperature dependence of the
damping rates in the bosonic case is strongly characterised by the presence of
the condensate and reflects the arising of its density. Moreover, we observed
that the thermalisation time strongly increases for values of the scattering
length around $-0.5R^{\star}$ (see Fig. 7) and thermalisation may not even be
achievable if the potential does not support any bound state. The former might
be related to the occurrence of a resonance as recently found experimentally
[15] and it will be a subject of future investigations. Furthermore, as
illustrated in Fig. 8, we found that a proper choice of the ion trap
parameters enables to cool the ion motion in a sodium gas to the quantum
regime, which affords prospects to quantum simulate impurity physics with
large atom-ion mass ratios. Finally, with the developed theoretical methods it
will be interesting to investigate how the gas quantum statistics affects the
Fock-state distribution of the ion motion, the coherence of ionic motional
superpositions, and to develop interferometric protocols for measuring the gas
temperature by reading out the thermal phonon distribution.

## Acknowledgements

This work is supported by the project “NE 1711/3-1” of the Deutsche
Forschungsgemeinschaft. We acknowledge H. Fürst, Z. Idziaszek, and K.
Jachymski for discussions. R.G. was supported by the Dutch Research Council
(Vrije Programma 680.92.18.05).

##  Appendix A Parameters of the regularised potential

To determine the parameters $b$ and $c$ of the regularised potential (4), we
follow the approach of Ref. [33]. Here, however, provide details that were not
discussed in that reference.

Since we have two free parameters, we need two physical conditions to
determine them. To this aim, we impose that:

  * (a)

The scattering length amplitude in first-order Born approximation (II) is
exactly equal to minus the atom-ion scattering length at zero-energy;

  * (b)

The potential (4) supports one bound state only.

The condition (b) is motivated by that fact that the energy separation between
bound states of the atom-ion polarisation potential (1) is rather large (order
of hundreds of $E^{\star}$), thus rendering very unlikely the population of
deeper bound states at typical atomic gases densities.

In the zero-energy limit the three-dimensional s-wave ion-atom scattering
length is indeed defined as

$\displaystyle a_{\mathrm{ia}}^{s}=-\lim_{q\rightarrow 0}f_{s}(q)$ (74)

with $f_{s}(q)$ being the full s-wave scattering amplitude at energy
$\hbar^{2}q^{2}/(2\mu)$, where $q=|\mathbf{k}-\mathbf{k}^{\prime}|$ is the
magnitude of the momentum transfer in the relative frame of reference. Hence,
the first aforementioned condition (a) reads

$\displaystyle a_{\mathrm{ia}}^{s}=-f(0).$ (75)

In the zero-energy limit $q\rightarrow 0$, and therefore, by expanding the
exponential functions in the last line of Eq. (II) to first order, we obtain

$\displaystyle f(0)=\pi(R^{\star})^{2}\frac{(b^{2}+2bc-c^{2})}{4b(b+c)^{2}}.$
(76)

Note that the scattering amplitude has the units of a length, which is
consistent with the definition (74).

The fulfillment of the second condition (b) is attained by determining the
s-wave scattering length as a function of either $b$ or $c$ by solving the
scattering problem at zero-energy. To this aim, we solve numerically the
radial time-independent Schrödinger equation

$\displaystyle\left[-\frac{\hbar^{2}}{2\mu}\frac{d^{2}}{dr^{2}}+V_{\mathrm{ai}}^{r}(r)\right]\psi(r)=0\qquad
r\in[0,+\infty),$ (77)

which in the $E^{\star}$ and $R^{\star}$ units reduces to

$\displaystyle\left[\frac{d^{2}}{dr^{2}}+\frac{r^{2}-c^{2}}{r^{2}+c^{2}}\frac{1}{(b^{2}+r^{2})^{2}}\right]\psi(r)=0.$
(78)

This differential equation is solved with boundary conditions $\psi(0)=0$ and
$\psi^{\prime}(0)=\epsilon$, where $\epsilon$ is a small number (e.g., 0.1).
We note, however, that the result does not rely on the particular numerical
value of $\epsilon$, as we have verified numerically. Thus, we fix the value
of the parameter $c$ (in units of $R^{\star}$) and we solve iteratively Eq.
(78) for different values of the parameter $b$ by evaluating the corresponding
scattering length, which becomes a function of $b$. We do the same for the
scattering amplitude (76) and thus search for the value of $b$ where
$a_{\mathrm{ia}}^{s}(b)$ and $-f(b)$ do cross, particularly where the first
zero-energy resonance occurs, which indicates that we have one bound state
only (see Fig. 9).

Let us now briefly explain some details about the numerical calculation of the
scattering length. First, we note that the solution to Eq. (78) behaves like
$r-a_{\mathrm{ia}}^{s}(b)$ at large distances, where the atom-ion interaction
vanishes. Hence, we have

$\displaystyle
a_{\mathrm{ai}}^{s}(b)=\lim_{r\rightarrow+\infty}\left[r-\frac{\psi(r)}{\psi^{\prime}(r)}\right].$
(79)

Numerically, we have noticed that a large grid size has to be chosen (a few
thousands of $R^{\star}$) such that the term on the right-hand side of the
limit (79) converges to a constant value, i.e. it is $r$-independent. An
example of such a calculation is given in Fig. 9.
For instance, by fixing $c=0.2239\,R^{\star}$, we find $b\simeq
0.078\,R^{\star}$, which yield $a_{\mathrm{ai}}^{s}\simeq 1.0054\,R^{\star}$.
Finally, let us remark that such a strategy relies crucially on the first-Born
approximation. Other strategies can be adopted in order to relax the latter
(see, e.g., Refs. [35, 36]). We chose, however, the method outlined above for
consistency, since we make use of the first-Born scattering amplitude in the
derivation of the master equation, as a consequence of the perturbative
description of the open system.

Figure 9:  (Color online). Atom-ion scattering length (blue line) computed via
Eq. (79) and minus the scattering amplitude (orange line) at zero-energy (76)
as a function of the $b$ parameter of the regularised atom-ion interaction
(4). Here, we have chosen $c=0.2039R^{\star}$. The two lines intersect at
$b\simeq 0.0770R^{\star}$, which yields a scattering length
$a_{\mathrm{ia}}^{s}\simeq 0.5385R^{\star}$.

##  Appendix B Thermal averages

The thermal averages of the double commutator (III.3) that yield non-zero
contributions are:

$\displaystyle\langle\tilde{\Gamma}^{\dagger}_{\mathbf{q}}(t)\tilde{\Gamma}_{\mathbf{q}^{\prime}}(t)\tilde{\Gamma}^{\dagger}_{\mathbf{k}}(t^{\prime})\tilde{\Gamma}_{\mathbf{k}^{\prime}}(t^{\prime})\rangle_{B_{0}}=e^{\frac{i}{\hbar}[(\varepsilon({\mathbf{q}})-\varepsilon({\mathbf{q}^{\prime}}))t+(\varepsilon({\mathbf{k}})-\varepsilon({\mathbf{k}^{\prime}}))t^{\prime}]}$
$\displaystyle\times[n_{\mathbf{q}}\delta_{\mathbf{q}^{\prime},\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}}+n_{\mathbf{q}}n_{\mathbf{q}^{\prime}}(\delta_{\mathbf{q}^{\prime},\mathbf{q}}\delta_{\mathbf{k},\mathbf{k}^{\prime}}+\delta_{\mathbf{q}^{\prime},\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}})]$
$\displaystyle\langle\tilde{\Gamma}^{\dagger}_{\mathbf{q}}(t)\tilde{\Gamma}_{\mathbf{q}^{\prime}}(t)\tilde{\Gamma}_{\mathbf{k}}(t^{\prime})\tilde{\Gamma}^{\dagger}_{\mathbf{k}^{\prime}}(t^{\prime})\rangle_{B_{0}}=e^{\frac{i}{\hbar}[(\varepsilon({\mathbf{q}})-\varepsilon({\mathbf{q}^{\prime}}))t+(\varepsilon({\mathbf{k}^{\prime}})-\varepsilon({\mathbf{k}}))t^{\prime}]}$
$\displaystyle\times
n_{\mathbf{q}}(1+n_{\mathbf{k}^{\prime}})(\delta_{\mathbf{q}^{\prime},\mathbf{q}}\delta_{\mathbf{k},\mathbf{k}^{\prime}}+\delta_{\mathbf{q}^{\prime},\mathbf{k}^{\prime}}\delta_{\mathbf{q},\mathbf{k}})$
$\displaystyle\langle\tilde{\Gamma}^{\dagger}_{\mathbf{q}}(t)\tilde{\Gamma}^{\dagger}_{\mathbf{q}^{\prime}}(t)\tilde{\Gamma}_{\mathbf{k}}(t^{\prime})\tilde{\Gamma}_{\mathbf{k}^{\prime}}(t^{\prime})\rangle_{B_{0}}=e^{\frac{i}{\hbar}[(\varepsilon({\mathbf{q}})+\varepsilon({\mathbf{q}^{\prime}}))t-(\varepsilon({\mathbf{k}^{\prime}})+\varepsilon({\mathbf{k}}))t^{\prime}]}$
$\displaystyle\times
n_{\mathbf{q}}n_{\mathbf{q}^{\prime}}(\delta_{\mathbf{q}^{\prime},\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}}+\delta_{\mathbf{q}^{\prime},\mathbf{k}^{\prime}}\delta_{\mathbf{q},\mathbf{k}})$
$\displaystyle\langle\tilde{\Gamma}_{\mathbf{q}}(t)\tilde{\Gamma}_{\mathbf{q}^{\prime}}(t)\tilde{\Gamma}^{\dagger}_{\mathbf{k}}(t^{\prime})\tilde{\Gamma}^{\dagger}_{\mathbf{k}^{\prime}}(t^{\prime})\rangle_{B_{0}}=e^{\frac{i}{\hbar}[(\varepsilon({\mathbf{k}})+\varepsilon({\mathbf{k}^{\prime}}))t^{\prime}-(\varepsilon({\mathbf{q}^{\prime}})+\varepsilon({\mathbf{q}}))t]}$
$\displaystyle\times(1+n_{\mathbf{q}}+n_{\mathbf{q}^{\prime}}+n_{\mathbf{k}}n_{\mathbf{k}^{\prime}})(\delta_{\mathbf{q}^{\prime},\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}}+\delta_{\mathbf{q}^{\prime},\mathbf{k}^{\prime}}\delta_{\mathbf{q},\mathbf{k}})$
$\displaystyle\langle\tilde{\Gamma}_{\mathbf{q}}(t)\tilde{\Gamma}^{\dagger}_{\mathbf{q}^{\prime}}(t)\tilde{\Gamma}^{\dagger}_{\mathbf{k}}(t^{\prime})\tilde{\Gamma}_{\mathbf{k}^{\prime}}(t^{\prime})\rangle_{B_{0}}=e^{\frac{i}{\hbar}[(\varepsilon({\mathbf{q}^{\prime}})-\varepsilon({\mathbf{q}}))t+(\varepsilon({\mathbf{k}})-\varepsilon({\mathbf{k}^{\prime}}))t^{\prime}]}$
$\displaystyle\times[n_{\mathbf{q}^{\prime}}\delta_{\mathbf{q},\mathbf{k}}\delta_{\mathbf{q}^{\prime},\mathbf{k}^{\prime}}+n_{\mathbf{k}}\delta_{\mathbf{q},\mathbf{q}^{\prime}}\delta_{\mathbf{k},\mathbf{k}^{\prime}}$
$\displaystyle+n_{\mathbf{q}^{\prime}}n_{\mathbf{k}}(\delta_{\mathbf{q}^{\prime},\mathbf{q}}\delta_{\mathbf{k},\mathbf{k}^{\prime}}+\delta_{\mathbf{q}^{\prime},\mathbf{k}^{\prime}}\delta_{\mathbf{q},\mathbf{k}})]$
$\displaystyle\langle\tilde{\Gamma}_{\mathbf{q}}(t)\tilde{\Gamma}^{\dagger}_{\mathbf{q}^{\prime}}(t)\tilde{\Gamma}_{\mathbf{k}}(t^{\prime})\tilde{\Gamma}^{\dagger}_{\mathbf{k}^{\prime}}(t^{\prime})\rangle_{B_{0}}=e^{\frac{i}{\hbar}[(\varepsilon({\mathbf{q}^{\prime}})-\varepsilon({\mathbf{q}}))t+(\varepsilon({\mathbf{k}^{\prime}})-\varepsilon({\mathbf{k}}))t^{\prime}]}$
$\displaystyle\times[n_{\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}}\delta_{\mathbf{q}^{\prime},\mathbf{k}}+(1+n_{\mathbf{k}}+n_{\mathbf{q}})\delta_{\mathbf{q},\mathbf{q}^{\prime}}\delta_{\mathbf{k},\mathbf{k}^{\prime}}$
$\displaystyle+n_{\mathbf{q}^{\prime}}n_{\mathbf{k}^{\prime}}(\delta_{\mathbf{q}^{\prime},\mathbf{q}}\delta_{\mathbf{k},\mathbf{k}^{\prime}}+\delta_{\mathbf{q}^{\prime},\mathbf{k}}\delta_{\mathbf{q},\mathbf{k}^{\prime}})].$
(80)

For these identities we used the relation (4.7) of Ref. [63].

##  Appendix C Ion motion in a Paul trap

Here, we provide details on the analytical solution of the ion motion in a
Paul trap using the notation of the review [18]. The goal is to provide the
relevant steps of its derivation such that the interested reader can
implemented them in numerics quickly and efficiently.

###  C.1 Classical solution of a charge in a Paul trap

Let us consider a particle of mass $M$ and charge $Z|e|$ in the quadruple
field

$\displaystyle\Phi(x,y,z,t)$ $\displaystyle=\frac{U}{2}(\alpha x^{2}+\beta
y^{2}+\gamma z^{2})$
$\displaystyle+\frac{\tilde{U}}{2}\cos(\Omega_{rf}t)(\alpha^{\prime}x^{2}+\beta^{\prime}y^{2}+\gamma^{\prime}z^{2}).$
(81)

For a linear Paul trap we have: $0<\gamma=-(\alpha+\beta)$,
$\alpha^{\prime}=-\beta^{\prime}$ and $\gamma^{\prime}=0$ and Poisson
equation, $\nabla^{2}\Phi=0$, is fulfilled. The Newton equation of motion
along the $x$-direction is given by (similarly for the other directions):

$\displaystyle\ddot{x}(t)=-\frac{|e|Z}{m}\frac{\partial\Phi}{\partial
x}=-\frac{|e|Z}{m}\left[U\alpha+\tilde{U}\alpha^{\prime}\cos(\Omega_{rf}t)\right]x.$

By introducing the dimensionless variable $\tau=\Omega_{rf}t/2$, the
corresponding rescaled equation reads

$\displaystyle\ddot{x}(\tau)+\left[a_{x}-2q_{x}\cos(2\tau)\right]x=0,$ (82)

where the newly introduced parameters are defined as: $a_{x}=4|e|Z\alpha
U/(M\Omega_{rf}^{2})$ and
$q_{x}=-2|e|Z\alpha^{\prime}\tilde{U}/(M\Omega_{rf}^{2})$. Note that for a
linear Paul trap we have: $q_{y}=-q_{x}\equiv q$ and $q_{z}=0$, where
$a_{y}=a_{x}=-a_{z}/2\equiv a$. The above outlined equation is solved by using
the following ansatz

$\displaystyle
x(\tau)=A_{x}e^{i\beta_{x}\tau}\\!\\!\\!\\!\sum_{n=-\infty}^{\infty}C_{2n}^{x}e^{i2n\tau}+B_{x}e^{-i\beta_{x}\tau}\\!\\!\\!\\!\sum_{n=-\infty}^{\infty}C_{2n}^{x}e^{-i2n\tau},$
(83)

where $A_{x},\,B_{x}$ are constants that depend on the initial conditions,
while the parameter $\beta_{x}$ and coefficients $C_{2n}^{x}$ have to be
determined recursively. Here, we use the same notation of Ref. [18], but we
note that the coefficients $C_{2n}^{x}$ could have been named $C_{n}^{x}$, as
we actually do in Eq. (IV.1). To this end, we insert the ansatz into Eq. (82)
and we obtain

$\displaystyle
C_{2n+2}^{x}-D_{2n}^{x}C_{2n}^{x}+C_{2n-2}^{x}=0\,\,\,\,\,\,D_{2n}^{x}=\frac{a_{x}-(\beta_{x}+2n)^{2}}{q_{x}}.$
(84)

Iterative application of the above identities yield the continued fraction
solutions

$\displaystyle\frac{C_{2n}^{x}}{C_{2n+2}^{x}}=\frac{1}{D_{2n}^{x}-\frac{1}{D_{2n-2}^{x}-\frac{1}{D_{2n-4}^{x}-\dots}}},$
$\displaystyle\frac{C_{2n}^{x}}{C_{2n-2}^{x}}=\frac{1}{D_{2n}^{x}-\frac{1}{D_{2n+2}^{x}-\frac{1}{D_{2n+4}^{x}-\dots}}}.$
(85)

With these expressions and Eq. (84), we obtain

$\displaystyle
D_{2n}^{x}=\frac{1}{D_{2n-2}^{x}-\frac{1}{D_{2n-4}^{x}\dots}}+\frac{1}{D_{2n+2}^{x}-\frac{1}{D_{2n+4}^{x}-\dots}}.$
(86)

Since $D_{0}^{x}=(a_{x}-\beta_{x}^{2})/q_{x}$ and $D_{2n}^{x}$ above, we have

$\displaystyle\beta_{x}^{2}=a_{x}-q_{x}\left[\frac{1}{D_{-2}^{x}-\frac{1}{D_{-4}^{x}\dots}}+\frac{1}{D_{2}^{x}-\frac{1}{D_{4}^{x}-\dots}}\right].$
(87)

Note that the expression in the brackets $[\dots]$ of Eq. (87) still depends
on $\beta_{x}$, $a_{x}$ and $q_{x}$. Nonetheless, given $a_{x}$ and $q_{x}$,
Eq. (87) can be solved with respect to $\beta_{x}$ after a few iterations of
the continued fraction.

For the numerical assessment of the coefficients $C_{2n}^{x}$ we proceed as
follows: First, we set the maximum of the $n$-index to some positive integer
$N_{F}$ and $C_{0}^{x}=1$ such that $C_{\pm 2n}^{x}=0$ $\forall n>N_{F}$.
Second, using the previously obtained value of $\beta_{x}$ and the formula
(84), we solve an inhomogeneous linear system of equations given by

$\displaystyle\left(\begin{array}[]{cccccccccccc}D_{2N_{F}}^{x}&-1&0&0&0&\dots&\dots&\dots&\dots&\dots&\dots&\dots\\\
-1&D_{2(N_{F}-1)}^{x}&-1&0&0&\dots&\dots&\dots&\dots&\dots&\dots&\dots\\\
0&-1&D_{2(N_{F}-2)}^{x}&-1&0&\dots&\dots&\dots&\dots&\dots&\dots&\dots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\vdots&\dots&\dots&\dots&0&-1&D_{2}^{x}&0&0&\dots&\dots&\dots\\\
\vdots&\dots&\dots&\dots&0&0&0&D_{-2}^{x}&-1&0&\dots&\dots\\\
\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots&\vdots\\\
\vdots&\dots&\dots&\dots&\dots&\dots&\dots&\dots&&\dots
0&-1&D_{-2N_{F}}^{x}\\\
\end{array}\right)\left(\begin{array}[]{c}C_{2N_{F}}^{x}\\\
C_{2(N_{F}-1)}^{x}\\\ C_{2(N_{F}-2)}^{x}\\\ \vdots\\\ C_{2}^{x}\\\
C_{-2}^{x}\\\ \vdots\\\ C_{-2N_{F}}^{x}\\\
\end{array}\right)=\left(\begin{array}[]{c}0\\\ 0\\\ 0\\\ \vdots\\\ 1\\\ 1\\\
\vdots\\\ 0\\\ \end{array}\right).$ (112)

The total number of numerically calculated coefficients is therefore $2N_{F}$.
Moreover, the normalisation condition

$\displaystyle\sum_{n=-\infty}^{\infty}C_{2n}^{x}\simeq\sum_{n=-N_{F}}^{N_{F}}C_{2n}^{x}=1$
(113)

has to be satisfied, from which we obtain the final normalised coefficients
$c_{2n}^{x}=C_{2n}^{x}/\sum_{j=-N_{F}}^{N_{F}}C_{2j}^{x}$. In this way we can
then immediately evaluate the reference harmonic oscillator frequency

$\displaystyle\nu_{x}$
$\displaystyle=\Omega_{rf}\sum_{n=-N_{F}}^{N_{F}}c_{2n}^{x}\left(\frac{\beta_{x}}{2}+n\right).$
(114)

Exactly the same procedure applies to the determination of the coefficients
$C_{2n}^{y,z},\,c_{2n}^{y,z}$ and frequencies $\nu_{y,z}$. We note that the
coefficients appearing in the ion solution (IV.1) are precisely the
coefficients $c_{2n}^{x,y,z}$.

Finally, the classical trajectory is obtained as $x(\tau=0)=A_{x}+B_{x}\equiv
x_{0}$ and

$\displaystyle\dot{x}(\tau=0)=i(A_{x}-B_{x})\left[\beta_{x}+2\sum_{n=-N_{F}}^{N_{F}}n\,c_{2n}^{x}\right]\equiv\frac{p_{0}}{M}.$

From these equalities we obtain

$\displaystyle
A_{x}=\frac{x_{0}}{2}-i\frac{p_{0}}{2M}\left[\beta_{x}+2\sum_{n=-N_{F}}^{N_{F}}n\,c_{2n}^{x}\right]^{-1},$
(115) $\displaystyle
B_{x}=\frac{x_{0}}{2}+i\frac{p_{0}}{2M}\left[\beta_{x}+2\sum_{n=-N_{F}}^{N_{F}}nc\,_{2n}^{x}\right]^{-1}.$
(116)

In the limit $a_{x}<|q_{x}|^{2}\ll 1$ and for $p_{0}=0$, an approximated
solution is given by

$\displaystyle
x_{\text{apx}}(\tau)=\frac{2x_{0}}{2-q_{x}}\cos(\beta_{x}\tau)\left[1-\frac{q_{x}}{2}\cos(2\tau)\right].$
(117)

This solution can be compared with the numerically obtained one from Eq. (83).
In Fig. 10 an example is shown, where we compare the solutions
$x_{\text{apx}}(\tau)$ (thin black line), Eq. (83) (blue slighter ticker line
than the black line), and the numerically solved Eq. (82) (tick yellow line).

Figure 10:  (Color online). Comparison of the numerically exact solution
obtained by solving the Newton equation of motion (82), the analytical
solution (83), and the approximated one (117). We have chosen the following
parameters: $a_{x}=0$, $q_{x}=0.28$, and $N_{F}=20$.

###  C.2 Quantum Hamiltonian

As it can be verified, Eq. (82) reproduces the motion of a parametric harmonic
oscillator with squared frequency

$\displaystyle
W_{\xi}(t)=\frac{\Omega_{rf}^{2}}{4}\left[a_{\xi}-2q_{\xi}\cos(\Omega_{rf}t)\right].$
(118)

As it can be shown formally with the Heisenberg equations of motion for
$\hat{r}_{\xi}$ and $\hat{p}_{\xi}$ [18], the following Hamiltonian in one
spatial direction

$\displaystyle\hat{H}_{I}^{trap}=\frac{\hat{p}_{\xi}^{2}}{2M}+\frac{M}{2}W_{\xi}(t)\hat{r}_{\xi}^{2}\qquad\xi=x,y,z,$
(119)

reproduces exactly the same equation of motion (82) for the operators
$\hat{r}_{\xi}$. Thus, Eq. (119) is the quantised version of the classical
Hamiltonian for an ion in a Paul trap, where we have promoted the ion position
and momentum variables to operators. Besides this, we note that the spatial
directions are uncoupled because of the form of the quadrupole field (C.1).

The eigenfunctions of $\hat{H}_{I}^{trap}$ are given by:

$\displaystyle\psi_{n}(r_{\xi},t)$
$\displaystyle=\left(\frac{M\nu_{\xi}}{\pi\hbar}\right)^{\frac{1}{4}}\frac{e^{-in\nu_{\xi}t+i\frac{M}{2\hbar}\frac{\dot{u}_{\xi}(t)}{u_{\xi}(t)}r_{\xi}^{2}}}{\sqrt{2^{n}n!u_{\xi}(t)}}$
$\displaystyle\times
H_{n}\left(\sqrt{\frac{M\nu_{\xi}}{\hbar|u_{\xi}(t)|^{2}}}r_{\xi}\right)$
(120)

where

$\displaystyle u_{\xi}(t)$
$\displaystyle=e^{i\beta_{\xi}\Omega_{rf}t/2}\sum_{n=-\infty}^{\infty}C_{2n}^{\xi}e^{in\Omega_{rf}t},$
(121)

with $u_{\xi}(0)=\sum_{n=-\infty}^{\infty}C_{2n}^{\xi}=1$,
$\dot{u}_{\xi}(0)=i\nu_{\xi}$, and $\nu_{\xi}$ given by Eq. (114).

Finally, as initial condition of the ion density matrix for the solution of
the master equation we have chosen

$\displaystyle\hat{\rho}_{\xi}=|\psi_{0}(t=T_{rf})\rangle\langle\psi_{0}(t=T_{rf})|,\qquad\forall\xi=x,y,z,$
(122)

since by starting from $|\psi_{0}(t=0)\rangle\langle\psi_{0}(t=0)|$ and by
computing the expectation value (68), one can numerically verify that the
minimum of the ion energy occurs at the time $t=2\pi/\Omega_{rf}$. This is the
energy minimum we assumed in our analyses and Eq. (122) as initial condition.
With that initial matrix we have computed the initial conditions for the
moments of Sec. V.

##  Appendix D Ion master equation details

In this section, we provide a few technical details on the calculation of the
double integral as a consequence of the double summation in momentum space in
Eq. (IV.1). The derivations are detailed for a bosonic bath only, while for a
fermionic one we simply provide the final result, as they are very similar.
Additionally, we provide details of the analytical calculation of the Cauchy
principal value within the Fröhlich model that yields the Lamb-shift.

###  D.1 Double integration in momentum space

In Eq. (IV.1) we have to evaluate terms of the kind

$\displaystyle\sum_{\mathbf{q},\mathbf{q}^{\prime}}n_{\mathbf{q}}(n_{\mathbf{q}^{\prime}}+1)(q_{\xi}^{\prime}-q_{\xi})^{2}|c_{\mathbf{q}^{\prime}-\mathbf{q}}|^{2}\delta(\omega_{\mathbf{q}^{\prime}}-\omega_{\mathbf{q}}\pm\omega_{s,\xi}),$
(123)

where $\omega_{s,\xi}=\Omega_{rf}(\beta_{\xi}/2+s)$ and
$\omega_{\mathbf{q}}\equiv\varepsilon(\mathbf{q})/\hbar$. To this end, we
first perform the centre-of-mass and relative coordinate transformation

$\displaystyle\mathbf{k}=\mathbf{q}^{\prime}-\mathbf{q},\qquad\qquad\mathbf{K}=\frac{\mathbf{q}^{\prime}+\mathbf{q}}{2},$
(124)

with

$\displaystyle\mathbf{q}^{\prime}$
$\displaystyle=\mathbf{K}+\frac{\mathbf{k}}{2},\qquad\qquad\mathbf{q}=\mathbf{K}-\frac{\mathbf{k}}{2},$
(125)

and therefore we have

$\displaystyle\omega_{\mathbf{q}^{\prime}}-\omega_{\mathbf{q}}=\frac{\hbar}{2m}\left[(\mathbf{q}^{\prime})^{2}-\mathbf{q}^{2}\right]=\frac{\hbar}{m}(\mathbf{k}\cdot\mathbf{K}).$
(126)

Given this, we can rewrite Eq. (123) as

$\displaystyle\sum_{\mathbf{k},\mathbf{K}}n_{\mathbf{k},\mathbf{K}}^{-}(n_{\mathbf{k},\mathbf{K}}^{+}+1)k_{\xi}^{2}|c_{\mathbf{k}}|^{2}\delta\left(\frac{\hbar}{m}\mathbf{k}\cdot\mathbf{K}\pm\omega_{s,\xi}\right),$
(127)

where

$\displaystyle n_{\mathbf{k},\mathbf{K}}^{-}$
$\displaystyle=\frac{1}{\exp\left\\{\beta_{\mathcal{T}}\left[\frac{\hbar^{2}}{2m}\left(\mathbf{K}^{2}+\frac{\mathbf{k}^{2}}{4}-\mathbf{k}\cdot\mathbf{K}\right)-\mu_{\mathrm{G}}\right]\right\\}-1},$
$\displaystyle n_{\mathbf{k},\mathbf{K}}^{+}$
$\displaystyle=\frac{1}{\exp\left\\{\beta_{\mathcal{T}}\left[\frac{\hbar^{2}}{2m}\left(\mathbf{K}^{2}+\frac{\mathbf{k}^{2}}{4}+\mathbf{k}\cdot\mathbf{K}\right)-\mu_{\mathrm{G}}\right]\right\\}-1}.$
(128)

We remind that $\mu_{\mathrm{G}}$ is the bosons’ chemical potential and
$\beta_{\mathcal{T}}=(k_{\mathrm{B}}\mathcal{T})^{-1}$. Thus, we perform the
continuum limit

$\displaystyle\sum_{\mathbf{k},\mathbf{K}}\rightarrow\frac{L^{6}}{(2\pi)^{6}}\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{k}\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{K},$
(129)

which transforms the double sum in Eq. (127) in the following two double
integrals

$\displaystyle\mathcal{I}_{1}+\mathcal{I}_{2}=\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{k}\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{K}\,n_{\mathbf{k},\mathbf{K}}^{-}k_{\xi}^{2}|c_{\mathbf{k}}|^{2}\delta\left(\frac{\hbar}{m}\mathbf{k}\cdot\mathbf{K}\pm\omega_{s,\xi}\right)$
$\displaystyle+\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{k}\int_{\mathbb{R}^{3}}\mathrm{d}\mathbf{K}\,n_{\mathbf{k},\mathbf{K}}^{-}n_{\mathbf{k},\mathbf{K}}^{+}k_{\xi}^{2}|c_{\mathbf{k}}|^{2}\delta\left(\frac{\hbar}{m}\mathbf{k}\cdot\mathbf{K}\pm\omega_{s,\xi}\right).$
(130)

Here, we have neglected the common factor $\left(\frac{L}{2\pi}\right)^{6}$.
In order to solve them, we first move to spherical coordinates

$\displaystyle k_{x}$ $\displaystyle=k\sin(\theta_{k})\cos(\varphi_{k}),$
$\displaystyle k_{y}$ $\displaystyle=k\sin(\theta_{k})\sin(\varphi_{k}),$
$\displaystyle k_{z}$ $\displaystyle=k\cos(\theta_{k}).$ (131)

Thus, the corresponding volume element is given by
$\mathrm{d}\mathbf{k}=\mathrm{d}k\mathrm{d}\theta_{k}\mathrm{d}\varphi_{k}k^{2}\sin(\theta_{k})$
with $k\equiv|\mathbf{k}|$. While a similar change of variables applies to the
centre-of-mass variable $\mathbf{K}$ as well, but with subscript $K$ for the
angular variables, we note that we choose the “z-axis” of the vector
$\mathbf{K}$ along the relative variable $\mathbf{k}$. In such a way the
scalar product appearing in the Dirac’s delta can be written as

$\displaystyle\mathbf{k}\cdot\mathbf{K}=kK\cos(\theta_{K}).$ (132)

Given this, the first integral becomes

$\displaystyle\mathcal{I}_{1}=2\pi^{2}(1+\delta_{z,\xi})\int_{0}^{\infty}\mathrm{d}kk^{4}\int_{0}^{\pi}\mathrm{d}\theta_{k}g_{\xi}^{2}(\theta_{k})\sin(\theta_{k})|c_{\mathbf{k}}|^{2}$
$\displaystyle\times\int_{0}^{\infty}\mathrm{d}KK^{2}\int_{0}^{\pi}\mathrm{d}\theta_{K}\sin(\theta_{K})n_{\mathbf{k},\mathbf{K}}^{-}$
$\displaystyle\times\delta\left(\frac{\hbar}{m}kK\cos(\theta_{K})\pm\omega_{s,\xi}\right),$
(133)

where we have performed the integrations of the variables $\varphi_{k}$ and
$\varphi_{K}$, since only $k_{\xi}$ relies on $\varphi_{k}$, while none of the
functions in the integrand depend on $\varphi_{K}$. Besides this, because of
$k_{\xi}$ we have introduced the angle function
$g_{\xi}(\theta_{k})=\delta_{z,\xi}\cos(\theta_{k})+(1-\delta_{z,\xi})\sin(\theta_{k})$
with $\delta_{z,\xi}$ being the Kronecker delta. Since neither
$|c_{\mathbf{k}}|^{2}$ nor $n_{\mathbf{k},\mathbf{K}}^{-}$ rely on
$\theta_{k}$, we can easily perform the integration

$\displaystyle\int_{0}^{\pi}\mathrm{d}\theta_{k}g_{\xi}^{2}(\theta_{k})\sin(\theta_{k})=\frac{2}{3}\left[\delta_{z,\xi}+2(1-\delta_{z,\xi})\right].$
(134)

Our next step is to integrate out the variable $K$. Towards this end, we first
rewrite the Dirac’s delta as

$\displaystyle\delta\left(\frac{\hbar}{m}kK\cos(\theta_{K})\pm\omega_{s,\xi}\right)=\frac{m}{\hbar
k|\cos(\theta_{K})|}$ $\displaystyle\times\delta\left(K\pm
K_{s,\xi}(\theta_{K},k)\right)\delta_{0,1\pm\mathrm{sgn}(\omega_{s,\xi}\cos(\theta_{K}))},$
(135)

where $\mathrm{sgn}(\cdot)$ is the sign function. The last Kronecker delta
ensures that $K_{s,\xi}(\theta_{K},k)=\frac{m\omega_{s,\xi}}{\hbar
k\cos(\theta_{K})}\geq 0$ in the minus case and $K_{s,\xi}(\theta_{K},k)\leq
0$ in the plus case, since $K\in\mathbb{R}^{+}$ and the integral over $K$
would be zero otherwise, and so would be $\mathcal{I}_{1}$. Hence, we obtain

$\displaystyle\mathcal{I}_{1}=\frac{8}{3}\pi^{2}\left(\frac{m}{\hbar}\right)^{3}\omega_{s,\xi}^{2}\int_{0}^{\infty}\mathrm{d}kk|c_{\mathbf{k}}|^{2}$
$\displaystyle\times\int_{0}^{\pi}\mathrm{d}\theta_{K}\frac{\tan(\theta_{K})n_{\mathbf{k},\mp\mathbf{K}_{s,\xi}(\theta_{K},k)}^{-}}{\cos(\theta_{K})|\cos(\theta_{K})|}\delta_{0,1\pm\mathrm{sgn}(\omega_{s,\xi}\cos(\theta_{K}))}.$
(136)

Finally, we perform the angular integral

$\displaystyle\int_{0}^{\pi}\mathrm{d}\theta_{K}\frac{\tan(\theta_{K})n_{\mathbf{k},\mp\mathbf{K}_{s,\xi}(\theta_{K},k)}^{-}}{\cos(\theta_{K})|\cos(\theta_{K})|}\delta_{0,1\pm\mathrm{sgn}(\omega_{s,\xi}\cos(\theta_{K}))}=$
$\displaystyle\int_{-1}^{1}\mathrm{d}u\frac{n_{\mathbf{k},\mp\mathbf{K}_{n,\xi}(u,k)}^{-}}{u^{2}|u|}\delta_{0,1\pm\mathrm{sgn}(\omega_{s,\xi}u)},$
(137)

where we performed the change of variable $u=\cos(\theta_{K})$. Hence,

$\displaystyle n_{\mathbf{k},\mathbf{K}_{s,\xi}(u,k)}^{-}$
$\displaystyle=\frac{1}{e^{\frac{\beta\hbar\omega_{s,\xi}}{2(k\ell_{s,\xi})^{2}}\left[\frac{1}{u^{2}}+\frac{(k\ell_{s,\xi})^{4}}{4}\pm(k\ell_{s,\xi})^{2}\right]-\beta\mu_{\mathrm{G}}}-1}$
(138)

with $\ell_{s,\xi}^{2}=\frac{\hbar}{m\omega_{s,\xi}}$ and (D.1) can be
rewritten as

$\displaystyle\mathcal{J}=\int_{0}^{1}\mathrm{d}u\frac{n_{\mathbf{k},\mathbf{K}_{s,\xi}(u,k)}^{-}}{u^{3}}=\int_{0}^{1}\frac{\mathrm{d}u}{u^{3}\left(e^{\frac{\alpha_{0}}{u^{2}}+\alpha_{\pm}}-1\right)},$
(139)

where

$\displaystyle\alpha_{0}$
$\displaystyle=\frac{\beta_{\mathcal{T}}\hbar\omega_{s,\xi}}{2(k\ell_{s,\xi})^{2}},$
$\displaystyle\alpha_{\pm}$
$\displaystyle=\alpha_{0}\left[\frac{(k\ell_{s,\xi})^{4}}{4}\pm(k\ell_{s,\xi})^{2}\right]-\beta_{\mathcal{T}}\mu_{\mathrm{G}}.$
(140)

In order to solve $\mathcal{J}$, we perform the change of variable:

$\displaystyle
z=\frac{\alpha_{0}}{u^{2}}+\alpha_{\pm}\,\,\Longrightarrow\,\,\mathrm{d}u=-\frac{\sqrt{\alpha_{0}}}{2}\frac{\mathrm{d}z}{(z-\alpha_{\pm})^{3/2}}.$
(141)

In such a way $\mathcal{J}$ is rewritten as

$\displaystyle\mathcal{J}=\frac{1}{2\alpha_{0}}\int_{\alpha_{0}+\alpha_{\pm}}^{+\infty}\frac{\mathrm{d}z}{e^{z}-1}=\frac{\alpha_{0}+\alpha_{\pm}-\ln(e^{\alpha_{0}+\alpha_{\pm}}-1)}{2\alpha_{0}},$

which holds as long as $\alpha_{0}>0$, as it is indeed the case. Hence, the
integral (D.1) is given by

$\displaystyle\left(\frac{L}{2\pi}\right)^{6}\mathcal{I}_{1}^{\pm}=\frac{1}{6}\left(\frac{m}{\pi\mu}\right)^{2}\frac{\hbar}{\beta_{\mathcal{T}}}\,\mathcal{F}_{s,\xi}^{(1),\pm}$
(142)

with

$\displaystyle\mathcal{F}_{s,\xi}^{(1),\pm}=\int_{0}^{\infty}\mathrm{d}k\,k^{3}\frac{\alpha_{0}(k)+\alpha_{\pm}(k)-\ln(e^{\alpha_{0}(k)+\alpha_{\pm}(k)}-1)}{|f(k)|^{-2}}.$
(143)

The integral over $k$ in Eq. (143) is computed numerically. We note that
$\alpha_{0}(k)\sim k^{-2}$ and therefore the exponential diverges for
$k\rightarrow 0$, which is not the case for $\alpha_{\pm}(k)\sim k^{2}$ that
tends to zero. Because of the logarithm, however, the exponent of the
exponential function compensates the -$\alpha_{0}(k)$ on the left-hand-side of
the logarithm so that the overall behaviour of the integrand is zero when
$k\rightarrow 0$. Instead, when $k\rightarrow\infty$, we have
$\alpha_{0}(k)\rightarrow 0$, while $\alpha_{\pm}(k)$ diverges. For the same
argument as before, the function of the integrand numerator tends to zero.
Therefore, the integral converges, even if $|f(k)|^{2}=1$. For the fermionic
bath, we get the expression

$\displaystyle\mathcal{F}_{s,\xi}^{(1),\pm}=\int_{0}^{\infty}\mathrm{d}k\,k^{3}\frac{\ln(1+e^{\alpha_{0}(k)+\alpha_{\pm}(k)})-\alpha_{0}(k)-\alpha_{\pm}(k)}{|f(k)|^{-2}}.$
(144)

The result is very similar to the bosonic case, but one has to remember that
the chemical potentials are different, especially for temperature below the
Fermi temperature and the critical temperature for condensation.

The integral $\mathcal{I}_{2}$ is almost the same, but $\mathcal{J}$ in Eq.
(139) is defined now as

$\displaystyle\mathcal{J}$
$\displaystyle=\int_{0}^{1}\mathrm{d}u\frac{n_{\mathbf{k},\mathbf{K}_{s,\xi}(u,k)}^{-}n_{\mathbf{k},\mathbf{K}_{s,\xi}(u,k)}^{+}}{u^{3}}$
$\displaystyle=\int_{0}^{1}\frac{\mathrm{d}u}{u^{3}\left(e^{\frac{\alpha_{0}}{u^{2}}+\alpha_{-}}-1\right)\left(e^{\frac{\alpha_{0}}{u^{2}}+\alpha_{+}}-1\right)}.$
(145)

To solve it, we first perform the change of variable (141), which yields

$\displaystyle\mathcal{J}$
$\displaystyle=\frac{1}{2\alpha_{0}}\int_{\alpha_{0}+\alpha_{-}}^{+\infty}\frac{\mathrm{d}z}{(e^{z}-1)(e^{z+\alpha_{+}-\alpha_{-}}-1)}.$
(146)

Thus, we perform the additional change of variable $y=e^{z}$ with d$z$=d$y/y$,
and we obtain

$\displaystyle\mathcal{J}$
$\displaystyle=\frac{1}{2\alpha_{0}e^{\alpha_{+}-\alpha_{-}}}\int_{a}^{+\infty}\frac{\mathrm{d}y}{y(y-1)(y-b)},$
(147)

with $a=e^{\alpha_{0}+\alpha_{-}}$ and $b=e^{\alpha_{-}-\alpha_{+}}$. The
integral in Eq. (147) can be solved analytically, which finally gives

$\displaystyle\mathcal{J}$
$\displaystyle=\frac{\ln(e^{\alpha_{0}+\alpha_{-}}-e^{\alpha_{-}-\alpha_{+}})-\alpha_{0}-\alpha_{-}}{2\alpha_{0}(1-e^{\alpha_{-}-\alpha_{+}})}$
$\displaystyle-\frac{e^{\alpha_{-}-\alpha_{+}}\Re\\{\mathrm{arctanh}(1-2e^{\alpha_{0}+\alpha_{-}})\\}}{\alpha_{0}(1-e^{\alpha_{-}-\alpha_{+}})}.$
(148)

Hence, $\mathcal{I}_{2}$ preserves the structure of $\mathcal{I}_{1}$, that
is,

$\displaystyle\left(\frac{L}{2\pi}\right)^{6}\mathcal{I}_{2}=\frac{1}{6}\left(\frac{m}{\pi\mu}\right)^{2}\frac{\hbar}{\beta_{\mathcal{T}}}\,\mathcal{F}_{s,\xi}^{(2)},$
(149)

where the radial integration in momentum space is given by

$\displaystyle\mathcal{F}_{s,\xi}^{(2)}=\int_{0}^{\infty}\mathrm{d}k\,k^{3}|f(k)|^{2}\,\left[\frac{\ln(e^{\alpha_{0}+\alpha_{-}}-e^{\alpha_{-}-\alpha_{+}})-\alpha_{0}-\alpha_{-}-2e^{\alpha_{-}-\alpha_{+}}\Re\\{\mathrm{arctanh}(1-2e^{\alpha_{0}+\alpha_{-}})\\}}{1-e^{\alpha_{-}-\alpha_{+}}}\right].$
(150)

For the fermions we obtain a similar expression

$\displaystyle\mathcal{F}_{s,\xi}^{(2)}=\int_{0}^{\infty}\mathrm{d}k\,k^{3}|f(k)|^{2}\,\left[\frac{e^{\alpha_{-}-\alpha_{+}}\ln(1+e^{\alpha_{0}+\alpha_{-}})-\ln(e^{\alpha_{-}-\alpha_{+}}+e^{\alpha_{0}+\alpha_{-}})}{e^{\alpha_{-}-\alpha_{+}}-1}-\alpha_{0}(k)-\alpha_{-}(k)\right].$
(151)

We note that in the numerical assessment of $\mathcal{F}_{s,\xi}^{(2)}$ we
found that both for the bosons and the fermions the integral is essentially
zero, as a consequence of the large numerical values taken by the exponents.
For high densities and high temperatures, however, $\mathcal{F}_{s,\xi}^{(2)}$
is not negligible anymore.

###  D.2 Cauchy principal value calculation

We need to compute the Cauchy principal value

$\displaystyle\mathcal{I}=\mathcal{P}\int_{\mathbb{R}^{3}}d\mathbf{q}\frac{f^{2}(q)q_{\xi}q_{\xi^{\prime}}}{\omega_{0}-\omega}\qquad\xi,\xi^{\prime}=1,2,3,$
(152)

where $\omega_{0}\equiv|\Omega_{rf}(\beta_{\xi}/2+s)|$,
$\omega\equiv\varepsilon(\mathbf{q})/\hbar$, and $f(q)$ given by Eq. (II).
Henceforth, we neglected the imaginary unit $i$ in Eq. (54). The above
integral can be rewritten as [$q_{s,\xi}$ is defined in Eq. (57)]:

$\displaystyle\mathcal{I}_{-}(q_{s,\xi})=\frac{2m}{\hbar}\mathcal{P}\int_{\mathbb{R}^{3}}d\mathbf{q}\,\frac{f^{2}(q)q_{\xi}q_{\xi^{\prime}}}{|q_{s,\xi}^{2}|-q^{2}}$
$\displaystyle=\frac{2m}{\hbar}\mathcal{P}\int_{0}^{\infty}dqq^{2}\int_{0}^{\pi}d\theta\sin(\theta)\int_{0}^{2\pi}d\varphi\frac{f^{2}(q)q_{\xi}q_{\xi^{\prime}}}{|q_{s,\xi}^{2}|-q^{2}}$
(153)

where we transformed the wave vector in spherical coordinates (D.1). The only
angular dependence in the integrand comes from $q_{\xi}q_{\xi^{\prime}}$, as
the other functions rely on $q$ only. Given this, one can verify that the
angular integral of $q_{\xi}q_{\xi^{\prime}}$ always vanishes for
$\xi\neq\xi^{\prime}$. Thus, we just need to compute the CPV for each
direction separately. The angular part of the integration is the same for all
directions, namely it yields $4\pi/3$. Thus, we need to assess the integral

$\displaystyle\mathcal{I}_{-}(q_{s,\xi})$
$\displaystyle=\frac{8m\pi}{3\hbar}\mathcal{P}\int_{0}^{\infty}dq\,\frac{f^{2}(q)q^{4}}{|q_{s,\xi}^{2}|-q^{2}}=\frac{8m\pi}{3\hbar}{\mathcal{J}}_{-}^{\prime}(q_{s,\xi}).$
(154)

In case the denominator of the integrand has a plus sign, we have an
integration without any singularity:

$\displaystyle\mathcal{I}_{+}(q_{s,\xi})$
$\displaystyle=-\frac{2m}{\hbar}\int_{\mathbb{R}^{3}}d\mathbf{q}\,\frac{f^{2}(q)q_{\xi}^{2}}{|q_{s,\xi}^{2}|+q^{2}}=-\frac{8m\pi}{3\hbar}{\mathcal{J}}_{+}^{\prime}(q_{s,\xi}).$
(155)

This integration can be in principle carried out analytically, but we refrain
to provide an expression, since it is quite involved and it includes Meijer
G-functions. A similar argument holds for the integrals involving the bosonic
occupation number, that is,

$\displaystyle\mathcal{I}_{\pm}^{n_{q}}(q_{s,\xi})$
$\displaystyle=\mp\frac{2m}{\hbar}\int_{\mathbb{R}^{3}}d\mathbf{q}\,\frac{f^{2}(q)q_{\xi}^{2}}{|q_{s,\xi}^{2}|\pm
q^{2}}\frac{1}{e^{\left[\beta_{\mathcal{T}}(\varepsilon(\mathbf{q})-\mu_{\text{G}})\right]}-1}$
$\displaystyle=\mp\frac{8m\pi}{3\hbar}{\mathcal{J}}_{\pm}^{n_{q}}(q_{s,\xi}).$
(156)

To compute them, we performed a numerical integration with Mathematica.

##  Appendix E Rescaled equations of motion

In our numerical simulations we solve the differential equations (V) in
rescaled units, namely we rescaled the time and space variables with respect
to $\nu_{\xi}$ and $l_{\xi}=\sqrt{\hbar/(M\nu_{\xi})}$, respectively, that is,
with respect to the reference harmonic oscillator frequency and length scale.
Hence we have: $\bar{r}_{\xi}^{2}=(\langle r_{\xi}\rangle/l_{\xi})^{2}$,
$\bar{p}_{\xi}^{2}=[\langle
p_{\xi}\rangle/(l_{\xi}\nu_{\xi}M)]^{2}=(l_{\xi}/\hbar)^{2}\langle
p_{\xi}^{2}\rangle$, and
$\bar{c}_{\xi}=c_{\xi}/(l_{\xi}^{2}\nu_{\xi}M)=c_{\xi}/\hbar$, where we
introduced a bar for indicating the dimensionless quantities. Given this,
equations (V) in rescaled units read

$\displaystyle\frac{d}{d\tau}\bar{r}_{\xi}^{2}=\bar{c}_{\xi}$
$\displaystyle\frac{d}{d\tau}\bar{p}_{\xi}^{2}=\left\\{\bar{\Gamma}_{\xi}\Im[\bar{\Phi}_{\xi}(\tau)]-\bar{{W}}_{\xi}^{\prime}(\tau)\right\\}\bar{c}_{\xi}-2\bar{\Gamma}_{\xi}\Im[\bar{\Lambda}_{\xi}(\tau)]\bar{p}_{\xi}^{2}$
$\displaystyle+\bar{\Gamma}_{\xi}\Re[\bar{\Phi}_{\xi}(\tau)]$
$\displaystyle\frac{d}{d\tau}\bar{c}_{\xi}=2\left\\{\bar{\Gamma}_{\xi}\Im[\bar{\Phi}_{\xi}(\tau)]-\bar{{W}}_{\xi}^{\prime}(\tau)\right\\}\bar{r}_{\xi}^{2}+2\bar{p}_{\xi}^{2}$
$\displaystyle-\bar{\Gamma}_{\xi}\left\\{\Im[\bar{\Lambda}_{\xi}(\tau)]\bar{c}_{\xi}-\Re[\bar{\Lambda}_{\xi}(\tau)]\right\\}.$
(157)

with $\tau=\nu_{\xi}t$, $\bar{\Omega}_{rf}^{\xi}=\Omega_{rf}/\nu_{\xi}$, and

$\displaystyle\bar{\Gamma}_{\xi}$
$\displaystyle=\frac{2\pi}{3}\frac{m}{M}\left(\frac{M}{\mu}\right)^{2}(n_{0}l_{\xi}^{3}),$
$\displaystyle\bar{{W}}_{\xi}^{\prime}(\tau)$
$\displaystyle=\frac{\bar{{W}}_{\xi}(\tau)+\delta\bar{{W}}_{\xi}(\tau)}{\nu_{\xi}^{2}}=\left(\frac{\bar{\Omega}_{rf}^{\xi}}{2}\right)^{2}$
$\displaystyle\times\left[a_{\xi}+\delta\bar{a}_{\xi}-2(q_{\xi}+\delta\bar{q}_{\xi})\cos(\bar{\Omega}_{rf}^{\xi}\tau)-2\delta\bar{q}_{\xi}^{\prime}\bar{g}_{\xi}(\tau)\right]$
$\displaystyle\bar{g}_{\xi}(\tau)$
$\displaystyle=\\!\\!\\!\sum_{s,s^{\prime}\notin
S_{\mathrm{i}}}\\!\\!\\!F_{s,s^{\prime}}^{\xi}\cos[(s-s^{\prime})\bar{\Omega}_{rf}^{\xi}\tau]\\!\\!\left[\bar{{\mathcal{J}}}_{+}^{\prime}(\bar{q}_{s,\xi})-\bar{{\mathcal{J}}}_{-}^{\prime}(\bar{q}_{s,\xi})\right]$
$\displaystyle\bar{Q}_{\xi}$
$\displaystyle=\frac{Q}{l_{\xi}}=\frac{32}{3}\frac{m}{M}\left(\frac{M}{\mu}\right)^{2}\frac{(n_{0}l_{\xi}^{3})}{(\bar{\Omega}_{rf}^{\xi})^{2}}$
(158)

and
$\bar{{\mathcal{J}}}_{\pm}^{\prime}(\bar{q}_{s,\xi})={\mathcal{J}}_{\pm}^{\prime}(q_{s,\xi})l_{\xi}$.
Here, for the sake of simple notation, we just refer to the bosonic case.
Moreover, the rescaled $\Phi$ and $\Lambda$ functions read:

$\begin{split}\bar{\Lambda}_{\xi}(\tau)=i&\sum_{s,s^{\prime}}C^{\xi}_{s}C^{\xi}_{s^{\prime}}\bigg{\\{}|\bar{q}_{s,\xi}|^{3}\bar{f}(\bar{q}_{s,\xi})^{2}(1-\phi)\\\\[4.30554pt]
&\bigg{[}i\sin\big{[}(s-s^{\prime})\bar{\Omega}_{rf}\tau\big{]}(1+2n_{\bar{q}_{s,\xi}})\\\\[4.30554pt]
&+\cos\big{[}(s-s^{\prime})\bar{\Omega}_{rf}\tau\big{]}\mathrm{sgn}\big{(}\beta_{\xi}/2+s\big{)}\bigg{]}\\\\[4.30554pt]
&-\bar{\eta}^{+}_{s,\xi}e^{-i(s-s^{\prime})\bar{\Omega}_{rf}\tau}+\bar{\eta}^{-}_{s,\xi}e^{i(s-s^{\prime})\bar{\Omega}_{rf}\tau}\bigg{\\}}\\\\[4.30554pt]
-\frac{2}{\pi}\bar{\Omega}_{rf}&\sum_{s,s^{\prime}}C^{\xi}_{s}C^{\xi}_{s^{\prime}}\bigg{\\{}\cos\big{[}(s-s^{\prime})\bar{\Omega}_{rf}\tau\big{]}\\\\[4.30554pt]
\Big{[}1-&2\Theta\big{(}\mathrm{sgn}(\beta_{\xi}/2+s)\big{)}\Big{]}\Big{[}\bar{\mathcal{J}}^{\prime}_{-}(\bar{q}_{s,\xi})+\bar{\mathcal{J}}^{\prime}_{+}(\bar{q}_{s,\xi})\\\\[4.30554pt]
&+2\big{(}\bar{\mathcal{J}}^{n_{q}}_{-}(\bar{q}_{s,\xi})+\bar{\mathcal{J}}^{n_{q}}_{+}(\bar{q}_{s,\xi})\Big{]}\\\\[4.30554pt]
-i\sin&\big{[}(s-s^{\prime})\bar{\Omega}_{rf}\tau\big{]}\Big{[}\bar{\mathcal{J}}^{\prime}_{-}(\bar{q}_{s,\xi})-\bar{\mathcal{J}}^{\prime}_{+}(\bar{q}_{s,\xi})\Big{]}\bigg{\\}}\end{split}$
(159)

$\begin{split}\bar{\Phi}_{\xi}(\tau)=\,\bar{\Omega}_{rf}&\sum_{s,s^{\prime}}C^{\xi}_{s}C^{\xi}_{s^{\prime}}(\beta_{\xi}/2+s^{\prime})\bigg{\\{}|\bar{q}_{s,\xi}|^{3}\bar{f}(\bar{q}_{s,\xi})^{2}\\\\[4.30554pt]
&(1-\phi)\bigg{[}\cos\big{[}(s-s^{\prime})\bar{\Omega}_{rf}\tau\big{]}(1+2n_{\bar{q}_{s,\xi}})\\\\[4.30554pt]
&+i\sin\big{[}(s-s^{\prime})\bar{\Omega}_{rf}\tau\big{]}\mathrm{sgn}\big{(}\beta_{\xi}/2+s\big{)}\bigg{]}\\\\[4.30554pt]
&+\bar{\eta}^{-}_{s,\xi}e^{-i(s-s^{\prime})\bar{\Omega}_{rf}\tau}+\bar{\eta}^{+}_{s,\xi}e^{i(s-s^{\prime})\bar{\Omega}_{rf}\tau}\bigg{\\}}\\\\[8.61108pt]