Quantum dynamics of Bose-polaron in a d-dimensional Bose Einstein condensate

Quantum dynamics of Bose-polaron in a d-dimensional Bose Einstein condensate M. Miskeen Khan,1, 2, 3 H. Ter¸cas,1, 3 J. T. Mendon¸ca,1, 3 J. Wehr,4 C. Charalambous,5, 2 M. Lewenstein,2, 6 and M. A. Garcia-March2, 7 1Instituto Superior T´ecnico, Universidade de Lisboa, Portugal 2ICFO – Institut de Ci`encies Fot`oniques, The Barcelona Institute of Science and Technology, Av. Carl Friedrich Gauss 3, 08860 Castelldefels (Barcelona), Spain 3Instituto de Plasmas e Fus˜ao Nuclear, Instituto Superior T´ecnico, Universidade de Lisboa, Portugal 4Department of Mathematics and Program in Applied Mathematics University of Arizona Tucson, AZ 85721-0089 USA 5Instituto de F´ısica Interdisciplinar y Sistemas Complejos IFISC (CSIC-UIB), 07122 Palma de Mallorca, Spain 6ICREA, Pg. Llu´ıs Companys 23, 08010 Barcelona, Spain 7Instituto Universitario de Matem´atica Pura y Aplicada, Universitat Polit`ecnica de Val`encia, E-46022 Val`encia, Spain We study the quantum motion of an impurity atom immersed in a Bose Einstein condensate in arbitrary dimension. The Bogoliubov excitations of the Bose Einstein condensate act as a bosonic bath for the impurity. We present a detailed derivation of the d-dimensional Langevin equations that describe the quantum dynamics of the system, and of the associated generalized tensor that describes the spectral density in the full generality. When the impurity is not trapped, we calculate the mean square displacement, showing that the motion is super diffusive. We obtain also explicit expressions for the super diffusive coefficient in the small and large temperature limits. We find that, in the latter case, the maximal value of this coefficient is the same in all dimensions. We study also the behaviour of the average energy and compare the results for various dimensions. In the trapped case, we study squeezing and find that the stronger position squeezing can be obtained in lower dimensions. We quantify the non-Markovianity of the particle’s motion, and find that it increases with dimensionality. I. INTRODUCTION − The concept of a quasiparticle plays a fundamental role in physics, allowing to greatly simplify the description of numerous complex phenomena. A paradigmatic classical problem, in which quasiparticles appear, is the study of an electron interacting with a surrounding dielectric crystal. Its dynamics can be approximated by a much simpler dynamics of an electron with a different mass, called polaron, traveling through free space. This classical theory (see historical note in [1]) keeps inspiring new developments in physics. In particular, it plays an important role in the recent studies of the Bose polaron – the quasiparticle associated with an impurity immersed in a Bose-Einstein condensate (BEC). Bose polarons were investigated in diverse experiments on impurities immersed in bosonic gases. To begin with, the quantum dynamics of impurities in Bose gases were examined in [2,3], while technical aspects of experiments with Cs impurities were studied in [4]. The phononic Lamb shift in the context of ultracold bosons was observed in [5]. In addition, in these first experiments, charged, ionic or fixed impurities and their dynamics were studied: a quantum spin of a localized neutral impurity [6], fermions in a Bose gas [7,8], ions embedded in a BEC [9,10]. Quantum dynamics of spin impurities and fermions immersed in a Bose gas in an optical lattice were studied in Refs. [11,12]. More recent experiments define the state-of-the-art of the field: in [13] existence of a well-defined quasiparticle state of an impurity interacting with a BEC was demonstrated, while in [14] the strong interacting regime, which is natural for polaron problems, was investigated. In [15], a Bose polaron was studied near criticality, which provided important insights into the physics of quasiparticles in the vicinity of quantum critical points, that are otherwise much more difficult to study in other physical systems. The polaron theory was first developed in the strong coupling limit, and later extended to the intermediate and weakly interacting regimes. In the context of the Bose polaron problem, a large part of the theoretical effort deals with the weak regime, described by the so called Fr¨ohlich Hamiltonian. This theoretical approach studies effective mass, quantum dynamics, [16–22], collision dynamics [23], the behaviour in a d-dimensional BEC near the critical temperature [24], and related aspects of the system. Some studies in the weak regime considered the impurity as a quantum Brownian particle in a BEC or in a so called Luttinger liquid [25–29]. Importantly, Monte Carlo studies, in some instances beyond the regime of validity of the Fr¨ochlich Hamiltonian, allow to benchmark the aforementioned theoretical results [30–33]. Other works focused on the intermediate and strong coupling regimes [34–42], and on the non-zero temperature systems [43–46]. Yet other works studied the quenched dynamics and orthogonality catastrophe using the multi-configuration time-dependent Hartree method, both in weak and strong interaction regimes [47–49]. Also, several papers investigated ejection and injection spectroscopy related to Bose polarons and the orthogonality catastrophe [50,51], as well as bound states [52–55]. A number of notable works study two polarons immersed in a BEC (Bose bi-polaron) [56,57] and the problem of an impurity in a arXiv:2007.04925v1 [quant-ph] 9 Jul 2020 2 two-component BEC [58,59]. There has recently been a renewed interest in the polaron problem in mathematical physics literature. See in particular: [60–62]. Finally a series of papers deal with applications of Bose polarons in quantum thermometry [63–66] and thermodynamics [67–69]. In the present paper we approach the Bose polaron problem from the quantum open systems perspective. This framework has been already used to understand several questions related to Bose polarons [25–29,57,59, 65,69]. The study of the Bose polaron from this perspective leads to quantum stochastic equations with inhomogeneous damping and multiplicative noise. This is not always possible in the case of non-Ohmic spectral densities [70–72], presenting a challenge for mathematical physics. Particular care is required in higher dimensions, where linearization of the inhomogeneous damping and diffusion is more difficult to control. In this paper, our first goal is to offer a detailed derivation of the quantum Langevin equations and the associated generalized d-dimensional spectral density, which describe from an open quantum system perspective the dynamics of an impurity immersed in a BEC. In one dimension, this was done in detail in [25], and some hints were given for two and three dimensions. Performing the calculation in two and three dimensions offers information, valuable for many experiments. When the impurity is not trapped, we calculate (in all dimensions) the mean square displacement and the average energy, together with analytic expressions valid in different limits, and discuss the effects of dimensionality on the dynamics. For the trapped impurity, we determine the stationary state, and analyze the position squeezing effect by employing the covariance matrix. We finally study the non-Markovian character of this dynamics in all dimensions. We perform our calculations with experimentally feasible parameters, making sure that the parameters used lie within the regimes of validity of this model. An important contribution of the paper is that it offers a comprehensive understanding of the quantum dynamics and non-Markovianity in different dimensions and experimentally realistic situations, not discussed in previous investigations. The paper is organized as follows: in section II we present the assumptions and derivations that permit to obtain the (linearized) quantum Brownian motion Hamiltonian (cf. Eq. (20)) from the initial second quantized one. In section III we obtain the generalized ddimensional spectral density. We study, in all dimensions, the quantum dynamics in the non-trapped and trapped case in section IV, and non-Markovianity in section V. We finally conclude in section VI. In appendices Aand Bwe include the detailed derivation of the generalized d-dimensional spectral density and vectorial quantum Langevin equations, as they are per se important results of this paper. In appendix Cwe present the expressions for the position and momentum variances of the generalized Langevin equations. Finally, in appendix D we discuss the the validity of our approximations. II. HAMILTONIAN AND BOGOLIUBOV MODES− We start by considering an impurity atom with mass mIimmersed in a d-dimensional ultracold gas of N bosons. The interaction between the bosons occurs through the scattering potential VB(r). We denote by Ψ(r) (Ψ†(r)) the annihilation (creation) field operator of the atoms at the position r, which fulfills canonical bosonic commutation relations [Ψ(r),Ψ†(r)] = δ(r−r0). The bosonic density therefore takes the form nB= Ψ†(r)Ψ(r). The total Hamiltonian is given by H=HI+HB+HBB +HIB.(1) Here, the four terms represent the Hamiltonians of the impurity being kept in an external potential Uext(r), bosons in a potential Vext(r), the boson-boson atomic interaction and the impurity-boson atomic interaction, respectively. Within the second quantization formalism, their explicit forms are [25] HI=P2 2mI +Uext(r),(2) HB=ZddrΨ†(r)P2 B 2mB +Vext(r)Ψ(r) =X k ka† kak,(3) HBB =gBZddrΨ†(r)Ψ†(r)Ψ(r)Ψ(r) =1 2VX k,k0,q VB(q)a† k0−qa† k+qak0ak,(4) HIB =gIBnB=1 VX k,q VIBρI(q)a† k−qak.(5) In the above expressions, Vext(r) denotes the external potential experienced by the Bosons which are contained in a (box of) volume Vof the hyperspace. From now on, we assume a homogenous BEC, that is, Vext(r) = 0 along the direction of the impurity motion. For the impurity, the external potential is Uext(r), and we will study two cases: a free or a parabolically trapped impurity. The bosonic operators ak(a† k) destroy (create) a boson of mass mB having wave vector kand energy k= (~k)2/(2mB)−µ, measured from its chemical potential µ. In addition, the quantities VBand VIB represent the Fourier transform Fq[.] of the impulsive (contact) boson-boson and impurity-boson interactions respectively. Their explicit expressions are: VB(q) = gBFq[δ(r−r0)],(6) VIB(q) = gIBFq[δ(r−r0)].(7) Here, the respective coupling strengths are gBand gIB. They are mainly determined by the corresponding scattering lengths and densities [73,74] and their explicit 3 expressions will be given later. We assume that the impurity density is low enough, which allows us to neglect the terms describing the interaction between impurities. The (dimensionless) density of the impurity in the momentum space is given by ρI(q) = Z∞ −∞ dr0e−iq.r0δ(r0−r).(8) Next, for the sake of completeness, we review how to construct the Fr¨ohlich Hamiltonian, which describes the linear interaction between the motional position quadrature of the impurity and the Bogoliubov bosonic modes of BEC. The goal is to show that such linear regime allows us to model the impurity as a quantum Brownian particle which experiences an effective environment formed by the Bogoliubov Bosonic modes of the BEC (as derived in [25]). Given that the Hamiltonian of the Bosonic interaction is not in bilinear form, we linearize it and replace the creation and annihilation operators by their average values pN0. Below a critical temperature, the atoms mainly occupy the ground state forming a BEC, however, we neglect terms proportional to Nk(k6= 0) i.e. the number of particles out of the ground state. In order to diagonalise the bath modes, we further apply the following Bogoliubov transformation ak=ukbk−vkb† −k, a−k=ukb−k−vkb† k.(9) The transformation coefficients are u2 k=1 2k+n0VB Ek + 1,(10) v2 k=1 2k+n0VB Ek−1,(11) where n0is the (constant) density of particles in the ground state of the homogeneous gas and the Bogoliubov energy spectrum is given by Ek=~ωk=~c|k|r1 + 1 2(ξk)2,(12) with ξ=~ √2gBmBn0 , c =~ √2mBξ,(13) representing the coherence length and the speed of sound, respectively. The effective bath Hamiltonian under such transformation reads [73] HB+HBB =X k6=0 Ekb† kbk.(14) Here, we have neglected the non-operator terms which simply shift the energy level of the atoms in BEC. We approximate the bosons-impurity interaction in a similar way. Further, we only keep the terms proportional to pN0, where the macroscopic occupation of the condensate holds as expressed by the condition Ni6=0 N0. By discarding the terms which might cause non-physical instabilities and those bilinear in pN0, we obtain the Hamiltonian HIB =n0VIB +rn0 VX k6=0 ρ(k)VIB(ak+a† k).(15) The first term represents simply the constant mean field energy and provides the shift of the energy of polaron. For the purposes here, it can be neglected. By further invoking the transformation from Eq. (9) into Eq. (15), one gets HIB =rn0 VX k6=0 ρ(k)VIB (uk−vk) (bk+b† −k) =rn0 VX k6=0 ρ(k)VIBrk Ek (bk+b† −k),(16) where once again we have discarded the non-operator terms. Since the density is dependent on the position of the impurity, we insert its expression into Eq. (16), which results in the interaction between impurity position and bath variables given by HIB =X k6=0 Vkeik·r(bk+b† −k).(17) Importantly, Vkcontains the impurity-Boson coupling coefficient, and takes the form Vk=gIBrn0 V(ξk)2 (ξk)2+ 21 4 .(18) The interaction in the Eq. (17) is the interaction part of the Fr¨ohlich Hamiltonian. Under the assumption that one restricts the calculation to the limit k·r1, the interaction reads HIB =X k6=0 Vk(1 + ik·r) (bk+b† −k).(19) We further simplify it by redefining the Bogoliubov modes operator bk→bk−vk/Ek1, to absorb terms proportional to identity operator. After all these simplifications, the final form of the Hamiltonian of the impurity in a BEC reads H=HI+X k6=0 Ekb† kbk+X k6=0 ~gk·rπk,(20) with gk=kVk/~, πk=ibk−b† k.(21) The Hamiltonian in Eq. (20) describes a linear interaction between the impurity center of mass motion and a 4 bath of the Bogoliubov modes of a BEC. It thus has a form of the QBM Hamiltonian, in which the impurity plays the role of a Brownian particle while the modes of BEC act as an effective Bosonic environment as represented by its (dimensionless) momenta πk. III. d-DIMENSIONAL SPECTRAL DENSITY− The Hamiltonian derived in the previous section allows us to study the quantum dynamics of an impurity, taking advantage of the analogy with the QBM model. To characterise the bath, we write its self-correlation function as C(τ) = X k6=0 ~gkhπk(τ)πk(0)i.(22) Here gk=gkgkTis the coupling tensor. The environment is made of bosons whose state at finite temperature Tfollows the Bose-Einstein statistics. Therefore, the mean number of bosons in each of the modes reads hb† kbki=1 exp(~ωk/kBT)−1.(23) In order to calculate the correlation, we invoke the expression for the dimensionless momenta and make use of Eq. (23) and Eq. (22) which results in C(τ) = X k6=0 ~gkcoth ~ωk 2kBTcos(ωkτ)−isin(ωkτ) ≡ν(τ)−iλ(τ),(24) where the real and imaginary part of the self-correlation function are given by ν(τ) = Z∞ 0 J(ω) coth ~ω 2kBTcos(ωτ)dω, (25) λ(τ) = Z∞ 0 J(ω) sin(ωτ)dω =−mI˙ Γ(τ).(26) Moreover, the damping kernel Γ(t) can be obtained from Γ(t) = (1/mI)Z∞ 0 dω(1/ω)J(ω) cos(ωt).(27) In the above expressions we have introduced the spectral density J(ω), which fully characterises the effects of the bath on the system. This information is contained in the coupling strengths of the various modes of the bath with the system. The spectral density is defined as J(ω) = X k6=0 ~gkδ(ω−ωk).(28) In the present case, the couplings of the impurity (system) and bosons (bath) interaction can be derived from first principles. It is therefore possible to obtain the exact expression for the spectral density. This scenario is in contrast to various complicated system-bath interactions, where it is hard to get an exact form of the spectral densities, such as in bulk mechanical structure akin to opto-mechanical setup [75–77]. While the case of spectral density in 1dhas been studied in [25], here we derive it systematically in d={1,2,3}dimensions of the quasi momentum space. In appendix A, we derive the expression for the spectral density tensor, which is given by Jd(ω) = d−1[Jd(ω)] Id×d,(29) where Id×dis the identity matrix and the scalar function Jd(ω) in d dimensions is given by Jd(ω) = Sd√2d(ηd)2(Λd)d+2 (2π)d!×        " mB [gB,d][d d+2 ]n0,d!qω2 (Λd)2+ 1 −1#(d+2 2) qω2 (Λd)2+ 1        .(30) For d= 1,2 and 3 we have S1= 2, S2= 2πand S3= 4πrespectively. Moreover, we have defined the d−dependent characteristic frequency Λd= (gB,dn0,d)/~ because the boson-boson coupling and the density differ in various dimensions. We also write the impurityboson coupling in the units of the boson-boson coupling as ηd= (gIB,d/gB,d). Such characterization allows us to study the long-time dynamics of the impurity in the following sense: one can identify two opposite limits in the above expression i.e. ωΛdand ωΛdin which Λd appears naturally as the characteristic cut-off frequency which distinguishes between the low and the high frequencies of the bath. The low-frequency behaviour is attributed to the linear part of the Bogoliubov spectrum [25]. From the Tauberian theorem [78], one can obtain the long-time behaviour of a function which is determined by the low frequency response of its Laplace transform. The above low-frequency choice is therefore a natural way of studying the dynamics perturbed by the bath that acts beyond the very short transient regime. Note that, for d= 1 the above expression reduces to the one dimensional spectral density used in [25]. To the lowest order of ω/Λd, the expression for the spectral density with the low frequency response of the bath is given by Jd(ω)≃Sd(ηd)2 2(2π)d mB [gB,d][d d+2 ]n0,d!(d+2 2) ×ωd+2.(31) This expression gives the scaling of the frequency for the spectral density function in all dimensions. We point out that due to the spherical symmetry of the bath, the spectral density tensor is a diagonal matrix. As a consequence, the noise and damping kernels given by Eq. (25) 5 and Eq. (27) are also diagonal. We now give further details about the other parameters involved. In dimension d, the coupling constant gB,dand boson density n0,dhave the form gB,d=Sd~2a3 mBp~/mBωd3−d, n0,d= (n0,1)d,(32) their units being J ·mdand m−drespectively. Here we have written these expressions in terms of the threedimensional scattering length a3and one-dimensional density n0,1. We have further assumed transverse confinement of the boson gas with a harmonic trap having a Gaussian ground state [79], which makes the cases d < 3 to be the quasi oneand two-dimensional. We emphasize that the parabolic potential is introduced only in the direction transverse to the direction under investigation. The dynamics we study is thus still confined to a box potential, making the homogeneity of Boson gas to be a valid approximation. Moreover, the zero point fluctuations of the condensate are characterised by the trapped frequencies ωd={ω1=ω⊥, ω2=ωz, ω3= 0}. For instance, when we consider one, twoand three-dimensional condensate to be confined in the xdirection, in the x−y plane or in the volume x−y−zrespectively, the explicit form of the potential may be given by Vext(r) =      (1/2)mBω2 ⊥y2+z2,for d= 1 (1/2)mBω2 zz2,for d= 2 0,for d= 3. (33) Note that for d= 3 there is no parabolic confinement and therefore the expressions are independent of the trapping frequency. It is then possible to define a characteristic time τdwhich is raised to the power din the expression for the spectral function, which is given by Jd(ω) = mI(τd)dωd+2 where (τd)d≡Sd(ηd)2 2(2π)dmI mB [gB,d][d d+2 ]n0,d!(d+2 2) .(34) It is evident that the spectral density has a super-ohmic dependence on the frequency in all dimensions. It is therefore expected that the bosonic bath would induce a non-Markovian dynamics of the impurity. Moreover, it can be shown that the increasing nature of the spectral density makes certain quantities, such as momentum dispersion, diverge. It is therefore customary to define the ultraviolet cut-off K(ω, Λd) in order to suppress the contribution of high frequencies. After this, the expression for the spectral density reads Jd(ω) = mI(τd)dωd+2K(ω, Λd).(35) In the following we study the impurity dynamics, varying the dimension and d-dependent cut-off function in the expression of the spectral density. IV. DYNAMICS AND CONTROL − To study the quantum dynamics of the impurity atom, we write down the equations of motion in the Heisenberg picture. The impurity, which is immersed in a bath of dimension d, is further trapped by a harmonic potential. In dimensions 1, 2 and 3, the potential is Uext(x) = (1/2)mIΩ2x2,Uext(x, y) = (1/2)mIΩ2x2+y2and Uext(x, y, z) = (1/2)mIΩ2x2+y2+z2respectively. Here, we have assumed equal trapping frequency in all the directions available to the impurity dynamics. The free QBM is therefore characterised by setting Ω = 0 in all the cases. We write the equation of motion in vectorial form as ˙ X(t) = i ~[H, X(t)] = ˙ P(t) mI ,(36) ˙ P(t) = i ~[H, P(t)] = −mIΩ2˙ X(t)−~X k gkπk(t),(37) ˙ bk(t) = i ~[H, bk(t)] = −iωkbk(t)−gkTX(t),(38) ˙ b† k(t) = i ~hH, b† k(t)i=−iωkb† k(t)−gkTX(t).(39) Here Hrepresents the Hamiltonian of the system given by Eq. (20). In general, the dimension of the vectors in the above equations is d, the dimension of the bath. In appendix B, we combined these equations to obtain an equation of motion for the impurity position vector: ¨ X(t)+Ω2X(t) + ∂tZt 0 Γ(t−s)X(s)ds = (1/mI)B(t).(40) Here, the quantum Brownian stochastic force B(t) stands for B(t) = X k i~gk(b† keiωkt−bke−iωkt).(41) In any given dimension, the diagonal damping kernel when equal weighting for all the directions is taken, will suffice to study the motion along any one of the coordinate axis. However for different dimensions, the spectral density will bring different level of super-ohmicity, as stated by Eq. (35). Therefore, the form of the noise and damping kernels will also differ according to the dimension involved. As a result, the impurity motion is different for different dimensions, despite being studied along one particular coordinate axis. One can then aim to study such a motion by constructing a unit vector (1,0,0) (i.e. along the x direction) and taking dot product with Eq. (40). This results in ¨x(t)+Ω2x(t) + ∂tZt 0 Γxx d(t−s)x(s)ds =1 mIBx(t). (42) 6 Note that the tensor components satisfy Γxy d= Γxz d= 0. Additionally, from the structure of the integral term in the Eq. (42), it is obvious that damping kernel is non-local in time. This implies that the dynamics of the impurity depends on its history. Therefore, in general the impurity motion displays memory effects. Only in the case of Ohmic spectral density (linear function of ω) the memory damping kernel reduces to a Dirac delta function and describes the time-local dynamics of the standard damped harmonic oscillator. The time local behaviour is violated in similar experimental configuration [77] and surge of non-Markovianity is addressed elsewhere [80–82]. The formal solution to Langevin-like Eq. (42) takes the form x(t) =G1,d(t)x(0) + G2,d(t) ˙x(0) + (1/mI)Zt 0 ds ×G2,d(t−s)Bx(s),(43) where the Green’s functions G1,dand G2,dare defined in terms of their Laplace transforms LS,d[G1,d(t)] = S S2+ Ω2+SLS,d[Γxx d(t)],(44) LS,d[G2,d(t)] = 1 S2+ Ω2+SLS,d[Γxx d(t)].(45) Moreover, they satisfy the following initial conditions G1,d(0) = 1,˙ G1,d(0) = 0,(46) G2,d(0) = 0,˙ G2,d(0) = 1.(47) From now on, the dynamics in the case when high frequencies are suppressed, will be analyzed by introducing a sharp cutoff. The latter is given by K= Θ(Λd−ω), where Θ is the Heaviside step function. We further use Eq. (27) to compute the damping kernel in any dimension d LS,d[Γxx d(t)] = (Λd)d+2 (τd)d2F11,d+2 2;d+4 2;−(Λd)2 S2 d(d+ 2)S, (48) with 2F1[.] denoting the hypergeometric function. A. Untrapped Case Let us first study the free QBM, that is the untrapped case Ω = 0. The quantities of interest in this case are the mean squared displacement MSDd(t) (see definition below, Eq. (53)) and the average kinetic energy Ed(t) of the impurity. The motion is fully characterised by the functions G1(t) and G2(t) which are the inverse Laplace transform of Eq. (44) and Eq. (45) respectively. Exact � �� �� �� �� ��� ��� ������ ������ ������ ������ ������ ������ ������ ������ Figure 1. (Color online) Dynamics of the propagator G2(t) in the untrapped case for various dimensions. The dark dashed lines show the asymptotic behaviour given by Eq. (51), while the light solid lines represent the solution obtained by employing the Zakian method. In the long time limit, both solutions match for all dimensions. The results refer to an impurity K with mass mI= 6.4924249 ×10−26kg, immersed in a gas of Rb with mass mB= 1.4192261 ×10−25kg. The one-dimensional boson density is n0,1= 7(µm)−1. We fix the three-dimensional scattering length a3= 100 a0, where a0is the Bohr radius. Here, the time axis is scaled with the one-dimensional characteristic frequency ω0≡(~n2 0,1/mI). analytical expressions for these functions are hard to obtain. Note however, that the Laplace transform of both of these functions are expressed in terms of the Laplace transform of the damping kernel. In the regime of interest |S|  Λd, which characterises the low frequency response, we have approximately LS,d[Γxx d(t)] = d−2(Λdτd)dS.(49) We therefore obtain the asymptotic expressions for the Laplace transforms of the position and momentum propagators LS,d[G1,d(t)] = 1 αdSand LS[G2,d(t)] = 1 αdS2, (50) where αd= 1 + d−2(Λdτd)d. Their time domain representations read G1,d(t)=1/αd, G2,d(t) = t/(αd).(51) We note that such expressions do not satisfy the boundary conditions stated in Eq. (46) and Eq. (47). However, this is justified since the above solution refers to long-time behaviour. Several algorithms exist for the numerical computation of the inverse Laplace transform of an arbitrary function. Here we employ the Zakian method [83] to compute the inverse Laplace transform of Eq. (44) and Eq. (45). This method approximates the inverse Laplace transform f(t) of a function F(S) through f(t)≃2 t N X j=1 <[kjF(Ξj/t)] ,(52) 7 with the values of the complex parameters kjand Ξj given in Ref. [83]. In order to check the equivalence between the asymptotic form of G2(t) and its Zakian approximation, we plot both in Fig. 1. It turns out that they agree in the long time limit. From here on we will be employing them interchangeably according to our computational convenience. It is evident from Eq. (43) that the function G2(t) is responsible for the propagation of the initial velocity of the impurity. From the results shown in Fig. 1, such a function follows a ballistic profile in any dimension, i.e. it is a linear function of time. 1. Mean Square Displacement In this section we discuss the mean squared displacement (MSD) of the impurity motion which is a measurable quantity in cold-atom experiments [2]. The MSD is defined as MSDx,d(t) = D[x(t)−x(0)]2E.(53) The expression of the MSD for the generalised Heisenberg-Langevin equations is given in the Appendix C. For the system at hand, we employ the asymptotic expressions of the Green’s functions to evaluate the MSD in different dimensions. Its dynamical part is given by MSDx,d(t) = t αd2 h˙x(0)2i+1 2(αdmI)2Zt 0 du Zt 0 dv ×(t−u)(t−v)h{Bx(u), Bx(v)}i.(54) Additionally, using the diagonal form of the noise tensor given in Eq. (25), one can obtain the relation between the correlation characterised by the positive commutator h{Bx(u), Bx(v)}i and the component νxx(t) of the noise kernel (fluctuation-dissipation relation) [27] . Explicitly, h{Bx(u), Bx(v)}i = 2~νxx(u−v).(55) Substituting the d-dimensional spectral function with a sharp cutoff into the noise component we get MSDx,d(t) = t αd2 h˙x(0)2i+~d−1(τd)d mI(αd)2Zt 0 du Zt 0 dv ZΛd 0 dω(t−u)(t−v) coth ~ω 2kBTcos [ω(u−v)]ωd+2. (56) By performing the two-dimensional integration over the time variables uand v, followed by an integration over the variable ω, we evaluate the expression for the low temperature regime, where coth (~ω/2kBT)→1 holds. In the long time limit, the resulting expression for the MSD is dominated by the terms proportional to t2and its explicit expression turns out to be MSDLT x,d(t) = "h˙x(0)2i+~(τd)d(Λd)d+1 mId(d+ 1) #t αd2 (57) Figure 2. (Color online) High temperature super diffusion coefficient in various dimensions as a function of the coupling strength. The solid, dashed and dotted curves represent the d= 1, d= 2 and d= 3 cases respectively. We set temperature T= 0.15µK, which fulfills the high temperature condition kBT⩾Max [~Λd] (see text). The rest of the parameters are the same as in Fig. 1. The vertical solid, dashed and dotted lines fix the critical coupling for the Fr¨ohlich Hamiltonian to be valid in one, two and three dimensions, respectively. In the regime which fulfills the conditions stated above, the MSD is proportional to the square of the time for all dimensions. In the normal diffusion scenario, the MSD shows a linear dependence on time. If, on the contrary, the MSD is non-linear in time, proportional to tα with an exponent higher than one, the diffusion is called anomalous and the motion is called superdiffusive. In the present case, superdiffusion is a consequence of the super-ohmic spectral density in every dimension. The coefficient in the second term is called the superdiffusion coefficient Dx,dand can be interpreted as the average of the square of the speed with which the impurity runs away. We thus have DLT x,d=~(τd)d(Λd)d+1 mId(d+ 1) (αd)2.(58) One can perform a similar analysis of the high temperature regime, which is followed by the approximation coth (~ω/2kBT)→(2kBT/~ω). We remark that the condition kBT⩾Max [~Λd] implies the high temperature regime in any dimension. Here Max[.] is the maximum of the cut-off frequencies of different dimensions. This means that all the Bogoliubov modes of the bath will be thermally populated in any of the considered dimensions. However, while the cut-off frequency in dimension dscales as Λd∼(n0,1)d, it also depends on the boson coupling constant in the corresponding dimension and therefore on the transverse confinement of the boson gas [cf. Eq. (32)]. Inserting the values of the parameters used in this article, we obtain Λ2>Λ3>Λ1. The high temperature regime holds as long as kBT⩾~Λ2. In this regime, the MSD again scales with the square of the time. The dimension-dependent superdiffusion coefficient takes 8 the form DHT x,d=2kBT(τd)d(Λd)d mId2(αd)2.(59) It is clear from this expression that the superdiffusion coefficient is proportional to the temperature of bath and inversely proportional to the mass of the impurity. One can further write the high temperature superdiffusion coefficient as an explicit function of the coupling parameter: DHT x,d=2kBT mI"βdη2 (βd)2d−2η4+ 2βdη2+d2#.(60) Here, we have defined the quantities βd≡(Λdτd,s)dwhere (τd,s)d≡η−2(τd)d.(61) In dimension d, the maximum of this function occurs at ηmax,d=d q(Λdτd,s)d.(62) and has the same maximal value in all dimensions: DHT,max x,d=1 mIkBT 2.(63) In Fig. 2, we plot the high temperature superdiffusion coefficient for a range of coupling strengths, covering the allowed critical coupling in every dimension (see appendix Dfor the validity of the Fr¨ohlich Hamiltonian). In all dimensions, for sufficiently weak coupling strengths, one observes a corresponding increase of the bath-induced momentum diffusions of the impurity as the coupling strength increases. However, above ηmax,d, the super diffusion coefficient is reduced as the coupling grows. Such a damped regime occurs only in the one-dimensional case within the Fr¨ohlich regime, where the impurity may exert both underdamped and overdamped motion. (In higher dimensions the over damped regime occurs past the vertical line that signals the value of the coupling, critical for the validity of the Fr¨ohlich Hamiltonian.) We finally remark that the condition for the occurrence of both of these characteristic motions within the Fr¨ohlich regime turns out to be ηmax,d< ηc,d. The latter are the critical couplings of Fr¨ohlich regime given in the appendix Dand correspond to the vertical lines shown in the Fig. 2. 2. Average Energy We now turn to the average kinetic energy E(t) of the impurity, corresponding to the observed coordinate. This can be computed from the variance of the corresponding momentum operator reading Ex,d(t) = hp2 x,d(t)i 2mI (64) The generalized expression for the variance of the momentum is given in the appendix C. Using the dimensiondependent asymptotic expressions for G1(t) and G2(t) we obtain Ex,d(t) = hp2 x,d(0)i 2mI(αd)2+~ 2mI(αd)2Zt 0 du Zt 0 dv νxx d(u−v). (65) For any arbitrary temperature, it is difficult to obtain from here an analytic expression. Here we are mainly interested in the ultracold regime. This means that all the bath modes are now in a collective vacuum state. Therefore, in the zero-temperature limit, the above expression further reduces to ELT x,d(t) = hp2(0)i 2mI(αd)2+~(Λd)d+1 (τd)d (αd)2d(d+ 1) + ~(Λd)d+1 (τd)d (αd)2d(d+ 1) 1F2d 2+1 2;1 2,d 2+3 2;−1 4(Λd)2t2. (66) Here, the first term represents the initial mean energy of the impurity determined by its initial momentum variance. Additionally, there is a rescaling of the mass of the impurity due to the interaction with the bath. The additional mass term depends on the dimensionality of the bath through Λdand τd. The second term is the steady state mean energy of the impurity which is determined by the impurity-bath coupling and density of the bath, again through the same parameters. The last term of the expression contains information about the energy variation in time. We plot the energy function in Fig. 3a-c) for different dimensions. In all dimensions, the energy oscillates in time. This clearly shows the energy exchange between the system and the bath. Where the energy increases, it is due to an energy absorbed from the bath. The back flow of energy is a manifestation of memory effects in the QBM [84]. Moreover, deep inside the weak coupling regime, the bath perturbs the system more strongly as the coupling strength increases. In any dimension, this results in the higher initial increasing peak for a larger coupling constant. The overall profiles of all the energy functions tend to approach their asymptotic steady state values. A similar analysis can be performed for the high temperature case, as was done in the previous section for the MSD. We treat the problem classically, meaning that the the symmetrised noise correlation function (cf. Eq. (55) and Appendix C) would act as the classical analogue in the present quantum formulation [85]. As for the dynamical part of the energy expression, the classical regime is further obtained by requiring the conditions t→ ∞, ~→0 and kBT~ω. In these limits, the energy Ecl,ss x,d becomes Ecl,ss x,d= 2kBT"βdη2 (βd)2d−2η4+ 2βdη2+d2#.(67) 9 � �� �� �� �� ��� ��� ���� ���� ���� ���� ���� �) � �� �� �� �� ��� ��� ���� ���� ���� ���� ���� ���� ���� �) � �� �� �� �� ��� ��� ���� ���� ���� ���� ���� ���� ���� �) Figure 3. (a-c): Dynamics of the average energy for 1d, 2dand 3d, respectively. The coupling strengths are shown in the legend for every case. The rest of the parameters are kept same as before. The above expression is once again maximised at ηmax,d giving an upper bound to the kinetic energy of the impurity reading Ecl,ss x,d,max =kBT/2.(68) Remarkably, this is the familiar equipartition theorem that holds in any dimension. It follows from these results that ηmax,dis the value of the system-bath coupling at which the impurity reaches thermal equilibrium with the Bogoliubov bath. B. Trapped Case In recent years, there has been an increased interest on trapped impurities within cold atomic media. For instance, the bound states of the trapped impurities provide a platform to test the existence of synthetic vacuum of the hosting medium by witnessing the induced Lambshifts [5]. Additionally, trapped impurities in BECs can serve as highly controlled phononic q-bits [7,86]. Hence, theoretical study of a trapped impurity as an open quantum system, as in this work, can be valuable in all of the aforementioned cases. In the trapped case, we confine the impurity into a harmonic trap with frequency Ω. We compute the functions G1(t) and G2(t) by employing the Zakian method. Their time dynamics is shown in Fig. 4(a-b). In all dimensions, both of these functions oscillate out of phase by π/2. This reflects the fact that position and momentum are the two quadratures of the impurity motion. Note that the information of the initial position and the momentum variances is carried by the functions G1(t) and G2(t) respectively (see appendix C). The decay of these functions provides insights into the system dynamics. First of all, such decay shows that the impurity dynamics is stable. In general, the stability analysis of the dynamics can be performed more rigorously, e.g. through the Routh–Hurwitz stability criterion [87]. However, given the absence of the analytical form of G1(t) and G2(t) (or their Laplace transforms), we rely on a numerical evaluation of their profiles. In fact, both of these functions approach to zero as t→ ∞. In effect, the system dynamics becomes independent of its initial conditions and its behaviour is completely determined by � �� �� �� �� ��� ��� ��� -��� ��� ��� ��� �) � �� �� �� �� ��� ��� ��� -��� ��� ��� �) Figure 4. (a-b) Dynamics of G1(t) and G2(t) in the trapped case for 1d, 2dand 3dwith the corresponding dimensions and couplings shown in the legend. In all the cases we have set Ω=4π×500 Hz. Rest of the parameters are same as in Fig. 1. the bath state. It can be seen that in the long time limit, each one of them collapses to a single curve for all the cases displayed in Fig. 4. On the contrary at initial short times, their amplitudes and phases are mismatched for different initial conditions. The differences coming from different coupling parameters and dimensions also vanish in the long time limit. This leads to the insight that the steady state regime is completely determined by the bath state and the system tends to equilibrate with the local state of the bath. We now turn to the study of the steady state dynamics of the impurity. Since the input Bogoliubov bath modes 16 In the above expression we have written the final equality in terms of the three dimensional scattering length a3 based on the harmonic confinement of the condensate in the transverse direction (see main text). From the second-to-last inequality one can get the expression for d= 3 gIB,3.2πc3(ξ3)2(D2) (where we put ~= 1). The same bound was reported in [94]. On the other hand, for d= 1, the last equality leads to the following scaling for the bound on the critical coupling: ηc,d∼√n0,1a1; where a1= (~/mBω⊥)/a3,(D3) stated also in [25,32]. Typically for a boson gas made of Rb87 atoms, the scattering length a3= 100a0with a0 denoting the Bohr radius [94]. Moreover, let the transverse frequencies be {ω⊥, ωz}= 2π×34 kHz as in the optical lattice [2]. This implies the following numerical bounds on the coupling, depending on the dimension: ηc,d∼     3.7,for d= 1, 4.4,for d= 2, 9.8,for d= 3, (D4) which are the values represented as vertical lines in the figures in Figs. 2and 6. [1] M. I. Dykman and E. I. Rashba, Physics Today 68, 10 (2015). [2] J. Catani, G. Lamporesi, D. Naik, M. Gring, M. Inguscio, F. Minardi, A. Kantian, and T. Giamarchi, Phys. Rev. A85, 023623 (2012). [3] N. Spethmann, F. Kindermann, S. John, C. Weber, D. Meschede, and A. Widera, Phys. Rev. Lett. 109, 235301 (2012). [4] M. Hohmann, F. Kindermann, B. G¨anger, T. Lausch, D. Mayer, F. Schmidt, and A. Widera, EPJ Quantum Technology 2, 23 (2015). [5] T. Rentrop, A. Trautmann, F. A. Olivares, F. Jendrzejewski, A. Komnik, and M. K. Oberthaler, Phys. Rev. X6, 041041 (2016). [6] F. Schmidt, D. Mayer, Q. Bouton, D. Adam, T. Lausch, N. Spethmann, and A. Widera, Phys. Rev. Lett. 121, 130403 (2018). [7] R. Scelle, T. Rentrop, A. Trautmann, T. Schuster, and M. K. Oberthaler, Phys. Rev. Lett. 111, 070401 (2013). [8] J. B. Balewski, A. T. Krupp, A. Gaj, D. Peter, H. P. B¨uchler, R. L¨ow, S. Hofferberth, and T. Pfau, Nature 502, 664 (2013). [9] C. Zipkes, S. Palzer, C. Sias, and M. K¨ohl, Nature 464, 388 (2010). [10] S. Schmid, A. H¨arter, and J. H. Denschlag, Phys. Rev. Lett. 105, 133202 (2010). [11] T. Fukuhara, A. Kantian, M. Endres, M. Cheneau, P. Schauß, S. Hild, D. Bellem, U. Schollw¨ock, T. Giamarchi, C. Gross, I. Bloch, and S. Kuhr, Nature Physics 9, 235 (2013). [12] S. Ospelkaus, C. Ospelkaus, O. Wille, M. Succo, P. Ernst, K. Sengstock, and K. Bongs, Phys. Rev. Lett. 96, 180403 (2006). [13] N. B. Jørgensen, L. Wacker, K. T. Skalmstang, M. M. Parish, J. Levinsen, R. S. Christensen, G. M. Bruun, and J. J. Arlt, Phys. Rev. Lett. 117, 055302 (2016). [14] M.-G. Hu, M. J. Van de Graaff, D. Kedar, J. P. Corson, E. A. Cornell, and D. S. Jin, Phys. Rev. Lett. 117, 055301 (2016). [15] Z. Z. Yan, Y. Ni, C. Robens, and M. W. Zwierlein, Science 368, 190 (2020). [16] J. Tempere, W. Casteels, M. K. Oberthaler, S. Knoop, E. Timmermans, and J. T. Devreese, Phys. Rev. B 80, 184504 (2009). [17] A. Shashi, F. Grusdt, D. A. Abanin, and E. Demler, Phys. Rev. A 89, 053617 (2014). [18] F. Grusdt, Y. E. Shchadilova, A. N. Rubtsov, and E. Demler, Sci. Rep. 5, 12124 (2015). [19] A. G. Volosniev, H.-W. Hammer, and N. T. Zinner, Phys. Rev. A 92, 023623 (2015). [20] F. Grusdt and E. Demler, Proc. of the Int. School of Physics ’Enrico Fermi’ course 191 arxiv:1510.04934 (2016). [21] Y. E. Shchadilova, F. Grusdt, A. N. Rubtsov, and E. Demler, Phys. Rev. A 93, 043606 (2016). [22] M. K. and R. Seiringer, (2020), arXiv:arXiv:2003.12371. [23] F. Lingua, L. Lepori, F. Minardi, V. Penna, and L. Salasnich, New Journal of Physics 20, 045001 (2018). [24] V. Pastukhov, Journal of Physics A: Mathematical and Theoretical 51, 195003 (2018). [25] A. Lampo, S. H. Lim, M. ´ A. Garc´ıa-March, and M. Lewenstein, Quantum 1, 30 (2017). [26] A. Lampo, C. Charalambous, M. A. Garc´ıa-March, and M. Lewenstein, Phys. Rev. A 98, 063630 (2018). [27] A. Lampo, M. Garcia-March, and M. Lewenstein, Quantum Brownian Motion Revisited: extensions and applications (Springer International Publishing, 2019). [28] J. Bonart and L. F. Cugliandolo, Phys. Rev. A 86, 023636 (2012). [29] J. Bonart and L. F. Cugliandolo, EPL (Europhysics Letters) 101, 16003 (2013). [30] L. A. P. Ardila and S. Giorgini, Phys. Rev. A 92, 033612 (2015). [31] L. A. P. Ardila and S. Giorgini, Phys. Rev. A 94, 063640 (2016). [32] F. Grusdt, G. E. Astrakharchik, and E. Demler, New Journal of Physics 19, 103035 (2017). [33] L. A. Pe˜na Ardila, N. B. Jørgensen, T. Pohl, S. Giorgini, G. M. Bruun, and J. J. Arlt, Phys. Rev. A 99, 063607 (2019). [34] F. M. Cucchietti and E. Timmermans, Phys. Rev. Lett. 96, 210401 (2006). [35] S. P. Rath and R. Schmidt, Phys. Rev. A 88, 053632 (2013). [36] D. Benjamin and E. Demler, Phys. Rev. A 89, 033615 (2014). 17 [37] R. S. Christensen, J. Levinsen, and G. M. Bruun, Phys. Rev. Lett. 115, 160401 (2015). [38] F. Grusdt and M. Fleischhauer, Phys. Rev. Lett. 116, 053602 (2016). [39] Y. E. Shchadilova, R. Schmidt, F. Grusdt, and E. Demler, Phys. Rev. Lett. 117, 113002 (2016). [40] F. Grusdt, K. Seetharam, Y. Shchadilova, and E. Demler, Phys. Rev. A 97, 033612 (2018). [41] R. Seiringer, (2019), arXiv:arXiv:1912.12509. [42] M. Drescher, M. Salmhofer, and T. Enss, (2020), arXiv:2003.01982. [43] J. Levinsen, M. M. Parish, R. S. Christensen, J. J. Arlt, and G. M. Bruun, Phys. Rev. A 96, 063622 (2017). [44] N.-E. Guenther, P. Massignan, M. Lewenstein, and G. M. Bruun, Phys. Rev. Lett. 120, 050405 (2018). [45] W. E. Liu, J. Levinsen, and M. M. Parish, Phys. Rev. Lett. 122, 205301 (2019). [46] D. Dzsotjan, R. Schmidt, and M. Fleischhauer, arXiv:1909.12856 (2019). [47] S. I. Mistakidis, F. Grusdt, G. M. Koutentakis, and P. Schmelcher, New Journal of Physics 21, 103026 (2019). [48] S. I. Mistakidis, G. C. Katsimiga, G. M. Koutentakis, T. Busch, and P. Schmelcher, Phys. Rev. Lett. 122, 183001 (2019). [49] S. I. Mistakidis, A. G. Volosniev, N. T. Zinner, and P. Schmelcher, Phys. Rev. A 100, 013619 (2019). [50] W. E. Liu, Z.-Y. Shi, M. M. Parish, and J. Levinsen, arXiv:2002.01211 (2020). [51] N.-E. Guenther, R. Schmidt, G. M. Bruun, V. Gurarie, and P. Massignan, (2020), arXiv:2004.07166. [52] J. Levinsen, M. M. Parish, and G. M. Bruun, Phys. Rev. Lett. 115, 125302 (2015). [53] M. Sun, H. Zhai, and X. Cui, Phys. Rev. Lett. 119, 013401 (2017). [54] S. M. Yoshida, S. Endo, J. Levinsen, and M. M. Parish, Phys. Rev. X 8, 011024 (2018). [55] Z.-Y. Shi, S. M. Yoshida, M. M. Parish, and J. Levinsen, Phys. Rev. Lett. 121, 243401 (2018). [56] A. Camacho-Guardian, L. A. Pe˜na Ardila, T. Pohl, and G. M. Bruun, Phys. Rev. Lett. 121, 013401 (2018). [57] C. Charalambous, M. Garcia-March, A. Lampo, M. Mehboudi, and M. Lewenstein, SciPost Phys. 6, 10 (2019). [58] Y. Ashida, R. Schmidt, L. Tarruell, and E. Demler, Phys. Rev. B 97, 060302 (2018). [59] C. Charalambous, M. ´ A. Garc´ıa-March, G. Mu˜noz-Gil, P. R. Grzybowski, and M. Lewenstein, Quantum 4, 232 (2020). [60] W. D. Roeck and J. Fr¨ohlich, Communications in Mathematical Physics 303, 613 (2011). [61] C. Mukherjee and S. Varadhan, arXiv preprint arXiv:1812.06927 (2018). [62] R. Seiringer, arXiv preprint arXiv:1912.12509 (2019). [63] M. Hohmann, F. Kindermann, T. Lausch, D. Mayer, F. Schmidt, and A. Widera, Phys. Rev. A 93, 043607 (2016). [64] L. A. Correa, M. Perarnau-Llobet, K. V. Hovhannisyan, S. Hern´andez-Santana, M. Mehboudi, and A. Sanpera, Phys. Rev. A 96, 062103 (2017). [65] M. Mehboudi, A. Lampo, C. Charalambous, L. A. Correa, M. A. Garc´ıa-March, and M. Lewenstein, Phys. Rev. Lett. 122, 030403 (2019). [66] Q. Bouton, J. Nettersheim, D. Adam, F. Schmidt, D. Mayer, T. Lausch, E. Tiemann, and A. Widera, Phys. Rev. X 10, 011018 (2020). [67] H. J. D. Miller and J. Anders, Nature Communications 9, 2203 (2018). [68] W. Niedenzu, I. Mazets, G. Kurizki, and F. Jendrzejewski, Quantum 3, 155 (2019). [69] C. Charalambous, M. A. Garcia-March, M. Mehboudi, and M. Lewenstein, New Journal of Physics 21, 083037 (2019). [70] D. Barik and D. S. Ray, Journal of Statistical Physics 120, 339 (2005). [71] S. H. Lim, J. Wehr, A. Lampo, M. ´ A. Garc´ıa-March, and M. Lewenstein, Journal of Statistical Physics 170, 351 (2018). [72] S. H. Lim, J. Wehr, and M. Lewenstein, Annales Henri Poincar´e 21, 1813 (2020). [73] L. Pitaevskii and S. Stringari, Bose-Einstein Condensation (Oxford University Press, 2003). [74] M. Lewenstein, A. Sanpera, and V. Ahufinger, Ultracold Atoms in Optical Lattices (Oxford University Press, 2012). [75] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt, Rev. Mod. Phys. 86, 1391 (2014). [76] M. M. Khan, M. J. Akram, M. Paternostro, and F. Saif, Phys. Rev. A 94, 063830 (2016). [77] S. Gr¨oblacher, A. Trubarov, N. Prigge, G. D. Cole, M. Aspelmeyer, and J. Eisert, Nature Communications 6(2015). [78] F. E. Nixon, Handbook of Laplace transformation: fundamentals, applications, tables, and examples (PrenticeHall, 1965). [79] T. J. Alexander, K. Heenan, M. Salerno, and E. A. Ostrovskaya, Phys. Rev. A 85, 063626 (2012). [80] Y. Minoguchi, P. Kirton, and P. Rabl, arXiv preprint arXiv:1904.02164 (2019). [81] M. Reitz, C. Sommer, B. Gurlek, V. Sandoghdar, D. Martin-Cano, and C. Genes, arXiv preprint arXiv:1912.02635 (2019). [82] W.-Z. Zhang, J. Cheng, W.-D. Li, and L. Zhou, Phys. Rev. A 93, 063853 (2016). [83] Q. Wang and H. Zhan, Advances in Water Resources 75, 80 (2015). [84] G. Guarnieri, C. Uchiyama, and B. Vacchini, Phys. Rev. A93, 012118 (2016). [85] A. A. Clerk, M. H. Devoret, S. M. Girvin, F. Marquardt, and R. J. Schoelkopf, Rev. Mod. Phys. 82, 1155 (2010). [86] M. I. Shaukat, E. V. Castro, and H. Ter¸cas, Phys. Rev. A95, 053618 (2017). [87] J. Morris, Aircraft Engineering and Aerospace Technology 34, 25 (1962). [88] U. Weiss, Quantum Dissipative Systems (WORLD SCIENTIFIC, 2008). [89] E. E. Wollman, C. U. Lei, A. J. Weinstein, J. Suh, A. Kronwald, F. Marquardt, A. A. Clerk, and K. C. Schwab, Science 349, 952 (2015). [90] H.-P. Breuer, E.-M. Laine, J. Piilo, and B. Vacchini, Rev. Mod. Phys. 88, 021002 (2016). [91] H.-P. Breuer, E.-M. Laine, and J. Piilo, Phys. Rev. Lett. 103, 210401 (2009). [92] R. Vasile, S. Maniscalco, M. G. A. Paris, H.-P. Breuer, and J. Piilo, Phys. Rev. A 84, 052118 (2011). 18 [93] A. Ferraro, S. Olivares, and M. G. Paris, arXiv preprint quant-ph/0503237 (2005). [94] M. Inguscio, W. Ketterle, and S. Stringari, Quantum Matter at Ultralow Temperatures, Vol. 191 (IOS Press, 2016). [95] M. Bruderer, A. Klein, S. R. Clark, and D. Jaksch, Phys. Rev. A 76, 011605 (2007).