Stochastic-field approach to the quench dynamics of the one-dimensional Bose polaron

We consider the dynamics of a quantum impurity after a sudden interaction quench into a one-dimensional degenerate Bose gas. We use the Keldysh path integral formalism to derive a truncated Wigner like approach that takes the back action of the impurity onto the condensate into account already on the mean-field level and further incorporates thermal and quantum effects up to one-loop accuracy. This framework enables us not only to calculate the real space trajectory of the impurity but also the absorption spectrum. We find that quantum corrections and thermal effects play a crucial role for the impurity momentum at weak to intermediate impurity-bath couplings.Furthermore, we see the broadening of the absorption spectrum with increasing temperature.


I. INTRODUCTION
The interaction of a mobile impurity with a surrounding many-body quantum system is one of the most prominent and oldest problems in condensed matter physics. The polaron, initially considered by Landau and Pekar [1,2] to describe an impurity electron interacting with the lattice vibrations of a solid is a prototypical scenario to study quasi-particle formation. In more recent years, neutral atoms immersed in a surrounding ultra-cold gas have attracted widespread attention due to great experimental controllability which enable the study of novel exotic regimes of the polaron. For example, the impuritybath coupling can be tuned via Feshbach resonances [3]. Here, the fundamental principles at work are the same as in the original problem; the impurity forms a polaron through interacting with the collective excitations of the surrounding superfluid. While the Fermi polaron has been subject to extensive experimental study [4][5][6][7][8][9][10][11][12], the Bose polaron has only been realised in a few experiment [13][14][15][16][17]. These experiments hint towards a delicate interplay between equilibrium and out-of-equilibrium effects.
There has been extensive work addressing the Bose polaron in equilibrium [18][19][20][21][22][23][24][25]. In this work, we focus, in contrast, on the out-of-equilibrium Bose polaron. More precisely, we consider quench dynamics, involving the abrupt immersion of an impurity into a homogeneous Bose gas, at finite temperature. The quench can be realised through a Feshbach resonant radiofrequency pulse [3] that switches on the impurity-condensate interaction nearly instantaneously. Such quench dynamics have been studied either at zero temperatures, or focused on the (extended) Fröhlich models (or very similar approximations) [26][27][28][29][30][31][32][33][34][35][36][37][38][39]. In much of the equilibrium and some of the out-of-equilibrium treatments, the Lee-Low-Pines transformation [40], which eliminates the impurity from the problem at the expense of adding an additional quartic vertex, has proven extremely useful. However, when considering finite temperature, this requires further approximations.
The vast majority of the existing literature has focused on three-dimensional systems, where the (extended) Fröhlich model is widely applicable. In 1D, the situation is different. In [21] it was shown that the ap-plicability of Fröhlich-type approximations are somewhat limited in 1D and while the mass balanced case is integrable for the Fermi-gas [41] no such limit exist in the case of a Bose gas. Instead, it is natural to incorporate the effects of the impurity on the condensate already at the mean-field level. This can be done efficiently for a single impurity at zero temperature using the LLP transformation, but this method does not extend to several impurities and is also not straightforward to generalise to finite temperatures. To circumvent those complications, several approaches do not eliminate the impurity from the problem and can address finite temperature [30,42], and for example, treat the impurity in a manner related to the coherent state representation of the impurity [43]. It has been shown in [33,44,45] that a product wave function within the tree-level approximation, yields inconsistent results with those obtained in the more accurate LLPframe in 1D. Therefore, it is perhaps more appropriate to treat the impurity in a position-momentum path-integral as it highlights the particle nature of the impurity. Using the coherent state path integral for the condensate has been shown to yield good results in 1D not only for the polaron but also the bipolaron problem [33,46,47]. This is conceptually close to the approach originally developed by Feynman [42] and applied to the Bose polaron in [19,24] with the main difference being that we do not expand the condensate around a homogeneous density, and our focus lies on out-of-equilibrium phenomena.
In this work, we develop a conceptually simple and numerically tractable approach to address quench dynamics at finite temperature in a general manner. Ultimately, this is achieved by mapping the dynamics to a set of deterministic differential equations with stochastic initial conditions. By averaging over the different trajectories, expectation values can be calculated within one-loop accuracy. We then use this methodology to study an impurity's evolution after a sudden interaction quench into the bath. We find that the impurity delocalises quickly for weak impurity-bath couplings and that observables like the impurity's velocity crucially depend on the incorporation of quantum and thermal effects. In the case of strong impurity-bath couplings, we observe self-trapping and quantum corrections and thermal effects to be e considerably less pronounced. This work is organised as follows. In Section II A we derive the truncated Wigner approximation from the Keldysh path integral representation of the problem at hand. In this section, we also show how to obtain the absorption in the language of semi-classical dynamics. We proceed by specifying the initial Wigner function in Section II B and show how to regularise the divergences that arise in the one-dimensional setting. We proceed by briefly outlining the numerical considerations in Section II C. To conclude, we discuss the results for an impurity at rest and finite momentum in Section III and outline further directions in Section IV.

