# Bose polaron as an instance of quantum Brownian motion

Aniello Lampo<EMAIL_ADDRESS>ICFO – Institut de Ciencies Fotoniques,
The Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain  Soon Hoe Lim  Department of Mathematics (Program in
Applied Mathematics), University of Arizona, Tucson, AZ 85721-0089, USA
Miguel Ángel García-March  ICFO – Institut de Ciencies Fotoniques, The
Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain  Maciej Lewenstein  ICFO – Institut de Ciencies Fotoniques,
The Barcelona Institute of Science and Technology, 08860 Castelldefels
(Barcelona), Spain  ICREA – Institució Catalana de Recerca i Estudis Avançats,
Lluis Companys 23, E-08010 Barcelona, Spain

###### Abstract

We study the dynamics of a quantum impurity immersed in a Bose-Einstein
condensate as an open quantum system in the framework of the quantum Brownian
motion model. We derive a generalized Langevin equation for the position of
the impurity. The Langevin equation is an integro-differential equation that
contains a memory kernel and is driven by a colored noise. These result from
considering the environment as given by the degrees of freedom of the quantum
gas, and thus depend on its parameters, e.g. interaction strength between the
bosons, temperature, etc. We study the role of the memory on the dynamics of
the impurity. When the impurity is untrapped, we find that it exhibits a
super-diffusive behavior at long times. We find that back-flow in energy
between the environment and the impurity occurs during evolution. When the
particle is trapped, we calculate the variance of the position and momentum to
determine how they compare with the Heisenberg limit. One important result of
this paper is that we find position squeezing for the trapped impurity at long
times. We determine the regime of validity of our model and the parameters in
which these effects can be observed in realistic experiments.

##  1 Introduction

The concept of polaron has been introduced by Landau and Pekar to describe the
behavior of an electron in a dielectric crystal [1]. The motion of the
electron distorts the spatial configuration of the surrounding ions, which let
their equilibrium positions to screen its charge. The movement of the ions is
associated to phonon excitations that dress the electron. The resulting
system, which consists of the electron and its surrounding phonon cloud, is
called a polaron. The concept of polaron has been extended to describe a
generic particle, the impurity, in a generic material, e.g. a conductor, a
semiconductor or a gas [2, 3]. One important example is that of an impurity
embedded in an ultracold gas. This system has been widely studied both
theoretically and experimentally, in the case of a ultracold Fermi [4, 5, 6,
7, 8, 9, 10] or Bose gas [11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36]. We focus on the latter
case, where the Bose gas hosting the impurity is a Bose-Einstein condensate
(BEC).

In this paper, we study the physics of the impurity as an open system in the
framework of quantum Brownian motion (QBM) model. In general, the QBM
describes a Brownian particle moving in a thermal bath consisting of a
collection of non-interacting harmonic oscillators satisfying the Bose-
Einstein statistics [37, 38, 39, 40, 41, 42, 43]. In our context, the impurity
plays the role of the Brownian particle, while the bath consists of the
Bogoliubov modes of the BEC. We assume that the bath is homogeneous, i.e. the
density of the BEC is space-independent.

The open quantum system approach has been used recently in the context of
ultracold quantum gases. For instance, the system of a bright soliton in a
superfluid in one dimension in [44], the system of a dark soliton in a one-
dimensional BEC coupled to a non-interacting Fermi gas in [45], the system of
the component of a moving superfluid [46], and the system of an impurity in a
Luttinger liquid in [47, 48].

There are a few formalisms within the open system approach that one can adopt
to study the dynamics of the system. One of them is the master equation which
describes the evolution of the reduced density matrix of the system. A recent
study of the master equation for QBM has been performed in [49], both for
linear and nonlinear coupling. There the master equation was derived using the
Born-Markov approximation. Therefore, one possible way to study the dynamics
of the impurity in both homogeneous and inhomogeneous BEC is to use such
master equations. However, the Born-Markov approximation is only valid at high
temperature and for sufficiently weak coupling between the environment and the
system. In fact, as we go beyond these restrictions we detect a breakdown of
the positivity of the solution, corresponding to a violation of the Heisenberg
uncertainty principle in these regimes. One remedy to fix these problems is to
cast the master equation in a Lindblad form by adding terms that are
negligible in these regimes to the Born-Markov master equation. In [50] a
Lindblad model for QBM has been studied, for both linear and nonlinear
coupling. This allows one to explore the low temperature regime, as well as to
consider values of the coupling between the environment and the system
stronger than those permitted by the Born-Markov treatment. The drawback to
this remedy is that the resulting Lindblad equation cannot be derived directly
from the microscopic Hamiltonian of the system. Therefore, it cannot be
adopted to investigate the dynamics of the impurity.

In view of the above remarks, we employ an alternative formalism to study our
system: the Heisenberg equations for observables. This formalism not only
permits us to go beyond the approximations underlying the works in [49, 50],
but is also simpler compared to the master equation formalism, since deriving
an exact master equation for QBM is more difficult analytically. Moreover, as
the main goal of the paper is to study the evolution of the observables of the
impurity, Heisenberg formalism is particularly well suited for this task.

We first show that the Hamiltonian of an impurity in an ultracold gas can be
cast as a Hamiltonian of QBM model. We find that the coupling between the
system and the environment in the resulting Hamiltonian is nonlinear. A key
assumption in this work and previous approaches [44, 47] is to approximate
this coupling by a linear one. Here we present a detailed study of the regimes
of parameters in which this assumption is valid, and evaluate them in view of
current experimental feasibility for this system. In particular, we find the
regimes of temperatures, interaction strength, experimental time and trapping
frequency for the impurity in which the model is valid. We derive the spectral
density (SD) of the system, which encapsulates the effect of the bath on the
central system. The spectral density is a super-ohmic one, corresponding to
the presence of memory effects in the dynamics of the system. We then derive
quantum Langevin equations to describe the evolution of relevant observables
associated to the impurity. We solve these equations analytically whenever
possible and numerically for both configurations of the impurity. In the
configuration where the impurity is untrapped, we calculate the mean square
displacement, a measurable quantity in experiments [15]. We find that the
particle exhibits a super-diffusive behavior, with a quadratic in time mean
square displacement at long times. We show that a back-flow of energy between
the system and the environment occurs during evolution. When the impurity is
trapped, we calculate its position and momentum variance, analyzing its
dependence on system parameters such as temperature and coupling strength. We
detect new effects concerning the structure of the uncertainties ellipse:
genuine position squeezing occurs in the system, corresponding to high spatial
localization of the impurity. This is a very relevant result of this paper. We
find the regime of experimentally feasible parameters to optimize such
squeezing.

The manuscript is organized as follows. In Sec. 2 we present the Hamiltonian
of an impurity in an ultracold gas, discuss the steps to derive a QBM model
for the system and justify the linear approximation for the coupling. In Sec.
3 we derive the SD of the system. In Sec. 4 we derive the quantum Langevin
equations, and solve them in Secs. 5 and 6 for the untrapped and trapped
configuration, respectively. In the Appendices, we provide the details of the
derivations, analysis of the feasibility of the approach and of the presented
results in the parameter regimes.

##  2 Hamiltonian

