Collision of impurities with Bose-Einstein condensates

Quantum dynamics of impurities in a bath of bosons is a long-standing problem of solid-state, plasma, and atomic physics. Recent experimental and theoretical investigations with ultracold atoms focused on this problem, studying atomic impurities immersed in a atomic Bose-Einstein condensate (BEC) and for various relative coupling strengths tuned by the Fano-Feshbach resonance technique. Here we report extensive numerical simulations on a closely related problem: the collision between a bosonic impurity made of few $^{41}$K atoms and a BEC made of $^{87}$Rb atoms in a quasi one-dimensional configuration and under a weak harmonic axial confinement. For small values of the interspecies interaction strength (no matter the sign of it), we find that the impurity, which starts from outside the BEC, simply oscillates back and forth the BEC cloud, but the frequency of oscillation depends on the interaction strength. For intermediate couplings, after a few cycles of oscillation the impurity is captured by the BEC and strongly changes its amplitude of oscillation. In the strong interaction regime, if the interspecies interaction is attractive, a local maximum (bright soliton) in the density of BEC occurs where the impurity is trapped; instead, if the interspecies interaction is repulsive, the impurity is not able to enter in the BEC cloud and the reflection coefficient is close to one. On the other hand, if the initial displacement of the impurity is increased, the impurity is able to penetrate in the cloud leading to the appearance of a moving hole (dark soliton) in the BEC.


Introduction
In 1933 Landau introduced the concept of polaron, an electron whose effective mass is affected by the coupling with the quantized lattice vibrations (phonons) of a crystal [1]. Later, Frölich derived a field-theoretical polaron Hamiltonian that covers all coupling strengths between the electron and the phonons [2]. The basic properties of polarons are now established (see, for instance, the review [3]) but the interest in the polaron dynamics and, more generally, in the dynamics of impurities interacting with a bosonic bath has recently gone through a vigorous revival, mainly in the context of ultracold atomic gases. In Ref. [4], localized bosonic impurities, made of few 41 K atoms, have been created in a one-dimensional (1D) configuration and their interactions with a Bose-Einstein condenstate (BEC) of 87 Rb atoms have been investigated by using a Fano-Feshbach resonance to tune the impurity-boson scattering length. More recently tunable BEC impurities, i.e. atomic impurity in a cloud of ultracold Bose-Einstein condensed atoms, have been obtained by other two experimental groups [5,6]. The Bose polaron problem has been addressed theoretically by using different techniques: quantum Langevin equation [4], mean-field theory with coupled Gross-Pitaevskii equations [7,8], time-dependent variational mean-field for lattice polarons [9,10], Feynman path integral and Jensen-Feynman variational principle [11], T-matrix [12] and perturbation [13] approaches, variational wavefunction [14,15], and quantum Monte Carlo schemes [16]. Recently, the transition of the impurity from the polaron to the soliton state has been studied in [17] by combining the Frölich Hamiltonian picture with the Landau-Brazovskii theory for first-order phase transitions. All these theories work quite well in the weakcoupling regime but show some deviations with respect to recent experiments [5,6] in the strong-coupling regime. However, a nonperturbative renormalization-group approach [18] seems able to give a reliable and unified picture of the Bose polaron problem from weak to strong coupling, also in one-dimensional configurations [19].
A common feature of the above investigations is that the impurity always remains inside the bosonic bath. In this paper we study a different but closely related nonequilibrium problem: the collision between a 41 K impurity and a 87 Rb BEC, where initially the impurity is outside the bosonic bath. In Section II we introduce the physical system, that is a quasi one-dimensional BEC of 300 87 Rb atoms located at the minimum of a weak harmonic axial trap and an impurity of 5 atoms that starts from the edge of the BEC cloud. In Section III we define the two coupled one-dimensional time-dependent Gross-Pitaevskii equations which are used to perform the numerical simulations. In Section IV we discuss our theoretical predictions which display a very rich phenomenology crucially depending on the impurity-BEC strength. Indeed our numerical simulations unequivocally show that the collision dynamics gives rise to a variety of highly-nonlinear effects such as dark and bright solitons, impurity trapping, and coupled oscillations in the confining trap, which, due to their macroscopic character, pave the way to forthcoming experiments. The paper is concluded by Section V.