II. METHODOLOGY
A. The Keldysh formalism for the Bose polaron In the section, we want to outline the derivation of the equations of motion and discuss the truncated Wigner approximation for the polaron problem. We start by con-sidering the Hamiltonian, where m (M ) denotes the mass of the bosons (impurity atom),φ(x) is the Bose field operator, g BB (g IB ) is the boson-boson (boson-impurity) interaction strength and X (P ) denotes the position (momentum) operator of the impurity. Moreover, we left the interaction potential between the impurities and the condensate general instead of directly assuming s-wave scattering. We will not employ a delta function here but a smoothed out version of it for numerical reasons, that will be discussed later. In the following, we are going to apply the Kedlysh formalism to (1). The expectation value of an arbitrary observable (see [48] for a detailed discussion) Ω(X,P ,φ † (x),φ(x), t) is given by where W is the Wigner function that depends on the initial density matrix and will be specified below. Ω W (X c (τ ), P c (τ ), φ * c (x, τ ), φ c (x, τ ), t) is the Weyl ordered operator of the observable one wants to calculate (again see [48] for more details). The subscript c denotes the classical field and the subscript q the quantum field, which describes the quantum fluctuation around the classical saddle point solution. Those two fields arise when mixing forward and backward contour in the Keldysh formalism. The D denotes the integration over all field configurations in space and time. In contrast to that the first integral in (2) can be understood as a normal integral in X 0 and P 0 . The one-loop approximation is now to drop all terms of order two and higher inh, which corresponds to an expansion up to second order in the quantum fields. We then find that there are, in fact, only terms linear in the quantum fields, and the action is given by It is now easy to see that all the quantum fields can be easily integrated out and yield functional delta distributions enforcing the classical equations of motion The only challenge remaining at this point is to integrate over the initial conditions weighted by the Wigner function. In practice this is done by sampling initial conditions according to the Wigner function, solving the classical equations of motion and then averaging the desired observable over the calculated trajectories. In this framework it is now straightforward to calculate the impurity dynamics. Before proceeding we would like to make some remarks about the classical equations of motion and make contact with the equilibrium case to show that even without taking the first-order correction into account those equations give satisfactory results in the limiting case of heavy impurities. For the equilibrium case we assume the impurity to travel at constant velocity d dτ X c (τ ) = v, implying d dτ P c (τ ) = 0, which in turn tells us that |φ c (x, τ )| 2 is symmetric around the impurity position. Together with the equilibrium assumption this directly leads to the conclusion that the bosonic field takes the following Putting it all together and defining x = x − X c (τ ) we find the equilibrium equation We can now compare this equation with the one obtained by performing a Lee-Low-Pines transformation and find that it has in fact the same form as the one found in [33,44,47,49], where it has been shown that quantities obtained like the effective polaron mass or the polaron energy are in excellent agreement with results obtained by quasi-exact quantum monte carlo methods. The only difference is that instead of the reduced mass the boson mass appears in front of ∂ 2 x , which can be traced back to the fact that in this derivation the effect of normal ordering is lost. However, this effect is unimportant for heavy impurities. To summarize, we showed that the equation obtained by employing a coherent state ansatz for the bosonic field paired with a position momentum representation for the impurity reduce to the correct mean-field equations in the equilibrium case if the impurity mass is large.
Another quantity of great interest is the impurity Green's function [18,27,31] from which the absorption spectrum can be calculated by taking the Fourier transformation. The impurity Green's function is defined by whereĤ 0 stands for (1) with V (x −X) = 0 andρ is the initial density matrix. We now observe that this has the same structure as the trace that is considered to derive the Keldysh path integral, with the only difference, that the forward and backward contour differ by an extra interaction term. We can therefore perform the same steps as when deriving (2). The only difference in S is the resulting impurity boson interaction where we have introduced the notation d n The resulting magnitude of the interaction is changed by a factor of 1/2 and a new purely classical term arises . Lastly, we note that there is also a quadratic term in φ q (x, τ ) and X q (τ ) now. If we want to keep the accuracy up to one loop order this term has to be considered, which somewhat complicates matters. An ad hoc approximation is to drop this term altogether and therefore staying in a strictly semi-classical regime. To see when this approximation is justified one can simply compare the arsing terms and their order of magnitude. We note, that the |φ q (x, τ )| 2 term directly competes with the |φ c (x, τ )| 2 term. As long as one is within the applicability region of a general c-field treatment, |φ q (x, τ )| 2 will be small compared with |φ c (x, τ )| 2 , whenever the condensate deformation is not large, which corresponds to small and intermediate η. For the corrections in the impurity degrees of freedom one has to com- We note that through partial integration the derivative terms can be brought onto bosonic field variables. One then realizes, that the second derivative of the fields is going to be small compared to the first derivative for weak coupling. Additionally, for large impurity masses M the magnitude of X q will stay small. Those two considerations show that one would expect the absorption spectrum yield reliable results for weak to intermediate impurity-boson coupling and potentially a wider range of couplings if the impurity is sufficiently heavy. We refer to section II D for more details on the validity the approach presented here.
We also note that this approximation is only made when calculating the absorption spectrum and the impurity Green's function and all dynamical results do not rely on this approximation. Henceforth, the additional stochastic term will be dropped. This leaves us with the expression where we denote the average with respect to the initial Wigner function as ... W . We are now in the position to calculate the real space trajectories of the impurities for finite temperature as well as the absorption spectrum.

