Dynamical variational approach to Bose polarons at finite temperatures
David Dzsotjan Richard Schmidt Michael Fleischhauer


	
	
	
	
	
	Department of Physics and Research Center OPTIMAS, University of Kaiserslautern, 67663 Kaiserslautern, Germany
	Wigner Research Center, Konkoly-Thege ut 29-33, 1121 Budapest, Hungary
	
	Max-Planck-Institute of Quantum Optics, Hans-Kopfermann-Strasse. 1, 85748 Garching, Germany
	Munich Center for Quantum Science and Technology (MCQST), Schellingstr. 4, 80799 München, Germany
	
	Department of Physics and Research Center OPTIMAS, University of Kaiserslautern, 67663 Kaiserslautern, Germany
	
	
	We discuss the interaction of a mobile quantum impurity with a Bose-Einstein condensate of atoms at finite temperature. To describe the resulting Bose polaron formation we extend the dynamical variational approach  of [Phys. Rev. Lett. 117, 11302 (2016)] to an initial thermal gas of Bogoliubov phonons. We study the  polaron formation after switching on the interaction, e.g., by  a radio-frequency (RF) pulse from a non-interacting to an interacting state. To treat also the strongly-interacting regime, interaction terms beyond the Fröhlich model are taken into account. We calculate the real-time impurity Green's function and discuss its temperature dependence. Furthermore we determine the RF absorption spectrum and find  good agreement with recent experimental observations. We predict temperature-induced shifts and a substantial broadening of spectral lines. The analysis of the real-time Green's function reveals a crossover to a linear temperature dependence of the thermal decay rate of Bose polarons as unitary interactions are approached. 
		
	
	
	
	
	
	Dynamical variational approach to Bose polarons at finite temperatures
    Michael Fleischhauer
    December 30, 2023
======================================================================




  
*Introduction.–

The interaction of a mobile impurity with a surrounding quantum bath is one of the paradigmatic models of many-body physics. The polaron introduced by Landau, Pekar and Fröhlich <cit.> to describe the motion of 
an electron in a lattice of ions, is formed by the dressing with lattice phonons and is a prime example of quasi-particle formation in condensed matter. More recently neutral atoms immersed in quantum degenerate gases of bosonic or fermionic atoms have attracted much attention since they are experimentally accessible
platforms allowing to study polaron physics with high precision and in novel regimes. Employing Feshbach resonances <cit.> it is possible  to tune the impurity-bath interaction from weak to strong coupling
and Rydberg states can be used to study impurities with
non-local interactions <cit.>.


	
	
The problem of a Fermi-polaron, i.e. an impurity interacting with a degenerate Fermi gas has been studied in 
a number of experiments in recent years <cit.>. This and related theoretical work <cit.> have led to a rather good understanding 
of this problem. In contrast, the description of impurities in a Bose-Einstein condensate (BEC), leading to the so-called Bose polaron is more involved <cit.>. The challenge for theory  is here directly related to the relatively large compressibility of the system, which allows for a much larger number of excitations that can be 
generated by the impurity. Also the experimental observation presented a major challenge due to three-body losses, and has only recently been achieved in experiments at JILA <cit.>, Aarhus <cit.>, and MIT 
<cit.>. Tuning through a Feshbach resonance  all regimes from weak to strong coupling were studied. While being in good general agreement with theoretical predictions, the Aarhus data showed deviations
for strong repulsive interactions, see Fig. <ref>, which were attributed to a nonzero temperature.  Following up on that, a recent extended T-matrix analysis predicted rather dramatic temperature effects <cit.>, most notably the appearance of new temperature-induced quasi-particle peaks <cit.>.  

	
	
These recent developments highlight the need of new theoretical approaches that take into account temperature as well as the creation of a large number of excitations in polaron formation. Indeed, while there exists by now a broad set of theoretical techniques to study Bose polarons at T=0,  extensions to T > 0 are not straight forward, so that, apart from first diagrammatic and functional determinant approaches <cit.>, theoretical progress remained highly limited so far.



We here discuss the Bose polaron at non-zero temperature extending the dynamic variational approach of <cit.> to an initial T
>  0 state of Bogoliubov phonons. The method is a non-equilibrium one and thus gives direct access to the dynamics of polaron formation after an excitation from a non-interacting state by an RF pulse, as it is observed in experiments. One of the key results is depicted in Fig. <ref>, where we show the calculated, trap averaged absorption spectrum A(ω) (color coding) for the Aarhus experiment. For 1/k_na_IB >0 a notable difference between the peak positions of T=0 calculations (full line) from the experimental values (green circles) 
was observed. In contrast, the results from our approach (red points) show good agreement. On the attractive side only uncorrected experimental data for the peak positions is available (blue circles, <cit.>). While this data exhibits a small deviation from our values [we define the mean peak response as ω̅=∫ dωωA(ω)], the onset of the polaron branch (purple circles, <cit.>) matches well the onset of the theoretical absorption spectra. Moreover, the calculated width of the absorption peaks (inset) agrees  well with the experiment. Finally while our approach predicts temperature-induced shifts, thermal quasi-particle broadening, and a temperature-dependent quasi-particle weight, in contrast to <cit.> we do not find evidence for a significant transfer of spectral weight to new quasi-particle peaks. 





  
*Model.– 