Properties of the system
We consider a setup very close to the one realized experimentally in Ref. [4]: a bosonic cloud of 87 Rb atoms and a bosonic impurity made of 41 K atoms in a 1D harmonic confinement. This quasi-1D configuration is obtained by a weak optical confinement along one direction and a strong optical confinement along the other two transverse directions. The atoms interact by intra-species and inter-species interactions. In particular, in this paper the inter-species interaction is controlled by the magnetic Feshbach resonance. On the contrary, the intra-species interactions are always repulsive and close to their background values in the range of magnetic field used to exploit the Feshbach resonance. In a realistic experiment, the variation of the inter-species s-wave scattering length a Rb,K between 87 Rb and 41 K atoms by a Feshbach resonance can induce also variations of intra-species scattering lengths a Rb and a K . However, these variations are quantitatively negligible [4]. In our simulations we use the following values for the 3D scattering lengths: a Rb = 100 a 0 [20] and a K = 63 a 0 [21] with a 0 the Bohr radius. We work with N Rb = 300 atoms for Bose-Einstein condensate and N K = 5 atoms for the impurity. Unlike the experiment of Ref. [4], in this work we assume zero temperature throughout. For 87 Rb atoms the frequencies of transverse and axial harmonic confinement are ω ⊥Rb = 2π × 34 × 10 3 Hz and ω Rb = 2π × 62 Hz. For 41 K atoms they are instead ω ⊥K = 2π × 45 × 10 3 Hz and ω K = 2π × 87 Hz.
Two lengths characteristic of the considered problem can be defined naturally for each species: the longitudinal and transverse harmonic oscillator lengths, respectively a ||,s = /(m s ω || s ) and a ⊥,s = /(m s ω ⊥,s ) with s =Rb, K. Due to the relations ω ⊥Rb ω Rb and ω ⊥K ω K between the confinement frequencies, the system behaves effectively as one-dimensional. The one-dimensional scattering lengths characterizing the intra-species interactions are obtained from the three-dimensional ones by the Olshanii formula [22]: a s,1D = −(a 2 ⊥ /a s )(1 − C(a s /a ⊥,s )) with C 1.4603/ √ 2. In this way, the 1D intra-species interaction strengths are given by with g 1 = 2.365 J·m and g 2 = 0.8598 g 1 ; while the 1D inter-species interaction strength between 87 Rb and 41 K atoms, g 12 , is numerically calculated following the analysis laid out in [23]. Following Ref. [24], we determine the 3D inter-atomic scattering length a Rb,K as a function of the magnetic field B near the Feshbach resonance at 78.2 Gauss. We then derive the corresponding 1D inter-atomic strength g 12 .

Theoretical approach
We describe the system under investigation in the frame of mean-field theory for both components. In the configuration of Ref. [4], the transverse confinement corresponds to a harmonic-oscillator length of approximately 1150a 0 for the majority component (Rb). As a consequence, the 1D gas parameter, i.e. the ratio of the healing length to interparticle distance, is ξ/d = a ⊥,Rb / √ 8a Rb d 1, justifying the use of mean-field theory for the 87 Rb cloud. In fact, 87 Rb atoms are in the 1D weak-coupling quasi-condensate regime, where the exact Lieb-Liniger theory of 1D bosons with repulsive contact interaction reduces to the 1D Gross-Pitaevskii model [25,27]. We complement the complex wavefunctions ψ Rb (x, t) with a similar wavefunction ψ K (x, t) for the impurity and impose that they satisfy the coupled Gross-Pitaevskii equations (GPEs) the atom numbers, which remain constant during the time evolution. In our model also the 41 K impurity is described by a 1D GPE, but the very small number of 41 K atoms implies that the nonlinear term proportional to g 2 is extremely small. We have verified that setting g 2 = 0 the numerical results of Section 4 are practically the same. The key parameter is instead g 12 , that is the strength of the density-density coupling between 87 Rb cloud and 41 K impurity. An accurate estimate of the axial density profile of the 87 Rb cloud (at the initial time and with good approximation also during the dynamics) can be achieved by the local density approximation, leading to where are the Thomas-Fermi radius of the Rb cloud and the corresponding chemical potential respectively. Inserting the numerical values for the various parameters one finds µ Rb = 75 ω Rb and R Rb = 12.2a Rb , where a Rb = /m Rb ω Rb = 1.37µm.

