Dynamics of polaron formation in 1D Bose gases in the strong-coupling regime

We discuss the dynamics of the formation of a Bose polaron when an impurity is injected into a weakly interacting one-dimensional Bose condensate. While for small impurity-boson couplings this process can be described within the Froehlich model as generation, emission and binding of Bogoliubov phonons, this is no longer adequate if the coupling becomes strong. To treat this regime we consider a mean-field approach beyond the Froehlich model which accounts for the backaction to the condensate, complemented with Truncated Wigner simulations to include quantum fluctuation. For the stationary polaron we find an energy-momentum relation that displays a smooth crossover from a convex to a concave dependence associated with a non-monotonous relation between impurity velocity and polaron momentum. For larger momenta the energy is a periodic function including regions of negative impurity velocity. Studying the polaron formation after turning on the impurity-boson coupling quasi adiabatically and in a sudden quench, we find a very rich scenario of dynamical regimes. Due to the build-up of an effective mass, the impurity is slowed down even if its initial velocity is below the Landau critical value. For larger initial velocities we find deceleration and even backscattering caused by emission of density waves or grey solitons and subsequent formation of stationary polaron states in different momentum sectors. In order to analyze the effect of quantum fluctuations we consider a trapped condensate to avoid 1D infrared divergencies. Using Truncated Wigner simulations in this case we show under what conditions the influence of quantum fluctuations is small.