B. The quench protocol and the initial Wigner function
In the following, we want to specify the quench protocol and the initial Wigner function. We start by briefly outlining the initial state and then specify the Wigner function for a 1D quasi condensate. Here, we will also discuss all the regularisation necessary to arrive at a divergence-free quasi condensate description. The quench protocol is the following: we start with a free impurity and an interacting superfluid at temperature T . At t = 0, the interaction between the superfluid and the impurity is turned on instantly. Experimentally this is realised through a Feshbach resonance [3]. We assume that the initial density matrix can be written as a direct product of the state of the impurity (which is assumed to be pure) and the thermal density matrix of the superfluid As a consequence, the Wigner function also factorises, and we can sample the initial conditions independently.
For the condensate, we employ a quasi condensate description; in 1D, this is best done by employing a density and phase representation. We note that this has been used before in [50], in the trapped gas case. Since we want to focus on a homogeneous gas in continuum here, we need to regularise the non-condensed part. In this representation, the condensate field operator can be written asφ The density operator and the phase operator can be expressed within the Bogoliubov approximation [51] aŝ where u k and v k are the usual Bogoliubov modes. For this treatment to be valid one has to be in the vicinity of a weakly interacting Bose gas, which can be characterised by the Tonks parameter γ = 1/(2n 2 0 ξ 2 ) [52, 53] which should be less than unity where ξ = 1 √ 2gBBn0 is the healing length. We refer to the section II D for a detailed discussion of the validity of the presented approach. In the path integral this corresponds to a shift of variables meaning that instead of integrating over φ k , corresponding to the operatorsâ † k . In the standard way, we can now write down the thermal Wigner function (within the coherent state representation) for the a k [48] with the Bogoliubov dispersion ω k . From this, it can be seen immediately, that the average condensate particle number is given by N 0 /L = ρ 0 . In order to account for the quantum and thermal depletion, we fix the total particle number N and then choose N 0 according to N 0 = N − N d (this is done for every realisation), where after proper regularisation (see [54,55]) with the single-particle dispersion e k = k 2 /(2m). Here, a first-order T-matrix approximation was employed, and in line with the Bogoliubov theory up to one loop, µ is the chemical potential within the Bogoliubov approximation and hence is not temperature-dependent. The 1/2 is needed to cancel the extra factor that stems from the symmetric ordering. After averaging, this reproduces the expected result for a thermal quasi condensate in 1D. It should be noted that µ has to be chosen consistently with the total particle number, which can be done by fixing one reference point, where the total particle number is known. Henceforth we assume a mean-field density at T = 0. This then fixes µ in the Bogoliubov approximation through µ = g BB n 0 (T = 0) We can now use (16) to determine the total particle number which remains fixed throughout the calculation. We can now sample individual realisations of the condensate, whose description is free of IR and UV divergences. Upon closer inspection, one might realise that even though the mean of the phase and density corrections are zero, the variance scales up to linearly with the system size L. A direct result of this is that the computational time needed to achieve convergence also scales with the system size. This computational challenge can be tackled by increasing the system size gradually until the results are independent of the system size and then validating certain data points with larger system size. Furthermore, this restricts the discretisation of space, as outlined in [56], which will be discussed in the next section. Because the effect of the impurity is local, we find relatively low dependence on the system size already for small systems. Lastly, we assume that the impurities are not entangled (namely can be represented by a product wave function) and are localised in space around x 0 or equivalently in momentum space localised around q at t = 0. It is therefore natural to choose a wave packet as their initial wave function Here a 0 is an external parameter that determines how localised the initial state is. The Wigner function in this setting is well known to be [48] W (X 0 , P 0 ) = 2 exp −2a 2