We here consider the interaction of a single impurity of mass M with a homogeneous Bose gas in d dimensions in a box of size L^d with periodic boundary conditions. Position and momentum operators of the impurity are  r̂⃗̂ and  p̂⃗̂.  We treat the BEC of condensate density n_0 in Bogoliubov approximation, i.e. in terms of non-interacting plane-wave excitations (phonons) of momentum k⃗, described by annihilation and creation operators  â_k⃗^(†). The condensate is characterized by the healing length ξ=1/√(2 g_BB n_0 m), where m is the mass of the BEC atoms and g_BB=2π a_BB /m describes their mutual interaction with a_BB the s-wave scattering length (ħ = 1). The Bogoliubov dispersion relation reads ω_= c k√(1+k^2ξ^2/2) with c= √(g_BB n_0/m) denoting the speed of sound and k=||. In order to formally decouple the impurity we
transform to  a co-moving frame  <cit.> using
Û= exp{ - i r̂⃗̂·_ph} where _ph=
∑_k⃗â_k⃗^†â_ is the total phonon momentum. Using Û^†  Û = - ∑_ â_^†â_ and Û^†â_k⃗Û
= â_k⃗  e^-i k⃗·r̂⃗̂, the Hamiltonian reads <cit.>


    H^LLP() =  1/2 M(-∑__k⃗_k⃗)^2 + ∑_ω__k⃗_k⃗ +  n_0 
             + g_IB/L^d/2∑_  n_0^1/2 W_ł_k⃗ + _-k⃗
             +/2 L^d∑_,^'[ V^+_,^'_k⃗_k⃗' + 1/2 V^-_,^'ł_k⃗_-k⃗' + _-k⃗_k⃗'],


where V^±_,^' = ł W_ W_'± W_^-1 W_'^-1$̊.now represents the total momentum of the system which is a constant of motion <cit.>. In the polaron frame the phonon dynamics attains a nonlinear term∼_ph^2which describes impurity-mediated phonon-phonon interactions that vanish in the limitM→∞. The impurity-BEC interaction strengthg_IB=2πa_IB/m_redis expressed in terms of thes-wave scattering lengtha_IBand  reduced massm_red= m M/(m+M), andW_=[k^2 ξ^2 / (2 + k^2 ξ^2)]^1/4. 
Crucially note that, due to thermal depletion, the condensate fractionn_0=n_0(T)is temperature dependent. For weak Bose-Bose interactions we haven_0(T) / n=1 - (T/T_c)^3/2at fixed total particle densityn.

Polaron  properties are encoded in the impurity Green's function S(t) = {e^i H_0 t e^- i H t ρ}, where the density matrixρdetermines the initial state of the system, andH_0andHare the Hamiltonian in absence and presence of the impurity bath interaction.S(t), also called `dynamical overlap', describes the dephasing dynamics of the system following a sudden quench ofg_IBat timet=0.  It can be measured using Ramsey spectroscopy as previously demonstrated in fermionic environments <cit.>. Fourier transformation ofS(t)in turn yields the (injection) absorption spectrumA(ω) = 2 Re∫_0^∞ dτ  e^iωτ  S(τ)in linear response, when the impurity is driven from a non-interacting state to a state with finiteg_IB<cit.>. 

 
 In contrast to previous studies of this problem <cit.>, we here consider finite temperature. This is accounted for by an initial density matrix ρ= ρ_T^ph ⊗|⟩⟨|where the phonon bath is in thermal equilibrium,ρ_T^ph= e^-β∑_ω___/Z(Z=[e^-βH_0]), for an impurity initially in an momentum eigenstate|⟩. 