Numerical results
For the 41 K impurity at the initial time of considered evolution, we assume for most of the considered cases a Gaussian wavefunction centered at a distance d 0 = 14.6a Rb = 20µm from the origin of the axial harmonic trap. The width σ is chosen to be equal to a Rb . Similarly, in the following (simulations included) all the lengths will be measured in unit of a Rb . Eqs. (3) and (4) are solved numerically by using a finite-difference predictorcorrector Crank-Nicolson algorithm [28]. The center of mass of the 87 Rb BEC remains practically constant while the center of mass of 41 K impurity, which is initially outside the 87 Rb BEC, evolves in time due to the axial harmonic potential and its dynamics strongly depends on interaction strength g 12 between the impurity and the 87 Rb condensate. As expected if g 12 = 0 the 41 K impurity crosses the 87 Rb cloud without perturbing it (see movie [30] for details). In this case the 41 K impurity simply oscillates back and forth with oscillation frequency Ω, that is exactly the frequency ω K of the axial harmonic confinement for 41 K atoms.
In Fig. 1 we plot the reflection coefficient R of the 41 K impurity for different values of the inter-atomic strength g 12 /g 1 at time t = π/ω ||K . The reflection coefficient is computed as where N K is the number of 41 K atoms. The coefficient R equals the unity when the 41 K impurity is entirely confined in the left side of the space domain (x < 0), while R vanishes when the impurity is completely in the right side of the domain (x > 0). As expected, at the time t = π/ω ||K , which corresponds to half period of oscillation, in the absence of inter-atomic interaction (g 12 = 0) the reflection coefficient is zero. For strongly repulsive (g 12 /g 1 1) interactions close to the resonance (i.e. for B 78.2 Gauss), Fig. 1 shows instead that the impurity is completely confined out the cloud (reflected back) leading to reflection coefficient R 1. For weakly attractive g 12 /g 1 , the coefficient R always lies at very small values close to zero. However, for strongly attracting values of g 12 /g 1 it is possible to notice a small increase in R. This occurs because, due to the strongly-attractive interaction, the 41 K wave-packet tends to be captured toward the center of mass of the 87 Rb cloud and, at the time t = π/ω ||K , a fraction of it is still in the left side (x < 0) of the space domain. For the same reason, the oscillation frequency of the strongly-attractive case reaches higher values with respect to ω K , as shown in the first panel of Fig. 2.

Weak-coupling and periodic motion of impurity
For small values of the inter-atomic strength, i.e. for |g 12 |/g 1 1 the 41 K the impurity simply oscillates back and forth inside the 87 Rb cloud: However, the frequency Ω of oscillation depends on g 12 .
To determine Ω we calculate the center-of-mass position x cm (t) of the 41 K impurity as a function of time t. Then we perform the Fourier transform of x cm (t) and plot its power spectrum |x cm (ω)| 2 vs ω/ω ||K in Fig. 2. The panels of this figure are obtained for different values of the inter-atomic strength g 12 between 41 K and 87 Rb atoms. The figure clearly shows that, as expected, for g 12 = 0 (middle panel) there is only one peak centered at ω/ω ||K = 1 and, consequently, the center of mass oscillates at the frequency Ω = ω ||K , that is the frequency of axial harmonic confinement of the k1 K atoms. Upper panels of Fig. 2 reveal that for small and negative (attractive) values of g 12 /g 1 the frequency Ω of oscillation of the 41 K impurity increases. Instead, for small and positive (repulsive) values of g 12 /g 1 (lower panels) a second mode appears at a lower frequency. This second mode becomes dominant growing g 12 /g 1 .