Introduction
The dynamics of a quantum impurity coupled to an interacting many-body environment is one of the most fundamental problems of many-body physics. Of particular interest is the dressing of the impurity with elementary excitations of the host systems leading to the formation of a quasiparticle. A paradigmatic model of such a quasiparticle in condensed matter physics is the polaron, introduced by Landau and Pekar [1,2] to describe the interaction of an electron with lattice vibrations in a solid, and which is key for understanding transport, response and induced interactions in many systems. In recent years ultra-cold quantum gases have become a versatile experimental testing ground for studying polaron physics with high precision and in novel regimes. For example, employing Feshbach resonances [3] for neutral atoms, the impurity-bath interaction can be tuned from weak to strong coupling. Furthermore many-body environments of different quantum statistics and with different interactions can be considered. While impurities in a degenerate Fermi gas, called Fermi-polarons, have been studied in a number of experiments only a small number of experiments exist on Bose polarons [4][5][6][7]. Here due to the large compressibility of the Bose gas a larger amount of excitations can be created by the impurity and interactions among the environment particles become increasingly important. , where M is the impurity mass, n the average boson density,c the speed of sound andξ the healing length (see text). The impurity emits density waves before converging into on of two stationary states, marked in red and blue. For large momenta (d) also solitons are created which can lead to a change in the direction of motion of the impurity. The evolution is shown for an impurity-Bose mass ratio M = 10m, Tonks parameter γ = 0.1 and impurity-Bose coupling constant g IB = gnξ.
While much theoretical work exists addressing the ground state properties of Bose polarons [8][9][10][11][12][13][14][15][16] there is still little understanding of finite temperature properties [7,[17][18][19] and even more so of its non-equilibrium dynamics [4,[20][21][22][23]. What happens when an impurity is injected into a weakly interacting Bose condensate? What is the dynamics of the formation of a polaron and under what conditions and on what time scales can a stable quasiparticle form at all? We will address these questions in the present paper considering a point impurity interacting with a weakly interacting, onedimensional Bose condensate in the full range of impurity-boson coupling strength, see Fig. 1. The limit of weak impurity-boson interaction can be well described by the generation and subsequent binding or emission of Bogoliubov phonons from the impurity [24] in terms of a Hamiltonian similar to that of the Froehlich model used in solid state systems [25]. Most existing studies of the non-equilibrium dynamics of Bose polarons is based on this model [8,[26][27][28][29][30][31]. It is however no longer well suited in the limit of strong impurity-boson coupling and thus we here follow a different approach. Starting from a full quantum description of the interaction of a mobile impurity with the condensate we employ a mean-field approach that takes the backaction of the impurity onto the condensate into account as in [32][33][34], but keeps the entanglement between impurity and BEC by working in a co-moving frame. This approach was shown to be very accurate for the prediction of ground state properties of Bose-polarons [35] and bi-polarons [36] even for very strong impurity-boson couplings as long as the Bose-Bose interaction is weak. Quantum effects are then taken into account by considering Bogoliubov excitations on top of the deformed condensate, which are here treated within a Truncated Wigner approximation [37][38][39]. The advantage of this approach as compared to the Froelich model and its extensions [12,14,28] is the substantially reduced number of deformed Bogoliubov phonons created by the impurity in such a description. As we will show the effect of these modified phonons can be neglected even in the non-equilibrium dynamics in many situations allowing for a comprehensive study of the polaron dynamics in terms of non-linear c-number differential equations.

Energy-momentum relation in a homogeneous 1D Bose gas
Before considering the time evolution of an impurity injected into a 1D Bose gas of neutral atoms, let us discuss the stationary properties of a polaron at finite momentum relative to the Bose gas. A widely used approach to describe Bose-polarons for weak boson-boson and impurity-boson couplings is to consider the interaction of the impurity with Bogoliubov excitations of the unperturbed condensate. The resulting model is reminiscent of the Froehlich model in condensed matter physics. Due to the large compressibility of the Bose condensate, this model is however no longer adequate if the impurity-boson coupling becomes large. In the latter case a growing number of phonons is generated at the location of the impurity and boson-boson interactions become relevant. For this reason we here start from the full quantum model and apply a different approximation scheme.

Model and modified mean-field approach
A single mobile impurity coupled to a Bose gas in a homogeneous one-dimensional system is described by the Hamiltonian for ℏ = 1. Here m (M ) is the boson (impurity) mass, g IB (g) the impurity-Bose (Bose-Bose) interaction constant,φ(x) the bosonic field operator andr (p) the impurity position (momentum) operator. In the following we consider the case of repulsive coupling between all particles, i.e. g, g IB > 0.
Since the system is homogeneous, its total momentump +P B is conserved, wherê is the momentum of the Bose gas. The infinite homogeneous system is treated here as the limit of a finite system of length L with periodic boundary conditions with L → ∞. In order to exploit the translational invariance we apply the Lee-Low-Pines (LLP) transformation [40]Û LLP = exp(−irP B ), leading to the transformed Hamiltonian Due to translation invariance,Ĥ LLP no longer depends onr, andp =Û † LLP (p + P B )Û LLP is the conserved total momentum in this frame and can be replaced by a c-number p. We note that due to the LLP transformation,φ ( †) (x) describes the creation/annihilation of a boson in a frame co-moving with the impurity, such that n(x) = ⟨φ † (x)φ(x)⟩ is the Bose gas density relative to the position of the impurity, or rather the relative impurity-boson density-density correlation function, to which we refer from now on for simplicity as the Bose gas density.
SinceĤ LLP is an interacting many-body Hamiltonian, a complete solution of the dynamics is difficult without further approximations. In [35,36] we have shown that for a weak boson-boson interaction, indicated by a small Tonks parameter γ = gm/n, the ground state properties of a single or a pair of Bose polarons at rest is very well captured by a mean-field approximation that takes the backacktion of the impurity to the condensate into account and goes beyond the standard Froehlich model, which fails if g IB ≫ gnξ. Hereξ = 1/ √ 2gnm is the rescaled healing length with the reduced massm = (1/m + 1/M ) −1 . This motivates us to use this mean-field approximation also for the case of a moving impurity. We will test its validity by taking into account quantum fluctuations within a Truncated Wigner approach in Section 4. The meanfield approximation amounts to replacing the field operatorφ(x) by a complex order parameter ϕ(x), whose time evolution is determined by the non-linear Schrödinger equation [32] i∂ t ϕ( Here v(t) = p − P B (t) /M is the impurity velocity.
Since we are interested in the formation dynamics of polarons either after a sudden quench or an adiabatic turn-on of the impurity-Bose coupling constant g IB , the Bose gas is assumed to be initially in its ground state ϕ(x, t = 0) = √ n at t = 0, where n is the average density of bosons. As the Bose gas carries no initial momentum, the conserved total momentum is equal to the initial impurity momentum p = M v(t = 0). We simulate the time evolution with periodic boundary conditions numerically using a Fourier split-step method [41]. If not stated otherwise we choose the system size L large enough such that signals are not able to reach the boundary for all times t considered, i.e. L ≫ct, wherec = gn/m is the rescaled speed of sound.

Stationary state and energy-momentum relation
We proceed by discussing the stationary properties of a polaron moving with non-zero total momentum p by characterizing the steady-state solutions of Eq. (3) [23]. Since the solution has a simpler form if the impurity velocity v is used as a parameter, rather than the conserved momentum p, it is useful to derive the stationary state as a function of v as in [32], and then calculate the corresponding momentum by . This is possible since v(t) = v is constant in the stationary state. The analytical expression of the state for a fixed v was derived in [16,32,35,42], see Appendix A for technical details and analytical expressions. As shown in [32], below a critical impurity velocity v c the system has two stationary states for every given value of v, see Fig. 2 b) for the corresponding density profiles. Above the critical value no stationary solution exists. Exactly at the critical velocity, the two states are equal two each other. v c depends on the coupling constant and can be found from the solution of [16,32] g IB gnξ It is plotted in Fig. 2 a). For small interaction g IB ≪ gnξ, it agrees with the prediction of the Froehlich model v c =c [24,43]. However in the limit of strong interactions g IB ≫ gnξ the condensate is strongly depleted at the impurity position, resulting in a vanishing critical velocity v c → 0. The polaron energy is given by the difference in energy of the full system with and without the impurity E pol = E(g IB , p) − E(g IB = 0, p = 0), where E(g IB , p) is the expectation value of the LLP Hamiltonian Eq. (2) in a coherent state with amplitude ϕ(x, t). E pol is shown as a function of the polaron momentum in Fig. 3 a), where the two stationary states are distinguished by solid and dashed lines. It becomes clear that the polaron state is unique for every momentum p and the two different states only refer to different parts of the energy-momentum relation. For the lower momentum state the energy-momentum relation was derived in [15]. The phase/density of the Bose gas is shown in b)/c) for g IB = gnξ and different total momenta p, below p/n = π indicated by the blue crosses in a). Other parameters are as in a) except a finite system size L = 30ξ. For νπ ≤ p/n ≤ (ν + 1)π and integer ν the condensate picks up a constant phase gradient of νπ over its total length.
From the existence of a critical velocity one might conclude that a Bose polaron exists only up to maximum polaron momentum p max depending on g IB , and the energy momentum relation E pol (p) would terminates at some value of momentum. This is not the case. Instead when increasing p further the solution smoothly crosses over into the second steady state with smaller kinetic velocity and larger condensate depletion. As can be seen from Fig. 3 c) the condensate depletion grows with increasing momentum. As a consequence the kinetic mass of the polaron, defined as the ratio of polaron momentum and velocity M * = p/v, increases with p as shown in Fig. 4 a). When crossing from the momentum regime of the first solution of Eq. (3) to the second solution, the increase of the mass with momentum becomes larger than linear. This leads to a non-monotonous relation between polaron velocity v and momentum p, plotted in Fig. 4 b). Note that the polaron velocity always stays below the weak-coupling critical value v c . When the momentum reaches the value p max = nπ, the state of the condensate is exactly equal to a dark soliton, such that the energy is E(p max ) = 4 3 nc and the impurity velocity goes to zero v(p max ) = 0. In this case the condensate phase winds by π over its entire length, which for periodic boundary conditions corresponds to half a flux quantum piercing through the ring. Since the density is fully depleted the kinetic polaron mass M * (p max ) diverges at this point.
We note that the relation between momentum, energy, and velocity of the polaron has already been found for momenta |p| < nπ [23]. A peculiar behavior of E(p) is however seen when the total momentum is increased further. The energy starts to decrease with increasing momentum, corresponding to a negative group velocity ∂E/∂p and consequently a negative impurity velocity. In fact, as can be seen from Fig. 3a, E(p) is a periodic function of p/n with period 2π. This is because the properties of the polaron are determined by collective excitations of the Bose gas, whose energymomentum relation in 1D is periodic with period 2πn due to Luttingers theorem. Here it can be interpreted as follows: If the total momentum is in the range (ν−1)π ≤ p/n ≤ νπ, with ν being an integer, the condensate picks up an integer winding of its phase over its whole length (period) in the stationary state. If ν is even, the momentum picked up by the background condensate exceeds the total momentum and the excess must be compensated by a relative motion of the impurity against the condensate, corresponding to negative values of v.
We will see in the following that the periodic behavior of the energy-momentum relation can give rise to negative asymptotic impurity velocities after injecting it into the condensate with large positive initial velocity. We note that the reversal of the impurity velocity has already been seen in experiments with a strongly interacting Bose gas (γ ≫ 1), where Bloch oscillations of an impurity subject to a constant force have been observed [20].