C. Numerical considerations
In this section, we will show that three quantities can describe the whole parameter space, and we will briefly outline the discretisation of space and all the subtleties involved. We now define the following scaled parameters where we choose τ s =ξ/c, x s = √ 2ξ =ξ, and where the speed of sound is given by c = gBBn0 m . Dropping the twiddle we then find the Hamiltonian +η where η remains unchanged and α = m M . From here the following equations of motion can be obtained In order to tame the extensive variance and ensure numerical stability, we have to choose the discretisation l = L/N grid as outlined in [56], namely l has to be large enough to satisfy ρ 0 l 1, while at the same time ensuring that the energy cut-off introduced by l does not alter the physics. This translates to l < ξ, λ, where ξ is the healing length that sets a natural length scale for our problem and λ is the thermal de Broglie wavelength. Lastly, we note that we make the following choice for the interaction potential which converges to the delta distribution as l → 0, but has the advantage of being smoother than δ x /l, where δ x is the Kronecker delta.

D. Validity of formalism
In this subsection we will address the regime of validity of the formalism used in this work. We start by giving some general arguments for the validity of the approach followed by a term-by-term discussion of the higher-order corrections and their order of magnitude. While for the absorption spectrum we have to restrict our considerations to weak or moderate boson-impurity couplings, we would like to stress that for all dynamical properties calculated (excluding the absorption spectrum) the result is strictly non perturbative in g BB and g IB . Hence we can safely say that the method presented is valid for γ ≤ 1, which is equivalent to the range for the Gross-Pitaevskii equation, which corresponds to the tree level approximation of (1). The same reasoning applies to g IB , leading to the conclusion that our results, at least for short times are valid across the whole range from weak to strong impurity-boson couplings. For longer times one has to consider the deformation of the condensate. Another interesting point to consider is that as soon as α is small (meaning the impurity is heavy) the accuracy of the presented approach will improve further, since the impurity will behave more classically. In fact α can be seen as a strict control parameter for the corrections arising due to higher order terms in X q (τ ). Combining those considerations with the known fact that in general the truncated Wigner approximation is exact for short time scales [48], we can conclude that for short to intermediate time scales our results are trustworthy for all values of boson-boson and boson-impurity coupling. For weak coupling it is ab initio reasonable to assume that the presented results hold for longer time scales, since higher order quantum corrections should accumulate slowly if at all. However, no such universal statement can be made for strong couplings. There are two contributions of orderh 2 or higher, which were neglected in our approach and would have to be added to (9) if one wanted to solve the problem exactly, the first one coming from the boson-boson interaction and the second one being due to the impurity-boson interaction. The first one takes the form We note that this term is closely related to the standard Bogoliubov approximation, with the main difference being that the classical field here is taken to be the deformed field. We note that in no point of our simulations the expectation value of |φ * c (x, τ )| 2 falls below the value of 3ξ, meaning that the it is safe, in the spirit of the well-established Bogoliubov approximation, to neglect higher-order terms in the quantum fields, which only scale linearly in φ c and are of O(φ c φ 3 q ) which is certainly small compared to O(|φ c (x, τ )| 3 φ q ), as long as the density of the condensate stays larger than the healing length. We note that in the case of very low density gases and strong boson-impurity coupling the above argument comes under question and it is a priori not clear whether the outlined method is reliable in this regime. This regime was not investigated in the present work. This argument is further underlined by the remarkable accuracy of pure c-field methods (which do not take firstorder quantum corrections into account) for the equilibrium polaron, see [35,44,46,49], where the c-field approximation was shown to also hold for low-density gases and strong coupling, and the fact that the boson-boson interaction to all orders in g BB and to first order inh are already taken into account in our calculations. A some-what more complicated expression is obtained when considering the impurity-boson interaction term. Here, one finds that already in (9) all orders of X q (τ ) are present and the higher order terms take the following form Again for approximating the impact of those terms it is convenient to bring the derivative terms to the bosonic field operators. It then becomes clear that all corrections can be understood as a gradient expansion in the bosonic field operators, whose impact is certainly small for short time scales and also for longer time scales as long as the coupling stays moderate. Besides the gradient terms we also note that each correction is accompanied by an increasing power of X q (τ ) terms, which are also small at short times and whose impact can be controlled by α.
To summarize, while the validity of this approach for very large couplings and long time scales can not be judged a priori, we note that for short time scales the results here hold regardless of impurity-boson coupling strength and also note that α serves as a control parameter for the approximation in the impurity degrees of freedom.