Intermediate coupling and impurity trapping
For intermediate couplings, i.e. around |g 12 |/g 1 0.5 no matter the sign of g 12 , after a few cycles of oscillation the 41 K impurity is captured by the 87 Rb cloud and it strongly changes its amplitude of oscillation. This phenomenon is shown in Fig. 3, where the contour plot of 41 K and 87 Rb density profiles is reported.
A difference between the attractive and the repulsive case is visible in the shape of the density profile of the impurity at fixed time. Indeed the repulsive case (lower panel) shows a mode of oscillation of the impurity featuring a density profile characterized by two main lobes for the distribution, with always an empty-region at the center of the oscillating wave packet. In other words, the impurity profile features, along the x-direction, two main maxima and one central minimum. On the other hand, the  attractive case (upper panel) displays an oscillation mode in which the wave packet features only a single lobe: a density profile featured by only one maximum. Moreover, in spite of the relevant magnitude of the attraction, no significant quantum reflection phenomena occurring for rapidly varying attractive potentials [31,32] are observed at the boundaries of the BEC condensate, neither the consequent vortex formation [33]. The latter absence, also holding for abrupt repulsive potential [34], is ascribable to the one-dimensionality of our simulations, as well as to the small size of the impurity. The full dynamics of the attractive and repulsive cases is well illustrated in Fig. 3  Fig. 1 this effect corresponds to a reflection coefficient R that becomes different from zero. Notably, when the inter-atomic interaction is strongly repulsive the impurity of 41 K is not able to enter the 87 Rb cloud, so that a barrier effect occurs (see movie [37] for details). The impurity then behaves as a classical object, similar to what observed with BEC solitons in presence of potential barriers much wider that their size [38]. In Fig. 1 this effect corresponds to a reflection coefficient R close to one. However, if the initial displacement of the impurity is increased so that to increase its initial potential energy, the 41 K cloud is able to penetrate the 87 Rb cloud. Due to the strong repulsive interaction a local minimum in the density of 87 Rb is observed in the correspondence of a sharp density peak of the 41 K impurity. The described situation is shown in the left panel of Fig. 4, where the normalized density profiles |Ψ Rb (x, τ )| 2 /N Rb and |Ψ K (x, τ )| 2 /N K are reported for two different values of g 12 /g 1 . The resulting moving hole in the 87 Rb BEC is a dark soliton created by the interaction with the 41 K impurity. In this regime, the hole-impurity pair is similar to the dark-bright solitons observed in the superfluid counterflow of miscible condensates [39]. Notice that in Fig. 4 the upper panels are contour plots of the normalized density profiles in the (x, y) plane. These contour plots are obtained adopting a Gaussian profile along the y axis with a width given by the characteristic length of transverse harmonic confinement. However, for the sake of visibility, the y direction is plotted not in scale. Movie [40] displays the full dynamics of this dark soliton.
full dynamics of this bright soliton is reported in the movie, see [41].

Capture and Localization effects
Movies [30,35,36,37,40,41] show a broad range of phenomenology. Part of the rich behavior highlighted in the movies can be efficiently summarized in the dynamics of the center of mass of the 41 K cloud. In Fig. 5 we plot the center of mass of the 41 K cloud as a function of time for several values of attractive (upper panel) and repulsive (lower panel) interaction g 12 /g 1 . Fig. 5 clearly shows how, for sufficiently strong repulsive values of g 12 /g 1 , the impurity cloud is completely blocked out (light-blue and orange line in the lower panel of Fig. 5), while in the opposite case, when the interaction is strongly attractive, the oscillations of 41 K impurity is completely captured and confined at the center of the trap. As pointed out in the previous section, for intermediate interaction strengths the system exhibit a more symmetric behavior. In both the attractive and the repulsive cases the impurity is always captured by the 87 Rb BEC. This can be clearly seen in Fig. 5 where for g 12 /g 1 ∈ [−0.5, −2.0] (upper panel) and g 12 /g 1 = 0.5 (lower panel) the evolution of the center of mass is damped in amplitude and confined inside the 87 Rb BEC (i.e. center of the trap). Interestingly, the "capture mechanism" seems to exhibit a set of common features among the different considered cases: (i) it occurs very quickly (about one oscillation cycle); (ii) smaller values of the interaction strength |g 12 | lead to a delay of the impurity capture; (iii) at least for intermediate interaction strengths, the impurity capture is always preceded by the appearance of an interference pattern in the density distribution of the impurity cloud. This sudden damping of the oscillation has been verified to be a rearrangement of internal energies, where part of the large initial kinetic energy of the impurity is handed over to the BEC cloud through inter-species interactions, hence producing a damping on the impurity oscillation.
In Fig. 6 we show two frames of the movies [35] and [36], for attractive g 12 /g 1 = −0.5 (left panel) and repulsive g 12 /g 1 = 0.5 (right panel) cases in which an interference pattern appears. The interaction of the wave packet Ψ K (x, t) with the Rb condensate produces a reflected counter-propagating wave that, by quantum interfering with the incoming packet, produces the interference pattern shown in Fig 6. This exclusively quantum phenomenon, previously encountered and studied in detail for instance in [33,39], proves to be driven by the interaction between the two atomic species g 12 , and arises in both the attractive and repulsive cases leaving a direct signature of the energy transfer between the impurity and the condensate. A difference in the spatial frequency of the interference pattern fringes is noticed between the attractive and repulsive cases with the former manifesting a generally higher spatial frequencies than the latter. However, the investigation of the true nature of such a difference goes beyond the goal of this paper and is left to future development.
A similar reflection-interference phenomenon appears at the boundary of the trap due to the reflection of the cloud between the parabolic walls. In that case, the intra- species interaction term g 2 provides that source of scattering between the 41 K particles capable of producing the reflected wave and the associated interference figure. We verified that if both g 2 = 0 and g 12 = 0 no interference patterns appear in the time evolution. In this case, the dynamical evolution of |Ψ K (x, t)| 2 is perfectly described by the typical coherent-state picture. The rich phenomenology described above reproduces the well-known localization effects characterizing mixtures both in the absence (see, e.g., [42]) and in the presence [43] of a superimposed optical lattice. In the repulsive case, two species, fully mixed for g 12 / √ g 1 g 2 < σ ≈ 1 (σ is determined in [42] and [43]), separates when g 12 is sufficiently larger than √ g 1 g 2 thus providing a spatial configuration where the density maximum of one species corresponds to density minimum of the other (the dark soliton in Fig. 4, lower left panel). A similar effect occurs if g 12 < 0: when |g 12 | is sufficiently larger than √ g 1 g 2 , a configuration crops up in which the density maxima of both species perfectly overlap (local supermixing) due to the attractive interaction (see Fig. 4, lower right panel). Such configurations clearly emerge in the oscillations of the impurity in the Rb cloud. In particular, after the capture of the K impurity by the Rb cloud, one can observe how the final part of the oscillations shown in [33], [34] feature a stable bond of the K density maximum with the Rb dark soliton (bright soliton) in the presence of a strong repulsive (attractive) interaction.