Mean-field description of polaron formation
We now discuss the dynamics when an impurity is injected into a homogeneous condensate with finite initial momentum. First we consider a quasi-adiabatic turn-on of the impurity-boson coupling and subsequently discuss a sudden quench.

Quasi-adiabatic evolution
Let us first investigate the dynamical properties of the system when the impurity-boson coupling constant g IB (t) is turned on slowly compared to the other time scales of the system. We will show that even though the energy spectrum of the full system is gapless in the thermodynamic limit, a local adiabatic following of the polaron ground state is possible if the initial velocity of the impurity is subsonic. In order to achieve a smooth turn-on protocol, we choose a time dependence of the coupling according to: g IB is the final coupling constant and T the turn-on timescale, which we chose large compared to the inverse chemical potential 1/gn. Since the critical velocity v c , below which stationary states exist, depends on g IB it is also time-dependent and the time evolution differs qualitatively whether the impurity momentum is below or above v c at any time. The time evolution of the condensate for the two cases is shown exemplarly in Fig. 6, where the blue and red lines indicate the two stationary states corresponding to the final velocity of the impurity. First, we focus on the case of a slow impurity v(t) < v c (t) for all t, such that Eq. (3) has a stationary solution for all v(t). The evolution of the impurity velocity is exemplary shown for different T in Fig. 5b). Even though the impurity is always subsonic it is decelerated to a finite value v(t f ) < v(0), which increases monotonously with T , see Fig. 5a). Notably also in the limit T → ∞ the impurity is still slowed down and v(t f ) does not converge to v(0). This is caused by the formation of the polaron. In the instantaneous ground state, the total conserved momentum p is the polaron momentum and is related to the impurity velocity by the effective mass m * off the polaron v(t) = p m * (t) .
Since the effective mass increases monotonously with the coupling constant g IB [35], the impurity must decelerate when the impurity-Bose coupling constant is turned on. For finite turn-on times T , density waves are created during the formation of the polaron leading to an additional friction force.
To quantify the quasi-adiabatic slow-down we fit an exponential α 1 e −α 2 /T + α 3 to the simulated impurity velocity (see Fig. 5a). From the fit the final impurity velocity in the adiabatic limit is determined by v ad = α 1 + α 3 . It is shown as a function of the conserved total momentum p = M v(0) in Fig. 5c), where it becomes apparent that the deceleration occurs for all p. Due to the gaplessness of the system to collective excitations the simulated final velocity is always slightly below the quasi-adiabatic value, following from Eq. (6) with m * replaced by the effective mass of the stationary solution, see Fig. 4. This is because the adiabatic theorem [44] strictly does not hold. Assuming strict adiabatic following we can derive the relation between v(t) and p as a function of the instantaneous coupling constant g IB (t) in the polaron ground state (see Appendix A). Fig. 5a) shows that the simulated time evolution follows the instantaneous ground state reasonably well for large T except for the small difference shown inFig. 5c), see Appendix B for details.
We now consider an impurity that is initially faster than the critical momentum v(0) >c. In that case, Eq. (3) does not have a stationary solution in the initial phase, which leads to the creation of density waves Fig. 6b) and a friction force acting on the impurity. Fig. 7a) shows that the impurity is quickly decelerated until its velocity is below the critical v c (g IB ). Afterward, the system again follows the instantaneous ground state quasi-adiabatically, resulting in a slower deceleration.
An important difference between these two processes is, that the second slow deceleration, related to the adiabatic formation of the polaron, is reversible, while the other one is not. This can be seen in Fig. 7b), where the impurity-Bose coupling constant is turned on and off again by g IB (t) = g IB sin 2 ( 1 2 π t/T ). Here impurities that are initially below the critical momentum are almost brought to a standstill when the coupling constant is at a maximum but accelerate again to the initial velocity when g IB is turned off again. In contrast, impurities starting above criticality are not reaching their initial momentum again after the sweep. The small deviation in the final velocity for subcritical trajectories is again caused by the system not being gaped in the thermodynamic limit.