A. Post-quench density profile
In this subsection, we focus on the density of the condensate for repulsive and attractive impurity-bath couplings at different times after the quench and supplement those findings with the variance of the position σ 2 X = (X − X ) 2 and the variance of the velocity We find a dynamically distinct behaviour for weak and strong couplings on the repulsive and the attractive side.
In Fig. 1, the condensate density at different times and the evolution of the variance of the position and the variance of the velocity for repulsive interactions are shown. Before discussing the results in more detail we notice that even for η = 50 the minimum of the condensate density is still larger than 3ξ, indicating that the approximations made later for the absorption spectrum are justified. We note that for weak coupling, the impurity delocalises faster than a free impurity would. For stronger couplings, the impurity stays localised much longer in time, indicating self-trapping. The velocity variance saturates after a finite time, and the time scale is inversely proportional to In all plots the parameters are α = 0.2, γ = 0.04, n0 = 5ξ and a0 = 3/ξ. As expected, we can clearly see that the impurity deforms the condensate over time. It also becomes apparent that the impurity delocalises over time and that this effect is slowed down with increasing η, which can be understood by a self-trapping argument given in the main text. In all plots the parameters are as follows α = 0.2, γ = 0.04, n0 = 5ξ and a0 = 3/ξ. We can see that the impurity attracts surrounding particles, which in turn depletes the condensate. As can be seen in (b), this appears as a repulsively interacting polaron.
the impurity boson coupling. We note that this can be explained by realising that two competing effects are determining the dynamics. Namely, the impurity tends to distribute the repulsion equally throughout the condensate, causing the impurity to delocalise and the opposing effect of self-trapping, where the impurity deforms the condensate and then self-traps in the deformation. It is intuitively clear that self-trapping will not occur for weak couplings, which explains the different behaviours seen in Fig. 1. It can also be seen that the variance of the impurity velocity can exceed the speed of sound, which is associated with the emission of non zero energy excitation, indicating an energy transfer from the impurity to the bath, which has also been observed in [57]. The temperature influences the value with which the position variance and velocity variance saturate; it does not significantly influence the timescale. For strong coupling, the temperature dependence becomes relatively weak, which can be explained by noting that the impurity-boson scat-tering length determines the relevant energy scale, which is much larger than the thermal length in this case. In Fig 2 the same situation for attractive couplings is shown. Here the difference between strongly attractive and moderate attractive couplings becomes obvious. We observe that in the case of moderately attractive couplings, the impurity not only diffuses but also forms a purely attractive polaron. In contrast, for strongly attractive couplings, an attractive polaron with repulsive interactions is observed, and the time scales of the polaron formation are prolonged. This difference is a dynamical effect and can be understood by noting that when the interaction is turned on, particles from the condensate start to accumulate around the impurity, which in turn depletes the condensate around the impurity. The superfluid is interacting itself, and therefore the depletion is filled by the particles around it, with the time scale being set by the boson-boson interaction. Meanwhile, the impurity-boson interaction strength determines the number of particles that can accumulate around the impurity before the boson-boson action prevents further accumulation. At the same time the impurity delocalises, which prevents the formation of a well-defined peak around the impurity and the interaction can thus seem repulsive for a long time scale. Ultimately this is of course only a metastable state. This process continues for a longer time when the impurity-boson interaction is large, resulting in a polaron that looks repulsive which can be observed in Fig 2(b).

