{'fraction_non_alphanumeric': 0.029926984425433237, 'fraction_numerical': 0.008429695089467582, 'mean_word_length': 4.507766655160511, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Real-space dynamics of attractive and repulsive polarons in Bose-Einstein condensates

We investigate the formation of a Bose polaron when a single impurity in a Bose-Einstein condensate is quenched from a non-interacting to an attractively interacting state in the vicinity of a Feshbach resonance. We use a beyond-Fr\"ohlich Hamiltonian to describe both sides of the resonance and a coherent-state variational ansatz to compute the time evolution of boson density profiles in position space. We find that on the repulsive side of the Feshbach resonance, the Bose polaron performs long-lived oscillations, which is surprising given that the two-body problem has only one bound state coupled to a continuum. They arise due to interference between multiply occupied bound states and therefore can be only found with many-body approaches such as the coherent-state ansatz. This is a distinguishing feature of the Bose polaron compared to the Fermi polaron where the bound state can be occupied only once. We derive an implicit equation for the frequency of these oscillations and show that it can be approximated by the energy of the two-body bound state. Finally, we consider an impurity introduced at non-zero velocity and find that, on the repulsive side, it is periodically slowed down or even arrested before speeding up again.


I. INTRODUCTION
The polaron is a general concept of many-body physics that naturally arises in different fields like solid state physics and the theory of ultracold gases. While it has long been used to describe electrons in a crystal lattice, only recent experimental advances allowed one to realize polarons in ultracold gases. Here, the Feshbach resonance allows for a high level of control and in particular for realizing the strong-coupling regime, which could not be done before. This gives access to interesting phenomena such as self-localization and bubble formation. Moreover, the interaction can be changed abruptly which allows for the investigation of the dynamics. Combined with the possibility of direct imaging, this allows to view polaron formation in position space, which is crucial for the physical intuition and interpretation of the time evolution.
The concept of polarons was originally invented by Landau [1]. He showed that an electron in a crystal lattice interacts with the surrounding atoms in such a way that it can be described as a quasiparticle with a higher effective mass, moving through free space. Describing the lattice deformations induced by the electron as phonons, the polaron can be imagined as an electron carrying a cloud of phonons around it. A very similar picture arises in ultracold bosonic gases: According to Bogoliubov theory, the elementary excitations of a BEC are phonons as well, so when an impurity is moving through the gas, the situation is analogous to that of an electron in a crystal. But in an ultracold gas, it is possible to tune the interaction between the particles via a Feshbach resonance and in particular to investigate the regime of strong coupling between impurity and host bosons.
A number of different theoretical approaches has been used to investigate different aspects of the Bose polaron. In 1954, Fröhlich introduced a Hamiltonian which is commonly used to study polarons [2]. It can be recovered from Bogoliubov theory with one further approximation [3]. This was first done for ultracold gases in [4], where the ground state properties were studied using a variational ansatz due to Feynman [5]. This ansatz works well for all couplings in the original case of electrons in a lattice, but in the ultracold gas, the regularization of the contact interaction leads to errors when the coupling becomes strong. This was discussed in [6] with Diagrammatic Monte Carlo calculations. These give access to the ground-state properties and are computationally intensive but numerically exact and provide valuable benchmarks for other methods. A coherent-state variational ansatz originally due to Lee, Low and Pines (LLP) [7] has been used to study dynamical properties [8,9]. It neglects entanglement in momentum space and is considered best for heavy impurities and weak couplings. More quantum fluctuations have been taken into account by a renormalization group technique for the ground state [10,11] and the dynamics [12] and by the correlated gaussian wave function ansatz [13], as well as a Hartree-Fock-Bogoliubov description [14]. There are some more works related to the Fröhlich Hamiltonian [15][16][17]; for a review, see [18]. The Bose polaron exhibits characteristic signatures also at finite temperature [19,20].
The interaction term in the Fröhlich Hamiltonian is, however, just an approximation in the case of ultracold gases and higher order terms become important in the regime of strong coupling. This was first observed in [21], where a T-matrix approximation was used. Subsequently, both the LLP coherent state ansatz [9] and the RG analysis [22] were applied to the fully interacting model, and recently the functional determinant approach [23]. In the one-dimensional case, some analytical results for heavy impurities have been obtained [24] and phononphonon interactions beyond Bogoliubov theory have been considered [25].
Approaches not based on Bogoliubov theory are more limited in number: Quantum Monte Carlo calculations [26] provide exact ground states for a limited number of parameters. Coupled Gross-Pitaevskii equations [27][28][29] can describe the spatial deformation of the BEC and the phenomena of self-localization and the bubble polaron but work on a mean-field level. A variational approach which treats the molecular state as an independent quasiparticle has been used to investigate three-body bound states [30,31].
Experimentally, Bose polarons in ultracold gases have been observed with a focus on absorption spectra and decoherence [32][33][34][35][36] for which some theoretical predictions have been made. Direct imaging experiments on the other hand are still in preparation and there have been few theoretical results concerning the real space dynamics of the bose polaron: the Monte Carlo calculations in [26] include density profiles but only statically for ground states while the Gross-Pitaevskii method in [29] considered a repulsive interaction. This is different from an attractive interaction with a positive scattering length in that it does not feature a bound state.
In this paper, we investigate the dynamics of polaron formation when an initially non-interacting impurity is quenched to an attractively interacting state. This situation has been studied before to compute radio-frequency absorption spectra [8,9] and, on the attractive side of the Feshbach resonance, polaron trajectories [12], as well as pre-thermalization dynamics [37]. Here, we focus on two new aspects: We compute the density profile of the BEC around the impurity as a function of time and thus view the formation of the polaron in position space. This can be directly measured with current imaging technologies and corresponding experiments are in preparation. On the other hand, we investigate the repulsive side of the Feshbach resonance where the scattering length is positive. Here, a two-body bound state exists and its interplay with the polaron leads to new effects, in particular characteristic oscillations and a depletion of the boson density in a halo around the impurity. These were inaccessible to many previous works based on the Fröhlich Hamiltonian, which depends only on the modulus of the scattering length and cannot describe the bound states. Our study, instead, uses the extended Hamiltonian including higher-order terms in the interaction, and we employ the LLP coherent-state ansatz to compute the dynamics.
The paper is organized as follows. In Sec. II, we review the construction of the Hamiltonian and the variational ansatz starting from Bogoliubov theory and discuss the stationary solution. Section III contains the results for the time evolution after a quench: we start with the case of an impurity initially at rest and compute boson density profiles as well as the total number of bosons gathering around the impurity. We then present a way to compute characteristic oscillation frequencies in the repulsive regime. Finally, we investigate the influence of a non-zero initial velocity and compute polaron trajectories.

II. MODEL
Our starting point is the Hamiltonian of a single impurity in a bath of bosons where m I and m B are the masses of impurity and bosons, p I andx I the impurity momentum and position operators and a ( †) k the bosonic creation and annihilation operators. n B (x) = a † x a x is the boson density and V BB and V IB are the boson-boson and impurity-boson interaction potentials. Our derivation follows Shchadilova et al. [9].
Since we are dealing with just one impurity, it is convenient to go to relative coordinates. This is achieved by the canonical transformation exp (iS) where It is known as the Lee-Low-Pines (LLP) transformation [7], see also [3]. Its effect on the operators is Note that formallyp =p I , but we have dropped the index after the transformation since the physical meaning is not the impurity but the total momentum. It is, of course, conserved and can be replaced by the initial impurity momentum p 0 such that the transformed Hamiltonian reads This transformation has simplified the interaction term and replaced the impurity momentum by the difference of total momentum and phonon momentum. Here, fourthorder terms in the phonon operators appear unless the impurity is taken to be infinitely heavy, i.e., stationary. A delocalized impurity thus induces effective interactions between the phonons.

Bogoliubov Theory
We use Bogoliubov theory which pre-supposes Bose-Einstein-Condensation in the k = 0 mode and approximates the low-temperature behaviour by discarding terms in 3rd and 4th order of boson operators with k = 0. The resulting bosonic part of the Hamiltonian is diagonalized by the Bogoliubov transformation b † k = cosh (ϕ k ) a † k − sinh (ϕ k ) a −k with exp (4ϕ k ) = ξ 2 k 2 /(2 + ξ²k²) and the healing length Up to a constant energy offset, with phonon dispersion a BB and a IB are the scattering lengths of the potentials V BB and V IB . Finally, n 0 is the condensate density, which is a free parameter in Bogoliubov theory. means that the sum runs over k = 0.

Contact interaction
In a dilute ultracold gas, the range of interactions is small compared to all other length scales. The effect of the interaction can therefore be described by a single number, the scattering length, while the precise shape of the potential does not matter and can be chosen arbitrarily. The most convenient choice is a zero-range pseudopotential. Taken literally, the Fourier transform of a delta potential would correspond to a potential in momentum space that is constant over an unbounded region, which does not make sense. Instead, one constructs it as a scaling limit, first cutting off all momentum sums at some large Λ, and then tuning the interaction strength in the limit Λ → ∞ such that the scattering length remains fixed at the desired value. The correctly regularized interaction strength is then given by But one must carefully note that a 3d delta distribution in real space is not the proper choice as it would not lead to a well-defined scattering length. In momentum space, this translates to the fact that when the potential is taken constant over an infinite region, it has to be multiplied by a properly chosen infinitesimal prefactor to yield the desired scattering length. Practically, this is done by cutting off all momentum sums at some large value Λ. The correctly regularized interaction strength is then given by with m red being the reduced mass of impurity and bosons, m red = m −1 I +m −1 B . Note that such a cutoff effectively corresponds to an interaction with a non-zero range of order 1/Λ. The cutoff Λ will be used implicitly in all sums and integrals throughout the paper. Also note that instead of a "hard" cutoff, one can also multiply the integrands by a decaying function such as exp −2k 2 /Λ 2 . This leads to smoother results when the cutoff is not large enough for perfectly converged behaviour.
In eq. (1), we still have to express the boson operators a k by phonon operators b k . The result is where W k = exp (ϕ k ). Here, N 0 = n 0 V is the number of condensed bosons and we have approximated 0|N |0 ≈ N 0 , i.e. neglected the ground state depletion, which gives a constant density shift [38] Inserting (2) into H Bog , we obtain the final Hamiltonian The first two lines of (3) correspond to a Fröhlich Hamiltonian, which is often used to study polarons. It has also been used for polarons in ultracold gases, even though for such systems the quadratic terms in the last two lines of eq. (3) are present. The so obtained results are still valid as long as the coupling between impurity and host atoms is sufficiently weak. One needs to take care however, that when using the Fröhlich Hamiltonian, the regularized contact interaction g IB may not be used and needs to be replaced by the result from the Born approximation g F r IB = 2πa IB /m red . Near the Feshbach resonance, the quadratic terms become important as pointed out in [9,21].

Coherent state ansatz
Also by Lee, Low and Pines [7], a variational ansatz for the Fröhlich model was suggested, which approximates the ground state by a coherent state. In the limit of infinitely heavy impurities, this ansatz becomes exact. In [8], a time dependent version of this ansatz has been applied to the Bose polaron described by a Fröhlich Hamiltonian. Subsequently, the same time-dependent ansatz has been applied to the full Hamiltonian (3) in [9].
Specifically, one considers wave functions of the form and projects the Schrödinger equation onto the submanifold spanned by these functions. Equivalent to this projection is the minimization of the functional with respect to the α k . Setting up the Euler-Lagrange where The initial value α k (0) = 0, i.e., |α(0) = |0 , corresponds to the situation of a quench. If the impurity is initially at rest, p 0 = 0, then p I [α] = 0 for all times due to spherical symmetry. In this case, the equation becomes R-linear and can be written in the form with a constant offset α (s) (the stationary solution, see below).

Stationary solution
Before turning to the dynamical solutions, it is instructive to look at the stationary solution obtained froṁ α k = 0 or equivalently ∂ α|H|α where C 1 and p I are determined by the implicit equations Note that these quantities are UV convergent: For large k, one has W 2 k = 1+O(k −2 ) and Ω k = k 2 2m red +O(1). The integrand is thus 1 Ω k + k·pI . The first term cancels with the divergence of g −1 IB , the second vanishes by antisymmetry. The momentum integrand is k and again, the first term is antisymmetric.
The integrals exist if and only if p I /m I < c, i.e., the stationary impurity velocity must always be below the speed of sound c = 1/ √ 2m B ξ. For too large initial momenta, no stationary solution exists (the same is true in the Fröhlich model where C 1 = 2πa IB √ n 0 /m red is a constant, see [8]).
In the special case p 0 = 0, one obtains p I = 0, In the repulsive case, the system keeps oscillating between the two limit curves that are shown. Calculations were carried out at a gas parameter of n0a 3 BB = 10 −5 , equal masses of bosons and impurity mB = mI and with the impurity initially at rest, p0 = 0. A soft momentum cutoff e −2k 2 /Λ 2 with Λ = 80ξ −1 was used.
Here, a + is one of two critical scattering lengths defined by which satisfy a − < 0 < a + . They were found in [22] and delimit three regions of different stability of the stationary solution. This will be reflected in the convergence behaviour of observables in our dynamical analysis: 1. a −1 IB < a −1 − : The stationary point is a minimum, observables converge. This is the region where the attractive Bose polaron is expected to form.
2. a −1 − < a −1 IB < a −1 + : The stationary point behaves like a saddle point, the system is dynamically unstable. This is well understood as coming from phase fluctuations growing without bounds [22]. This behaviour seems unphysical and means that the approach cannot cover very strong couplings, because one of the two approximations involved fails: The Bogoliubov approximation, which neglects interactions between phonons, or the coherent state ansatz, which does not include all quan-tum fluctuations.
3. a −1 + < a −1 IB : The stationary point behaves like a maximum, observables are oscillating. In this region, the stationary solution is usually interpreted as a repulsive polaron due to its positive energy, while at negative energies, a molecular state is expected. We will discuss to what extent the system can be described as oscillating between these two states.
In Fig. 1 we show how the boundaries of the three regimes change with the reduced mass. For the case of light impurities, m red → 0, the unstable region grows. Here, quantum fluctuations become especially important as the impurity is delocalized. One may thus suspect that the instability arises from the failure of the LLP ansatz to describe these fluctuations correctly. However, even in the case of an infinitely heavy impurity, m red → m B , then unstable region remains. This suggests that it is rather the Bogoliubov approximation that fails to describe the strongly interacting region correctly. It must be noted however, that even for m I → ∞, the coherent state ansatz is an approximation since we are not dealing with the Fröhlich but with the fully interacting Hamiltonian.

A. Time evolution of density profiles
Fourier transforming the numerical solution of (4) back to position space, we can compute the boson density at distance r from the impurity. More precisely, since the impurity is itself a quantum particle, the quantity to consider is correlation function Expressing the boson density by phonon operators and applying the variational ansatz, n(x) takes the following form in terms of the coefficients α k : where denotes the transformation to position space. Figure 2 shows the results for the three different regimes: 1. Attractive regime: Bosons are gathering around the impurity and the profile quickly converges to form the attractive Bose polaron. The final shape matches precisely that of the stationary solution. 2. Unstable regime: The impurity keeps pulling in more and more bosons. As mentioned before, this unphysical behaviour reflects the failure of the Bogoliubov approximation when interactions are too strong. Including phonon interactions, i.e., higherorder terms in the bosonic operators a k , might prevent this.
3. Repulsive regime: Close to the impurity, the boson density is strongly increased but there is a halo of reduced density around it. There is no convergence to a ground state profile, but instead, the solution keeps oscillating between two states of the coupled system of impurity and surrounding condensate: At some times, the bath is completely depleted at a certain distance while at other times, there is still about half the original density left.
Comparing these results with the quantum Monte Carlo calculations of the ground state profile in [26], the results are qualitatively similar, even though quantitatively slightly different (our parameters correspond to a IB /a BB = −20.94, ∞ and 20.94).
The complete depletion is a feature that is present even in the stationary solution and for all scattering lengths above a + : Since α (s) is real, the last term in (7) vanishes and one can always find an r so that the first term vanishes as well. The length scale on which the depletion takes place is given by the scattering length, as shown in Fig. 3. This is not surprising since this is the scale of the two-body bound state. The return to the condensate density then happens on the order of the healing length (not shown in the figure).

B. Boson Number
From the momentum space coefficients, we can compute the total change in the number of bosons. This is not zero because the Bogoliubov theory does not preserve particle number. It can be seen as a measure of how many particles the impurity attracts in total. The formula is Results are shown in Fig. 4a. The characteristics of the three regimes -convergence, unbounded growth and oscillations -are clearly visible. Note that in the repulsive case, the maxima of the boson number correspond to the more extreme density profiles, i.e., those with full depletion. Doing the same computations for many different scattering lengths, we arrive at Fig. 4b. As the critical scattering lengths a + and a − are approached from the attractive or repulsive regime, the total boson number grows rapidly as well as the time to convergence. In the figure, this is indicated by the fact that in these areas, the curves are still washed out, therefore not yet converged.

C. Oscillation Frequencies
In the case of an impurity initially at rest, the frequencies of the oscillations in the repulsive regime can be predicted by making an ansatz for the long-time solution. The coefficients C 1 and C 2 from the differential equation (4) show the same qualitative behaviour as the other observables: convergence, divergence or oscillations, ac-cording to the regime. We therefore make an asymptotic ansatz where the coefficients λ can take finitely many complex values with Re λ ≥ 0. This covers all of the three cases: convergence if only λ = 0 is present, exponential growth if a λ > 0 exists and oscillations for imaginary λ. The case Re λ < 0 would be interesting as well to describe the speed of convergence, but the restriction to Re λ ≥ 0 will prove necessary for the calculation. Since C 1 and C 2 are real, we must have A λ = A λ and B λ = B λ (the bar denotes complex conjugation). Our aim is to derive conditions on λ to be able to predict the exponential growth rate or oscillation frequency of the physical observables. We thus insert the ansatz into the differential equation (4) with p I [α] = 0 and find the solution and unknown s k , which depend on the full history of the time evolution.
Re-inserting the expression for α into the definitions of . Predictions for exponential growth rate Re λ and oscillation frequency Im λ in the long-time limit for the unstable and repulsive regime. The oscillation frequency is well approximated by 1/2m red a 2 IB , the energy of the two-body bound state. As the small figures show, this remains true for different mass ratios (note the difference in vertical scale). Parameters are as in Fig. 2 with a hard momentum cutoff at Λ = 1000ξ −1 .
These equations should be regarded only as determining the solution asymptotically because the integrals over e −iΩ k t will decay while no finite sum of exponentials e λt with all Re λ ≥ 0 can ever be decaying. But since we are interested in the long-time limit, we can ignore the oscillatory integrals. The need to drop these terms simply reflects the fact that a system never exactly reaches its asymptotic solution but only comes arbitrarily close.
We want to use the linear independence of e λt with different λ, but first, the Re and Im need to be expanded as 2Re λ b kλ e λt = λ (b kλ e λt + b kλ eλ t ) = λ (b kλ + b kλ )e λt and similarly for Im . We then find λ = 0 : Recall that A λ and B λ need not be real even though the sums in (9) are. Also note that these equations would not be true for Re λ < 0 because the oscillating integrals could not be ignored. Using the expressions for b kλ and the relations we arrive at Written in this way, all integrals are UV convergent. Multiplying the last two equations finally leads to an equation for λ 2 , independent of A λ and B λ : In the above expressions, many of the integrals do in fact not exist if λ is purely imaginary since the integrands have a pole at k = k c in this case. What does exist, however, are the Cauchy principal value (PV) integrals  Such integrals are not invariant under coordinate transformations because the way in which the pole is approached is crucial, so it is not immediately clear how to make sense of (10) for negative λ 2 . This becomes clearer if, instead of making an ansatz for an asymptotic solution, one considers the time evolution operator in (5).
In fact, the product H (2) H (1) determines the dynamics completely, so it is sufficient to compute its spectrum. One obtains, once again, equation (10), where now λ 2 are the eigenvalues of H (2) H (1) . But in the case λ 2 < 0, one finds that (10) must hold with the integrals replaced by PV integrals in any choice of coordinates, as long as the same is used in all three integrals. Therefore, a coordinate transformation will change the value of the individual integrals, but when it is applied to all of them, the equation must stay true. Solving (10) numerically for different parameters, we find that it has 1. no solution in the attractive regime, so only λ = 0 is possible here; 2. one solution for positive real λ 2 in the unstable regime; 3. one solution for negative real λ 2 in the repulsive regime.
Using PV integrals in a particular choice of coordinates, one may find an additional negative real solution in cases 2 and 3, but they are not valid because they change when different coordinates are used.
The correct values, on the other hand, give predictions of the exponential growth rate or frequency, respectively, which are in perfect agreement with the simulations. Figure 5 shows λ over a range of different scattering lengths. The oscillatory behaviour in the observables is reminiscent of oscillations between two quantum states with a frequency corresponding to their energy difference. Here, Im λ is very well matched by the energy of the two-body bound state, 1/2m red a 2 IB . For weak interactions, this is also the energy difference between the repulsive polaron and the molecular state since the former approaches zero. To describe the repulsive polaron at stronger coupling, we consider a variational state proposed in [9] where a single phonon is added to the stationary solution ψ (s) , i.e., k γ k b † k ψ (s) . We find that it has an energy very close to E (s) − 1/2m red a 2 IB . This suggests that the additional particle does indeed form a molecule with the impurity and that neither the interaction with the other phonons nor the difference between a phonon and a boson plays a big role. We thus believe that the oscillations occur between two quantum states, one of which we interpret as the repulsive polaron state and the other one as an interacting molecular state. Note, however, that such an effective picture applies only to observables that integrate over momentum space and not to the full wave function itself because we discarded the highly oscillatory terms in the coefficients α.

D. Moving Impurity
We now turn to the case where the impurity has an initial velocity v 0 = p 0 /m I . This case has been investigated for the Fröhlich Hamiltonian in [8] using the LLP variational ansatz and in [12] with a time-dependent renormalization group method. In the latter reference, also the full Hamiltonian was investigated on the attractive side of the Feshbach resonance and it was argued that in this case, the second order terms lead only to a shift of the inverse scattering length. On the repulsive side, this is not true since the Fröhlich Hamiltonian cannot describe molecule formation.
In Fig. 6 we show the time evolution of the impurity velocity v 0 v I (t )dt (according to Ehrenfest's theorem). In the attractive regime, the behaviour is simple: If the total momentum is not too high, it converges to a non-zero final value, which agrees with the stationary solution. This matches the picture of a polaron with an increased effective mass. If the total momentum is too large such that no stationary solution exists, the velocity converges to the speed of sound. Note, however, that close to the resonance, it has been predicted that quantum fluctuations beyond the LLP ansatz lead to an enhanced damping or even recoil effects, cf. [12].
On the repulsive side, the behaviour is different: the velocity is oscillating with the same frequency as the density profile and boson number. We imagine this as being due to the fact that the bound state has a different effective mass from the polaron state such that oscillations are visible in the velocity. The effect is most striking close to the critical scattering length a −1 + : Here the impurity velocity quickly reaches zero but has periodic revivals. This is in accordance with the fact that the polaron state has a very high effective mass in contrast to the bound state.
On both sides of the resonance, the initial velocity does not matter much as long as it is below or close to the speed of sound: it leads only to a rescaling of the velocity at later times. This is also reflected in the density profiles (Fig. 7), which are still symmetric around the impurity. Above the speed of sound, however, the number of attracted bosons is increased, leading to a faster decay of the velocity. The density profiles now become asymmetric with some depletion in front of the impurity.

IV. DISCUSSION
We investigated the dynamics of polaron formation in a BEC after a quench, focussing on real space density profiles and the behaviour for positive scattering lengths. These could not be investigated in previous works that used the Fröhlich Hamiltonian.
We found that three regions of qualitatively different behaviour exist, where the strong-coupling region is unstable, as expected from the stationary analysis in [22]. The fact that the instability persists even in the limit of heavy impurities is a hint that Bogoliubov theory is not adequate to investigate the strong-coupling regime: the deformation of the BEC is too important for the Bogoliubov approximation to hold.
For positive scattering lengths, oscillations can be observed in the expectation values of many observables and we presented a way to compute their frequency. It is well approximated by 1/2m red a 2 IB , which is the energy of the two-body bound state. We thus expect that the system can be effectively described as oscillating between two quantum states, one of which might have an interpretation as a repulsive polaron and the other one as a bound state. This can be an interesting topic for future work. The challenge here is in particular to find two suitable extremal states and a physical interpretation for them.
Remarkably, these oscillations are present even in the impurity velocity, leading to striking "stop-and-go" polaron trajectories. The effect is most pronounced for strong coupling when oscillations are slow compared to the velocity relaxation: Here the position is advanced in steps. It will be interesting to see if this can be observed in experiments.
In position space, the positive scattering length leads to a halo of reduced condensate density around the im-purity, whose size corresponds to a scattering length and is independent of the mass. This is a version of a bubble polaron, where the impurity has, however, still a core of increased density around it and the profile oscillates in intensity. At certain times, the depletion is even perfect, which was not visible in ground state calculations.
It will be interesting for future work to investigate how the system behaves for strong coupling where the coherent state ansatz becomes unstable. However, suitable techniques still have to be developed. Within Bogoliubov theory, the so-called correlated gaussian wave functions are a promising way since they should become exact in the limit of heavy impurities. On the other hand, it will be important to find out in which region Bogoliubov theory is not reliable any more and how the system can be described in this region.