Quench
Here we examine the evolution of the system when the coupling constant is abruptly quenched at t = 0.
The time evolution of the impurity velocity after a sudden turn-on of the interaction with the condensate is shown in Fig. 8a) and b) for different initial velocities and different ratios of bare impurity mass M to boson mass m. The quench leads to radiation of density waves until the system reaches a steady state. For all initial conditions and parameters, a friction force [23] is exerted on the impurity and slows it down until the final velocity is reached smaller than the critical v c . This agrees again with the analytic prediction, that the system has a stationary state only below v c . However, for a large initial momentum, the impurity velocity is non-monotonic, which cannot be explained by a frictional force alone.
If the impurity is heavy a rather unexpected behaviour is found for sufficiently large  initial momentum, see Fig. 8 b). First, as expected the deceleration is slower in the case of heavier impurities. However for a sufficiently large mass, the impurity is not only slowed down, but the velocity can change its direction before converging to a constant velocity. In this case, the background condensate attains an additional momentum by building up a finite phase gradient away from the impurity locally approaching a stationary polaron solution with negative impurity velocity as discussed in the previous section.
This effect is examined in more detail in Fig. 9 a). It shows the final impurity velocity as a function of the conserved momentum p = M v(0) of a heavy impurity M = 10m. Again the final velocity is for all initial conditions in the interval [−v c , v c ]. However, as mentioned in Section 2.2 the system has two stationary states for each velocity and we examine which of these states are populated. For this we first focus on Fig. 1 b) -d), showing the evolution of the Bose gas density for different total momenta. The red and blue lines are the analytically calculated stationary states, where we used the final impurity velocity from the simulation as a parameter, rather than the total momentum. Depending on the initial conditions the system converges in either of these two states. In order to quantify the overlap of the final polaron state with either of the two stationary states we determine the generalized contrast Here |ϕ(t)⟩ is the simulated state of the system and |ϕ i ⟩ the steady states, evaluated at the simulated final velocity. The contrast is 1 if the system converges into the first state and −1 for the second. This definition differs from the standard expression of contrast by the term ⟨ϕ 1 |ϕ 2 ⟩ in the denominator, which we need since the stationary states are not orthogonal. We evaluate the scalar products over an interval l which we choose such that the time evolution is properly converged within it. The contrast is depicted in Fig. 9 a) by the color code and it becomes apparent that the polaron changes its state when its velocity intersects with the critical value v c . This is explained by the states being equal at the critical velocity.
The explanation for the non-monotionic impurity velocity and the change in the movement direction at larger total momenta can be seen in Fig. 1 d). At large momenta, the impurity carries enough energy that a grey soliton is created besides density waves which carry additional momentum. This is visible as the local depletion moving away from the impurity at a slower speed as the initially created density waves in Fig. 1 d). Note, however, that this interpretation is based on numerical evidence only where a non-trivial behaviour of the impurity velocity always coincided with soliton emission.