In order to unambiguously identify the role of bath temperature and to allow direct comparison with previous studies <cit.> we focus in the following on=0.
While forT=0the initial state is invariant under the
Lee-Low-Pines transformation, this is no longer the case at finiteT. Introducing a projectorΠ̂_on eigenstates of_phwith eigenvalueone findsU^†ρU
=  ∫ d^3Q   Π̂_  ρ_T^ph ⊗|+⟩⟨+| .
We represent the thermal state of phonons as a Gaussian average over coherent states|ξ_⟩<cit.>:ρ_T^ph =∏_ ∫d^2ξ_k⃗    e^-|ξ_k⃗ |^2/n̅_k⃗  /πn̅_k⃗     |ξ_k⃗⟩⟨ξ_k⃗|.Heren̅_k⃗= 1/(e^βω_-1)is the average phonon number in modek⃗.
Thus we find 


    S(t)    =   ∫ d^3  Q  (⟨ψ_0(t)| Π̂_ |ψ(t) ⟩.

Here|ψ(t)⟩= e^-iH^LLP()t |ξ_⟩and|ψ_0(t)⟩=e^-iH^LLP_0()t |ξ_⟩describe the time-evolution of the initial states|ξ_⟩underH^LLP()andH_0^LLP()(Eq. (<ref>) forg_IB=0), respectively, where the impurity-momentum operatoris replaced by the c-number.
The overbar denotes the average over theξ_'s
and  we have used thatH_0^LLPcommutes with the projectorΠ̂_.  



  
*Dynamical variational ansatz.–

We calculateS(t)using wave functions|ψ(t)⟩and|ψ_0(t)⟩in a variational submanifold of Hilbert space constructed by time-dependent multi-mode coherent states|β(t)⟩=∏_e^β_â^†_-h.c.|0⟩=∏_|β_(t)⟩including a time-dependent phase|ψ(t) ⟩= e^-iϕ(t) |β(t)⟩. With|ψ_0(t) ⟩= e^-iϕ_0(t) |β^0(t)⟩, and definingΔϕ≡ϕ-ϕ_0one findsS(t)  = ∫ d^3  Q   ( e^-iΔϕ(t)  ⟨β^0(t) |  Π̂_  |β(t)⟩.The minimization of the LagrangianL=⟨ψ(t)|i∂_t -H^LLP|ψ(t)⟩, and similarlyL_0,  gives the Euler-Lagrange equationsd/d t (∂L/∂β̇_)-∂L/∂β_=0, 


    id/dtβ_(t) = (ω_ +^2/2M -k⃗/M·( - _ph)) β_(t) 
      + √(n_0)/L^d/2 W_
              + /L^d∑_(W_ W_ Re[β_(t)]+i   W_^-1W_^-1 Im[β_(t)]).

and similarly forβ_^0(t)(where=0), with (random) initial valuesβ_(0) = β_^0(0)=ξ_. The total phonon momentum in state|ψ(t)⟩is given by_ph=∑_|β_(t)|^2. The time-dependent Schrödinger equation impliesL=L_0=constand thus


    id Δϕ/dt  = (n_0 +√(n_0/L^d)∑_ W_β_^'(t) )-(_ph -^0_ph)^2/2M,

withβ_^'= Re[β_]. Note that different fromT=0, also the phonon momentum^0_phwithout interactions enters.









  
*Infinitely heavy impurity.–





We first discuss the limit of an infinitely heavy impurity,M→∞. In this caseH_0^LLP()andH^LLP()become independent ofand Eq. (<ref>) becomes


    S(t)  = ( e^-iΔϕ(t)⟨β^0(t) | β(t)⟩.

Moreover, the equations of motion (EOM) forβ_(t),β_^0(t), andΔϕ(t)become linear. This allows one to express the overlap⟨β^0(t)|  β(t)⟩as a matrix-Gaussian function in terms of the random initial variablesξ_, and the thermal average can be carried out analytically (see Supplementary Materials). Forg_BB=0our coherent state approach becomes exact and we have verified that our results match those from a functional determinant approach <cit.>. Eq. (<ref>) allows one to determine the temperature dependence ofS(t)and from its Fourier transform the absorption spectrum as shown in Fig. <ref> for a fixed impurity-bath interaction strength. One notices a substantial broadening with increasing temperature accompanied with a small shift of the peak position. 






  
*Finite impurity mass.–
	


For a finite impurity mass the variational states are-dependent and the corresponding integration in Eq. (<ref>) cannot be carried out upfront. Furthermore the EOM Eqs. (<ref>) become nonlinear  due to the presence of the total phonon momentum_ph.   ForT=0,_phis proportional to the conserved polaron momentum and thus vanishes in the case of an impurity initially at rest. This does not hold, however, at finite temperatures where_phalso contains the random initial amplitudesξ_. For this reasons the case of a finite impurity mass is substantially more involved compared to zero temperature even within the coherent-state variational approach and one has to resort to approximations. 

First, the projector on total-momentum eigenstates can be written asΠ̂_Q=1/(2π)^3∫ d^3z   e^iz⃗·(_ph-)where the action of the operatore^iz⃗·_ph=∏_e^iz⃗·  â_^†â_can be absorbed in a phase shift of the coherent amplitudesβ_^0(t)→β_^0(t) e^-iz⃗·. Thus evaluatingS(t)for a finite impurity mass is formally analogous to the infinite-mass case, however, demanding two additional integrations∫ d^3zand∫ d^3Q, which presents a numerical challenge. To address this issue, we here replace the phonon-momentum operator inΠ̂_by its expectation value with respect to the variational wavefunction, followed by an average over the random thermal amplitudes, i.e.e^iz⃗·_ph →e^iz⃗·_ph.
For an impurity initially at rest one has_ph=0andΠ̂_Qbecomes the unity operatorΠ̂_Q→δ^(3)(-_ph)=δ^(3)(). Theintegration is then trivial andS(t)obeys again Eq. (<ref>), whereβ_(t),β^0_(t), andΔϕ(t)now, however, follow equations (<ref>) and (<ref>) for a finite-mass impurity.










Accordingly, the EOM ofβ_(t)andβ_^0(t), which are given by Eq. (<ref>) with=0, are still nonlinear. As outlined in the Supplementary material, the nonlinear terms in Eqs. (<ref>) and (<ref>) can be approximated by a mean-field ansatz, where the quantities_ph β_(t)and(_ph-^0_ph)^2are effectively replaced by   n_  β_(t), and 2 ∑_k^2   n_(|β_(t)|^2-|ξ_|^2)respectively.
To this end we note that att=0theβ_(t)are  Gaussian random variables given byξ_and we can assume that they remain Gaussian for all times.
This finally  renders the EOM forβ_(t)in a linear form that is amenable to an analytical solution. Note that the equation for the total phase remains nonlinear but can be readily integrated. As a result the dynamical overlap can be expressed as matrix-Gaussian functions in terms ofξ_and  the thermal averaging can be carried out analytically.

 


The dynamical overlapS(t)calculated within the mean-field approximation is shown in Fig. <ref>(a) as a function of time. Results are shown for the attractive polaron for increasing temperatures in units of the critical temperatureT_c= (2π/m) (n/ζ(3/2))^2/3of a non-interacting gas in a box . 
While for zero temperatureS(t)approaches a finite value at large times given by the zero-temperature quasiparticle weight, it decays exponentially forT> 0with an asymptotic behavior|S(t)|∼Z(T)   e^-γ(T)  t.  The decay ratesγ(T)and thermal weightsZ(T)obtained from fits of the asymptotic tails are plotted in Fig. <ref>(b) and (c) as function ofT/T_cfor different impurity-Boson interaction strengths (for more details on the analysis of emerging, subleading quasiparticle branches see the Supplementary Materials). 
One recognizes an asymptotic power-law scaling of both quantities as function ofT/T_c. While this is reminiscent to the fermionic case <cit.>, where the exponents are given by the scattering phase shift at the Fermi momentum, it remains an open question to find analytical expressions for the exponents in the case of Bose polarons where no such a special finite momentum exists and scattering should predominantly take place at small momenta, or a momentum scale∼√(k_BT)determined by the thermal de Broglie wave length (at sufficiently large T). Remarkably, we find that close to unitary interactions the broadening of the quasiparticle peak, determined byγ(T), shows a crossover to a linear temperature dependence, which may be attributed to quantum critical behavior of impurities  in a Bose gas <cit.>.




	




From the Fourier-transform ofS(t), such as shown in Fig. <ref>(a), we have calculated the absorption spectrum for parameters of the experiment of Jø rgensen et al.<cit.> and <cit.>. The result, shown in Fig. <ref>, is in good agreement with the experiment. The peak positions extracted from the numerical simulations (red points) on the repulsive side coincide with the experimental values (green points) 
determined by Gaussian fits. In order to take into account the inhomogeneous density distribution in the experiment and the finite resolution of the spectrometer we have made a trap average and included a Gaussian broadening using the experimental parameter.  On the attractive side the experimental values for the lowest polaron branch, extracted from the onset of the measured absorption spectrum (purple points), agree with the onset of the absorption spectrum obtained from our numerical simulations.




  
*Summary.–

We discussed the physics of a single, mobile quantum impurity interacting with a BEC of atoms at finite temperature. Extending the dynamical variational approach 
of <cit.> to the case of an initial thermal state of Bogoliubov phonons, we showed how thermal effects  enter the Hamiltonian in the Lee-Low-Pines frame that is used to decouple the impurity from the phonon dynamics. To describe  polaron formation we calculated the real-time polaron Green's functionS(t)and from it the absorption spectrumA(ω)for a transition of the impurity from a state  non-interacting to a state interacting with the environment.  Strong  impurity-BEC interactions are accounted for by the inclusion of two-phonon terms in the 
 Hamiltonian <cit.>. Within the proposed variational approach one restricts the dynamics to a submanifold of coherent-state wavefunctions that are thermally averaged with Gaussian, random initial amplitudes, with weights determined by the temperature. In this approach observables, such as the the Ramsey signal given by the polaron Green's function are calculated from the solution of nonlinear EOMs for the coherent amplitudes that are subsequently averaged over the thermal distribution. While in the limit of an infinitely heavy impurity the EOMs become linear and the thermal average can be performed analytically, for a finite impurity mass  a mean-field
approximation is required to allow for an analytical thermal average. We calculated the temperature dependence  ofS(t)for different interaction strengths,  and found an asymptotic exponential decayS(t)∼Z(T)   e^-γ(T)  t. The  extracted decay ratesγ(T)and thermal weightsZ(T)show a  power-law dependence onT/T_c. Close to unitarity,1/k_na=0, the inverse polaron quasiparticle lifetime shows a linear dependence on temperature that is indicative of non-Fermi liquid behavior  in vicinity of the underlying quantum critical point <cit.>. The comparison of the theoretical absorption spectra with a recent experiment <cit.> shows that the inclusion of finite temperature corrections leads to excellent agreement between theory and experiment  on the repulsive side that was lacking in previous comparisons withT=0calculations. Our results are also in excellent agreement with the measured lowest polaron branch as well as the widths of the absorption spectra. In contrast to recentT-matrix calculations <cit.> we do not find a splitting or separate temperature-induced
quasi-particle peaks of substantial spectral weight.
 


 §.§ Acknowledgments

The authors like to thank N. Jørgensen, J. Arlt for providing the experimental data from <cit.>. We also
thank G. Bruun, E. Demler, Y. Shchadilova, and M. Zwierlein for useful comments and discussions. The work of D. D. and M. F. has been supported by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) within the SFB-TR 185 – 277625399.  R. S. is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy – EXC-2111 – 390814868. 

	

Supplementary Material





In the following we provide  details about the calculations presented in the main text.






§ I. INFINITE IMPURITY MASS



  §.§.§ A. Dynamical overlap and solution of equations of motion


In the limit ofM→∞the multi-mode coherent state dynamical variational ansatz yields the dynamical overlap  given by
Eq. (<ref>) of the main text. Introducing real and imaginary parts according toβ_= β_^'+ iβ_^''andξ_= ξ_^'+ iξ_^''we find


    S(t)   =   ( e^-iΔϕ(t)⟨β^0(t) | β(t)⟩
        =   ∏_𝐤∫dξ^'_∫dξ^''_1/πn_ e^-ξ_^'^2 
     + ξ_^''^2/n_e^-iΔϕ(t) e^-1/2∑_𝐪|β_(t) - β^0_(t)|^2 
     e^i∑_𝐪Im{β_(t)β^0_(t)^*}.

Hereβ_(t),β^0_(t)andΔϕ(t)depend onξ_^'andξ_^''through the initial conditions of the equations of motion. For convenience, we introduce the variablesβ̃_= √(L^d) β_andξ̃_= √(L^d) ξ_, and omit the tilde from here on. Using the rescaled variables, the EOMs read



    d/dtβ_^'    = ω_β_^'' + g_IB/L^dW_^-1∑_𝐪 W_^-1β_^'',        β_^'(0) = ξ_^', 
    d/dtβ_^''    = -ω_β_^' - g_IB/L^dW_∑_𝐪 W_β_^' - g_IB√(n_0) W_,        β_^''(0) = ξ_^'',

and


    d/dtΔϕ 
      =  g_IB√(n_0)/L^d∑_𝐤 W_β_^' + g_IBn_0,               Δϕ(0) = 0.

These equations can be expressed in matrix notation,

    β̇⃗̇(t)=K β⃗(t) + f⃗,

whereKis a coefficient matrix,f⃗is the vector containing the inhomogeneous terms, and the vectorβ⃗contains the dynamical variables as


    β⃗(t) = ([  β⃗_^'(t); β⃗_^''(t) ]).

Being an ordinary differential equation, the solution is easily found to be:


    β⃗(t) = e^Kt(ξ⃗ + K^-1f⃗) - K^-1f⃗,

where the initial conditionβ⃗(0)=ξ⃗was used. For the following calculations the coefficient matrixKhas to be diagonalized. SinceKis non-hermitian this requires the determination of  the right and left eigenvectors which can be done numerically, to obtain the diagonal matrixΛ= U_L K U_R, whereU_LandU_Rcontain the left and right eigenvectors ofK, respectively. It is convenient to reexpress Eq. (<ref>) in the compact form


    β⃗(t)= A(t)ξ⃗+(A(t)-1)w⃗,

where we introducedA(t) = e^K tandw⃗ = K^-1f⃗. This expression is used in Eq. (<ref>) which after integration over time yields


    Δϕ(t)    = g_IB√(n_0)∑_𝐪  W_[B(t)ξ⃗+(B(t)-1· t)w⃗]_ + c_0(t),

whereB(t) = U_R (e^Λt - 1 )Λ^-1U_L andc_0(t)=g_IB n_0 t.

The matrices and vectors can be written in separate blocks,

    A(t)   =([ C(t) D(t); E(t) F(t) ]),        B(t)=([ G(t) H(t); I(t) J(t) ]),     w⃗=([ x⃗_; y⃗_ ]),

which allows one to derive the formal solution for the real and imaginary parts ofβ_:


    β⃗^'(t)   =C(t)ξ⃗^'+D(t)ξ⃗^''+(C(t)-1)x⃗+D(t)y⃗, 
    β⃗^''(t)   =E(t)ξ⃗^'+F(t)ξ⃗^''+E(t)x⃗+(F(t)-1)y⃗, 
    Δϕ(t)   =g_IB√(n_0)∑_𝐪 W_[G(t)ξ⃗^'+H(t)ξ⃗^''+(G(t)-1· t)x⃗+H(t)y⃗]_ + c_0(t).



  §.§.§ B. Thermal average


Having derived the explicit solutions of the EOMs, we can now calculate the dynamical overlap at finite temperature. To this end the thermal average, encoded in the Gaussian integrations in Eq. (<ref>), has to be performed.

Since we assumep→0for the impurity momentum, the problem becomes spherically symmetric. We can thus perform an angular average ofβ_(t)andξ_(t)and only the dependence on the radial component ofremains. Discretizing the radial momentum component ofask_qwithq=1...N_one can write


    1/L^3∑_𝐤→∑_q=1^N_kΔ k/2π^2 k_q^2.

HereΔkdenotes the difference between the discretized radialkvalues. Taking this into account, and substituting the solutions (<ref>) for the rescaledβ_^',β_^''and forΔϕ,  the exponents in Eq.(<ref>) can be written in a vectorised form. Grouping them according to terms quadratic, linear and constant inξ, one  obtains


    -iΔϕ(t)    = -iρ⃗^I ξ⃗ + -iκ, 
    
      -1/21/L^3∑_k⃗ |β_(t) - β_^0(t)|^2    = - 1/2ξ⃗^T M^II_1 ξ⃗ + ρ⃗^IIξ⃗ - 1/2w⃗^T M^II_2 w⃗, 
    i1/L^3∑_k⃗Im{β_(t) β_^0(t)}   = iξ⃗^T M^III_1 ξ⃗  + ρ⃗^IIIξ⃗.

The detailed form of the  matricesMfollows directly from the substitution of the solutions (<ref>). Also the Gaussian probability distribution can be written in a vectorised form,


    ∏_𝐤∫dξ^'_∫dξ^''_1/πn_k e^-ξ_^'^2 
     + ξ_^''^2/n_k = ∫d^2N_kξ  √(Γ/(2π)^2N_k)   e^-1/2ξ⃗^T Γξ⃗,

whereΓis a2N_k ×2N_kdiagonal matrix, containing the average phonon occupation numbers,n_kas 

    Γ = Δ k/π^2( [ [     k_1^2/ n_1                              ;                             ⋱               ;                               k_N_k^2/ n_N_k ]                                                                                                                                            0;                                                                                                                                            0 [     k_1^2/ n_1                              ;                             ⋱               ;                               k_N_k^2/ n_N_k ] ]).

Rewriting Eq.(<ref>) in terms of these vectorized expressions, one obtains a simple,2N_k- dimensional Gaussian integral:


    S(t) = ∫d^2N_ξ  √(Γ/(2π)^2N_)   e^-1/2ξ⃗^T (Γ + M_T)ξ⃗   e^ρ⃗·ξ⃗   e^ -1/2w⃗^TM_0w⃗    e^-iκ,

whereM_T  = M^II_1 + i (M^III + (M^III)^⊤),M_0 = M^II_1,ρ⃗ = ρ⃗^II + i(ρ⃗^III - ρ⃗^I), andκis a phase term originating fromΔϕ(t).
The integral can be performed analytically. Importantly, however, the form Eq. (<ref>) allows for an  instructive physical interpretation in terms of the different scattering processes mediated by the impurity that contribute to the many-body dynamics of the system. In fact, the exponent quadratic inξ⃗represents scattering of excitations from the thermal reservoir back into the thermal reservoir. The linear-ξ⃗exponent contains the contribution of scattering of excitations from the BEC into the thermal part and vice versa, while the third exponent expresses the scattering processes from the BEC back to the BEC. At zero temperature, only the latter term survives while above the critical condensation temperatureT_c, only the first term is present. This shows that our model indeed encompasses all possible channels of scattering processes in the system which on its own is an important finding.  In particular, it is the second term that allows the conversion between atoms in the condensate fraction and the thermal contribution in BECs at finite temperatures. In previous approaches using functional determinants to, e.g., describe Rydberg Bose polarons in finite temperature BECs <cit.> this contribution was not accounted for and the present approach shows not only how it can be derived explicitly, but it also provides a form amenable to numerical evaluation.

Performing the Gaussian integrals, we obtain a closed, analytic expression for the dynamical overlap:


    S(t) = e^1/2ρ⃗^T(Γ + M_T)^-1ρ⃗/√((1 + Γ^-1 M_T))  e^ -1/2w⃗^TM_0w⃗    e^-iκ.

Besides the advantage that this result enables one to forego any numerical sampling and averaging, it also has a very clear structure: the finite-temperature dynamical overlap can be written as the product of aT=0part and a finite-temperature factor. Hence the absorption spectrum, obtained by Fourier transformation becomes a simple convolution integral. AsT →0, due to Bose-statistics,  the occupation of finite-momentum states goes to zero. As a result,Γ^-1 →0, thus the finite-temperature part disappears, leaving 


    lim_T→ 0S(t) = e^ -1/2w⃗^TM_0w⃗    e^-iκ,

which means we are indeed left with the impurity-mediated scattering processes from the BEC back into the BEC, and a result that recovers the previous findings at zero temperature <cit.>.







§ II. FINITE IMPURITY MASS AND MEAN-FIELD APPROXIMATION


For a finite impurity massMandT>0, the Lee-Low-Pines transformation acting on the initial density matrix makes the initial state of the impurity effectively dependent on the total phonon momentum (according to the relationU^†ρU
=  ∫ d^3Q   Π̂_  ρ_T^ph ⊗|+⟩⟨+| stated in the main text), which has nonzero fluctuations in a thermal state. Consequently, the variational wavefunctions must be evaluated for all finite valuesand the dynamical overlap is subsequently obtained as the average over(see Eq. (<ref>) in the main text),


    S(t)  = ∫ d^3  Q  ( e^-iΔϕ(t)⟨β^0(t) | Π̂_ |β(t)⟩
    = 1/(2π)^3∫ d^3  Q ∫ d^3 z  ( e^-iΔϕ(t)⟨β^0(t) |  e^-i(_ph-)·z⃗ |β(t)⟩.

Making use ofe^i_ph·z⃗|β_^0(t)⟩= |e^i·z⃗ β_^0(t)⟩≡|β̃^0(z,t)⟩this expression can again be expressed in terms of coherent-state overlaps


    S(t)  = 1/(2π)^3∫ d^3  Q ∫ d^3 z   e^i·z⃗ ( e^-iΔϕ(t)⟨β̃^0(z,t) 
     |β(t)⟩.


To deal with the challenge that the double integration in this expression prohibits a direct analytical approach to the problem,  we employ here a mean-field approximation which replaces the phonon momentum operator in the expression e^-i_ph·z⃗by its expectation value in the variational wavefunctions, averaged over the random initial values, i.e. e^-i_ph·z⃗. Since_phis proportional to the conserved total momentum, for an impurity initially at rest it holds_ph=0. With this theandz⃗integrations can be performed explicitly
and we obtain an expression that is formally identical to the dynamical overlap in theM→∞case, i.e.


    S(t)= ( e^-iΔϕ(t)⟨β^0(t) | β(t)⟩.

Note, however, that unlike the infinite mass case,  the EOMs determining the time-evolution of the coherent amplitudes and phases are now nonlinear equations


    id/dtβ_(t) = (ω_ +^2/2M +k⃗/M·_ph) β_(t) 
      + √(n_0)/L^d/2 W_
     + /L^d∑_(W_ W_ Re[β_(t)]+i   W_^-1W_^-1 Im[β_(t)]),



    id/dtΔϕ(t) =  n_0 +√(n_0/L^d)∑_ W_Re[β_(t)] 
    -1/2 M(_ph -^0_ph)^2,

 owing to the fact that they still contain the total phonon momentum_ph=∑_     β_^*(t)β_(t)which is a random variable due to the initial values ofβ_andβ_^0. 

In order to treat these nonlinearities in the EOMs we make further approximations: First one notes that for an impurity initially at rest, the total phonon momentum vanishes when thermally averaged over the stochastic initial values,_ph=0. Consequently, when  replacing the phonon momentum by its average over the stochastic initial values, the nonlinear terms in Eqs. (<ref>) and (<ref>) would disappear. This would however be a too strong approximation and we proceed instead as follows:



Att=0the coherent amplitudesβ_(t=0)=ξ_are indeed Gaussian random variables.  Assuming that they remain approximately Gaussian with respect to the average over theξ_variables one can apply a mean-field factorization (δx ≡x- x)


    δ x δ y δ z   =   δ xδ y δ z+δ xδ z δ y+δ yδ z δ x
    δ w δ x δ y δ z    =   δ wδ x δ y δ z + δ w δ x δ yδ z -  δ wδ x δ yδ z
        +   δ wδ y δ x δ z + δ w δ y δ xδ z -  δ wδ y δ xδ z
       +    δ wδ z δ x δ y + δ w δ z δ xδ y -  δ wδ z δ xδ y.

In the present case, for an impurity initially at rest one has∑_β_^*β_ =0. Moreover, there is only a coupling between coherent amplitudes of different modes if the two-phonon terms are significant. Neglecting this cross coupling one finds 


    β_^*(t)β_(t) -β_^*(t) β_(t)   =   ξ_^*ξ_ -ξ_^* ξ_ =
    δ_n_,
    β_(t)β_(t) -β_(t) β_(t)   =   ξ_ξ_ -ξ_ ξ_ =0,

as well as


    β_^0 = 0,    β_^0(t)β_^0(t)  =0,    β_^0*(t)β_^0(t) = ξ_^*ξ_ =
    δ_n_.

Furthermore, recognizing that 


    β_(t)β_^0(t)  =0,    β_^*(t)β_^0(t) = ξ_^*ξ_ =
    δ_n_,


	
one is finally led to the following mean-field approximation of the nonlinear terms in the EOMs


    _ph β_(t)=∑_β_^*(t)β_(t) β_(t) ≈∑_(β_^*β_ β_ +(β_^*β_-β_^* β_)  β_ +(β_β_-β_ β_) β_^*)≈ n_ β_(t),

and 


    (_ph-_ph^0)^2 =∑_   (|β_(t)|^2 
    - |β_^0(t)|^2 )(|β_(t)|^2 
    - |β_^0(t)|^2 )
            ≈ 2 ∑_   ( ( β_(t)^*β_(t) -β_(t)^* β_(t) )
     β_(t)β_^*(t)
    +β_(t)^0*β_^0(t) β_^0(t)β_^0*(t)
    -β_(t)^0*β_(t) β_^0(t)β_^*(t)
    -β_(t)^*β_^0(t) β_(t)β_^0*(t) )
            = 2 ∑_^2 n_ ( |β_(t)|^2 - |ξ_|^2).

	
These expressions give, after insertion into the EOM forβ_the linear equations


    id/dtβ_(t)    =   (ω_ +^2/2M(1+2 n_)) β_(t) 
      + √(n_0)/L^d/2 W_
     + /L^d∑_(W_ W_ Re[β_(t)]+i   W_^-1W_^-1 Im[β_(t)]),


 
 
    id/dtΔϕ(t)    =    n_0 +1/M∑_^2 n_ [|β_(t)|^2-|ξ_|^2]
    +√(n_0/L^d)∑_ W_Re[β_(t)].

When two-phonon terms are not important, we can further replace|β_(t)|^2in the second term of the latter equation by the solution of Eqs.(<ref>) without mode-mixing two-phonon terms


    d/dtβ_(t) ≈ -i (ω_ +^2/2M(1+2 n_)) β_(t) 
      -i √(n_0)/L^d/2 W_
    =
     -iΩ_β_(t) -ig_IB√(n_0/L^d) W_k,

which reads


    β_(t) = -i g_IB√(n_0/L^d) W_k/Ω_k(1-e^-iΩ_k t) +ξ_ e^-iΩ_k t= β_(t) +ξ_ e^-iΩ_k t.

With this the solution of the EOMs forβ_(t)andΔϕ(t)have the same form as the corresponding equations in theM→∞limit and hence all further calculations can be performed in a similar way as outlined for the case of an  infinitely heavy impurity.  

In order to illustrate the impact of a finite impurity mass and the related  recoil the impurity experiences in collisions with bath excitations in Fig. <ref> we show the dependence of the impurity Green's functionsS(t)on the mass ratioM/m_Bat fixed temperatureT/T_c=0.2and interaction strength1/k_na_IB=-1. In Fig. <ref> we give examples for the dynamical overlap and the absorption spectra for different temperatures for the cases of attractive, strongly attractive and repulsive interactions. TheT=0results agree with those from Ref.<cit.>.








§ III. SPECTRAL ANALYSIS OF THE TIME-DEPENDENT GREEN'S FUNCTION S(T)

	



The real-time Green's functionS(t)and the corresponding spectral functionA(ω), such as shown in Fig. <ref>, contain  information about possible quasi-particle excitations (here at zero momentum), including their energy, lifetime, and quasiparticle weights. In a regime where quasiparticle excitations dominate the spectral response,  the full time-resolved Green's functionS(t)may be approximated as a sum over such quasiparticle contributions

    S(t)≈ G_QP(t)≡∑_i^N_QP Z_i e^-i E_i t e^-γ_i t e^-i ϕ_i.

HereE_iand1/γ_irepresent  quasiparticle energies and  lifetimes, respectively, and the quasiparticle weights are given by0≤Z_i≤1. Theϕ_iaccount for a possible phase arising due to overlapping spectral weights. 




In Ref. <cit.> is was predicted using a diagrammatic approximation that at finite temperatureTthe attractive polaron peak splits into two quasiparticle excitations which share the weight of the original (T=0) polaron peak. Our approach does not predict such a splitting. By fitting the approximation (<ref>) to the full numerical data ofS(t), we analyse the frequency-resolved properties of the impurity Green's function and we find that, while the polaron peak shifts energetically, it does pertain and even increase its quasiparticle weight as temperature increases,  contrary to the prediction in Ref. <cit.>.  





Importantly, the analysis of the numerically obtained spectral function (see Fig. <ref>) shows that, for mobile impurities, weight is accumulated in emergent, additional  quasiparticle-like spectral features that accompany polaron excitations present atT=0.  As an example we show in Fig. <ref> the results of fitting Eq. (<ref>) withN_QP=3to the numerical dataS(t)in the strongly interacting regimek_n a_IB = -1. In each row of Fig. <ref> we show from the left to right as solid curves the magnitude ofS(t), its real and imaginary part, and the spectral functionA(ω)resulting from Fourier transformation. The columns in turn correspond to increasing temperature fromT/T_c = 0toT/T_c = 0.4. For finiteTthe quasiparticle fit functionG_QP(t)is applied to the data and the resulting fits are shown as dashed lines. 

It is evident that not only the absolute value|S(t)|(left column) is extremely well reproduced by the fit but even the full complex functionS(t)from short to long times (middle column). Consequently, also the spectra are fit well by the simple approximation (<ref>) of the full impurity Green's function. Moreover, we find that the fits usingG_QP(t)yield lifetimes, quasiparticle weights and energies that are fully consistent with the simpler single quasiparticle fit at long times used in the analysis  discussed in the main text (see corresponding Fig. <ref>), supporting the robustness of our analysis. 


	

The analysis in Fig. <ref> also not only reveals  how the polaron energyE_1 = E_polshifts withT/T_c(Fig. <ref>, c.f. also the vertical, dashed lines in the right column in Fig. <ref>), but also helps to identify the temperature-dependence of the two additional spectral features atE_2,3accompanying the dominant polaron peak atE_1. Indeed we find that while the features at energiesE_2,3do not absorb weight from the dominant polaron peak, they do become increasingly well defined for increasing temperature and show a clear energetic separation from the main peak at energyE_1(Fig. <ref>). 

We emphasize that the Fourier transform ofS(t)yields the frequency resolved (retarded) Green's functionG(ω). Hence, our analysis also gives insight into the approximate pole structure ofG(ω), which can, to a good approximation, be described by a small number of poles in the complex frequency plane with increasingly larger residues as temperature increases.  It remains an interesting and open question whether these features are robust under more advanced variational wavefunction manifolds that, e.g., include squeezing terms beyond the LLP and Bogoliubov unitary transformations included in our approach, and whether the emergent spectral features can be efficiently occupied in an adiabatic preparation in experiments while avoiding three-body losses.