B. Impurity velocity
Another quantity that is of great interest is the impurity velocity. The out-of-equilibrium case gives insight into the polaron formation and the time scales at work. It is also of great importance for the equilibrium case since it can be used to calculate the effective mass of the polaron [58].
Here, the impurity is not at rest when the quench occurs but instead carries some finite momentum. The sudden quench of the impurity-boson interaction leads to a momentum transfer from the impurity to the surrounding particles, and we expect a slow down of the impurity. The time evolution for the velocity of the impurity is depicted in Fig. 3. Here it becomes apparent that quantum corrections have a noticeable impact on the evolution of the velocity at weak to intermediate coupling. This can easily be understood by noting that the impurity is treated as a point particle on the mean-field level, and as observed in Fig 1, which is not a valid approximation for weak couplings. This also explains why the MF velocity is lower than the corrected solution. We also note that the steady-state velocity increases with the temperature, which can be understood by noting that the surrounding gas has a higher average squared velocity for increasing temperature, and therefore the momentum transfer will be smaller. In contrast, for strong couplings, the impurity stays localised and approximating it as a point particle is, therefore, less of a simplification. The same behaviour can be observed for the temperature dependence, which is more critical for weak and intermediate couplings, which is again due to the scattering length being larger than the thermal length. We also note that the impurity transfers some of its momentum to the Bose-gas and then relatively quickly reaches an equilibrium velocity for not too strong interactions. For stronger interactions, we observe different behaviour. After an initial abrupt slow down of the impurity, the impurity velocity changes sign, which directly results from the back action from the condensate. We note that a similar abrupt slowdown has been observed in the three dimensional setting in [59]. The velocity then performs a damped oscillation around its final velocity.