Truncated-Wigner approximation of a harmonically trapped polaron
The above mean-field analysis has neglected quantum fluctuations. Although it has been shown in [35,36] that the ground state properties of Bose polarons and bipolarons are well described by the modified mean-field approach even in the limit of strong boson-impurity coupling g IB ≫ gnξ, provided the Bose-Bose interaction is weak, i.e. if γ ≪ 1, it is not clear if this still holds in the non-equilibrium case. For this reason we now consider the effect of small quantum fluctuations using a   7). The evolution of the Bose gas densities at the marked positions is shown in Fig. 1 b)-d). The colored dashed line shows the velocity-momentum relation of the stationary state. It agrees well with the simulated result for sub-critical initial velocities v(0) < v c . For higher total momenta emission of density waves or solitons leads to initial friction forces, causing the deviation.
truncated Wigner approach [37][38][39]. This approach fully captures the influence of (deformed) Bogoliubov phonons on top of the condensate in quadratic (Bogoliubov) approximation. However, as stated by the Mermin-Wagner-Hohenberg theorem [45,46] there is no true Bose condensation in homogeneous 1D gases, which manifests itself by infrared divergencies when considering lowest-order quantum fluctuations. The latter also holds for finite systems with periodic boundary conditions. Thus in order to describe quantum fluctuations we can no longer approximate the one-dimensional gas as being homogeneous and have to take into account the presence of a harmonic trapping potential. (The different regimes of quantum degeneracy in trapped 1D Bose gases are discussed e.g. in [47].) This complicates the theoretical description as the Lee-Low-Pines transformation, conveniently used in homogeneous systems, no longer leads to a decoupling of the total momentum of the polaron. We will show however that the total momentum obeys a simple equation of motion if also the impurity is trapped.

Lee-Low-Pines Hamiltonian of a trapped 1D Bose gas
We start by adding a harmonic potential with frequency ω for bosons and Ω for the impurity to the Hamiltonian Eq. (1) to avoid infrared divergencies in the Bogoliubov theory of boson-boson interactions. Since these potentials break the translation invariance of the system the total momentum of the system is no longer conserved. Therefore, the transformation into a relative and a center of mass coordinate by the LLP transformation does not eliminate the impurity operators. Nevertheless, it is useful to apply the transformation, sincep andr only appear up to quadratic order in the LLP To simulate the time evolution of the system we derive the Heisenberg equation of motion (x, t). The advantage of the LLP transformation even in the case of harmonic trapping becomes clear here since the equations forp(t) andr(t) are formally solvable and we get Here N is the particle number of bosons andX B (t) = dx xφ † (x, t)φ(x, t)/N their center of mass position. In the case of equal trapping frequencies ω = Ω, the last term in the solution of the total momentump(t), Eq. (9), vanishes, and its time evolution corresponds to that of an uncoupled harmonic oscillator p(t) =p(0) cos Ωt + 1 Ωṗ (0) sin Ωt, for Ω = ω.
The remaining equation for the bosonic fieldφ(x, t) then reads:

Truncated Wigner simulation of the boson field
We now solve the system of Eqs. (9), (10) and (12) in a limit where quantum fluctuation of the impurity position and momentum are negligible, but include fluctuation of the Bose field using a truncated Wigner phase space approach (TWA). For reviews on the TWA methods see [38,39]. To be able to apply a TWA we have to first treat the quantum evolution ofr(t) andp(t). For this we apply another approximation and replace the total momentum and position operator of the Bose field in the dynamical equations of the impurity by expectation valueŝ This approximation is reasonable since the number of bosons is large N ≫ 1 such that fluctuations of their center of mass coordinate are small. From Eqs. (9) and (10) then follows, that the fluctuation ofp andr do not grow in time. It is possible to prepare the impurity in a state where fluctuation are small as long as its harmonic oscillator length scale l I = 1/ √ ΩM is small when compared to all other length scales of the system, especially the healing length of the condensate ξ = 1/ √ 2gnm, where n is the peak density of the ground state of the trapped Bose gas without an impurity. The semiclassical treatment of impurity position and total momentum is therefore justified only in the regime where gnm ΩM ≪ 1. Under this condition the operatorŝ r → ⟨r⟩ := r, andp → ⟨p⟩ := p (14) are replaceable by expectation values. In order to calculate the time evolution of the Bose field in a Wigner phase space description, an expression for the initial ground state of the Bose gas including quantum fluctuation is needed. Since we consider a weakly interacting Bose gas, it suffices to do this by replacing the field operatorsφ(x, t = 0) by the mean-field ground state ϕ 0 (x) of a trapped Bose gas and add quantum fluctuation within Bogoliubov-de Gennes (BdG) approximationφ Here u n (x) and v n (x) are the BdG coefficients of the trapped Bose gas andb n the phonon operators of the respective modes, see Appendix C for more details. The position r appears in this expression since the Bose gas ground state is transformed into the LLP frame, corresponding to the shift by r. In the Wigner phase-space description, the phonon operators are replaced by stochastic c-numbersb where all stochasticity is in the initial state. Since this state is the phonon vacuum they are set to Gaussian random variables with mean and variance given by ⟨β n (0)⟩ = 0 and ⟨β n (0)β * m (0)⟩ = 1 2 δ n,m , corresponding to a virtual occupation of half a phonon per mode on average. By symmetric ordering of the Heisenberg equation Eq. (12) and replacing the operators by c-numbersφ → ϕ we get the c-number equation of motion Here n w (x) is the virtually occupied particle density due to the Wigner description and is given by Note that we have truncated the number of modes taken into account, which is necessary in TWA, since if all BdG modes are included n w (x) would diverge n w (x) = 1 2 δ(0). The truncation of higher modes is commonly used in TWA simulations of trapped gases [37,38] and is physically justified as quantum fluctuations of high frequency modes can be neglected. We simulate Eq. (16) multiple times for different initial conditions and average about the different realization. In order to obtain expectation values all operators need to be symmetrically ordered first, e.g. the Bose gas density is given by where r(t) appears, such that the expression describes the density in the laboratory and not LLP frame. Figure 10. TWA simulation of a trapped system for different impurity-Bose coupling constants. Here trapping frequencies for bosons and impurity are equal Ω = ω = 0.3gn, the Tonks parameter is γ = 0.1, the mass ratio is M/m = 10 and the initial impurity momentum is p = 0.75M c. The upper panels show the time evolution up to t = 25/gn = 7.5/ω, where the dashed line marks the position of the impurity. The lower panels show the final density (blue) and compare it to the stationary state of the homogeneous system (dashed) modulated by the initial mean-field density of the trapped gas. The TWA simulation is averaged over 5000 noise realizations.
The time evolution of the trapped system is illustrated in Fig. 10. A problem of this approach is that the system is finite, so that density waves created by the impurity oscillate in the trap and return back to the impurity. The system does therefore only reach locally a stationary state in very shallow traps with very small trapping frequency ω ≪ gn. This however conflicts with the condition of a classical impurity gnm ΩM ≪ 1 for reasonably heavy impurities and equal trapping potential. However, although the system is not in a stationary state the stationary solution described in Section 2.2 agrees reasonably well with the observed density distribution, see Fig. 10 when applying a local density approximation.
Next, in order to test the validity of the mean-field approach used in Sections 2 and 3, we compare the time evolution of the trapped system obtained from TWA calculations to a mean-field simulation. For the latter, we set the initial virtual particle occupation to zero β n = 0. The impurity velocity v(t) = p(t) − P B (t) /M and position r(t) obtained in that way are shown in Fig. 11. For a small coupling constant g IB ≪ gnξ the evolution remains sinusoidal, however, at large coupling it deviates strongly from an harmonic motion. The agreement between TWA and mean-field is reasonably well, given that we consider a quite strongly interacting gas with γ = 0.1. The impurity position Fig. 11b) gets in some cases an overall shift between TWA and mean-field, the motion is however qualitatively similar. Especially the deviation in the impurity velocity is small, see Fig. 11a). From this, we conclude that the mean-field simulation is sufficient to predict the time evolution at least qualitatively.  Figure 11. a) Impurity velocity and b) impurity position in a trapped system with equal trapping frequencies Ω = ω. Solid lines show the TWA simulation, dashed lines the classical mean-field approximation. The TWA simulation is averaged over 5000 noise realizations.