Conclusions
In this paper we have analyzed the behavior of a system characterized by quasi onedimensional Bose-Einstein condensate made of three hundreds 87 Rb atoms interacting with a bosonic impurity made of five 41 K atoms, which starts from outside the Bose condensate and collides on it. Despite the specific physical system under investigation, the results we obtain are quite general and, depending on the range and on the sign of the inter-species interaction, different regimes can be identified. Among them, the full reflection of the impurity, the impurity trapping, and also the emergence of dark and bright solitons in the Bose condensate. Preliminary experiments have measured the reflection coefficient as a function of interatomic strength in the system of Ref.
[4] at finite temperature. Thus, our zero-temperature theoretical predictions provide a useful benchmark for forthcoming experimental and theoretical investigations on impurity-BEC collisions. We think that our mean-field simulations, based on coupled Gross-Pitaevskii equations, are quite reliable but they could be improved adopting more sophisticated time-dependent approaches where quantum depletion is taken into account. In particular, when the 87 Rb -41 K scattering length is very large, the meanfield density-density interaction is questionable and beyond-mean-field effects could be relevant. Finally, it is important to stress that a deeper connection between impurity in bosonic atomic clouds and the solid-state polaron of Landau and Frölich can be obtained with a bosonic lattice polaron [9], i.e. a single-impurity atom confined to a optical lattice and immersed in a homogeneous Bose-Einstein condensate. We are now planning to investigate this difficult but stimulating problem by using both mean-field and beyond-mean-field techniques.
Equations (3) and (4) are solved numerically by using a finite-difference predictorcorrector Crank-Nicolson algorithm [28]. Numerical discretization is performed on a fixed mesh-grid with constant spatial spacing dx/a Rb = 1.50 × 10 −3 and constant temporal spacing dt/ω Rb = 1.25 × 10 −3 . A single computation run, performed on a laptop with an Intel i7 2.90 GHz processor and 16 GB RAM, can take up to 40 min with the above discretization.
To test the numerical accuracy of our results we checked the validity of the conservation laws for our solutions, namely, the conservation of energy E and the conservation of total number of particles N Rb and N K . Energy conservation is tested by computing the energy functional at any instant of time t E(t) = Ψ Rb (x, t)Ĥ Rb Ψ Rb (x, t) dx + Ψ K (x, t)Ĥ K Ψ K (x, t) dx (A.1) whereĤ Rb andĤ K are the Hamiltonian operators on the right-hand side of equations (3) and (4) respectively. Similarly, the test of the conservation of particle numbers is performed by checking the correct normalization of the wave-functions Ψ Rb (x, t) and Ψ K (x, t), verifying that equations (5) and (6) hold for any time t. We verified that the normalization of the wave-functions is always almost perfectly conserved, and observed that, with the given discretization also the energy is conserved (> 98%). Only for strong enough attractive interactions, the strong non-linearity of the soliton-solution leads to a decrease of energy conservation. However, we verified that, by increasing time and space discretization, this problem is easily overcome at the expense of simulation-time, with no qualitative change in our results.