Collision of impurities with Bose–Einstein condensates

In 1933, Landau introduced the concept of the polaron, an electron whose effective mass is affected by 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 [4], localized bosonic impurities, consisting of a few ⁴¹K atoms, have been created in a one-dimensional (1D) configuration and their interactions with a Bose–Einstein condenstate (BEC) of ⁸⁷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 impurities in a cloud of ultracold Bose–Einstein condensed atoms—have been obtained by two other experimental groups [5, 6]. The Bose polaron problem has been addressed theoretically using various 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 weak-coupling 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, including 1D 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 non-equilibrium problem: the collision between a ⁴¹K impurity and a ⁸⁷Rb BEC, where the impurity is initially outside the bosonic bath. In section 2, we introduce the physical system: a quasi-1D BEC of 300 ⁸⁷Rb atoms located at the minimum of a weak harmonic axial trap and an impurity of five atoms that starts from the edge of the BEC cloud. In section 3, we define the two coupled 1D time-dependent Gross−Pitaevskii equations which are used to perform the numerical simulations. In section 4, 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 in section 5.

We consider a set-up very close to that realized experimentally in [4]: a bosonic cloud of ⁸⁷Rb atoms and a bosonic impurity consisting of ⁴¹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 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. In contrast, 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 ⁸⁷Rb and ⁴¹K atoms by a Feshbach resonance can also induce 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₀ [20] and a_K = 63 a₀ [21] with a₀ 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 [4], in this work we assume zero temperature throughout. For ⁸⁷Rb atoms, the frequencies of transverse and axial harmonic confinement are ω_⊥Rb = 2π × 34 × 10³ Hz and ω_∥Rb = 2π × 62 Hz. For ⁴¹K atoms, they are instead ω_⊥K = 2π × 45 × 10³ Hz and ω_∥K = 2π × 87 Hz.

Two lengths characteristic of the problem under consideration can be defined naturally for each species: the longitudinal and transverse harmonic oscillator lengths, respectively ${a}_{| | ,s}=\sqrt{{\hslash }/({m}_{s}\,{\omega }_{| | s})}$ and ${a}_{\perp ,s}=\sqrt{{\hslash }/({m}_{s}\,{\omega }_{\perp ,s})}$ with s = Rb, K. Due to the relations ω_⊥Rb ≫ ω_∥Rb and ${\omega }_{\perp {\rm{K}}}\gg {\omega }_{\parallel {\rm{K}}}$ between the confinement frequencies, the system effectively behaves 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}_{\perp }^{2}/{a}_{s})(1-C({a}_{s}/{a}_{\perp ,s}))$ with $C\simeq 1.4603/\sqrt{2}$. In this way, the 1D intra-species interaction strengths are given by

$$\begin{eqnarray}&&{g}_{1}=-\displaystyle \frac{2\,{{\hslash }}^{2}}{{m}_{{\rm{Rb}}}\,{a}_{{\rm{Rb}},1D}},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{g}_{2}=-\displaystyle \frac{2\,{{\hslash }}^{2}}{{m}_{{\rm{K}}}\,{a}_{{\rm{K}},1D}},\end{eqnarray} \tag{ 2 }$$

with g₁ = 2.365 J m and g₂ = 0.8598 g₁; while the 1D inter-species interaction strength between ⁸⁷Rb and ⁴¹K atoms, g₁₂, is numerically calculated following the analysis laid out in [23]. Following [24], we determine the 3D inter-atomic scattering length ${a}_{{\rm{Rb}},{\rm{K}}}$ as a function of the magnetic field B near the Feshbach resonance at 78.2 Gs. We then derive the corresponding 1D inter-atomic strength g₁₂.

We describe the system under investigation in the context of mean-field theory for both components. In the configuration of [4], the transverse confinement corresponds to a harmonic oscillator length of approximately 1150a₀ for the majority component (Rb). As a consequence, the 1D gas parameter, i.e. the ratio of the healing length to interparticle distance, is $\xi /d={a}_{\perp ,{\rm{Rb}}}/\sqrt{8{a}_{{\rm{Rb}}}d}\lesssim 1$, justifying the use of mean-field theory for the ⁸⁷Rb cloud. In fact, ⁸⁷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 wave-functions ψ_Rb(x, t) with a similar wavefunction ψ_K(x, t) for the impurity, and impose that they satisfy the coupled Gross−Pitaevskii equations (GPEs)

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial {{\rm{\Psi }}}_{{\rm{Rb}}}}{\partial t}=\left[-\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{{\rm{Rb}}}}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}+\displaystyle \frac{1}{2}{m}_{{\rm{Rb}}}\,{\omega }_{| | {\rm{Rb}}}^{2}\,{x}^{2}+{g}_{1}| {{\rm{\Psi }}}_{{\rm{Rb}}}{| }^{2}+{g}_{12}| {{\rm{\Psi }}}_{K}{| }^{2}\right]{{\rm{\Psi }}}_{{\rm{Rb}}},\end{eqnarray} \tag{ 3 }$$