Summary
In the present paper we have discussed the dynamics of the formation of a Bose polaron when an impurity is injected into a 1D weakly-interacting Bose gas by performing time-dependent simulations of mean-field equations of the condensate amplitude complemented by truncated Wigner simulations to include quantum fluctuations. In a homogeneous gas with periodic boundary conditions the total momentum p of the system is conserved and can be used as independent parameter to characterize different dynamical regimes. Analyzing steady state solutions of the mean-field equations first, we showed that stationary solutions exist only for impurity velocities below a critical value, which in the limit of weak impurity-boson couplings g IB agrees with the Landau critical velocity, as predicted by the Froehlich model, but monotonously decreases with increasing interaction and eventually approaches zero. This is because with growing values of g IB the condensate is more and more depleted in the vicinity of the impurity, which leads to a reduced local speed of sound. While, as first shown in [32] for a given velocity of the impurity below the critical value, there are always two stationary solutions of the condensate equations, the solution is unique when fixing the total momentum. For momentum values with a convex energy-momentum relation, ∂ 2 E/∂p 2 > 0, one of the two solutions applies and in regions with ∂ 2 E/∂p 2 < 0 the other solution holds. We showed moreover that the stationary energy-momentum relation is periodic in p as the background condensate away from the impurity can pick up additional quantized amounts of momentum corresponding to integer windings of the condensate phase over its length. As a consequence the relation between impurity velocity v and polaron momentum p is also periodic and includes regions of momentum where the impurity velocity is negative. In these regions the Bose gas stabilizes only a steady state with momentum exceeding the total momentum which must then be compensated by an opposite motion of the impurity. While a direct measurement of the energy-momentum relation of the polaron is challenging, its non-monotonous form can have interesting experimental consequences. E.g. injecting impurities with finite velocity into a small ring condensate can induce a finite circular current corresponding to a finite number of enclosed flux quanta.
To study the formation of the polaron we considered two cases, a slow, quasiadiabatic turning on of the Bose-impurity coupling and a sudden quench. If in the quasi-adiabatic situation the initial velocity is chosen small enough such that it stays below the critical value at all times, the impurity is decelerated only due to the increase in its effective mass, associated with the formation of the polaron. As the system evolves quasi-adiabatic this reduction in velocity is reversible. We find that the polaron quasiparticle is formed on the timescale of the inverse chemical potential 1/gn (see Fig. 8), which for parameters of a recent experiment in 1D gases [4] is on the order of 60µs. If in the quasi-adiabatic scheme the initial velocity is above the critical value the impurity emits density waves irrespective how slowly the interaction is turned on leading to irreversible friction. Switching on the impurity-boson interaction suddenly a rich scenario of dynamical regimes is observed. Depending on the mass ratio of particles, the total momentum and the impurity-boson coupling strength the impurity is slowed down by emission of density waves or grey solitons. The latter happens for large momenta and large impurity masses and is specific for the regime of strong impurityboson coupling. In this case asymptotic states can form where the impurity velocity changes its sign, i.e. backscattering occurs, which cannot occur in the weak coupling regime dominated by Cherenkov radiation of phonons. While an in-situ measurement of the impurity motion is difficult in an experiment, the emission of grey solitons can be directly observed by density measurements of the condensate. We here considered impurities in one-dimensional condensates. The modified mean-field approach including the backaction of the impurity to the condensate can however be applied also to higher dimensions. Theories predicting the polaron dynamics based on the Froehlich model, using e.g. a coherent variational ansatzes [28][29][30] or master equation [24,31] are capable of capturing the evolution as long as the condensate deformation is not substantial. From a straightforward dimensional analysis of the mean-field equation in D dimension we estimate that the condensate deformation becomes significant for g IB /g ≳ nξ D . Some of the predicted effects are expected to carry over from one to higher dimensions. E.g. the reversible slowdown of a sub-sonic impurity due to the formation of a polaron and the friction forces experienced by a super-sonic impurity due to emission of density waves will be very similar. The emission of grey solitons and the dragging of the impurity towards the grey solitons, on the other hand, is an effect specific to one-dimensional gases. In two dimensions a heavy, supersonic impurity might instead emit vortex antivortex pairs.
To justify the validity of the mean-field approximation we performed truncated Wigner simulations of the full quantum problem in a trapped gas. The TWA accounts for quantum fluctuations due to Bogoliubov phonons on the deformed condensate background up to quadratic order. To avoid infrared divergencies related to the onedimensional setup, enforced by the Mermin-Wagner-Hohenberg theorem, we considered a harmonically trapped gas. Although the total momentum is no longer conserved it follows a simple equation of motion, which we solve in semiclassical approximation. The TWA simulations show that the mean-field description of the dynamics of polaron formation is well justified as long as the Tonks parameter of the Bose gas is small, i.e. for a weakly interacting gas. The case of strong boson-boson interactions requires different analytical and numerical tools and will be discussed elsewhere.
as long as the system size L is large compared to the rescaled healing lengthξ = 1/ √ 2gnm. Here a = v/c, b = √ 1 − a 2 and the parameters φ 1 and φ 2 are chosen such that the phase of the solution is continuous at x = 0 and fulfills periodic boundary condition. The parameter d shifts the grey soliton wave function such that the boundary condition generated by the delta distribution in Eq. (3) is fulfilled. From this, it can be deduced that tanh d must be the solution of a cubic equation The three solutions of the equation are shown in Fig. A1 b). The one which is real for all parameters is always less than or equal −1, such that d is not a real number corresponding to a nonphysical state. The other two solutions are real and between 0 and 1, if the impurity velocity v is below the critical velocity v c = a cc , see Fig. A1 a).
Here a c can be determined by solving which is equivalent to a cubic equation in a 2 c . For a small coupling constant g IB ≪ gnξ the critical velocity isc and agrees with the prediction of the Froehlich model [24,43]. However for strong repulsion g IB ≫ gnξ the critical velocity converges to zero. Substituting the two physical solutions of Eq. (A.4) into Eq. (A.2) yields the two stationary states mentioned in the main part of this work. The two solutions are equal at the critical momentum, explaining why the stationary states merge at criticality.
Next, since p and not v is conserved under time evolution it is important to derive an expression relating the parameters for the stationary solutions. It follows from Eqs. (A.1) and (A.2) and is given by In order to compare this analytic expression to the time-dependent simulation we solve it numerically for v. The polaron energy is given by E pol = E(g IB , p)−E(g IB = 0, v = 0), where E(g IB , p) is the expectation value of the LLP Hamiltonian Eq. (2) (A.7) In the approach described so far Eq. (A.6) results only in momentum values with −πn ≤ p ≤ πn. In order to reach higher momenta the stationary solution Eq. (A.2) must be modified by an additional phase gradient ϕ(x, t) = e ix 2π L ν ϕ(x, t) with ν ∈ Z.
(A.8) picks up an additional term, which explains why the observables in Fig. 3a) and Fig. 4b) are periodic in p, with a period length of 2πn.