C. The absorption spectrum
Next we turn our attention to the (injection) absorption spectrum A(ω) = 2Re ∞ 0 G(t)e iωt dt [10,12], which can be measured using Ramsey spectroscopy [9,10,14]. It can be calculated by taking the Fourier transform of the impurity green's function, which characterises the dephasing of the system and is closely related to the Lochschmidt echo [60]. The absorption spectrum gives essential information about the polaron formation and can be used to estimate the polaron energy and lifetime [61,62]. At this point, we also want to stress that pure ∞ 0 G(t)e iωt dt for different temperatures T calculated using the truncated Wigner approximation for α = 0.2, γ = 0.04, n0 = 5ξ and a0 = 3/ξ mean-field (MF) calculations in position momentum basis are not sufficient to calculate the impurity Green's function. This can be seen by noting that there is no averaging for a classical calculation, which means there is no dephasing between different trajectories, and therefore one always obtains |G| = 1 in purely classical calculations.
Our results are depicted in Fig. 4, it can be observed that the quasi-particle peak widens with increasing temperature, which is also reported in 3D [31]. However, in contrast to the 3D case, [27], where the extended Fröhlich model was considered, we do not find several peaks on the repulsive side. We also note that the overall amplitude decreases with η, and the quasi-particle peak gets washed out with increasing η, which is a direct consequence of the orthogonality catastrophe [63][64][65], hence the emergence of a clear quasi-particle peak comes under question. We note that the absorption spectrum shows a functional dependence associated with quasi-particle behaviour, and we do not find an infrared dominated regime as observed in other one dimensional systems [61]. Those findings are supplemented by the overlap G(t). Here we can see another vastly different feature of the one-dimensional case, compared to the three-dimensional case [27,31], namely that |G|, approaches zero even for moderate couplings, signalling the onset of the orthogonality catastrophe. As expected, the dephasing becomes more rapid with increased temperature and increased g IB .

IV. SUMMARY & OUTLOOK
In summary, by leveraging the Keldysh formalism, we derived a truncated Wigner approach to study dynamical properties of the bose polaron in 1D. This allowed us to reduce the problem to simulating semi-classical equations of motion with stochastic initial conditions. We showed how to adequately account for temperature effects of the surrounding bath by sampling the phase and density of the condensate. We also discussed how to regularise the arising divergences that typically occur in such one dimensional systems. The method presented here takes the back action of the impurity onto the condensate into account and is therefore applicable from weak to strong impurity bath couplings.
We then used this framework to calculate the dynamics of an impurity after sudden immersion in a surrounding bath and the absorption spectrum. By considering the condensate density and the position/velocity variance, we showed that there is a distinct dynamical behaviour associated with the strong and weak coupling regime, namely our results indicate self-trapping of the impurity for strongly repulsive interactions, and we also find a repulsive polaron on the attractive side. We also investigated the temperature dependence of the polaron formation and found a substantial influence of quantum corrections on dynamical properties like the velocity of the impurity, showing the necessity to go beyond pure meanfield considerations. Lastly, we considered the absorption spectrum and the impurity Green's function. Here, we observed a clear quasi-particle peak for weak to intermediate couplings. In contrast to that, we see that the quasi-particle peak is washed out for strong couplings and that temperature effects widen the quasi-particle peak. In contrast to the higher dimensional case, the impurity Green's function approaches zero even for weak couplings.
At this point, we want to stress that our approach is neither limited to 1D nor a single impurity. It could therefore serve as an exciting starting point to explore higher dimensional systems, as well as the interplay of several impurities. While the generalization to several impurities is quite straightforward we want to stress that, as pointed out for example in [66][67][68], the generalisation to higher dimensions is highly non-trivial in general. First we note that in higher dimensions it is not possible to use bare contact interactions for the boson-boson interaction and the boson-impurity interaction simultaneously when employing the approach. One has to resort to using more realistic interaction potentials for at least one of them as has, for example, been done in the three dimensional context in [35,64]. Another major challenge is a purely numerical one, as it becomes more costly to sample the Bose fields in higher dimensions. Nevertheless we expect that the presented method may be paired with some small approximations to address higher dimensional systems.