$$\begin{eqnarray}&&{\rm{i}}{\hslash }\displaystyle \frac{\partial {{\rm{\Psi }}}_{{\rm{K}}}}{\partial t}=\left[-\displaystyle \frac{{{\hslash }}^{2}}{2{m}_{{\rm{K}}}}\displaystyle \frac{{\partial }^{2}}{\partial {x}^{2}}+\displaystyle \frac{1}{2}{m}_{{\rm{K}}}\,{\omega }_{| | {\rm{K}}}^{2}\,{x}^{2}+{g}_{2}| {{\rm{\Psi }}}_{{\rm{K}}}{| }^{2}+{g}_{12}| {{\rm{\Psi }}}_{{Rb}}{| }^{2}\right]{{\rm{\Psi }}}_{{\rm{K}}},\end{eqnarray} \tag{ 4 }$$

where

$$\begin{eqnarray}&&{N}_{{\rm{Rb}}}=\int | {{\rm{\Psi }}}_{{\rm{Rb}}}(x,t){| }^{2}\ {\rm{d}}{x},\end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray}&&{N}_{{\rm{K}}}=\int | {{\rm{\Psi }}}_{{\rm{K}}}(x,t){| }^{2}\ {\rm{d}}{x}\ \end{eqnarray} \tag{ 6 }$$

indicate the number of atoms of each type, and remain constant during the time evolution. In our model, the ⁴¹K impurity is also described by a 1D GPE, but the very small number of ⁴¹K atoms implies that the nonlinear term proportional to g₂ is extremely small. We have verified that the numerical results of section 4 are practically the same on setting g₂ = 0. The key parameter is instead g₁₂—that is, the strength of the density−density coupling between ⁸⁷Rb cloud and ⁴¹K impurity. An accurate estimate of the axial density profile of the ⁸⁷Rb cloud (at the initial time and, to a good approximation, also during the dynamics) can be achieved by the local density approximation, leading to

$$\begin{eqnarray}{\rho }_{{\rm{Rb}}}(x,0)=| {{\rm{\Psi }}}_{{\rm{Rb}}}(x,0){| }^{2}\simeq \left\{\begin{array}{ll}\tfrac{1}{{g}_{1}}\left({\mu }_{{\rm{Rb}}}-\tfrac{1}{2}\,{m}_{{\rm{Rb}}}\,{\omega }_{\parallel {\rm{Rb}}}^{2}\,{x}^{2}\right) & {\rm{if}}\ | x| \leqslant {R}_{{\rm{Rb}}}\\ 0 & {\rm{o}}{\rm{t}}{\rm{h}}{\rm{e}}{\rm{r}}{\rm{w}}{\rm{i}}{\rm{s}}{\rm{e}}\end{array}\right.,\end{eqnarray} \tag{ 7 }$$

where

$$\begin{eqnarray}&&{R}_{{\rm{Rb}}}=\sqrt{(2\,{\mu }_{{\rm{Rb}}})/({m}_{{\rm{Rb}}}\,{\omega }_{\parallel {\rm{Rb}}}^{2})},\end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray}&&{\mu }_{{\rm{Rb}}}={((3{g}_{{\rm{Rb}}}{N}_{{\rm{Rb}}}\sqrt{{m}_{{\rm{Rb}}}{\omega }_{\parallel {\rm{Rb}}}^{2}})/(4\sqrt{2}))}^{2/3}\end{eqnarray} \tag{ 9 }$$

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 ${\mu }_{{\rm{Rb}}}=75{\hslash }{\omega }_{\parallel {\rm{Rb}}}$ and R_Rb = 12.2a_∥Rb, where ${a}_{\parallel {\rm{Rb}}}=\sqrt{{\hslash }/{m}_{{\rm{Rb}}}{\omega }_{\parallel {\rm{Rb}}}}$ = 1.37 μm.