Appendix B. Gapless adiabaticity
In Section 3.1 we showed that the system evolves quasi-adiabatic if the impurity-Bose coupling constant is turned on slowly compared to the other timescales. However, there always remains a small but finite difference to the instantaneous stationary state in Fig. 5c). In this section we show, that this difference originates from the system not being energetically gaped in the thermodynamic limit. To this end, the time evolution of a large system L ≫cT is compared to a small one L ≪cT , with otherwise equal parameters. The small system is gaped due to finite-size effects, such that the adiabatic theorem strictly holds. This is shown in Fig. B1 a)-c), where the impurity momentum as well as the density and phase of the Bose gas agree with the instantaneous ground state. In contrast in the large system. Here in particular the phase disagrees at a large distance from the impurity |x| ≫ξ, see Fig. B1 d), explaining the small discrepancy of the impurity momentum Fig. B1a). In the large system the stationary state is not reached globally, but only locally at the position of the impurity. This is however sufficient for the system to evolve quasi-adiabatic. Figure B1.
Comparison of the quasi-adiabatic (large system) to true-adiabatic (small system) evolution for v(0) = p/M = 0.1c, g IB = gnξ, γ = 0.1 and M = 3m. a) Evolution of the impurity velocity for two turn-on timescales T . Dashed lines represent the instantaneous stationary state. b) -d) State at the end of the evolution at t gn = 2 · 10 4 for T gn = 3 · 10 3 . For a small system L = 100ξ, both phase b) and density c) agree with the instantaneous stationary state (dashed). For a large system L = 5 · 10 4ξ the phase d) disagrees at a large distance from the center.

Appendix C. Bogoliubov-de Gennes in a trap
In order to express the initial ground state of the trapped Bose gas, before the interaction with the impurity, we diagonalize the Bose gas Hamiltonian approximately using a Bogoliubov-de Gennes (BdG) approach [48]. In the first step, the mean-field ground state is determined by the Gross-Pitaevskii equation (GPE) − ∂ 2 x 2m + 1 2 mω 2 x 2 + g|ϕ 0 (x)| 2 − µ ϕ 0 (x) = 0, (C.2) which we solve numerically using imaginary time evolution. Here µ is the mean-field chemical potential. In case of a weakly interacting Bose gas, it is sufficient to only include small fluctuation on top of the mean-field solution, which is done by expressing the bosonic field operators bŷ ϕ(x) = ϕ 0 (x) + n u n (x)b n + v n (x) * b † n , (C. 3) and only keep terms up to quadratic order in the operatorsb where ϵ n are the eigenenergies of the corresponding BdG modes. To solve this equation, we expand it in a finite number of eigenfunctions of the free harmonic oscillator and diagonalize the resulting matrix numerically.