We consider a single impurity atom with mass $m_{\rm I}$ immersed in a
$d-$dimensional gas of $N$ bosons with mass $m_{\rm B}$. The interactions
among the bosons are described by a potential $V_{\rm B}(\textbf{x})$. Let
$\Psi(\textbf{x})$, $\Psi^{\dagger}(\textbf{x})$ denote the annihilation and
creation operators of atoms at the position x. They fulfill the canonical
bosonic commutation relation,
$[\Psi(\textbf{x}),\Psi^{\dagger}(\textbf{x'})]=\delta(\textbf{x}-\textbf{x'})$.
Bosonic density thus reads $n_{\rm
B}=\Psi^{\dagger}(\textbf{x})\Psi(\textbf{x})$. The Hamiltonian of the system
is

$H=H_{\rm I}+H_{\rm B}+H_{\rm BB}+H_{\rm IB},$ (1)

where:

$\displaystyle H_{\rm I}=\frac{{\textbf{p}}^{2}}{2m_{\rm I}}+U(\textbf{x}),$
(2) $\displaystyle H_{\rm B}=\int
d^{d}\textbf{x}\,\Psi^{{\dagger}}(\textbf{x})\left(\frac{\textbf{p}_{\rm
B}^{2}}{2m_{\rm B}}+V(\textbf{x})\right)\Psi(\textbf{x})$ (3)
$\displaystyle=\sum_{\textbf{k}}\epsilon_{\rm\textbf{k}}{a}^{\dagger}_{\textbf{k}}{a}_{\textbf{k}},$
(4) $\displaystyle H_{\rm BB}=g_{\rm B}\int
d^{d}\textbf{x}\,\Psi^{\dagger}(\textbf{x})\Psi^{\dagger}(\textbf{x})\Psi(\textbf{x})\Psi(\textbf{x})$
(5) $\displaystyle=\frac{1}{2V}\sum_{\textbf{k},\textbf{k'},\textbf{q}}V_{\rm
B}(\textbf{q}){a}^{\dagger}_{\textbf{k'}-\textbf{q}}{a}^{\dagger}_{\textbf{k}+\textbf{q}}{a}_{\textbf{k'}}{a}_{\textbf{k}},$
(6) $\displaystyle H_{\rm IB}=g_{\rm IB}n_{\rm B}(\textbf{x})$ (7)
$\displaystyle=\frac{1}{V}\sum_{\textbf{k},\textbf{q}}V_{\rm
IB}(\textbf{k}){\rho}_{\rm
I}(\textbf{q}){a}^{\dagger}_{\textbf{k}-\textbf{q}}{a}_{\textbf{k}},$ (8)

in which the operator $a_{\textbf{k}}$ ($a^{\dagger}_{\textbf{k}}$) destroys
(creates) a boson of mass $m_{\rm B}$, wave vector k, and energy
$\epsilon_{\textbf{k}}=(\hbar k)^{2}/(2m_{\rm B})-\mu$, measured from its
chemical potential, $\mu$. Equation (2) is the Hamiltonian for the free
impurity. Equation (4) is the Hamiltonian of non-interacting bosons in a
potential $V(\textbf{x})$. We consider in the following that the potential is
homogeneous and the system is enclosed in a box of volume $V$. The last two
terms in Eqs. (6) and (8) are the interaction among the atoms of the gas and
between them and the impurity, respectively. The quantities

$\displaystyle V_{\rm B}(\textbf{q})=\mathcal{F}_{\textbf{q}}\left[g_{\rm
B}\delta(\textbf{x}-\textbf{x'})\right],$ (9) $\displaystyle V_{\rm
IB}(\textbf{q})=\mathcal{F}_{\textbf{q}}\left[g_{\rm
IB}\delta(\textbf{x}-\textbf{x'})\right],$ (10)

represent respectively the Fourier transforms of the boson-boson and impurity-
boson interaction, with

$g_{\rm B}=4\pi\hbar^{2}a_{\rm B}/m_{\rm B},\quad g_{\rm
IB}=2\pi\hbar^{2}a_{\rm IB}/m_{\rm R}.$ (11)

Here $a_{\rm B}$ is the scattering length between two identical bosons and
$a_{\rm IB}$ represents that between the impurity and the BEC bosons. The
reduced mass is $m_{\rm R}=m_{\rm B}m_{\rm I}/(m_{\rm B}+m_{\rm I})$ and the
(dimensionless) density of the impurity in the momentum domain is

${\rho}_{\rm
I}(\textbf{q})=\int^{+\infty}_{-\infty}e^{-i\textbf{q}\cdot\tilde{\textbf{x}}}\delta\left(\tilde{\textbf{x}}-\textbf{x}\right)d^{3}\tilde{\textbf{x}}.$
(12)

In experiments, we usually have more than one impurity in the gas, so one can
include a term modeling the interaction between several impurities in Eq.
$\left(\ref{Ham1}\right)$. Here we consider that the impurities concentration
is low enough as to neglect such an additional interacting term. It is
important to realize that the Hamiltonian is positively defined, and as such
cannot lead to instabilities – this is clearly seen from the form of the
various parts of the Hamiltonian in the position representation.

To obtain the QBM form of the Hamiltonian we first replace the
creation/annihilation operator in the fundamental state by its average value
$\sqrt{N_{0}}$. For bosons below the critical temperature the atoms mainly
occupy the ground mode, with negligible fluctuations to other modes, thus
forming a BEC. Consequently, we neglect the terms proportional to
$N_{\textbf{k}}$ $(\textbf{k}\neq\textbf{0})$, i.e. the number of particles
out of the ground state. Then, we apply the Bogoliubov transformation

${a}_{\textbf{k}}=u_{\textbf{k}}{b}_{\textbf{k}}-v_{\textbf{k}}{b}_{-\textbf{k}}^{\dagger},\quad{a}_{-\textbf{k}}=u_{\textbf{k}}{b}_{-\textbf{k}}-v_{\textbf{k}}{b}_{\textbf{k}}^{\dagger},$
(13)

with

$\displaystyle
u^{2}_{\textbf{k}}=\frac{1}{2}\left(\frac{\epsilon_{\textbf{k}}+n_{0}V_{\rm
B}}{E_{\textbf{k}}}+1\right),$ (14) $\displaystyle
v^{2}_{\textbf{k}}=\frac{1}{2}\left(\frac{\epsilon_{\textbf{k}}+n_{0}V_{\rm
B}}{E_{\textbf{k}}}-1\right),$ (15)

in which

$E_{\textbf{k}}=\hbar
c|\textbf{k}|\sqrt{1+\frac{1}{2}\left(\xi\textbf{k}\right)^{2}}\equiv\hbar\omega_{\textbf{k}}$
(16)

is the Bogoliubov spectrum and $n_{\rm 0}$ is the density of particles in the
ground state. Since we consider a homogeneous gas, the density $n_{\rm 0}$ is
constant. In Eq. (16) the quantities

$\xi=\frac{\hbar}{\sqrt{2g_{\rm B}m_{\rm B}n_{\rm 0}}},\quad
c=\sqrt{\frac{g_{\rm B}n_{\rm 0}}{m_{\rm B}}}=\frac{\hbar}{\sqrt{2}m_{\rm
B}\xi}$ (17)

represent respectively the coherence length and the speed of sound.

The transformations (13) diagonalize the terms describing the condensed atoms
(see e.g. [51])

${H}_{\rm B}+H_{\rm
BB}=\sum_{\textbf{k}\neq\textbf{0}}E_{\textbf{k}}{b}^{\dagger}_{\textbf{k}}{b}_{\textbf{k}},$
(18)

apart of a few non-operator terms, simply providing a shift of the energy
levels of the atoms of the BEC.

We treat in the same manner the bosons-impurity interaction. We are going now
to keep only the terms proportional to $\sqrt{N_{0}}$ – this is in principle a
well motivated approximation, since the condensate is macroscopically occupied
and $N_{i\neq 0}\ll N_{0}$. Unfortunately, this approximation is dangerous.
Our model generically has an ultraviolet divergence, like most of models of
the non-relativistic quantum field theory. We do all the calculations with a
physical cut-off, so that the ultraviolet divergences do not really affect our
theory. Still, at large values of the cut-off we might expect large negative
shifts of the impurity energy, that might cause unphysical instability. After
dropping out the terms bilinear in $N_{i\neq 0}$ we obtain

${H}_{\rm IB}=n_{0}V_{\rm IB}+\sqrt{\frac{n_{\rm
0}}{V}}\sum_{\textbf{k}\neq\textbf{0}}{\rho}_{\rm I}(\textbf{k})V_{\rm
IB}\left({a}_{\textbf{k}}+{a}^{\dagger}_{-\textbf{k}}\right).$ (19)

The first term on the right hand-side represents the mean field energy. It is
a constant just providing a shift of the energy of the polaron, so we will
neglect it in what follows. Inserting Eq. (13) into Eq. (19), one gets

$\displaystyle{H}_{\rm IB}=$ $\displaystyle\sqrt{\frac{n_{\rm
0}}{V}}\sum_{\textbf{k}\neq\textbf{0}}{\rho}_{\rm I}(\textbf{k})V_{\rm
IB}\left(u_{\textbf{k}}-v_{\textbf{k}}\right)\left({b}_{\textbf{k}}+{b}^{\dagger}_{-\textbf{k}}\right)$
(20) $\displaystyle=$ $\displaystyle\sqrt{\frac{n_{\rm
0}}{V}}\sum_{\textbf{k}\neq\textbf{0}}{\rho}_{\rm I}(\textbf{k})V_{\rm
IB}\sqrt{\frac{\epsilon_{\textbf{k}}}{E_{\textbf{k}}}}\left({b}_{\textbf{k}}+{b}^{\dagger}_{-\textbf{k}}\right),$

where we again neglected the non-operator terms.
By some algebra, the expression in Eq. (20) reads

${H}_{\rm
IB}=\sum_{\textbf{k}\neq\textbf{0}}V_{\textbf{k}}e^{i\textbf{k}\cdot{\textbf{x}}}\left({b}_{\textbf{k}}+{b}^{\dagger}_{-\textbf{k}}\right),$
(21)

in which

$V_{\textbf{k}}=g_{\rm IB}\sqrt{\frac{n_{\rm 0}}{V}}\left[\frac{(\xi
k)^{2}}{(\xi k)^{2}+2}\right]^{\frac{1}{4}}.$ (22)

The expression in Eq. (21) is known in the literature as the interaction term
of the Fröhlich Hamiltonian. We restrict to the limit

$\textbf{k}\cdot\textbf{x}\ll 1.$ (23)

Accordingly Eq. (21) assumes the form:

${H}_{\rm
IB}=\sum_{\textbf{k}\neq\textbf{0}}V_{\textbf{k}}\left(\mathbb{I}+i\textbf{k}\cdot\textbf{x}\right)\left({b}_{\textbf{k}}+{b}^{\dagger}_{-\textbf{k}}\right).$
(24)

Equation $\left(\ref{mainAssumption}\right)$ is a crucial assumption of this
paper. It allows for a linear coupling between the impurities and the bosons,
which will play the role of an environment. This is crucial when considering
this system from an open quantum system perspective. In Appendix C we justify
when this assumption is valid as a function of the relevant physical
parameters in the problem. A similar assumption is considered in [44] to study
the physics of a bright soliton in a superfluid, and in [47] where QBM has
been employed to treat the dynamics of an impurity in a Luttinger Liquid.

The resulting Hamiltonian of the impurity in a BEC is

${H}=H_{\rm
I}+\sum_{\textbf{k}\neq\textbf{0}}E_{\textbf{k}}{b}^{\dagger}_{\textbf{k}}{b}_{\textbf{k}}+\sum_{\textbf{k}\neq\textbf{0}}\hbar
g_{\textbf{k}}\pi_{\textbf{k}}\textbf{x},$ (25)

with

$g_{\textbf{k}}=\textbf{k}V_{\textbf{k}}/\hbar,\quad{\pi}_{\textbf{k}}=i\left({b}_{\textbf{k}}-{b}^{\dagger}_{\textbf{k}}\right).$
(26)

To get Eq. (25) we redefined the Bogoliubov modes operators as
${b}_{\textbf{k}}\rightarrow{b}_{\textbf{k}}-V_{\textbf{k}}/E_{\textbf{k}}\mathbb{I},$
in order to get rid of the term in Eq. (21) proportional to the identity
operator. This operation yields a non-operator term which has been neglected
in agreement with the procedure above.

The Hamiltonian in Eq. (25) describes an impurity coupled to a bath of
Bogoliubov modes through an interaction term linearly dependent on the
impurity position. This is exactly the same situation of the QBM. Here the
impurity plays the role of the Brownian particle, while the Bogoliubov modes
represents the environment. The Hamiltonian in Eq. (25) is in fact almost the
same of that of the QBM model. The only difference lies in the dependence of
the interaction term on the Bogoliubov modes operator: while in the QBM model
it depends on their positions, for the present system it depends on their
(dimensionless) momenta, $\pi_{\textbf{k}}$. We will see in the following that
the two situations are equivalent, and the theory of the QBM can be exploited
to investigate the impurity problem.

Unfortunately, the Hamiltonian in Eq. (25) is not positively defined (note
that at this point we match the situation described in [52], where
translational symmetry is broken). We could repair this by taking into account
bilinear terms in ${a}_{\textbf{k}}$ and ${a}^{\dagger}_{\textbf{k}}$ in Eq.
(19), or equivalently ${b}_{\textbf{k}}$ and ${b}^{\dagger}_{\textbf{k}}$, as
in [53, 17, 54, 29, 30]. One way of curing it is to include these terms in the
theory in an exact manner, which is possible, but technically complex. It is
much easier to use the same trick as Caldeira and Leggett used in their
seminal paper – i.e. complete the Hamiltonian to a positively defined one by
writing

${H}=H_{\rm
I}+\sum_{\textbf{k}\neq\textbf{0}}E_{\textbf{k}}\left({b}^{\dagger}_{\textbf{k}}+i\hbar
g_{\textbf{k}}\textbf{x}/E_{\textbf{k}}\right)\left({b}_{\textbf{k}}-i\hbar
g_{\textbf{k}}\textbf{x}/E_{\textbf{k}}\right).$ (27)

Clearly, Caldeira-Leggett remedy leads directly to the trapping harmonic
potential for the impurity that cancels the negative harmonic frequency shift
that appears in the absence of the compensation term. In Appendix A it is
shown how this negative contribution arises by evaluating its effect directly
in the equations of motion derived from Hamiltonian (25). The equations of
motion are presented in Sec. 4.

The Hamiltonian in Eq. (25) is the first step to put the Bose polaron problem
in the framework of quantum Brownian motion. This result has been obtained by
considering the Hamiltonian in Eq. (1), which is a conventional choice in the
context of polaron physics and it is largely used in the literature, for
example in [55, 19, 56]. Nevertheless, this Hamiltonian is not fully
appropriate if we push our analysis towards the strong coupling regime. Here,
the interaction term in Eq. (8) needs to be generalized by including a
quadratic dependence on the operators of the impurity. Accordingly, when one
introduces Bogoliubov operators, the interaction term in Eq. (21) includes
additional terms manifesting a quadratic, rather than linear, dependence on
these operators. In this way one goes beyond the Fröhlich paradigm. Such a
generalization is widely studied nowadays, for instance in [53, 17, 54, 29,
30]. In particular there is still an open debate concerning the validity
regime of the Fröhlich Hamiltonian, i.e. for which values of the system
parameters the quadratic Bogoliubov operators terms can be dropped out. In the
present paper we neglect these terms, looking to the traditional Fröhlich
interaction term in Eq. (21). Of course this choice does not allow us to
explore the strong coupling regime, but we shall explain in the following that
it is appropriate for the values of the system parameters we consider.

##  3 Spectral Density

The last term of Hamiltonian (25) is a coupling between the impurity position
and the momenta of the bath oscillators, $\pi_{\textbf{k}}$. This plays the
role of the interaction Hamiltonian between the system and the bath
oscillators, when compared with the QBM Hamiltonian. The difference is that
the coupling occurs with the momenta of the oscillators rather than on their
positions. In the following, we show that both situations are equivalent. In
fact, the interaction with the Bogoliubov modes enters in the dynamics of the
system only through the self-correlation function of the environment, defined
as

$\mathcal{C}(\tau)=\sum_{\textbf{k}\neq\textbf{0}}\hbar
g^{2}_{\textbf{k}}\langle{\pi}_{\textbf{k}}(\tau){\pi}_{\textbf{k}}(0)\rangle.$
(28)

Replacing the expression for the dimensionless momenta $\pi_{\textbf{k}}$ in
Eq. (26) and recalling that the environment is made up by bosons,

$\langle{b}^{\dagger}_{\textbf{k}}{b}_{\textbf{k}}\rangle=\frac{1}{\exp\left(\hbar\omega_{\textbf{k}}/k_{\rm
B}T\right)-1},$ (29)

we find

$\displaystyle\mathcal{C}(\tau)=$
$\displaystyle\sum_{\textbf{k}\neq\textbf{0}}\hbar
g^{2}_{\textbf{k}}\left[\coth\left(\frac{\hbar\omega_{\textbf{k}}}{2k_{\rm
B}T}\right)\cos\left(\omega_{\textbf{k}}\tau\right)-i\sin\left(\omega_{\textbf{k}}\tau\right)\right]$
(30) $\displaystyle\equiv$ $\displaystyle\nu(\tau)-i\lambda(\tau),$ (31)

where

$\displaystyle\nu(\tau)=\int^{\infty}_{0}J(\omega)\coth\left(\frac{\hbar\omega}{2k_{\rm
B}T}\right)\cos\left(\omega\tau\right)d\omega,$ (32)
$\displaystyle\lambda(\tau)=\int^{\infty}_{0}J(\omega)\sin\left(\omega\tau\right)d\omega=-m_{\rm
I}\dot{\Gamma}(t),$ (33)

are respectively the noise and dissipation kernel, representing the real and
imaginary part of the self-correlation function. The latter is related to the
damping kernel, which is defined as

$\Gamma(\tau)=\frac{1}{m_{\rm
I}}\int^{\infty}_{0}\frac{J(\omega)}{\omega}\cos(\omega\tau)d\omega.$ (34)

A crucial function in the expressions above is the spectral density  (SD)

$J(\omega)=\sum_{\textbf{k}\neq\textbf{0}}\hbar
g^{2}_{\textbf{k}}\delta\left(\omega-\omega_{\textbf{k}}\right).$ (35)

The SD plays a fundamental role on characterizing the dynamics of QBM. It
contains the information about the environment after averaging over its
degrees of freedom. Note that this quantity depends on the square of the
coupling constant. This is the reason why, as we will see in the following,
our theory does not depend on the sign of the interaction.

The rest of this section is devoted to the derivation of the SD. We first
assume that the environment is large, that is, the large number of oscillators
within it allows to switch from a discrete to a continuous distribution of
Bogoliubov modes in the frequency domain. Then, in the definition of the SD,
Eq. (35), we turn the discrete sum into an integral,

$\sum_{\textbf{k}}\rightarrow\int\frac{V}{(2\pi)^{\rm d}}d^{\rm d}k.$ (36)

Using the relation

$\delta\left(\omega-\omega_{\textbf{k}}\right)=\frac{1}{\partial_{\textbf{k}}\omega_{\textbf{k}}|_{\textbf{k}=\textbf{k}_{\omega}}}\delta\left(\textbf{k}-\textbf{k}_{\omega}\right),$
(37)

one finds

$J(\omega)\\!=\\!\frac{n_{\rm 0}g^{2}_{\rm IB}}{\hbar}\frac{S_{\rm
d}}{(2\pi)^{\text{d}}}\int dkk^{\rm
d+1}\sqrt{\frac{(\xi\textbf{k})^{2}}{(\xi\textbf{k})^{2}+2}}\frac{\delta\left(\textbf{k}-\textbf{k}_{\omega}\right)}{\partial_{\textbf{k}}\omega_{\textbf{k}}|_{\textbf{k}=\textbf{k}_{\omega}}},$
(38)

where

$k_{\omega}=\xi^{-1}\sqrt{\sqrt{1+2\left(\xi\omega/c\right)^{2}}-1},$ (39)

is the inverse of the Bogoliubov spectrum in Eq. (16). The quantity $S_{\rm
d}$ is the surface of the hypersphere in the momentum space with radius $k$ in
$\rm{d}-$dimensions. For $\rm{d}=1,2,3$ it reads

$S_{1}=2,\quad S_{2}=2\pi,\quad S_{3}=4\pi^{2}.$ (40)

Hereafter, we focus on $\text{d}=1$, but the generalization to higher
dimensions is conceptually immediate. Thus, in one dimension the SD is

$\displaystyle J_{\rm 1d}(\omega)=m_{\rm I}\tilde{\tau}\omega^{3}\chi_{\rm
1d}(\omega),$ (41) $\displaystyle\tilde{\tau}=\frac{\eta^{2}}{2\pi m_{\rm
I}}\left(\frac{m_{\rm B}}{n_{\rm 0}g^{\frac{1}{3}}_{\rm
B}}\right)^{\frac{3}{2}},$ (42) $\displaystyle\chi_{\rm
1d}(\omega)=2\sqrt{2}\left(\frac{\Lambda}{\omega}\right)^{3}\frac{\left[\sqrt{1+\frac{\omega^{2}}{\Lambda^{2}}}-1\right]^{\frac{3}{2}}}{\sqrt{1+\frac{\omega^{2}}{\Lambda^{2}}}},$
(43)

where $\tilde{\tau}$ represents a relaxation time scale. We introduced the
interaction strength

$\eta=g_{\rm IB}/g_{\rm B},$ (44)

which expresses the strength of the impurity-bosons interaction in units of
the bosons-bosons one. The majority of our results are expressed as a function
of such a parameter because in experiments with ultracold gases it can be
tuned [15]. This parameter is the crucial one to define the regime of validity
of our theory. In fact, in Sec. 2 we precised that the form of the Hamiltonian
we use, showing a linear, rather than quadratic, dependence on the Bose
operators, does not work for strong coupling. The parameter in Eq. (44) allows
to quantify to which extent the coupling has to be weak. In particular in [27,
57] it has been shown that the standard Fröhlich Hamiltonian in Eq. (21)
without quadratic terms in the Bose operators holds if

$\eta\lesssim\eta_{c}\equiv\pi\sqrt{n_{\rm 0}\left|a_{\rm
B}\right|}=\pi\sqrt{\frac{2n_{\rm 0}}{m_{\rm B}g_{\rm B}}}.$ (45)

This equation has been derived for a gas in one dimension, but in [27] the
corresponding $3\text{D}$ generalization is presented. We consider the same
situation of [15], i.e. a gas of Rb atoms with

$n_{\rm 0}=7\left(\mu\text{m}\right)^{-1},\quad g_{\rm B}=2.36\cdot
10^{-37}\text{J}\cdot\text{m}.$ (46)

Accordingly we obtain

$\eta_{c}\approx 7,$ (47)

providing an upper bound for the acceptable values for the coupling strength.

The quantity

$\Lambda=n_{0}g_{\rm B}/\hbar$ (48)

is a characteristic frequency distinguishing the high-frequency domain from
the low-frequency one. The SD is ohmic when it depends linearly on the
frequency of the oscillators of the environment. This is not the case of the
physical system we are dealing with. For $\omega\ll\Lambda$ we find

$J_{1\rm d}\sim\omega^{3},$ (49)

namely SD is proportional to the third power of the frequency of the
Bogoliubov modes. In our case, the SD is super-ohmic. In higher dimensions,
the SD is also super-ohmic, as in general

$J_{\rm d}\sim\omega^{2+\text{d}}.$ (50)

The super-ohmic dependence in Eq. (49) has been found in [47, 58] for an
impurity immersed in a Luttinger liquid, and in [44] for a bright soliton in a
superfluid. This behavior is associated to the the linear part of the
Bogoliubov spectrum. Hereafter we shall focus on the former considering

$J(\omega)=m_{\rm I}\tilde{\tau}\omega^{3}\theta(\omega-\Lambda),$ (51)

namely we introduce a sharp ultraviolet cut-off to regularize the SD at high-
frequency. We emphasize that the results we will present are independent on
this ultraviolet regularization. Indeed, we reproduced the calculation by
introducing an exponential, rather than a sharp cut-off, and we obtained the
same results (see Appendix B).

However, apart of its analytical form, the cut-off plays a crucial role. The
results we shall present in the following depend on its presence, because it
allows us to focus on the low frequency portion of the spectral density,
forgetting about the high frequency one. Such a way to proceed is absolutely
fitting if one is interested in the long-time dynamics. In this context the
cut-off frequency in Eq. (48) is not an artificial quantity, but it arises in
a natural manner: it is the characteristic frequency distinguishing the
phononic linear part of the Bogoliubov spectrum from the quadratic one. This
is more clear if we express it in terms of the traditional parameters of the
Bogoliubov spectrum:

$\Lambda=c/\left(\sqrt{2}\xi\right)$ (52)

Precisely, considering $\omega/\Lambda\ll 1$, the inverse of the Bogoliubov
spectrum in Eq. (39) takes the following form:

$k_{\omega}\sim\omega,$ (53)

namely the Bogoliubov dispersion relation gets linear. Replacing it in Eq.
(38) we obtain the cubic behavior in Eq. (41).

Finally, the negligibility of the high frequency part of the SD, i.e. the
ultraviolet cut-off, arises naturally if one looks to the phonon linear part
of the Bogoliubov spectrum, which is reasonable if we want to study to long-
time dynamics. It is important to highlight that in such a regime our results
do not depend on whether the cut-off is present or not, as we proved in the
end of Appendix B. The SD takes thus the polynomial cubic structure in Eq.
(41). Developing the theory of QBM for such a cubic SD is a central part of
the current work. In the following section we shall show that such a super-
ohmic behavior is related to memory effects.

##  4 Heisenberg Equations

In this section we derive the Heisenberg equations of motion. For the sake of
simplicity, we focus on the case where the BEC is confined in one dimension.
We consider that the impurity is trapped in a harmonic potential, i.e.,

$U(x)=m_{\rm I}\Omega^{2}x^{2}/2,$ (54)

which is a common set-up in ultracold atom systems. The Heisenberg equations
are

$\displaystyle\dot{x}(t)=\frac{i}{\hbar}\left[{H},{x}(t)\right]=\frac{{p}(t)}{m_{\rm
I}},$ (55)
$\displaystyle\dot{p}(t)\\!=\\!\frac{i}{\hbar}\left[{H},{p}(t)\right]=-m_{\rm
I}\Omega^{2}{x}(t)-\hbar\sum_{k}g_{k}\pi_{k}(t),$ (56)
$\displaystyle\dot{b}_{k}(t)=\frac{i}{\hbar}\left[{H},{b}_{k}(t)\right]=-i\omega_{k}{b}_{k}(t)-g_{k}{x}(t),$
(57)
$\displaystyle\dot{{b}}^{\dagger}_{k}(t)=\frac{i}{\hbar}\left[{H},{b}^{\dagger}_{k}(t)\right]=+i\omega_{k}{b}^{\dagger}_{k}(t)-g_{k}{x}(t),$
(58)

where the Hamiltonian is given by in Eq. (25). Equations (55)-(58) can be
combined to obtain

$\ddot{x}(t)+\Omega^{2}x(t)+\frac{\partial}{\partial
t}\int^{t}_{0}\Gamma(t-s)x(s)ds=\frac{B(t)}{m_{\rm I}},$ (59)

as detailed in Appendix A. This is the equation of motion for the impurity
position. It may be viewed as the quantum analogue of a classical stochastic
differential equation, where

${B}(t)=\sum_{k}i\hbar g_{\rm k}({b}^{\dagger}_{\rm k}e^{i\omega_{\rm
k}t}-{b}_{\rm k}e^{-i\omega_{\rm k}t})$ (60)

plays the role of a stochastic force, and $\Gamma(t)$ is the damping kernel
introduced in Eq. (34), which expression may be obtained through an
integration by parts:

$\Gamma(t)=\frac{\tilde{\tau}}{t^{3}}\left[2\Lambda t\cos\left(\Lambda
t\right)-(2-\Lambda^{2}t^{2})\sin\left(\Lambda t\right)\right].$ (61)

An important feature of such an equation is that it is non-local in time,
namely the temporal evolution of $x$ at time $t$ depends on its past history,
i.e. $x(s)$ with $s<t$. We can state that in general the dynamics of an
impurity in a homogeneous BEC in one dimension, described by Eq. (59), carries
a certain amount of memory effects. Actually, there is one special situation
where this equation reduces to a local one, and in particular to a traditional
damped harmonic oscillator with a stochastic force. This is the case
constituted by a damping kernel proportional to a Dirac delta. It is immediate
to prove that this expression of the damping kernel results by a ohmic SD,
i.e. an SD depending linearly on the frequency. As we showed in Sec. 3 this is
not the case of the present system, where the SD is super-ohmic. In conclusion
we can state that the dynamics of an impurity in a BEC always carries a
certain amount of memory effects.

Equation (59) is a second-order linear non-homogeneous differential equation.
Its solution is the sum of a particular one and that of the related
homogeneous equation. For the latter we can proceed by applying the Laplace
transform operator, in order to obtain the solution in terms of the initial
position and velocity. The particular solution, instead, is the convolution
product of the Green function and the position operator. Therefore, the
solution for the impurity position in the Heisenberg picture takes the
following form

$x(t)\\!=\\!G_{1}(t)x(0)+G_{2}(t)\dot{x}(0)+\frac{1}{m_{\rm
I}}\int^{t}_{0}G_{2}(t-s)B(s)ds,$ (62)

where the functions $G_{1}$ and $G_{2}$ are defined through their Laplace
transforms

$\displaystyle\mathcal{L}_{z}[G_{1}(t)]=\frac{z+\mathcal{L}_{z}[\Gamma(t)]}{z^{2}+\Omega^{2}+z\mathcal{L}_{z}[\Gamma(t)]},$
(63)
$\displaystyle\mathcal{L}_{z}[G_{2}(t)]=\frac{1}{z^{2}+\Omega^{2}+z\mathcal{L}_{z}[\Gamma(t)]},$
(64)

and satisfy

$\displaystyle G_{1}(0)=1,\quad\dot{G}_{1}(0)=0,$ (65) $\displaystyle
G_{2}(0)=0,\quad\dot{G}_{2}(0)=1.$ (66)

Equations (63) and (64) depend on the Laplace transform of the damping kernel
which is calculated in Appendix B. Here we present only the final result

$\mathcal{L}_{z}\left[\Gamma(t)\right]=z\tilde{\tau}\left[\Lambda-z\arctan\left(\Lambda/z\right)\right].$
(67)

In order to find the expression of the position impurity as a function of the
time, and so to characterize its motion, we have to invert the Laplace
transforms in Eqs. (63) and (64). In the following two sections, we will
invert this transformations to obtain the complete solution both analytically
and numerically, for the case of untrapped (Sec. 5 ) and trapped impurity
(Sec. 6).

##  5 Untrapped Impurity

Here we consider an untrapped impurity, that is $\Omega=0$ in Eq. (54). In
this case the impurity position in the Heisenberg picture is described by Eq.
(62) with

$\displaystyle\mathcal{L}_{z}[G_{1}(t)]=\frac{z+\mathcal{L}_{z}[\Gamma(t)]}{z^{2}+z\mathcal{L}_{z}[\Gamma(t)]}=1/z,$
(68)
$\displaystyle\mathcal{L}_{z}[G_{2}(t)]=\frac{1}{z^{2}+z\mathcal{L}_{z}[\Gamma(t)]}.$
(69)

To completely characterize the motion of the impurity one has to invert these
Laplace transforms. For Eq. (68) this calculation can be performed
straightforwardly to get

$G_{1}(t)=1.$ (70)

Thus the contribution of the initial position in Eq. (62) is constant in time.
We highlight that such a result may be found even without knowledge of the
form of the damping kernel, and therefore about the form of the SD. Thus, Eq.
(70) holds for any untrapped impurity, regardless of the details of the
coupling.

It is more difficult to handle with the Laplace transform in Eq. (69), due to
the arctangent in the Laplace transform of the damping kernel in Eq. (67). We
are interested in the long-time regime, which corresponds to $z\ll\Lambda$.
The asymptotic expansion of Eq. (67) in this limit reads

$\mathcal{L}_{z}\left[\Gamma(t)\right]=z\tilde{\tau}\Lambda+o(z^{2}),$ (71)

and therefore

$\mathcal{L}_{z}[G_{2}(t)]=\frac{1}{(1+\Lambda\tilde{\tau})z^{2}}\equiv\frac{1}{\alpha
z^{2}},$ (72)

which can be inverted providing

$G_{2}(t)=t/\alpha.$ (73)

We note that such an expression does not fulfill the boundary conditions in
Eq. (66), in particular $\dot{G}_{2}(0)\neq 1$, but this is justified by the
fact that Eq. (73) refers to a long-time behavior.

Equation (73) has been obtained by expanding the Laplace transform of the
damping kernel at the first order in $z/\Lambda$. In general one may perform
such an expansion to the $n^{\text{th}}$ order, and the Laplace transform in
Eq. (69) takes the form of the inverse of a polynomial of $(n+1)^{\text{th}}$
degree. The Laplace transform can be henceforth inverted by computing the
roots of such a polynomial, say $\bar{z}_{\rm k}$, with $k=1,...,n+1$, and
finally one gets

$G_{2}(t)=\sum^{n+1}_{k=1}a_{\rm k}e^{\bar{z}_{\rm k}t},$ (74)

where the $a_{\rm k}$ are constants depending on the roots $\bar{z}_{\rm k}$.
However, one has to notice that some roots have a positive real part,
corresponding to divergent run-away solutions [59]. These divergent roots have
to be dropped out. In fact it is possible to see that they do not satisfy the
condition $z\ll\Lambda$ that we assumed in the beginning of the calculation.
Then, they are negligible in the long-time limit. One may invert the Laplace
transform in Eq. (69) through the general approach of the Bromwich integral.
Here the run-away roots remain out of the integration contour, and therefore
the related residuous do not contribute. We conclude that Eq. (73) represents
the complete expression for the Laplace transform in Eq. (69) in the long-time
limit.

In order to go beyond any asymptotic expansion of Eq. (67), we have to perform
numerically the inversion of the Laplace transform in Eq. (69). There exist
several algorithms aimed to deal with this problem, as carefully discussed in
[60]. Here, we consider the Zakian method, where the inverse Laplace transform
$f(t)$ of a function $F(z)$ is approximated as

$\tilde{f}(t)=\frac{2}{t}\sum^{N}_{j=1}\text{Re}\left[k_{\rm
j}F\left(\alpha_{\rm j}/t\right)\right],$ (75)

with $\alpha_{\rm j}$ and $k_{\rm j}$ constants that can be either complex or
reals. In the limit $N\rightarrow\infty$ it turns that
$\tilde{f}(t)\rightarrow f(t)$. Many studies show that the error in the
approximation is negligible in several situations already for $N=5$, and the
parameters $\alpha_{\rm j}$ and $k_{\rm j}$ are listed in Table 1 of [60].
Applying this method to our problem we invert the Laplace transform in Eq.
(69) without performing any asymptotic approximation. The result is presented
in Fig. 1.

Figure 1: Inverse Laplace transform of Eq. (69). The plot refers to a K
impurity embedded in a gas of Rb atoms with a density $n_{\rm 0}=7\cdot
10^{6}\text{m}^{-1}$ and a coupling constant $g_{\rm B}=2.36\cdot
10^{-37}\text{J}\cdot\text{m}$. Here $\omega_{0}=\hbar n^{2}_{0}/m_{\rm I}$
represents the characteristic frequency. Solid lines refer to result of the
approximation in Eq. (75) for $N=5$, while dashed ones represent analytic
long-time predictions from Eq. (73).

We note that in the long-time limit the numerical solution matches the
analytic one presented in Eq. (73). This linear divergence is approached
through damped oscillations, which characterize the short- and middle-time
regimes. The same behavior is detected by reproducing the calculation by means
of other algorithms in [60], such as the Fourier and Week ones.

In conclusion, the position operator in the Heisenberg picture for an
untrapped impurity is described by the expression in Eq. (62), where $G_{1}$
is given by Eq. (70), while $G_{2}$ is represented in Fig. 1. Only in the
long-time limit it is possible to exhibit an analytic expression for the
second function, Eq. (73). Such a term shows a ballistic form, namely it is
proportional to time. It means that, as time flows the position impurity
becomes larger and larger. The untrapped impurity does not approach the
equilibrium, but runs away from its initial position. This is a reasonable
behavior when we remove trapping. It is natural to characterize quantitatively
the motion of the untrapped impurity with the mean square displacement (MSD)

$\mbox{MSD}(t)=\langle\left[x(t)-x(0)\right]^{2}\rangle,$ (76)

which provides the deviation between the position at time $t$ and the initial
one. In experiments dealing with ultracold gases such a quantity can be
measured [15]. In the long-time limit it is possible write

$\displaystyle\mbox{MSD}(t)=\left(\frac{t}{\alpha}\right)^{2}\langle\dot{x}(0)^{2}\rangle$
$\displaystyle+$ $\displaystyle\frac{1}{2\left(\alpha m_{\rm
I}\right)^{2}}\\!\int^{t}_{0}\\!ds\int^{t}_{0}\\!d\sigma(t-s)(t-\sigma)\langle\\{B(s),B(\sigma)\\}\rangle,$
(77)

where we considered a factorizing initial state
$\rho(t)=\rho_{S}(0)\otimes\rho_{B}$. The initial conditions of the impurity
and bath oscillators are then uncorrelated. Then, averages of the form
$\langle\dot{x}(0)B(s)\rangle$ vanish. The second term on the right-hand side
of Eq. (5) can be treated noting that:

$\langle\\{B(s),B(\sigma)\\}\rangle=2\hbar\nu(s-\sigma)$ (78)

and remembering the definition of the noise kernel in Eq. (32) it turns:

$\displaystyle\nu(s-\sigma)=m_{\rm
I}\tilde{\tau}\int^{\Lambda}_{0}d\omega\cos\left[(s-\sigma)\omega\right]\omega^{3},$
(79)

where the hyperbolic cotangent in the noise kernel in Eq. (32) has been
approximated to one assuming low-temperature. This is an important assumption.
It is possible to check that for realistic values of the physical quantities
such an approximation for the hyperbolic cotangent is reasonable. Replacing
Eq. (79) in the second term of the right-hand side of Eq. (5), and integrating
with respect of time and frequency, it turns

$\mbox{MSD}(t)=\left[\langle\dot{x}(0)^{2}\rangle\\!+\\!\frac{\hbar\tilde{\tau}\Lambda^{2}}{2m_{\rm
I}}\right]\left(\frac{t}{\alpha}\right)^{2}.$ (80)

When this quantity shows a linear dependence on time, the impurity experiences
normal diffusion [38]. Conversely, the MSD is proportional to the square of
time in Eq. (80). Such a behavior is termed super-diffusion and provides a key
signature of the motion of the impurity.

In this context super-diffusion is a consequence of the presence of memory
effects. In [38] a similar calculation is performed considering an ohmic SD,
associated to the absence of memory effects, and leads to a diffusive
behavior. Super-diffusion in Eq. (80) arises due to the super-ohmic character
of the SD. Therefore, it represents a witness of memory effects for a
measurable observable.

Apart from the position, the result we presented above permits to obtain an
expression for the momentum of the impurity in the Heisenberg picture. In
fact, by inserting within Eq. (55) the expression in Eq. (62) with the $G_{2}$
function in Eq. (73), one infers

$p(t)=\frac{1}{\alpha}\left[p(0)+\int^{t}_{0}B(s)ds\right].$ (81)

Equation (81) can be used to compute the average energy of an untrapped
impurity

$E(t)=\frac{\langle p^{2}(t)\rangle}{2m_{\rm I}}.$ (82)

Proceeding as in the derivation of the MSD we find in the low temperature
limit

$\displaystyle E^{LT}(t)=$ $\displaystyle\eta g_{\rm B}n_{\rm
0}+\frac{E(0)}{\alpha^{2}}+\frac{\hbar\tilde{\tau}}{2\alpha^{2}}\Lambda^{2}$
$\displaystyle-$ $\displaystyle\frac{\hbar\tilde{\tau}}{(\alpha
t)^{2}}\left[\cos\left(\Lambda t\right)+\Lambda t\sin\left(\Lambda
t\right)-1\right],$ (83)

where we recover the mean-field energy term, represented now by the first term
on the right hand side.

Figure 2: Average of the energy in the low temperature limit in Eq. (5) as a
function of the time for different values of the coupling strength. The plot
refers to the situation where the initial average energy E(0), as well as the
mean-field term, have been set to zero.

Such a quantity is plotted in Fig. 2, where it is shown that the average of
the energy oscillates initially and after a long time approaches a constant
value. These oscillations represent a non-monotonic behavior of the energy.
Importantly, the increasing parts correspond to a flow of energy directed from
the environment towards the impurity. In [61] such a backflow energy has been
studied for QBM and indicated as a witness of memory effects.

We note that after a long time, in the weak coupling limit, the only term
surviving in Eq. (5) is the mean-field one. In fact the terms in the second
line vanish after a long time. The second term on the right hand-side vanish
in the weak coupling limit, since it is proportional to $\eta^{2}$. In the
end, in the weak coupling limit, the asymptotic value of the energy approaches
the mean-field one.

In the high-temperature limit super-diffusion equally emerges, though with a
different coefficient. Also the energy shows qualitatively the same behavior
as in the low-temperature limit. Nevertheless, the linear approximation in Eq.
(23) fails, as discussed in Appendix C.

##  6 Trapped Impurity

In this section we restore the presence of a harmonic trap, i.e. we look to
the situation where $\Omega>0$ in Eq. (54). To investigate the dynamics of the
impurity we proceed as in the previous section, namely we invert the Laplace
transforms in Eqs. (63) and (64), in order to characterize the expression of
the position operator in the Heisenberg picture, given by in Eq. (62). Here,
an important difference with Sec. 5 is that the expression for $G_{\rm 1}(t)$
cannot be obtained regardless of any information about the analytic structure
of the damping kernel. Conversely, for a trapped impurity the calculation of
$G_{1}$ requires the same attention of $G_{2}$.

We approach the problem from a numerical point of view, with the Zakian method
presented in [60], briefly described in Sec. 5.

Figure 3: Inverse Laplace transforms of Eqs. (63) (Up) and (64) (Down)
obtained through the Zakian method. The figures refer to an impurity of K with
$\Omega=2\pi\cdot 500$ Hz in a gas made up by Rb with a density of $n_{\rm
0}=7(\mu\text{m})^{-1}$ and a coupling strength $g_{\rm B}=2.36\cdot
10^{-37}$J$\cdot$m.

The results are presented in Fig. 3. Both functions exhibit oscillations which
get damped in the long-time regime. In particular, in the long-time limit one
gets

$\lim_{t\rightarrow\infty}G_{1}(t)=0,\quad\lim_{t\rightarrow\infty}G_{2}(t)=0.$
(84)

The boundary conditions in Eqs. (65) and (66) are satisfied. Employing
alternative numerical methods described in [60] to invert the Laplace
transforms we obtain the same result. From the physical point of view it means
that, after a long-time, the contributions of the initial position and
velocity vanish.

For a trapped impurity we are also interested in the long-time behavior of the
particle. We characterize it by means of the position variance, which is a
measurable quantity. Taking into account the behavior of $G_{1}$ and $G_{2}$
functions showed in Fig. 3, the expression for the position variance can be
easily found starting from Eq. (62)

$\langle x^{2}(t)\rangle=\int^{t}_{0}ds\int^{t}_{0}d\sigma
G_{2}(t-s)G_{2}(t-\sigma)\frac{\nu(s-\sigma)}{m^{2}_{\rm I}\hbar^{-1}},$ (85)

where we used the first line in Eq. (79). Replacing the expression for the
noise kernel one gets

$\displaystyle\langle x^{2}(t)\rangle$
$\displaystyle\\!=\\!\frac{\hbar}{m^{2}_{\rm
I}}\int^{\Lambda}_{0}J(\omega)\coth\left(\hbar\omega/2k_{\rm
B}T\right)d\omega$ $\displaystyle\\!\times\int^{t}_{0}ds\int^{t}_{0}d\sigma
G_{2}(t-s)G_{2}(t-\sigma)\cos[\omega(s-\sigma)]$
$\displaystyle\equiv\frac{\hbar}{m^{2}_{\rm
I}}\int^{\Lambda}_{0}J(\omega)\coth\left(\hbar\omega/2k_{\rm
B}T\right)\phi(t,\omega)d\omega.$ (86)

We stress that once the long-time conditions in Eq. (84) are fulfilled, it is
not necessary to suppose that the state of the global system is initially non-
correlated, as we did in Sec. 5. Equation (86) turns into

$\displaystyle\phi(t,\omega)=$
$\displaystyle\frac{1}{2}\int^{t}_{0}ds\int^{t}_{0}d\sigma
G_{2}(t-s)G_{2}(t-\sigma)$ $\displaystyle\times$
$\displaystyle\left[e^{i\omega s}e^{-i\omega\sigma}+c.c.\right]$
$\displaystyle=$
$\displaystyle\frac{1}{2}\int^{t}_{0}d\tilde{s}e^{-i\omega\tilde{s}}G_{2}(\tilde{s})$
(87) $\displaystyle\times$
$\displaystyle\int^{t}_{0}d\tilde{\sigma}e^{i\omega\tilde{\sigma}}G_{2}(\tilde{\sigma})+c.c.,$
(88)

where we introduced

$\tilde{s}=t-s,\quad\tilde{\sigma}=t-\sigma.$ (89)

We are interested in the long-time limit, $t\rightarrow\infty$. In this limit,
one gets

$\phi(t,\omega)=\mathcal{L}_{-i\omega}\left[G_{2}(t)\right]\mathcal{L}_{+i\omega}\left[G_{2}(t)\right].$
(90)

Replacing Eq. (64) into Eq. (90), we obtain the final expression for the
position variance:

$\langle
x^{2}\rangle=\frac{\hbar}{2\pi}\int^{+\Lambda}_{-\Lambda}d\omega\coth\left(\hbar\omega/2k_{\rm
B}T\right)\tilde{\chi}^{\prime\prime}(\omega),$ (91)

where

$\tilde{\chi}^{\prime\prime}(\omega)=\frac{1}{m_{\rm
I}}\frac{\zeta(\omega)\omega}{\left[\omega\zeta(\omega)\right]^{2}+\left[\Omega^{2}-\omega^{2}+\omega\theta(\omega)\right]^{2}},$
(92)

and

$\displaystyle\zeta\left(\omega\right)=\text{Re}\\{\mathcal{L}_{\tilde{z}}\left[\Gamma(t)\right]\\}=\frac{\pi\tilde{\tau}}{2}\omega^{2}+o\left(\frac{\omega}{\Lambda}\right)^{5},$
(93)
$\displaystyle\theta\left(\omega\right)=\text{Im}\\{\mathcal{L}_{\tilde{z}}\left[\Gamma(t)\right]\\}=-\tilde{\tau}\Lambda\omega+\frac{\tilde{\tau}}{\Lambda}\omega^{3}+o\left(\frac{\omega}{\Lambda}\right)^{5}.$
(94)

with $\tilde{z}=-i\omega+0^{+}$. The expression in Eq. (91), endowed by Eqs.
(92) and (93), completely determines the position variance for an impurity
trapped in a harmonic potential. We emphasize that such an expression has been
obtained just by considering the long-time limit of the solution of the
Heisenberg equations in Eq. (62). It is possible to note, however, that it
corresponds to that achieved in the context of the linear response theory by
means of the fluctuation-dissipation theorem [38], as discussed in detail in
Appendix D. Indeed, $\tilde{\chi}^{\prime\prime}$ can be seen as the imaginary
part of the Fourier transform of the linear response to an external force
applied to the system, at the equilibrium.

In conclusion, in presence of a harmonic trap the impurity approaches the
equilibrium in the long-time limit. We describe such a state through position
and momentum variance (a similar expression to Eq. (91) is also found for the
momentum). In particular we look to

$\delta_{x}=\sqrt{\frac{2m_{\rm I}\Omega\langle
x^{2}\rangle}{\hbar}},\quad\delta_{p}=\sqrt{\frac{2\langle
p^{2}\rangle}{m_{\rm I}\hbar\Omega}},$ (95)

which represents the position and momentum variances regularized in order to
be dimensionless. In these units the Heisenberg principle reads as
$\delta_{x}\delta_{p}\geq 1$, so the Heisenberg threshold is set to be equal
to one. These quantities do not depend on time, because they refer to an
equilibrium stationary state. We study the dependence of $\delta_{x}$ on the
system parameters, such as temperature and interaction strength, that can be
tuned in experiments.

In Fig. 4 we show the position variance for a trapped impurity as a function
of the temperature, for several values of the coupling strength.

Figure 4: Dimensionless position variance in Eq. (95) as a function of the
temperature, for several values of the interaction strength. The red dashed
line represents the function $\sqrt{2T}$, associated to the equipartition
theorem. The figure refers to an impurity: of K with $\Omega=2\pi\cdot 500$ Hz
in a gas made up by Rb with a density of $n_{\rm 0}=7(\mu\text{m})^{-1}$ and a
coupling strength $g_{\rm B}=2.36\cdot 10^{-37}$J$\cdot$m.

Such a result follows from both a numerical and analytic integration. In the
second case, one may proceed by noting that the integrand function rapidly
vanishes as $\omega$ increases, so we approximate Eq. (91) with an integral
over the whole real axis (also note that the integrand is an even function).
Therefore it is possible to apply the Residue theorem. We also stress that Eq.
(91) refers to the asymptotic expansion in Eq. (93), justified in the long-
time limit. Such an expansion has been performed till the fifth order in
$\omega/\Lambda$, but even going to higher orders we recover the same result
as in Fig. 4.

In physical grounds we note from Fig. 4 that for $k_{\rm
B}T/\hbar\Omega\gtrsim 0.5$, the position variance grows as the square root of
the temperature. Indeed, we detect the behavior provided by
$\delta_{x}\sim\sqrt{2T}$ (red dashed line), in agreement with the
equipartition theorem. We note that as the temperature increases the value of
the position variance turns to be less dependent of the interaction strength.

In the other limit, that is when one approaches the zero-temperature limit
(for $k_{\rm B}T/\hbar\Omega\lesssim 0.5$), we find that $\delta_{x}<1$. This
means that the position variance is smaller than the value related to the
Heisenberg threshold in these units. This important effect is named genuine
position squeezing, and corresponds to high spatial localization of the
impurity.

Figure 5: Dimensionless position variance in Eq. (95) as a the interaction
strength, for several values of the temperature, satisfying $k_{\rm
B}T/\hbar\Omega\ll 0.5$. The figure refers to an impurity: of K with
$\Omega=2\pi\cdot 500$ Hz in a gas made up by Rb with a density of $n_{\rm
0}=7(\mu\text{m})^{-1}$ and a coupling strength $g_{\rm B}=2.36\cdot
10^{-37}$J$\cdot$m.

This is an important resource in quantum technologies, and therefore we would
like to find under which parameters such an effect is enhanced. In Fig. 5 we
plot the value of the dimensionless position variance in Eq. (95) as a
function of the interaction strength, obtained by exploring values of the
temperature as low as possible. We point out that genuine squeezing is
maximized, i.e. $\delta_{x}$ gets smaller, if the temperature tends to zero,
and the interaction becomes stronger. For instance, tuning $\eta\approx 5$ and
setting $T=2$nk it is possible to reach $\delta_{x}\approx 0.7$. This degree
of squeezing may be enhanced by increasing the value of $\eta$. In this case,
our results would lie at the border of the validity regime of our Hamiltonian
model. Note that in experiments the value of $\eta$ may be pushed towards very
high values, for instance $\eta=30$ [15]. Of course, the present theory cannot
be employed to investigate such a regime. However, we remark that for
$\eta\lesssim 7$, where the theory is well defined, we already find a good
degree of genuine position squeezing. We underline, anyway, that even for
$\eta\lesssim 7$, the investigation of such a limit would be not possible
through the tools developed in [12], basically due to the violations of the
Heisenberg principle at low temperature. Here, instead, Heisenberg principle
is fulfilled.

Finally, we would like to underline that the genuine position squeezing
appears both for attractive and repulsive interactions. This is in agreement
with the results presented in [15] (see Fig. 4), where it has been shown that
the position variance of the impurity does not depend on the sign of the
interaction. To understand this, let us first note that in the presence of a
trap confining the BEC, the density of the condensate, will be spatially
changing and peaked in the center of the trap. If impurity-boson interactions
are repulsive, one should expect that the impurity will be pushed away from
the center and localized around the distance $D$ resulting from am interplay
between the force trapping the impurity and the force resulting from the mean
field interactions. While due to the parity symmetry $\langle x\rangle=0$, one
should expect that $\langle x^{2}\rangle\simeq D^{2}$ and, unless $D$ is very
small, there will be no squeezing of the position. In contrast, for attractive
interactions the impurity will be localized in the center of the trap and
squeezing will be possible. In Ref. [62] it is considered the model case when
the bath is coupled to the square of the position of the impurity. This is the
first step towards the study of the full problem, in which impurity couples to
a confined condensate with spatially dependent density and complicated
spatially dependent Bogoliubov-de Gennes modes.

In contrast, in the present case we consider a spatially homogeneous
condensate with constant density and Bogoliubov-de Gennes modes, which are
plane waves. Moreover, we linearize the spatial dependence of modes, assuming
self-consistently that the impurity is localized in the region of $x$ allowing
for such a linearization. The above arguments do not apply in thus case, and
our results: i) do not depend on the sign of impurity-bath interaction; ii)
squeezing is thus possible for both repulsive and attractive interactions.

##  7 Conclusions

We presented a discussion concerning the physics of an impurity embedded in a
homogeneous BEC, adopting an open quantum system point of view, in particular
recalling the paradigmatic quantum Brownian motion model. In this framework
the impurity may be treated as a Brownian particle interacting with a bath of
Bogoliubov modes. Such an approach is suitable to describe the dynamics of the
impurity, namely to characterize its motion. We pursued this task employing
Heisenberg equations, finding an equation for the impurity position (see Eq.
(59)) which can be seen as the quantum analogue of the stochastic Langevin
equation. Equation (59) is non-local in time, suggesting that the dynamics of
an impurity in a homogeneous BEC carries a certain amount of memory effects.
The presence of memory effects is embodied in the structure of the spectral
density, which we calculated explicitly in Sec. (3). The spectral density
shows a super-ohmic behavior.

The solution of the equation for the impurity position, together with the
spectral density, permits to qualify the motion of the impurity. We
distinguished two situations: the case where the impurity is trapped in a
harmonic potential (Sec. 6) and that in which there is not any trap (Sec. 5).
When the impurity is untrapped it does not approach the equilibrium, and runs
away from the point where was initially localized. We characterized
quantitatively such a behavior by computing the mean square displacement,
providing information about the portion of space explored by the impurity
during its random motion. This quantity can be measured in experiments [15].
Such a quantity shows a super-diffusive form, namely it is proportional to the
square of time. The super-diffusive behavior is a consequence of the presence
of memory effects. Therefore, this result constitutes a witness of non-
Markovianity on a measurable observable.

In [61] the presence of memory effects has also been related to a back-flow of
energy, directed from the environment to the impurity. We evaluated the
average of energy for an untrapped impurity as a function of the time, and we
found non-monotonic behavior: the impurity does not just release energy to the
environment (dissipation), but also acquires energy from it. It is an open
question to understand if such a back-flow of energy may be employed as a
resource for quantum technologies.

In presence of a trap the impurity approaches the equilibrium after a long-
time evolution. It is also interesting to study the position and momentum
variance of such a final stable state. We study the behavior of both
quantities as a function of the temperature and of the interaction strength.
We find that as the coupling gets stronger, and as the temperature tends to
zero, the impurity experiences genuine position squeezing, namely its position
variance takes values smaller than the Heisenberg principle. This feature
corresponds to high localization in space of the impurity. Such an effect may
be optimized by increasing the strength of the coupling between the impurity
and the BEC. Here, one has to take into account the constraint imposed on the
interaction in order to ensure the validity of the Fröhlich Hamiltonian model,
consisting of an upper bound on the quantity $\eta$. Also within the
limitations defined by this condition we find a good degree of genuine
position squeezing.

Finally, we would like to highlight that the method we developed in this work
can be exported to others ultracold systems. For instance, in [44] it has been
proved that the dynamics of a bright soliton in a superfluid in one dimension
is described by an equation showing the same form of that in Eq. (59). The
spectral density for the system is in some circumstances proportional to the
third power of the frequency of the Bogoliubov modes, in agreement with the
current work.

Insightful discussion with Victor Galitski, Jan Wehr, Pietro Massignan, and
Nils-Eric Guenther are gratefully acknowledged. This work has been funded by a
scholarship from the Programa Màsters d’Excel-lència of the Fundació
Catalunya-La Pedrera, ERC Advanced Grant OSYRIS, EU IP SIQS, EU PRO QUIC, EU
STREP EQuaM (FP7/2007-2013, No. 323714), Fundació Cellex, the Spanish MINECO
(SEVERO OCHOA GRANT SEV-2015-0522, FOQUS FIS2013-46768, FISICATEAMO
FIS2016-79508-P), and the Generalitat de Catalunya (SGR 874 and
CERCA/Program).

##  Appendix A Derivation of the equation for the impurity position

In this Appendix we show in detail the calculation leading to Eq. (59). The
starting point is constituted by the Eqs. (55)-(58). The first step is to
derive both sides of Eq. (55) and to replace the result in Eq. (56), thus
obtaining an equation just for the position of the impurity

$\ddot{x}{(t)}+\Omega^{2}x(t)=-i\sum_{k}\frac{\hbar g_{\rm k}}{m_{\rm
I}}\left[b_{\rm k}(t)-b^{\dagger}_{\rm k}(t)\right].$ (96)

The time dependence of the Bogoliubov modes can be extracted from Eqs. (57)
and (58). They are linear but non-homogeneous first-order differential
equations. Therefore their solution is the sum of that of the related
homogeneous one and a particular integral. The former can be easily obtained
since it is just that of a harmonic oscillator,

$b_{\rm k}(t)=b_{\rm k}e^{-i\omega_{\rm k}t}+h^{-}_{\rm k}(t),\quad
b^{\dagger}_{\rm k}(t)=b^{\dagger}_{\rm k}e^{+i\omega_{\rm k}t}+h^{+}_{\rm
k}(t).$ (97)

The quantities $h^{-}_{\rm k}$ and $h^{+}_{\rm k}$ represent the particular
solutions of Eqs. (57) and (58), that may be expressed as convolution product
of the unknown function $x(t)$ and the Green function

$h^{\pm}_{\rm k}(t)=\int^{t}_{0}G^{\pm}_{\rm k}(t-s)x(s)ds.$ (98)

Then, the problem of solving the Heisenberg equations for the Bogoliubov modes
reduces to that of finding the Green function of Eqs. (57) and (58). The Green
function is defined by the equation

$\dot{G}^{\pm}_{\rm k}(t)\mp i\omega_{\rm k}{G}^{\pm}_{\rm k}(t)=-g_{\rm
k}\delta(t),$ (99)

which applying Fourier transform turns into

$\mathcal{F}_{\tilde{\omega}}\left[G^{\pm}_{\rm k}(t)\right]=-\frac{ig_{\rm
k}}{\tilde{\omega}\pm\omega}.$ (100)

Thus, the problem of determining the Green function is solved if one performs
the inversion of the Fourier transform in Eq. (100). Namely, one has to
calculate the following integral

$G^{\pm}_{\rm k}(t)=-\frac{ig_{\rm
k}}{2\pi}\mathcal{P}\int^{+\infty}_{-\infty}\frac{e^{i\tilde{\omega}t}}{\tilde{\omega}\pm\omega_{\rm
k}}d\tilde{\omega},$ (101)

where we introduced the principal Cauchy part $\mathcal{P}$ because otherwise
the integral is not well-defined in $\tilde{\omega}=\pm\omega_{\rm k}$.

Figure 6:  Path in complex plane to solve the integral in Eq. (101).

The integral in Eq. (101) can be solved recalling the Jordan Lemma, selecting
the path in Fig. 6. It follows

$G^{\pm}_{\rm k}(t)=\frac{g_{\rm k}}{2}\exp\left(\mp i\omega_{\rm k}t\right).$
(102)

In conclusion, Eq. (96) takes the form

$\displaystyle\ddot{x}(t)+\Omega^{2}x(t)-\frac{\hbar}{m_{\rm
I}}\sum_{k}g^{2}_{\rm k}\int^{t}_{0}x(s)\sin\left[\omega_{\rm
k}(t-s)\right]ds$ (103) $\displaystyle=$ $\displaystyle\frac{B(t)}{m_{\rm
I}},$

where

${B}(t)\equiv-\sum_{k}\hbar g_{\rm k}{\pi}_{\rm k}(t)=\sum_{k}i\hbar g_{\rm
k}({b}^{\dagger}_{\rm k}e^{i\omega_{\rm k}t}-{b}_{\rm k}e^{-i\omega_{\rm
k}t}).$ (104)

Equation (103) can be expressed in terms of the dissipation kernel, Eq. (33)

$\ddot{x}(t)+\Omega^{2}x(t)-\frac{1}{m_{\rm
I}}\int^{t}_{0}\lambda(t-s)x(s)ds=\frac{B(t)}{m_{\rm I}}.$ (105)

One can also introduce the damping kernel in Eq. (34). The third term in the
first hand-side of Eq. (103) so writes as

$\displaystyle-\frac{1}{m_{\rm
I}}\int^{t}_{0}\lambda(t-s)x(s)ds=\int^{t}_{0}\dot{\Gamma}(t-s)x(s)ds=$ (106)
$\displaystyle=\frac{\partial}{\partial
t}\int^{t}_{0}\Gamma(t-s)x(s)ds-\Gamma(0)x(t).$

Accordingly it is possible to express Eq. (105) as

$\ddot{x}(t)+\tilde{\Omega}^{2}x(t)+\frac{\partial}{\partial
t}\int^{t}_{0}\Gamma(t-s)x(s)ds=\frac{B(t)}{m_{\rm I}},$ (107)

in which we introduced the renormalized frequency

$\tilde{\Omega}^{2}=\Omega^{2}-\Gamma(0).$ (108)

Hereafter we neglect the contribution to the frequency provided by
$\Gamma(0)$. This term grows as the interaction strength increase and could
lead to a negative renormalized frequency. In this context the impurity
experiences instability. By deriving the equations of motion directly from
Hamiltonian Eq. (27) one avoids the instability problem, at the cost of
introducing and additional counter term ad hoc. With the procedure here
presented we can identify in the equations of motion the effect of the absence
of such a term.

##  Appendix B Laplace tranform of the damping kernel

In this Appendix we derive the expression of the Laplace transform of the
damping kernel presented in Eq. (67). We perform the calculation for a cubic
SD with a general ultraviolet regularization

$J(\omega)=m_{\rm I}\tilde{\tau}\omega^{3}\Theta(\omega,\Lambda),$ (109)

with $\Theta(\omega,\Lambda)>0$ specifying the dependence on the cut-off of
the SD. The SD in Eq. (67) corresponds to the particular case when

$\Theta(\omega,\Lambda)=\theta(\omega-\Lambda),$ (110)

with $\theta(\cdot)$ the Heaviside step function. We will also compute the
Laplace transform considering an exponential cut-off

$\Theta(\omega,\Lambda)=\exp\left(-\Lambda\omega\right),$ (111)

showing that a much complicated expression turns out.

In general, we assume that $\Theta(\omega,\Lambda)$ decays fast enough so that

$\int_{0}^{\infty}e^{-zt}\left(\int_{0}^{\infty}|\omega^{2}\Theta(\omega,\Lambda)\cos(\omega
t)|d\omega\right)dt<\infty.$ (112)

Note that such a behavior covers both the sharp and the exponential dependence
on the cut-off. Therefore, by Fubini-Tonelli theorem one can interchange the
integrals in the following. For $z>0$, it results in

$\displaystyle\mathcal{L}\left[\Gamma\left(t\right)\right]_{z}$
$\displaystyle\\!=\\!\tilde{\tau}\int_{0}^{\infty}e^{-zt}\left(\int_{0}^{\infty}\omega^{2}\Theta(\omega,\Lambda)\cos(\omega
t)d\omega\right)dt$ (113)
$\displaystyle=\tilde{\tau}z\int_{0}^{\infty}\frac{\omega^{2}}{\omega^{2}+z^{2}}\Theta(\omega,\Lambda)d\omega,$
(114)

where we used the expression

$\int^{\infty}_{0}\exp\left(-zt\right)\cos(\omega
t)dt=\frac{z}{\omega^{2}+z^{2}},$ (115)

that is the result of an integration by parts.

Choosing the cut-off function in Eq. (110), we obtain

$\mathcal{L}\left[\Gamma\left(t\right)\right]_{z}=z\tilde{\tau}\left[\Lambda-z\arctan\left(\Lambda/z\right)\right],$
(116)

as we presented in Eq. (67). If we consider the cut-off function in Eq. (111),
the Laplace transform of the damping kernel results to be expressed in terms
of the Mejer function. Such a function is not suitable to handle in an
analytic calculation. For the environment given by a BEC, the analytical
behavior together with the form of the Bogoliubov spectrum justify the
consideration of a sharp cut-off.

Let us finally note that, regardless of its analytic form, the presence of the
cut-off introduced in Eq. (51) does not affect our results. This is because
our results refer to the long-time limit, which is not influenced by the high-
frequency part of the SD. To prove this, one should compare the Laplace
transform of the damping kernel induced by the SD in Eq. (51) (shown in Eq.
(116)) with that induced by the SD in Eq. (41), where there is not any cut-
off. This is enough because such a Laplace transform is the only object
carrying information about the SD along all the theory.

However, the calculation of the Laplace transform of the damping kernel
induced by the “uncutted" SD in Eq. (41) is not an easy task, due to the very
complicated form at high-frequency (see Eq. (43)). It is given by:

$\mathcal{L}\left[\Gamma\left(t\right)\right]_{z}=\tilde{\tau}z\int_{0}^{\infty}\frac{\omega^{2}}{\omega^{2}+z^{2}}\chi_{1d}(\omega,\Lambda)d\omega.$
(117)

We may evaluate it at the first-order in $z$, which is the relevant one at
long-times,

$\mathcal{L}\left[\Gamma\left(t\right)\right]_{z}=\tilde{\tau}z\int_{0}^{\infty}\chi_{1d}(\omega,\Lambda)d\omega+o(z/\Lambda)^{2}.$
(118)

Since

$\int_{0}^{\infty}\chi_{1d}(\omega,\Lambda)d\omega=\Lambda,$ (119)

we recover the expression in Eq. (116), derived with a sharp cut-off.

In the end, we conclude that at the long-times our results do not depend on
whether the SD is cutted or not, i.e. on the existence of the cut-off. This
could also be inferred by recalling the Tauberian theorem [63, 64] according
to which the long time behavior of a function (in the time domain) is
determined by the low frequency behavior of its Laplace transform (in the
frequency domain). Since the low frequency behavior of the Laplace transform
of the damping kernels for both spectral densities coincides, the long time
dynamics of the impurity do not depend on whether there is a cutoff or not.

##  Appendix C Validity of the linear approximation

In Sec. 2 we proved that the Hamiltonian of an impurity in a gas may be
expressed in the form of that of the QBM. A crucial step to perform this task
consists of a linear expansion of the exponential appearing in Eq. (21). The
present Appendix is devoted to discuss the validity of such an approximation,
namely we wonder for which values of the system parameters the condition

$kx\ll 1,$ (120)

holds, allowing the expansion in Eq. (24). Here, $k$ represents the wave
number of the Bogoliubov modes, depending on their frequency as showed in Eq.
(39). This function increases monotonically as

$\displaystyle k\approx\omega/c\quad\text{for}\quad\omega\ll\Lambda,$ (121)
$\displaystyle k\sim\sqrt{\omega}\quad\text{for}\quad\omega\gg\Lambda,$ (122)

hence we can minimize the left-hand side of Eq. (120) by looking to the small
frequency regime in Eq. (121). Note that this is in agreement with our
treatment, because all the results presented here refer to the phonon linear
part of the Bogoliubov dispersion relation. In fact such a portion of the
Bogoliubov spectrum is embodied in the super-ohmic form of the SD we have
considered (see Sec. 3). Moreover, the ultraviolet sharp cut-off, $\Lambda$,
permits to get rid of the contribution due to the non-linear part in Eq.
(122). The condition in Eq. (120) is then

$\frac{\omega}{c}x\ll 1,$ (123)

and recalling the expression

$c=\sqrt{2}\xi\Lambda,$ (124)

it turns into

$\frac{\omega}{\Lambda}\frac{x}{\xi}\ll\sqrt{2}.$ (125)

We point out that, in order to fulfill the condition, the position of the
impurity may acquire large values, provided we consider a frequency much
smaller than $\Lambda$.

The value of the frequency depends in general on the temperature. For a system
of bosons the energy at a given temperature $T$ is

$E(\omega)=\frac{\hbar\omega}{e^{\frac{\hbar\omega}{k_{\rm B}T}}-1}\leq k_{\rm
B}T,$ (126)

accordingly one can evaluate Eq. (125) as

$\frac{k_{\rm B}T}{\hbar c}x\ll 1.$ (127)

Note that, because of the inequality in Eq. (126), the condition in Eq. (127)
provides an upper bound for Eq. (120).

We have to evaluate now the position of the impurity. The quantity appearing
in Eq. (120) is an operator. The impurity is modeled as a harmonic oscillator,
so its state is Gaussian. Therefore it ensues

$x=\langle x\rangle+\Delta_{x}=\Delta_{x},$ (128)

where we considered that, for a harmonic oscillator, the average value of the
position is zero. Finally, we can state that the linear approximation
underlying our analysis is fulfilled if it is satisfied that

$\frac{k_{\rm B}T}{\hbar c}\Delta_{x}\ll 1.$ (129)

To evaluate $\Delta_{x}$, and thus Eq. (129), we have to distinguish again the
case of trapped and untrapped impurity.

For the trapped case we consider in Eq. (129) the fluctuation $\Delta_{x}$ as
that given by $\delta_{x}$ in Eq. (95),

$\chi^{(\text{Tr})}\equiv\frac{k_{\rm B}T}{\hbar c}\delta_{x}(T,\Omega).$
(130)

If $\chi^{(\text{Tr})}<1$ it is possible to state that the linear
approximation, Eq. (23), is fulfilled.

Figure 7: Behavior of $\chi^{(\text{Tr})}$ [see Eq. (130)] as a function of
the temperature and of the frequency of the trap. For $\chi^{(\text{Tr})}<1$
the linear approximation, Eq. (23), is appropriate for the trapped impurity.
The figure refers to an impurity of K in a gas made up by Rb with a density of
$n_{\rm 0}=7(\mu\text{m})^{-1}$ and a coupling strength $g_{\rm B}=2.36\cdot
10^{-37}$J$\cdot$m. The interaction strength is $\eta=1$.

We plot in Fig. 7 the quantity in Eq. (130) as a function of the temperature
and the trap frequency. Note that $\delta_{x}$ depends also on the interaction
strength, but we focus on the case where $\eta=1$. In general, as the
temperature grows the position variance approaches the value predicted by the
equipartition theorem, so it gets coupling-independent. At low temperature,
instead, the value of the position variance approaches one at weak coupling,
i.e. $\eta\approx 1$, while in general it is smaller (the impurity experiences
genuine position squeezing). Therefore, the value of the position variance at
$\eta=1$ represents an upper bound. Therefore, if the linear approximation is
fulfilled for $\eta=1$, it also holds for other, smaller or larger, values of
$\eta$.

Figure 7 shows that the linear approximation for a trapped impurity is
fulfilled at low temperature. Moreover, it is very well satisfied as the trap
potential gets more and more steep, as pointed out also in [47]. In
particular, we see that for the value we selected to detect genuine position
squeezing, i.e. $\Omega=2\pi\cdot 500$ and $T\geq 2$nm the threshold in Eq.
(130) takes a value smaller than $0.1$, suggesting that the condition allowing
the linear approximation is very well satisfied.

Next we consider the untrapped impurity. In this case $\Delta_{x}$ may be
evaluated by means of the MSD in Eq. (80). Considering the case in which
$\langle\dot{x}^{2}(0)\rangle=0$, the left-hand side in Eq. (127) is

$\chi^{(\text{Un})}\equiv\frac{k_{\rm B}T}{\hbar
c}\sqrt{\frac{\hbar\tilde{\tau}}{2m_{\rm I}}}\frac{t\Lambda}{\alpha(\eta)}.$
(131)

Again, in order to state that our linear approximation is satisfied the
quantity in Eq. (131) has to be smaller than one.

Figure 8: Behavior of $\chi^{(\text{Un})}$ [see Eq. (131)] as a function of
the temperature and of the time. For $\chi^{(\text{Un})}<1$ the linear
approximation, Eq. (23), is appropriate for the untrapped case. The figure
refers to an impurity of K in a gas made up by Rb with a density of $n_{\rm
0}=7(\mu\text{m})^{-1}$ and a coupling strength $g_{\rm B}=2.36\cdot
10^{-37}$J$\cdot$m. The interaction strength is $\eta=1$.

In Fig. 8 we plotted the value of the quantity $\chi^{(\text{Un})}$ introduced
in Eq. (131). Even in this case we find that the linear approximation
underlying our treatment is very well fulfilled at low-temperature. For an
untrapped impurity the validity of such an approximation holds until a certain
value of the time

$t^{(\text{cr})}=\frac{\hbar c}{k_{\rm B}T}\sqrt{\frac{2m_{\rm
I}}{\hbar\tilde{\tau}}}\frac{\alpha(\eta)}{\Lambda},$ (132)

corresponding to the instant in which the quantity in Eq. (7) is equal to one,
i.e. we have

$\chi^{(\text{Un})}(t^{(\text{cr})})=1.$ (133)

Figure 9: Critical time in Eq. (132) as a function of the temperature for
different values of the interaction strength. Below the curves the linear
approximation, Eq. (23), is fulfilled. The figure refers to an impurity of K
in a gas made up by Rb with a density of $n_{\rm 0}=7(\mu\text{m})^{-1}$ and a
coupling strength $g_{\rm B}=2.36\cdot 10^{-37}$J$\cdot$m.

The critical time in Eq. (132) is plotted in Fig. 9 as a function of the
temperature, for different values of the interaction strength. Each line
corresponds to a different value of $\eta$. Once one fixes $\eta$, the
corresponding line determines the value of time and temperature for which the
threshold in Eq. (131) is equal to one. Accordingly, above the line we have a
forbidden area. Again, we find that this forbidden area is enlarged as the
temperature grows, suggesting that our approximation works at low temperature.

##  Appendix D Derivation of the position variance by means of the
fluctuation-dissipation theorem

In this Appendix we present an alternative way to derive the position variance
for a trapped impurity, employing the linear response theory. For this goal we
introduce the response function,
$\chi(t)=\frac{i}{\hbar}\theta(t)\langle[x(t),x(0)]\rangle$, describing the
linear response of the system to an external force at the equilibrium. Here
$\theta(t)$ is the step function, specifying causality. The response function
defined here is also known as generalized susceptibility. Taking the
commutator with $x(0)$ on both sides of (59) and then averaging (note the
initial slip term and the stochastic forcing term are absent once a commutator
of $x(0)$ is applied to them), we obtain the following c-valued equation for
$\chi(t)$:

$\ddot{\chi}(t)+\Omega^{2}\chi(t)+\int_{0}^{t}\Gamma(t-s)\dot{\chi}(s)ds=0$
(134)

with $\chi(0)=0$, $\dot{\chi}(0)=1/m_{\rm I}$.

Next, consider the symmetric position autocorrelation function, defined as

$S(t)=\frac{1}{2}\langle\\{x(t),x(0)\\}\rangle.$ (135)

We recall the fluctuation-dissipation relation of Callen-Welton [38], which
relates the equilibrium fluctuation to the response function in the frequency
domain as follows:

$\tilde{S}(\omega)=\hbar\coth\left(\frac{\hbar\omega}{2k_{\rm
B}T}\right)\tilde{\chi}^{\prime\prime}(\omega),$ (136)

where

$\tilde{S}(\omega)=\int^{\infty}_{-\infty}dt\cos(\omega t)S(t),$ (137)

and

$\tilde{\chi}^{\prime\prime}(\omega)=\text{Im}\\{\mathcal{F}_{\omega}\left[\chi(t)\right]\\}$
(138)

Using this, we can obtain the symmetric position autocorrelation function:

$S(t)=\int^{\infty}_{-\infty}d\omega\frac{\hbar}{2\pi}\coth\left(\frac{\hbar\omega}{2k_{\rm
B}T}\right)\tilde{\chi}^{\prime\prime}(\omega)\cos(\omega t).$ (139)

from which the position variance at the equilibrium follows:

$\langle
x^{2}\rangle=S(0)=\int^{\infty}_{-\infty}d\omega\frac{\hbar}{2\pi}\coth\left(\frac{\hbar\omega}{2k_{\rm
B}T}\right)\tilde{\chi}^{\prime\prime}(\omega)$ (140)

Note that
$\tilde{\chi}^{\prime\prime}(w)=\text{Im}\\{\mathcal{L}_{\bar{z}}[\chi(t)]\\}$,
in which

$\displaystyle\mathcal{L}_{\bar{z}}[\chi(t)]$ $\displaystyle=\frac{1}{m_{\rm
I}}\mathcal{L}_{\bar{z}}[G_{2}(t)]$ (141) $\displaystyle=\frac{1}{m_{\rm
I}}\frac{1}{(\Omega^{2}-\omega^{2})-i\omega\mathcal{L}_{\bar{z}}[\Gamma(t)]},$
(142)

where we used the expression in Eq. (64) for the Laplace transform of $G_{2}$,
and $\bar{z}=-\omega+0^{+}$ as defined in Sec. 6. Therefore, we have

$\tilde{\chi}^{\prime\prime}(\omega)=\frac{1}{m_{\rm
I}}\frac{\omega\xi(\omega)}{\left[\Omega^{2}-\omega^{2}-\omega\theta(\omega)\right]^{2}+\left[\omega\xi(\omega)\right]^{2}}.$
(143)

This expression has the same form as Eq. (92). Accordingly, the position
variance in Eq. (140), obtained recalling the fluctuation-dissipation theorem,
corresponds to that in Eq. (91), calculated by solving Heisenberg equation and
adopting the SD in Eq. (51). Therefore, we prove that the two methods lead to
the same result.

One can proceed in a similar way to evaluate the MSD for an untrapped impurity
in the context of linear response theory, rather than solving Heisenberg
equation. This method leads however to a very complicated expression for the
MSD, that does not provide a more convenient alternative to Eq.(5). This topic
is discussed in detail in [65].

## References

  * Landau and Pekar [1948] L. D. Landau and S. I. Pekar,  _Effective mass of a polaron_ , Zh. Eksp. Teor. Fiz. (1948).
  * Fröhlich [1954] H. Fröhlich,  _Electrons in lattice fields_ , Advances In Physics 3(11):325 (1954).
  * Alexandrov and Devreese [2009] A. Alexandrov and J. Devreese,  _Advances in Polaron Physics_, Springer Series in Solid-State Sciences (Springer, 2009).
  * Schirotzek  _et al._ [2009] A. Schirotzek, C.-H. Wu, A. Sommer, and M. W. Zwierlein,  _Observation of Fermi Polarons in a Tunable Fermi Liquid of Ultracold Atoms_ , Phys. Rev. Lett. 102, 230402 (2009).
  * Kohstall  _et al._ [2012] C. Kohstall, M. Zaccanti, M. Jag, A. Trenkwalder, P. Massignan, G. M. Bruun, F. Schreck, and R. Grimm,  _Metastability and coherence of repulsive polarons in a strongly interacting Fermi mixture_ , Nature 485, 615 (2012).
  * Koschorreck  _et al._ [2012] M. Koschorreck, D. Pertot, E. Vogt, B. Fröhlich, M. Feld, and M. Köhl,  _Attractive and repulsive Fermi polarons in two dimensions_ , Nature 485, 619 (2012).
  * Massignan  _et al._ [2014] P. Massignan, M. Zaccanti, and G. M. Bruun,  _Polarons, dressed molecules and itinerant ferromagnetism in ultracold Fermi gases_ , Reports on Progress in Physics 77, 034401 (2014).
  * Lan and Lobo [2014] Z. Lan and C. Lobo,  _A single impurity in an ideal atomic Fermi gas: current understanding and some open problems_ , J. Indian I. Sci. 94, 179 (2014).
  * Levinsen and Parish [2014] J. Levinsen and M. M. Parish,  _Strongly interacting two-dimensional Fermi gases_ , (2014), arXiv:1408.2737 .
  * Schmidt  _et al._ [2012] R. Schmidt, T. Enss, V. Pietilä, and E. Demler,  _Fermi polarons in two dimensions_ , Phys. Rev. A 85, 021602 (2012).
  * Côté  _et al._ [2002] R. Côté, V. Kharchenko, and M. D. Lukin,  _Mesoscopic Molecular Ions in Bose-Einstein Condensates_ , Phys. Rev. Lett. 89, 093001 (2002).
  * Massignan  _et al._ [2005] P. Massignan, C. J. Pethick, and H. Smith,  _Static properties of positive ions in atomic Bose-Einstein condensates_ , Phys. Rev. A 71, 023606 (2005).
  * Cucchietti and Timmermans [2006] F. M. Cucchietti and E. Timmermans,  _Strong-Coupling Polarons in Dilute Gas Bose-Einstein Condensates_ , Phys. Rev. Lett. 96, 210401 (2006).
  * Palzer  _et al._ [2009] S. Palzer, C. Zipkes, C. Sias, and M. Köhl,  _Quantum Transport through a Tonks-Girardeau Gas_ , Phys. Rev. Lett. 103, 150601 (2009).
  * Catani  _et al._ [2012] J. Catani, G. Lamporesi, D. Naik, M. Gring, M. Inguscio, F. Minardi, A. Kantian, and T. Giamarchi,  _Quantum dynamics of impurities in a one-dimensional Bose gas_ , Phys. Rev. A 85, 023623 (2012).
  * Spethmann  _et al._ [2012] N. Spethmann, F. Kindermann, S. John, C. Weber, D. Meschede, and A. Widera,  _Dynamics of Single Neutral Impurity Atoms Immersed in an Ultracold Gas_ , Phys. Rev. Lett. 109, 235301 (2012).
  * Rath and Schmidt [2013] S. P. Rath and R. Schmidt,  _Field-theoretical study of the Bose polaron_ , Phys. Rev. A 88, 053632 (2013).
  * Fukuhara  _et al._ [2013] T. Fukuhara, A. Kantian, M. Endres, M. Cheneau, P. Schauß, S. Hild, D. Bellem, U. Schollwöck, T. Giamarchi, C. Gross, I. Bloch, and S. Kuhr,  _Quantum dynamics of a mobile spin impurity_ , Nature Physics 9, 235 (2013).
  * Shashi  _et al._ [2014] A. Shashi, F. Grusdt, D. A. Abanin, and E. Demler,  _Radio-frequency spectroscopy of polarons in ultracold Bose gases_ , Phys. Rev. A 89, 053617 (2014).
  * Benjamin and Demler [2014] D. Benjamin and E. Demler,  _Variational polaron method for Bose-Bose mixtures_ , Phys. Rev. A 89, 033615 (2014).
  * Grusdt  _et al._ [2014a] F. Grusdt, A. Shashi, D. Abanin, and E. Demler,  _Bloch oscillations of bosonic lattice polarons_ , (2014a), arXiv:1410.1513 .
  * Grusdt  _et al._ [2014b] F. Grusdt, Y. E. Shchadilova, A. N. Rubtsov, and E. Demler,  _Renormalization group approach to the Fröhlich polaron model: application to impurity-BEC problem_ , (2014b), arXiv:1410.2203 .
  * Christensen  _et al._ [2015a] R. S. Christensen, J. Levinsen, and G. M. Bruun,  _Quasiparticle Properties of a Mobile Impurity in a Bose-Einstein Condensate_ , Phys. Rev. Lett. 115, 160401 (2015a).
  * Levinsen  _et al._ [2015] J. Levinsen, M. M. Parish, and G. M. Bruun,  _Impurity in a Bose-Einstein Condensate and the Efimov Effect_ , Phys. Rev. Lett. 115, 125302 (2015).
  * Ardila and Giorgini [2015] L. A. P. Ardila and S. Giorgini,  _Impurity in a Bose-Einstein condensate: Study of the attractive and repulsive branch using quantum Monte Carlo methods_ , Phys. Rev. A 92, 033612 (2015).
  * Volosniev  _et al._ [2015] A. G. Volosniev, H.-W. Hammer, and N. T. Zinner,  _Real-time dynamics of an impurity in an ideal Bose gas in a trap_ , Phys. Rev. A 92, 023623 (2015).
  * Grusdt and Demler [2016] F. Grusdt and E. Demler,  _New theoretical approaches to Bose polarons_ , arxiv 1510.04934 (2016).
  * Grusdt and Fleischhauer [2016] F. Grusdt and M. Fleischhauer,  _Tunable Polarons of Slow-Light Polaritons in a Two-Dimensional Bose-Einstein Condensate_ , Phys. Rev. Lett. 116, 053602 (2016).
  * Shchadilova  _et al._ [2016a] Y. E. Shchadilova, R. Schmidt, F. Grusdt, and E. Demler,  _Quantum Dynamics of Ultracold Bose Polarons_ , Phys. Rev. Lett. 117, 113002 (2016a).
  * Shchadilova  _et al._ [2016b] Y. E. Shchadilova, F. Grusdt, A. N. Rubtsov, and E. Demler,  _Polaronic mass renormalization of impurities in Bose-Einstein condensates: Correlated Gaussian-wave-function approach_ , Phys. Rev. A 93, 043606 (2016b).
  * Castelnovo  _et al._ [2016] C. Castelnovo, J.-S. Caux, and S. H. Simon,  _Driven impurity in an ultracold one-dimensional Bose gas with intermediate interaction strength_ , Phys. Rev. A 93, 013613 (2016).
  * Ardila and Giorgini [2016] L. A. P. Ardila and S. Giorgini,  _Bose polaron problem: Effect of mass imbalance on binding energy_ , Phys. Rev. A 94, 063640 (2016).
  * Robinson  _et al._ [2016] N. J. Robinson, J.-S. Caux, and R. M. Konik,  _Motion of a Distinguishable Impurity in the Bose Gas: Arrested Expansion Without a Lattice and Impurity Snaking_ , Phys. Rev. Lett. 116, 145302 (2016).
  * Jørgensen  _et al._ [2016] N. B. Jørgensen, L. Wacker, K. T. Skalmstang, M. M. Parish, J. Levinsen, R. S. Christensen, G. M. Bruun, and J. J. Arlt,  _Observation of Attractive and Repulsive Polarons in a Bose-Einstein Condensate_ , Phys. Rev. Lett. 117, 055302 (2016).
  * Hu  _et al._ [2016] M.-G. Hu, M. J. Van de Graaff, D. Kedar, J. P. Corson, E. A. Cornell, and D. S. Jin,  _Bose Polarons in the Strongly Interacting Regime_ , Phys. Rev. Lett. 117, 055301 (2016).
  * Rentrop  _et al._ [2016] T. Rentrop, A. Trautmann, F. A. Olivares, F. Jendrzejewski, A. Komnik, and M. K. Oberthaler,  _Observation of the Phononic Lamb Shift with a Synthetic Vacuum_ , Phys. Rev. X 6, 041041 (2016).
  * Gardiner and Zoller [2004] C. Gardiner and P. Zoller,  _Quantum Noise: A Handbook of Markovian and Non-Markovian Quantum Stochastic Methods with Applications to Quantum Optics_ , Springer Series in Synergetics (Springer, Berlin, 2004).
  * Breuer and Petruccione [2007] H. Breuer and F. Petruccione,  _The Theory of Open Quantum Systems_ (OUP, Oxford, 2007).
  * Schlosshauer [2007] M. Schlosshauer,  _Decoherence and the Quantum-To-Classical Transition_, The Frontiers Collection (Springer, 2007).
  * Schlosshauer [2005] M. Schlosshauer,  _Decoherence, the measurement problem, and interpretations of quantum mechanics_ , Rev. Mod. Phys. 76, 1267 (2005).
  * Zurek [2003] W. H. Zurek,  _Decoherence, einselection, and the quantum origins of the classical_ , Rev. Mod. Phys. 75, 715 (2003).
  * Caldeira and Leggett [1983a] A. Caldeira and A. Leggett,  _Path integral approach to quantum Brownian motion_ , Physica A: Statistical Mechanics and its Applications 121, 587 (1983a).
  * Caldeira and Leggett [1983b] A. Caldeira and A. Leggett,  _Quantum tunnelling in a dissipative system_ , Annals of Physics 149, 374 (1983b).
  * Efimkin  _et al._ [2016] D. K. Efimkin, J. Hofmann, and V. Galitski,  _Non-Markovian Quantum Friction of Bright Solitons in Superfluids_ , Phys. Rev. Lett. 116, 225301 (2016).
  * Hilary M. Hurst [2016] I. B. S. V. G. Hilary M. Hurst, Dmitry K. Efimkin,  _Kinetic theory of dark solitons with tunable friction_ , arxiv (2016).
  * Keser and Galitski [2016] A. C. Keser and V. Galitski,  _Analogue Stochastic Gravity in Strongly-Interacting Bose-Einstein Condensates_ , arXiv:1612.08980 (2016).
  * Bonart and Cugliandolo [2012] J. Bonart and L. F. Cugliandolo,  _From nonequilibrium quantum Brownian motion to impurity dynamics in one-dimensional quantum liquids_ , Phys. Rev. A 86, 023636 (2012).
  * Bonart and Cugliandolo [2013] J. Bonart and L. F. Cugliandolo,  _Effective potential and polaronic mass shift in a trapped dynamical impurity Luttinger liquid system_ , EPL (Europhysics Letters) 101, 16003 (2013).
  * Massignan  _et al._ [2015] P. Massignan, A. Lampo, J. Wehr, and M. Lewenstein,  _Quantum Brownian motion with inhomogeneous damping and diffusion_ , Phys. Rev. A 91, 033627 (2015).
  * Lampo  _et al._ [2016] A. Lampo, S. H. Lim, J. Wehr, P. Massignan, and M. Lewenstein,  _Lindblad model of quantum Brownian motion_ , Phys. Rev. A 94, 042123 (2016).
  * Pitaevskii and Stringari [2003] L. Pitaevskii and S. Stringari,  _Bose-Einstein Condensation_ (Oxford University Press, Oxford, 2003).
  * Canizares and Sols [1994] J. S. Canizares and F. Sols,  _Translational symmetry and microscopic preparation in oscillator models of quantum dissipation_ , Physica A: Statistical Mechanics and its Applications 212, 181 (1994).
  * Bruderer  _et al._ [2007] M. Bruderer, A. Klein, S. R. Clark, and D. Jaksch,  _Polaron physics in optical lattices_ , Phys. Rev. A 76, 011605 (2007).
  * Christensen  _et al._ [2015b] R. S. Christensen, J. Levinsen, and G. M. Bruun,  _Quasiparticle Properties of a Mobile Impurity in a Bose-Einstein Condensate_ , Phys. Rev. Lett. 115, 160401 (2015b).
  * Tempere  _et al._ [2009] J. Tempere, W. Casteels, M. Oberthaler, S. Knoop, E. Timmermans, and J. Devreese,  _Feynman path-integral treatment of the BEC-impurity polaron_ , Phys. Rev. B 80, 184504 (2009).
  * Casteels and Wouters [2014] W. Casteels and M. Wouters,  _Polaron formation in the vicinity of a narrow Feshbach resonance_ , Phys. Rev. A 90, 043602 (2014).
  * Grusdt  _et al._ [2017] F. Grusdt, G. E. Astrakharchik, and E. A. Demler,  _Bose polarons in ultracold atoms in one dimension: beyond the Fröhlich paradigm_ , arXiv:1704.02606 (2017).
  * Peotta  _et al._ [2013] S. Peotta, D. Rossini, M. Polini, F. Minardi, and R. Fazio,  _Quantum Breathing of an Impurity in a One-Dimensional Bath of Interacting Bosons_ , Phys. Rev. Lett. 110, 015302 (2013).
  * Ford and O’Connell [1991] G. W. Ford and R. F. O’Connell,  _Radiation reaction in electrodynamics and the elimination of runaway solutions_ , Physics Letters A 157, 217 (1991).
  * Wang and Zhan [2015] Q. Wang and H. Zhan,  _On different numerical inverse Laplace methods for solute transport problems_ , Advances in Water Resources 75, 80 (2015).
  * Guarnieri  _et al._ [2016] G. Guarnieri, C. Uchiyama, and B. Vacchini,  _Energy backflow and non-Markovian dynamics_ , Phys. Rev. A 93, 012118 (2016).
  * Lim  _et al._ [2017] S. H. Lim, J. Wehr, A. Lampo, M. A. García-March, and M. Lewenstein,  _On the Small Mass Limit of Quantum Brownian Motion with Inhomogeneous Damping and Diffusion_ , arXiv:1704.02606 (2017).
  * Nixon [1965] F. E. Nixon,  _Handbook of Laplace transformation: fundamentals, applications, tables, and examples_ (Prentice-Hall, 1965).
  * Feller [1971] W. Feller,  _An Introduction to Probability Theory and Its Applications_ (Wiley, 1971).
  * Grabert  _et al._ [1987] H. Grabert, P. Schramm, and G.-L. Ingold,  _Localization and anomalous diffusion of a damped quantum particle_ , Phys. Rev. Lett. 58, 1285 (1987).