For the ⁴¹K impurity at the initial time of considered evolution, we assume for most of the cases considered a Gaussian wavefunction

$$\begin{eqnarray}&&{{\rm{\Psi }}}_{{\rm{K}}}(x,0)=\displaystyle \frac{{N}_{k}}{{\pi }^{1/4}{\sigma }^{1/2}}\ {{\rm{e}}}^{-{(x+{d}_{0})}^{2}/(2{\sigma }^{2})}\end{eqnarray} \tag{ 10 }$$

centered at a distance d₀ = 14.6a_∥Rb = 20 μm from the origin of the axial harmonic trap. The width σ is chosen to be equal to ${a}_{\parallel {\rm{Rb}}}$. Similarly, in the following (simulations included) all the lengths will be measured in units of ${a}_{\parallel {\rm{Rb}}}$.

Equations (3) and (4) are solved numerically by using a finite-difference predictor-corrector Crank−Nicolson algorithm [28, 29]. The center of mass of the ⁸⁷Rb BEC remains practically constant, while the center of mass of the ⁴¹K impurity, which is initially outside the ⁸⁷Rb BEC, evolves in time due to the axial harmonic potential, and its dynamics strongly depends on interaction strength g₁₂ between the impurity and the ⁸⁷Rb condensate. As expected, if g₁₂ = 0 the ⁴¹K impurity crosses the ⁸⁷Rb cloud without perturbing it (see movie [30] for details). In this case, the ⁴¹K impurity simply oscillates back and forth with oscillation frequency Ω, which is exactly the frequency ω_∥K of the axial harmonic confinement for ⁴¹K atoms.

In figure 1, we plot the reflection coefficient R of the ⁴¹K impurity against the inter-atomic strength g₁₂/g₁ at time $t=\pi /{\omega }_{| | {\rm{K}}}$. The reflection coefficient is computed as

$$\begin{eqnarray}&&R(t)=\displaystyle \frac{1}{{N}_{{\rm{K}}}}{\int }_{-\infty }^{0}| {{\rm{\Psi }}}_{{\rm{K}}}(x,t){| }^{2}\ {\rm{d}}{x},\end{eqnarray} \tag{ 11 }$$

where N_K is the number of ⁴¹K atoms. The coefficient R equals unity when the ⁴¹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=\pi /{\omega }_{| | {\rm{K}}}$, which corresponds to half the period of oscillation, in the absence of inter-atomic interaction (g₁₂ = 0) the reflection coefficient is zero. For strongly repulsive (g₁₂/g₁ ≫ 1) interactions close to the resonance (i.e. for $B\simeq 78.2$ Gs), figure 1 shows instead that the impurity is completely confined outside the cloud (reflected back), leading to reflection coefficient $R\simeq 1$. For weakly attractive g₁₂/g₁, the coefficient R always lies at very small values close to zero. However, for strongly attractive 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 ⁴¹K wave packet tends to be pulled toward the center of mass of the ⁸⁷Rb cloud and, at time $t=\pi /{\omega }_{| | {\rm{K}}}$, a fraction of it is still in the left side (x < 0) of the space domain. For the same reason, the oscillation frequency in the strongly attractive case reaches higher values with respect to ω_∥K, as shown in the first panel of figure 2.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

4.1. Weak-coupling and periodic motion of impurity

For small values of the inter-atomic strength, i.e. for $| {g}_{12}| /{g}_{1}\lesssim 1$ the ⁴¹K the impurity simply oscillates back and forth inside the ⁸⁷Rb cloud: however, the frequency Ω of oscillation depends on g₁₂.

To determine Ω, we calculate the center-of-mass position x_cm(t) of the ⁴¹K impurity as a function of time t. Then we perform the Fourier transform

$$\begin{eqnarray}&&{\tilde{x}}_{{\rm{cm}}}(\omega )=\int {x}_{{\rm{cm}}}(t)\ {{\rm{e}}}^{{\rm{i}}\omega t}\ {\rm{d}}{t}\ \end{eqnarray} \tag{ 12 }$$

of x_cm(t), and plot its power spectrum $| {\tilde{x}}_{{\rm{cm}}}(\omega ){| }^{2}$ versus $\omega /{\omega }_{| | {\rm{K}}}$—see figure 2. The panels of this figure are each obtained for different values of the inter-atomic strength g₁₂ between ⁴¹K and ⁸⁷Rb atoms. The figure clearly shows that, as expected, for g₁₂ = 0 (middle panel) there is only one peak centered at $\omega /{\omega }_{| | {\rm{K}}}=1$ and, consequently, the center of mass oscillates at the frequency ${\rm{\Omega }}={\omega }_{| | {\rm{K}}}$, that is the frequency of axial harmonic confinement of the ⁴¹K atoms.

The upper panels of figure 2 reveal that for small negative (attractive) values of g₁₂/g₁ the frequency Ω of oscillation of the ⁴¹K impurity increases. In contrast, for small positive (repulsive) values of g₁₂/g₁ (lower panels), a second mode appears at a lower frequency. This second mode becomes dominant as g₁₂/g₁ grows.

4.2. Intermediate coupling and impurity trapping

For intermediate couplings, i.e. around $| {g}_{12}| /{g}_{1}\simeq 0.5$ regardless of the sign of g₁₂, after a few cycles of oscillation the ⁴¹K impurity is captured by the ⁸⁷Rb cloud, and strongly changes its amplitude of oscillation. This phenomenon is illustrated in figure 3, where the contour plot of ⁴¹K and ⁸⁷Rb density profiles is reported.

Zoom In Zoom Out Reset image size

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, always with 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 characterized by only one maximum. Moreover, in spite of the relevant magnitude of the attraction, no significant quantum reflection phenomenon occurring for rapidly varying attractive potentials [31, 32] is observed at the boundaries of the BEC condensate, nor is there any consequent vortex formation [33]. The latter absence, also holding for abrupt repulsive potential [34], is attributable 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 figure 3, and can also be seen in the movies [35] and [36], respectively.

4.3. Strong coupling and solitary waves

4.3.1. Repulsive inter-atomic strength

When the inter-atomic interaction is repulsive and sufficiently strong, only a minor part of the ⁴¹K impurity ends up on the other side of the ⁸⁷Rb cloud. In figure 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 ⁴¹K is not able to enter the ⁸⁷Rb cloud, so that a barrier effect occurs (see movie [37] for details). The impurity then behaves as a classical object, similar to what is observed with BEC solitons in the presence of potential barriers much wider that their size [38]. In figure 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 ⁴¹K cloud is able to penetrate the ⁸⁷Rb cloud. Due to the strong repulsive interaction, a local minimum in the density of ⁸⁷Rb is observed in the correspondence of a sharp density peak of the ⁴¹K impurity. The situation described is shown in the left panel of figure 4, where the normalized density profiles $| {{\rm{\Psi }}}_{{\rm{Rb}}}(x,\tau ){| }^{2}/{N}_{{\rm{Rb}}}$ and $| {{\rm{\Psi }}}_{{\rm{K}}}(x,\tau ){| }^{2}/{N}_{{\rm{K}}}$ are reported for two different values of g₁₂/g₁. The resulting moving hole in the ⁸⁷Rb BEC is a dark soliton created by the interaction with the ⁴¹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 figure 4 the upper panels are contour plots of the normalized density profiles in the (x,y) plane. These contour plots are obtained by 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 not plotted to scale. Movie [40] displays the full dynamics of this dark soliton.

Zoom In Zoom Out Reset image size

4.3.2. Attractive inter-atomic strength

When the sign of g₁₂ is taken to be negative, the ⁴¹K impurity enters and oscillates in the ⁸⁷Rb BEC; at each oscillation a part of it is captured around the center, so that after a certain time all the impurity gets confined around this point. Increasing the magnitude of the negative g₁₂, the number of oscillations prior to complete trapping decreases. Moreover, when full confinement is reached, a local maximum in the ⁸⁷Rb density occurs, corresponding to the peak in the ⁴¹K impurity. This moving peak is a bright soliton, shown in the right panel of figure 4. The full dynamics of this bright soliton is reported in movie [41].

4.4. Capture and Localization effects

Movies [30, 35–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 ⁴¹K cloud. In figure 5, we plot the center of mass of the ⁴¹K cloud as a function of time for several values of attractive (upper panel) and repulsive (lower panel) interaction g₁₂/g₁. Figure 5 clearly shows how, for sufficiently strong repulsive values of g₁₂/g₁, the impurity cloud is completely blocked out (light blue and orange lines in the lower panel of figure 5), while in the opposite case, when the interaction is strongly attractive, the ⁴¹K impurity is completely captured and confined at the center of the trap.

Zoom In Zoom Out Reset image size

As pointed out in the previous section, for intermediate interaction strengths the system exhibits a more symmetric behavior. In both attractive and repulsive cases, the impurity is always captured by the ⁸⁷Rb BEC. This can be clearly seen in figure 5 where for ${g}_{12}/{g}_{1}\in [-0.5,-2.0]$ (upper panel) and g₁₂/g₁ = 0.5 (lower panel) the evolution of the center of mass is damped in amplitude and confined inside the ⁸⁷Rb BEC (i.e. at the center of the trap). Interestingly, the 'capture mechanism' seems to exhibit a set of common features among the various cases considered: (i) it occurs very quickly (about one oscillation cycle); (ii) smaller values of the interaction strength $| {g}_{12}|$ lead to delayed 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 of the impurity oscillation.

In figure 6 we show two frames of the movies [35] and [36], for attractive ${g}_{12}/{g}_{1}=-0.5$ (left panel) and repulsive g₁₂/g₁ = 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 figure 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₁₂, 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 scope of this paper, and is left to future development.

Zoom In Zoom Out Reset image size

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₂ provides that source of scattering between the ⁴¹K particles capable of producing the reflected wave and the associated interference figure. We verified that if both ${g}_{2}=0$ and g₁₂ = 0, no interference patterns appear in the time evolution. In this case, the dynamical evolution of $| {{\rm{\Psi }}}_{{\rm{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}/\sqrt{{g}_{1}{g}_{2}}\lt \sigma \approx 1$ (σ is determined in [42] and [43]), separates when g₁₂ is sufficiently larger than $\sqrt{{g}_{1}{g}_{2}}$, thus providing a spatial configuration where the density maximum of one species corresponds to the density minimum of the other (the dark soliton in figure 4, lower left panel). A similar effect occurs if g₁₂ < 0: when $| {g}_{12}|$ is sufficiently larger than $\sqrt{{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 figure 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] features a stable bond of the K density maximum with the Rb dark soliton (bright soliton) in the presence of a strong repulsive (attractive) interaction.

In this paper we have analyzed the behavior of a system consisting of a quasi-1D Bose−Einstein condensate made of 300 ⁸⁷Rb atoms interacting with a bosonic impurity made of five ⁴¹K atoms, which starts from outside the Bose condensate and collides with it. Despite the specific physical system under investigation, the results we obtain are quite general, and, different regimes can be identified depending on the range and on the sign of the inter-species interaction. These include the full reflection of the impurity, trappping of the impurity, and also the emergence of dark and bright solitons in the Bose condensate. Preliminary experiments have measured the reflection coefficient as a function of inter-atomic strength in the system of [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 by adopting more sophisticated time-dependent approaches where quantum depletion is taken into account. In particular, when the ⁸⁷Rb—⁴¹K scattering length is very large, the mean-field density−density interaction is questionable, and beyond-mean-field effects could be relevant. Finally, it is important to stress that a deeper connection between impurities 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 an optical lattice and immersed in a homogeneous Bose–Einstein condensate. We are now planning to investigate this difficult but stimulating problem using both mean-field and beyond-mean-field techniques.

We thank Giacomo Lamporesi for useful discussions. FM acknowledges funding from FP7 Cooperation STREP Project EQuaM (Grant n. 323714). LS acknowledges Project BIRD164754 of University of Padova for partial support.

Equations (3) and (4) are solved numerically by using a finite-difference predictor-corrector Crank−Nicolson algorithm [28]. Numerical discretization is performed on a fixed mesh-grid with constant spatial spacing dx/a_∥Rb = 1.50 × 10⁻³ and constant temporal spacing dt/ω_∥Rb = 1.25 × 10⁻³. 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 was tested by computing the energy functional at any instant of time t:

$$\begin{eqnarray}&&E(t)=\int {{\rm{\Psi }}}_{{\rm{Rb}}}^{\star }(x,t){\hat{H}}_{{\rm{Rb}}}{{\rm{\Psi }}}_{{\rm{Rb}}}(x,t){\rm{d}}{x}+\int {{\rm{\Psi }}}_{{\rm{K}}}^{\star }(x,t){\hat{H}}_{{\rm{K}}}{{\rm{\Psi }}}_{{\rm{K}}}(x,t){\rm{d}}{x},\end{eqnarray} \tag{ A.1 }$$

where ${\hat{H}}_{{\rm{Rb}}}$ and ${\hat{H}}_{{\rm{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 was 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 the energy is also conserved (>98%). Only for strong enough attractive interactions does the strong non-linearity of the soliton-solution lead 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.