Many-body collisional dynamics of impurities injected into a double-well trapped Bose-Einstein condensate

We unravel the many-body dynamics of a harmonically trapped impurity colliding with a bosonic medium confined in a double-well upon quenching the initially displaced harmonic trap to the center of the double-well. We reveal that the emerging correlation dynamics crucially depends on the impurity-medium interaction strength allowing for a classification into different dynamical response regimes. For strong attractive impurity-medium couplings the impurity is bound to the bosonic bath, while for intermediate attractions it undergoes an effective tunneling. In the case of weak attractive or repulsive couplings the impurity penetrates the bosonic bath and performs a dissipative oscillatory motion. Further increasing the impurity-bath repulsion results in the pinning of the impurity between the density peaks of the bosonic medium, a phenomenon that is associated with a strong impurity-medium entanglement. For strong repulsions the impurity is totally reflected by the bosonic medium. To unravel the underlying microscopic excitation processes accompanying the dynamics we employ an effective potential picture. We extend our results to the case of two bosonic impurities and demonstrate the existence of a qualitatively similar impurity dynamics.


I. INTRODUCTION
Due to their extraordinary controllability ultracold atoms have been used to study various properties of many-body quantum systems. Indeed, they can be confined in arbitrary trapping geometries and dimensions [1][2][3][4], the underlying interatomic interactions are tunable via Feshbach resonances [5][6][7][8][9] while mixtures of quantum gases namely Bose-Bose [10,11], Bose-Fermi [12,13] and Fermi-Fermi ones [14][15][16][17] can be realized. Recently, major attention has been placed on strongly particle imbalanced mixtures where for instance a single impurity is immersed in a manybody environment. Here, the concept of a polaron [18], which has been exhaustively studied in solid state physics, can be recovered where the impurity plays the role of an effective particle dressed by the excitations of its surrounding. In this context, the existence and characteristics of Fermi [19][20][21][22][23] and Bose polarons [24][25][26][27][28] have been unveiled, mainly focusing on their stationary properties [29][30][31][32] and more recently on the dynamics [33][34][35] of these quasi-particles. Another important aspect of such impurity settings concerns their transport properties through the environment. Here, especially the involved confining potentials have a major impact on the dynamical behavior of the impurity. For instance it has been shown that the impuritymedium interaction quench dynamics in a harmonic trap leads to oscillatory [36], dipole-like and dissipative impurity [37] motion depending on the impuritybath coupling strength or to temporal orthogonality catastrophe events for strong repulsion [28]. Also the tunneling dynamics of impurities confined in a doublewell and coupled to a lattice trapped medium have been studied in the context of an effective potential [38,39]. Additionally, dephasing and clustering processes [40,41] as well as the tunneling dynamics [42] were observed for impurities confined in lattice poten-tials. Furthermore, the collisional mean-field dynamics of impurities with a harmonically confined Bose-Einstein condensate (BEC) has been investigated [43] unveiling the complete reflection of the impurities at the BEC, their trapping within the bath as well as generation of dark and bright solitons.
However, the collisional dynamics of impurities with a BEC background is possible to exhibit scattering processes caused exclusively by the presence of impurity-medium correlations. In this sense besides the mean-field collision channels most importantly new ones can be triggered including, for instance, dephasing dynamics of the impurities associated with enhanced energy redistribution processes [37,44] or the formation of impurity-medium bound states. Specifically we will investigate here the dynamical response of impurities colliding with a bosonic environment confined in a double-well. Of immediate interest will be the collisional properties of the impurity and the BEC background for different impurity-bath coupling strengths including the induced scattering and transmission processes of the impurity. Since the doublewell, being the minimal lattice model, allows for a rich dynamical behavior we shall address the questions whether a pinning of the impurity within the bath or an impurity-medium bound state formation for attractive interactions are feasible and whether the relevant scattering events are accompanied by phase separation at repulsive couplings. Moreover, the backaction of the impurity on the bosonic background where the former is expected to trigger a tunneling of the medium [45] is certainly of interest. Due to the collision of the impurities with their environment strong impurity-medium correlations are expected to emerge giving rise to beyond-mean-field collisional channels. To trace the non-equilibrium quantum dynamics, we employ the Multi-Layer Multi-Configuration Time-Dependent Hartree method for atomic mixtures (ML-arXiv:2009.12147v1 [cond-mat.quant-gas] 25 Sep 2020 MCTDHX) [46][47][48] which is capable of capturing all relevant inter-and intraspecies correlations.
To address these aspects we consider a harmonically trapped impurity which is coupled via a contact interaction potential to a bosonic environment confined in a double-well. The dynamics is induced by quenching the initially displaced harmonic confinement of the impurity to the center of the double-well. By steering the impurity-medium interaction strength from attractive to strongly repulsive values we are able to identify five dynamical response regimes in the case of a single impurity [49]. These regimes range from a bound state formation between the impurity and its environment for strong attractive impurity-bath interaction strengths to its dissipative oscillatory motion [37,50] within the bosonic background at weak attractive and repulsive couplings and finally its total reflection from the medium for strong repulsive interactions. For intermediate attractive or repulsive interaction strengths the impurity effectively tunnels between the sites of the double-well potential or it is pinned between the later, respectively [38]. In all of the above-mentioned cases we reveal the build up of a significant impurity-medium entanglement [28] which is mostly pronounced in the dissipative oscillation and the pinning regimes. To unravel the microscopic processes participating in the dynamics we construct an effective potential [35,37,42,51]. This picture enables us to understand the dynamical behavior of the impurity in all response regimes and, in particular, uncover hidden excitations in the pinning regime. Extending our results to the two-impurity case we identify five qualitatively similar response regimes as compared to the single-impurity scenario. Moreover, we explicate the involvement of single-and two-particle excitation processes of the impurities within the effective potential [52].
This work is structured as follows. In section II we introduce the system under investigation and specify the used quench protocol. The employed variational method to trace the many-body dynamics is outlined in section III. Section IV provides a detailed classification and analysis of the dynamical response regimes in dependence of the impurity-medium interaction strength. We extend our results to two impurities in section V and conclude this work in section VI providing a summary and an outlook of possible future research directions. In the appendices we further elaborate on the features of the identified dynamical response regimes discussing energy redistribution processes (appendix A), the impurity-medium two-body correlation dynamics at strong attractions (appendix B), the effective mass of the impurity (appendix C) and the exposure of hidden excitations revealed for repulsive impurity-medium couplings (appendix D). Appendix E demonstrates the collisional dynamics when considering a bosonic bath in a triple-well.

II. SETUP AND QUENCH PROTOCOL
Our setup consists of two different species of bosons B and I, also referred to as the medium and impurity species, respectively. For the two species we consider N B and N I particles of mass m B and m I , respectively. We operate in the ultracold regime and thus s-wave scattering is the dominant process allowing us to model the interaction between the atoms with a contact interaction potential [53]. Therefore, we employ for the impurity-bath interaction a contact interaction potential of strength g BI , while particles of the same species interact with a contact interaction potential among each other with strengths g BB for the environment and g II for the impurity species. Each species is confined in a different one-dimensional optical potential V σ at zero temperature. This can be easily achieved experimentally [11,[54][55][56] especially for the mass-imbalanced case under consideration of, i.e. 87 Rb atoms for the bath and 133 Cs atoms for the impurity species. The resulting Hamiltonian of the system readsĤ =Ĥ B +Ĥ I +Ĥ int wherê is the Hamiltonian of species σ ∈ {B, I}. The bosonic medium and impurity species are coupled viaĤ int = g BI . The impurities are confined in a harmonic oscillator potential where ω I is the trapping frequency with x I 0 being the spatial displacement of the trap. The environment is trapped in a double- which is constructed by superimposing a harmonic oscillator potential of frequency ω B and a Gaussian of width w B and height h B [57]. We consider for the bosonic medium N B = 20 87 Rb atoms with mass m B = 1, and for the impurity species N I = 1, 2 133 Cs atoms with mass m I = 133/87 [58].
Experimentally, a one-dimensional potential can be realized by employing a strong harmonic confinement along the transverse direction in order to freeze out the relevant degrees of freedom [53,59]. Subsequently, we provide the energy of the HamiltonianĤ in terms ofẼ =hω, whereω is the frequency of the perpendicular confinement. The length and time scales are then expressed in units ofx = h/(m Bω ) and ω −1 =hẼ −1 , respectively. Regarding the frequency of the harmonic oscillator potential of the impurities we use ω I /ω = 0.2. For the harmonic contribution to the double-well potential we employ a frequency of ω B /ω = 0.15 and for the Gaussian a width of w B /x = 0.8 and a height of h B /Ẽx = 2.0 leading to a central barrier of the double-well below which six energetically lowest eigenstates of the corresponding one-body Hamiltonian are located. We prepare the system in its ground state with the harmonic trap of the impurities displaced by x I 0 . The spatial overlap between the species in the noninteracting case is of about 3.5 % and it increases (decreases) for attractive (repulsive) interactions. After a g BI -dependent ground state is found the dynamics is induced by quenching, at t = 0, the trap center of the impurity's harmonic potential to the center of the bosonic environment, i.e. setting x I 0 = 0 (see Figure 1 for t > 0). Thereby, a collision of the initially displaced impurities with the bosonic medium is triggered, a process that strongly depends on the impurity-medium coupling strength g BI as we shall demonstrate below.

III. MANY-BODY WAVE FUNCTION ANSATZ
To calculate the quantum dynamical behavior of the binary system we employ the ab initio Multilayer Multi-Configuration Time-Dependent Hartree method for atomic mixtures (ML-MCTDHX) [46][47][48]. Within this approach we express the many-body wave function |Ψ MB (t) of the binary mixture using the Schmidt decomposition [60,61] |Ψ For our purposes we expand each species in M = 6 species functions |Ψ σ i with σ ∈ {B, I}. Moreover, the species functions are weighted with the timedependent Schmidt-coefficients λ i which contain information about the entanglement between the two species. For instance, in the case of only one nonvanishing Schmidt-coefficient the species are considered to be not entangled since the system can be described by a single product state (species meanfield ansatz) [61,62]. Next, each species function is expanded in a set of time-dependent number states | n σ i (t) with time-dependent coefficients C σ i, n . Each number state | n σ i (t) determines the configurational occupa-tion of N σ particles on d σ single-particle functions (SPFs) where, at the same time, the number of occupied SPFs must add up to the total particle number N σ (indicated by n|N σ ). In this work we employed d B = d I = 6 SPFs. Eventually, the single-particle functions are represented in a time-independent discrete variable representation (DVR) [63]. The propagation in time is performed by employing the Dirac-Frenkel variational principle [64,65] leading to a set of equations of motion for the system (see for more details [48,66]). The advantage of this method is its underlying multi-layering architecture of the total wave function combined with its time-dependent basis set [Eqs. (2) and (3)]. Especially, with the latter a co-moving basis set is utilized leading to a significant reduction of required basis functions compared e.g. to an exact diagonalization approach. On the other hand, the multi-layering structure provides access to all relevant inter-and intraspecies correlations of the system in an efficient manner.

IV. DYNAMICAL RESPONSE REGIMES OF A SINGLE IMPURITY
In the following we analyze the collisional dynamics of a single impurity (N I = 1) trapped in a harmonic oscillator and interacting with a bosonic medium confined in a double-well. Initially, the system is prepared in its ground state with the impurity's harmonic trap being spatially shifted by x I 0 /x = 8 with respect to the center of the double-well. The dynamics is induced by quenching the harmonic oscillator of the impurity to x I 0 /x = 0. Varying the impurity-medium interaction strength g BI from the strongly attractive to the strongly repulsive regime we discuss the emergent dynamical response of the impurity and its backaction to the environment. The intraspecies interaction strength between the medium particles is fixed to g BB /Ẽx = 0.5.
To reveal the overall dynamical response of the system we initially inspect the one-body density ρ is the bosonic field operator of the corresponding species. The spectral decomposition of the one-body density [67,68] reads where n σ j denote the natural populations and Φ σ,j the natural orbitals of species σ. Figure 2 presents the time-evolution of ρ (1) σ (x, t) for different impurity-medium interaction strengths g BI ranging from strongly attractive to strongly repulsive values. In Figures 2(a1) and (a2) we monitor the onebody density of the impurity ρ  B (x, t) in the course of the evolution for g BI /Ẽx = −0.9 namely for strong attractive couplings. As a result of the quench the impurity starts to oscillate with a small amplitude which decays during the time-evolution. Thereby, the spatial maximum of ρ (1) I (x, t) remains in the vicinity of the left site of the double-well potential. Also ρ (1) B (x, t) exhibits a maximum at the same location [cf. Figure 2(a2)], a phenomenon that is attributed to the strong impuritybath attraction [69]. Due to this behavior, i.e. the enhanced spatial localization tendency of the impurity and the medium, and the fact that the impurity dominantly occupies a state with negative eigenenergy [see for details Figure 5(a1)] we refer to this dynamical response regime as the steady bound state regime.
For intermediate attractive impurity-medium interactions, g BI /Ẽx = −0.7, corresponding to the onebody densities shown in Figures 2(b1) and (b2) the impurity does not localize exclusively on one site of the double-well anymore as in the above discussed case. Rather, we observe a decay of ρ (1) I (x, t) on the left site and a simultaneous increase at the right site of the double-well. This response of the impurity is reminiscent of the tunneling dynamics of a single particle confined in a double-well [70]. Later on, we will show that the effective potential encountered by the impurity resembles a double-well since it accounts for the effects of the attractive impurity-medium coupling [cf. Figure 5(b1)]. In this sense we label this response region as the (effective) tunneling regime. Note that the back-action of the impurity on the bosonic medium leads to a shift of the maximum of ρ (1) B (x, t) following the impurity's tunneling behavior.
The next response regime, which we will refer to as the dissipative oscillation regime, emerges at weakly attractive or repulsive impurity-medium interaction strengths, e.g. for g BI /Ẽx = −0.2, 0.3. Here, the impurity-medium coupling is sufficiently small such that the impurity is able to completely penetrate its environment. Consequently, the impurity initiates a tunneling of ρ (1) B (x, t) from one site of the double-well to the other which decays in the course of time. The resulting impurity dynamics turns out to be a decaying oscillatory motion with an initial amplitude as large as the spatial extent of ρ (c2) and (d1), (d2). We attribute this decay process to a continuous energy transfer from the impurity to the medium [see Appendix A for details]. A more detailed analysis of this dissipative behavior estimating also the effective mass of the emergent quasiparticle can be found in Appendix C.
In the case of intermediate repulsive impurity-bath couplings, e.g. g BI /Ẽx = 1.0, we observe a spatial localization tendency of the impurity at the trap center accompanied by vanishing oscillations [cf. Figure  2(e1)]. Therefore, we refer to this response regime as the pinning regime. During the impurity's localization process, an intra-well dynamics is induced on the bosonic background. Here, the central barrier of the double-well is effectively enlarged due to the material barrier created by the impurity and leading to an oscillatory motion of the bath cloud in each site of the double-well. Further increasing g BI to strong impurity-medium interaction strengths, e.g. g BI /Ẽx = 2.0, the underlying repulsion between the impurity and the bath particles in the left site of the double-well becomes sufficiently large such that the impurity is totally reflected [see Figure 2(f1)]. Thereby, the bosonic environment experiences a population imbalance in the double-well potential [cf. Figure 2(f2)] and we shall address this behavior as the total reflection regime.
In conclusion, we have captured five different dynamical response regimes of the impurity depending on the impurity-medium interaction strength. Remarkably, by tuning the impurity-bath coupling g BI it is possible to control the location of the impurity. It is also important to note that in all of the abovedescribed response regimes, besides the total reflection one, the impurity is trapped within the bosonic medium and thus can be dressed by the excitations of the latter giving rise to a Bose polaron [28,71,72]. However, in the total reflection regime the impurity phase separates from its environment after the first collision and, therefore, the polaron, even it is formed for very short evolution times, decays. Similar manifestations of a decaying polaron formation at strong repulsive impurity-medium couplings in the course of the evolution have been already reported in the literature being referred to as temporal orthogonality catastrophe events [28,72,73]. In order to further classify the above-discussed response regimes we invoke the mean position of the impurity x I (t) , see Figure 3(a), for different impuritymedium interactions corresponding to the aforementioned dynamical regimes. Evidently, x I (t) exhibits individual characteristics in each regime allowing for their clear distinction. For instance, in the case of weak impurity-bath couplings, e.g. g BI /Ẽx = −0.2, an oscillatory behavior of x I (t) takes place as expected from ρ (1) Figure 2(c1)]. Turning to the total reflection regime, e.g. for g BI /Ẽx = 2.0, x I (t) captures the irregular behavior of the impurity on the left edge of the double-well, thus indicating its total reflection from the bosonic environment. The long-time evolution of x I (t) for g BI /Ẽx = −0.2 is illustrated in the inset of Figure 3(a). As can be seen, the decreasing amplitude of x I (t) becomes evident which is caused by the continuous energy transfer from the impurity to the bosonic medium, see also Appendix A.
To provide the complete response phase diagram of the impurity we show the behavior of x I (t) in dependence of the impurity-medium interaction strength g BI , thus capturing the dynamical crossover between the aforementioned regimes, see Figure 3(b). For convenience the five identified response regimes are labeled from I to V. Regime I corresponds to the steady bound state formation, see the small amplitude oscillations of ρ (1) I (x, t) in the vicinity of the left site of the double-well. In regime II we find the expected behavior of x I (t) represented by its low frequency oscillations around the trap center as shown in Figure 3(a). For weak impurity-medium interaction strengths, corresponding to regime III, the dissipative oscillatory motion of x I (t) occurs characterized by a relatively large amplitude of the underlying oscillations. Increasing g BI to intermediate repulsive values we reach the pinning regime (cf. regime III) where the mean position saturates towards x I = 0. For larger g BI the impurity is not able to penetrate the bath anymore and it is totally reflected at the edge of the latter. This regime corresponds to the total reflection one and it is denoted by V in Figure 3(b). Finally, we comment on the response regimes in which the mean position obtains values close to zero, viz. the pinning and the tunneling regime. Even though the impurity's mean position within these regimes is well distinguishable [see Figure 3(a)] in a corresponding experiment a clear distinction might be challenging. To ensure a clear distinction between these regimes one can use the experimentally accessible position variance s = x 2 − x 2 [74,75]. Since in the tunneling regime the impurity is distributed over the double-well the respective variance is larger than the one in the pinning regime where the impurity is localized at the trap center (not shown here).
In order to expose the robustness of the impurity dynamics with respect to parametric variations we present in Figure 4 x I (t) for a wide range of system parameters. As we emphasized previously, it is possible to distinguish between the dynamical response regimes by inspecting the behavior of x I (t) . Therefore, we choose the impurity-bath coupling strength g BI such that we obtain a behavior of x I (t) which can be in turn associated with a specific dynamical response regime. Figures 4(a) and (b) show x I (t) obtained for x I 0 /x = 5 and x I 0 /x = 10, respectively. Here, each mean position exhibits the same behavior as the corresponding one depicted in Figure 3(a) for x I 0 /x = 8. As expected within the dissipative oscillation regime, an amplification of the oscillation amplitude occurs as the initial displacement increases.
Varying the intraspecies coupling strength g BB between the bath particles we are again able to realize the respective dynamical response regimes but for shifted g BI [see . We attribute this effect to the broadened spatial extension of the bosonic medium for larger g BB , whereas smaller values of g BB lead to a less wider cloud of the bath in the double-well. Therefore, it is easier for the impurity to overcome the bosonic medium at the left site of the double-well in the case of a larger g BB , i.e. for a broadened background, than in the case of a smaller g BB , e.g. compare the mean positions of g BI /Ẽx = 2.0 in Figures 4(c) and (d). Hence, we conclude that the impurity dynamics is robust with respect to the variation of the initial displacement x I 0 and the intraspecies interaction strength g BB of the bath.
Subsequently, we aim to quantify the associated impurity-medium entanglement by monitoring the von-Neumann entropy [76], which reads This expression possesses an upper bound for maximal entanglement between the species, viz. λ i = 1/M leading to S V N max = ln M = 1.79 in our case (M = 6), and vanishes when no entanglement is present, e.g. λ i = 1 with λ i>1 = 0. In Figure 3(c) we provide the long-time evolution of the von-Neumann entropy for different values of g BI corresponding to the five dynamical response regimes of the impurity. We find in all regimes a finite impurity-medium entanglement [28,37] which tends to saturate for larger times (t/ω −1 > 1500) besides the tunneling regime where S V N (t) performs an oscillatory motion. Among the investigated regimes the steady bound state and the total reflection regime appearing at large attractive and repulsive g BI couplings experience the smallest amount of entanglement. Indeed, S V N (t) is maximized within the pinning (g BI /Ẽx = 1.0) and the dissipative oscillation regime (g BI /Ẽx = −0.2). Additionally, in the latter response regime the entanglement increases with time and reaches a plateau at around t/ω −1 = 300. During this time-interval the impurity penetrates the bosonic medium twenty times, thereby, enhancing the entanglement at each penetration. In the pinning regime, the system becomes maximally entangled after the impurity penetrates its environment a single time and, subsequently, becomes pinned between the effective barriers raised by the bosonic medium.
To obtain a better understanding of the underlying microscopic mechanisms appearing in the respective response regimes we analyze, in the following, the impurity dynamics with respect to an effective potential [35,37,38], which reads Here, we choose x I 0 = 0 for the harmonic confinement V I (x I ). This effective potential is based on the assumption of a product state ansatz |Ψ MB (t) = |Ψ B (t) ⊗ |Ψ I (t) where the degrees of freedom of the bosonic medium are integrated out. Thus, the effective potential is the superposition of the harmonic confinement V I (x I ) after the quench and the one-body density of the environment weighted with the particle number N B and the impurity-medium interaction strength g BI [38]. Since the bosonic medium remains to a certain extend well localized at the sites of the double-well in the course of the evolution the averaging of the effective potential over the total propagation time T /ω −1 = 200 is justified. Note that even though we considered a product state ansatz for the construction of V eff I (x I ), beyond mean-field effects are included in the one-body density ρ   This picture changes drastically for repulsive impurity-medium interactions. Here, the bosonic medium imprints a potential barrier with two maxima located at the two sites of its actual doublewell. Therefore, in this case, e.g. for g BI /Ẽx = 1.0, the effective potential V eff I (x I ) obtains the shape of a deformed harmonic oscillator having an additional prominent dip at the trap center [ Figure 5(c1)]. As g BI increases to g BI /Ẽx = 2.0 the aforementioned two density maxima of ρ As a next step, we construct for each V eff I (x I ) the single-particle HamiltonianĤ (1) =h 2 2m I ∂ 2 ∂x 2 +V eff I and calculate its first six energetically lowest-lying eigenfunctions ψ eff i . The corresponding absolute squares of the eigenfunctions shifted by their eigenenergies are shown in Figures 5(a1)-(d1). In the following, theses sets of eigenfunctions ψ eff i are taken as basis sets in order to analyze the underlying microscopic mechanisms in the course of the impurity dynamics. The time-dependent probability for the impurity to occupy the state ψ eff i reads where {φ B j } is an arbitrary basis set covering the whole subspace of the bosonic medium. By summing over all basis states of the bath we single out the probability to find the impurity in ψ eff i . Note, that |Ψ MB (t) is the full many-body wave function defined via ML-MCTDX [Eq. (2)], while ψ eff i serves as a basis set to unravel the underlying participating dynamical processes. The occupation probabilities for the respective basis sets [ Figures 5(a1)-(d1)] are presented in Figures 5(a2)-(d2). In order to justify the quality of the basis we sum up all non-vanishing occupation probabilities P 1 i (t) and determine the time at which the sum exceeds and subsequently remains above 0.94 (dashed gray lines).
For g BI /Ẽx = −0.9 the impurity predominantly populates the energetically lowest eigenstate on the left site of the tilted effective double-well V eff I (x I ), while the probability to occupy energetically higherlying states is strongly suppressed as time evolves [see Figure 5(a2)]. Based on this behavior of P 1 i (t), i.e. the spatial localization of the impurity, and the fact that the eigenenergies are negative we associate the bound state formation with the energetically lowest eigenstate of the effective potential V eff I (x I ) [37]. A further analysis of this steady bound state is pro-vided in Appendix B in terms of the involved twobody density. For intermediate attractive impuritymedium couplings corresponding to the tunneling regime V eff I (x I ) has, in contrast to the steady bound state regime, the shape of a nearly symmetric doublewell. Accordingly, the impurity dynamics is mainly determined by the superposition of the two energetically lowest eigenstates of V eff I (x I ) [see Figure 5(b2)]. In the case of intermediate repulsive interactions a pinning of the impurity between the density maxima of the bosonic bath is realized [cf. Figure 5(c2) with g BI /Ẽx = 1.0]. Here, the occupation probabilities start to saturate after t/ω −1 = 50 which in turn leads to the energetically lowest eigenstate of V eff I (x I ) being predominantly populated. However, we observe that the probability to find the impurity in an energetically higher-lying eigenstate is approximately 40 %. In this sense, the broadening of ρ (1) I (x, t) around the trap center [see Figure 2(e1)] can be interpreted as impurity excitations with respect to the effective potential. We refer to those excitations as hidden excitations since they can only be identified by such a microscopic analysis. It is worth noticing that the spatial structure of the first three species functions of the impurity |Ψ I i are in a good agreement with the three energetically lowest eigenfunctions of V eff I (x I ) (see Appendix D). For strong repulsive couplings, e.g. g BI /Ẽx = 2.0, the effective potential exhibits the shape of a triplewell where the two potential barriers stem from the bosonic medium being localized at the sites of the double-well potential. Since these potential barriers are comparatively large, the impurity is totally reflected by the left barrier and predominantly occupies the eigenstates located in the left site of the triple-well [see Figure 5(d2)].
In summary, the analysis of the effective potential enables us to unravel the underlying microscopic mechanisms of the impurity dynamics. In particular, it allows for a deep understanding of the steady bound state formation and proves to be crucial in order to identify the hidden excitations in the pinning regime. Note that for an analogous analysis of the dissipative oscillation regime a much larger set of eigenfunctions has to be taken into account to achieve a comparable quality of the employed single-particle basis (similar to a coherent state in a harmonic oscillator). In Appendix C we provide a discussion of the impurity dynamics in the dissipative oscillation regime where we compare the motion of its mean position to a damped harmonic oscillator.

V. DYNAMICAL RESPONSE REGIMES OF TWO IMPURITIES
To generalize our findings, in the following, we consider two non-interacting impurities (g II = 0) coupled to the bosonic environment. Therefore, all interactions between the impurities are induced by their coupling to the bosonic medium [27,35,77]. The quench protocol is the same as described in section II.
In order to characterize the dynamics of the two bosonic impurities we monitor the time-evolution of their mean position, i.e. their center of mass position, for different impurity-bath coupling strengths g BI [ Figure 6(a)]. Analogously to the one-impurity case we identify five dynamical response regimes depending on g BI and appearing for similar interaction strengths as for N I = 1, see also Figure 3(b). However, we find in the response regime II [ Figure  6(a)] that the mean position is almost constant with x I = 0, whereas in the single-impurity case x I (t) oscillates with a small amplitude around zero. By inspecting the impurity's one-body density ρ (1) I (x, t) for the corresponding regime II (not shown here) it is observed that ρ (1) I (x, t) is almost equally distributed over the effective double-well potential (mediated by the bosonic bath) yielding a mean position close to zero. Moreover, the discrepancy between the tunneling and the pinning regime regarding the impurities' mean position can be revoked by inspecting the variance s analogously to the single-impurity case [see also Figure 3(b)]. Additionally, we observe a broadening of the transition from the pinning to the total reflection regime for increasing g BI with respect to the singleimpurity case where we identified a sharper transition [see Figure 6(a)]. As a case example we shall focus on the steady bound state regime and unravel the microscopic processes in the case of two impurities.
To provide a qualitative understanding of the dynamics in the steady bound state response regime we will describe the impurity dynamics within an effective potential picture, similarly to the previously discussed single-impurity case. For this purpose, we calculate the effective potential V eff I (x I ) [see Eq. (6)] for g BI /Ẽx = −0.9 which we present in Figure 6(b) 1 . Additionally, we compute the effective single-particle eigenstates ψ eff i shifted by their eigenenergies E i [see Figure 6(b)]. A comparison with the effective potential in Figure 5(a1) reveals that the left site of V eff I (x I ) in the two-impurity case is much deeper. Therefore, the two energetically lowest eigenfunctions ψ eff 1 and ψ eff 2 in the two-impurity scenario are located at the left site of V eff I (x I ). Next, we unravel the interplay between the impurities by studying the conditional probabilities to occupy specific eigenstates of their effective potential. In particular, we define P 2 ij (t) as the probability for one impurity to occupy the effective eigenstate ψ eff i while at the same time the other impurity populates the eigenstate ψ eff j . In Figure 6(c) we demonstrate P 2 ij (t) with respect to the five energetically lowest-lying eigenstates of the effective potential for g BI /Ẽx = −0.9 [cf. Figure  6(b)] 2 . Interestingly, we can infer that the two impurities predominantly occupy simultaneously the same energetically lowest eigenstate, which is in accordance with the observations made in the single-impurity case [ Figure 5(a2)]. However, also single-particle excitations in the second eigenstate ψ eff 2 [P 2 12 (t)] as well as two particle excitations contribute to the many-body wave function of the impurities [38,42]. For instance, we find a small but non-vanishing probability for observing two impurities in the second eigenstate P 2 22 (t). Notice that an analogous analysis for the other response regimes, apart from the dissipative oscillation one, leads to similar observations where one or two impurities occupy simultaneously the same or different excited eigenstates of V eff I (x I ) (not shown here). In summary, we have deduced that the dynamical response regimes of two non-interacting impurities are similar to the single-impurity case with a small modification regarding the tunneling regime. Furthermore, we can gain insights into the time-dependent microscopic configuration of the two impurities by investigating the associated conditional probability to find the impurities in two particular eigenstates of their effective potential. Particularly, we exemplified that for the steady bound state regime the impurities predominantly occupy simultaneously the lowest-lying eigenstate of V eff I (x I ). However, we also observed the occurrence of single-and two-particle excitations in higher-lying eigenstates.

VI. SUMMARY AND OUTLOOK
We have investigated the dynamical behavior of bosonic impurities colliding with a BEC trapped in a double-well. The impurities are initially confined in a harmonic oscillator which is spatially displaced with respect to the double-well of the bosonic medium. Upon quenching the harmonic potential to the trap center of the double-well the quantum dynamics is induced such that the impurities collide with the bosonic environment. The correlated non-equilibrium dynamics is tracked with the variational ML-MCTDHX method which enables us to access the full many-body wave function of the system, thereby, including all relevant inter-and intraspecies correlations.
By varying the impurity-medium interaction strength g BI from strongly attractive to repulsive values we are able to control the collisional dynamics of the impurity and identify five distinct dynamical response regimes by inspecting the associated one-body density evolution. These response regimes correspond to the steady bound state regime, the tunneling regime, the dissipative oscillation motion, the pinning regime and the total reflection regime. We demonstrate that they can be easily identified by monitoring the mean position of the impurity. Moreover, by calculating a cross-over phase diagram of the impurity's mean position with respect to the impurity-medium coupling strength we obtain an overview of the emergent dynamical response regimes as a function of g BI and identify smooth transitions between two consecutive ones. Additionally, we explicate the robustness of the response regimes for different parametric system variations, i.e. the intraspecies interaction strength of the bath and the initial displacement of the impurity's harmonic trap.
To provide a better understanding of the involved microscopic mechanisms we employ a time-averaged effective potential picture. By projecting the total many-body wave function onto the eigenstates associated with this effective potential allows us to gain insights into the underlying excitation processes for different interactions. In particular, we find that the impurity is bound to the bosonic medium for strong attractive impurity-bath interaction strengths corresponding to the steady bound state regime and unveil hidden excitations in the pinning regime occurring for intermediate repulsive g BI . We extend our study to the two-impurity case where we showcase the emergence of similar dynamical response regimes as in the single-impurity scenario. Furthermore, we unravel the underlying microscopic mechanisms of the impurities' dynamics analogously to the single-impurity case.
Here, for the steady bound state regime the participation of single-as well as two-particle excitations into energetically higher-lying states of the effective potential is demonstrated.
The results of this work are beneficial for future ultracold atom experiments of impurity-medium scattering for investigating the corresponding collisional channels caused exclusively by presence of the impurity-bath entanglement. Furthermore, this setup can be extended by implementing an additional spindegree-of-freedom for two non-interacting or weaklyinteracting impurities. In this case it would be interesting to identify the individual spin configurations and related spin-mixing processes in dependence of the impurity-medium coupling and whether the impurities evolve as a 'Cooper-pair'. Certainly the generalization of our results to higher dimensions is an intriguing perspective.
Appendix A: Impurity-medium energy transfer processes Due to the initial quench of the harmonic oscillator potential the impurities collide with their bosonic background and, thereby, a transfer of energy to the latter is triggered [28,49,72]. The total energy of the system can be written as E tot = E B tot + E I tot + E int where E σ tot = Ψ MB |Ĥ σ |Ψ MB represents the total energy of species σ ∈ {B, I} and E int = Ψ MB |Ĥ int |Ψ MB the interaction energy between the species.
In order to capture the quench-induced impuritymedium energy transfer we present in Figure 7 the relative energy of species σ defined as E σ rel (t) = E σ tot (t) − E σ tot (0). In each dynamical response regime we observe an energy transfer from the impurity to  its environment. The smallest energy transfer occurs for attractive impurity-medium interaction strengths where, due to the attraction to the bath, the impurity resides closer to the trap center leading to a smaller initial total energy E I tot (0) than in the case of repulsive g BI [compare with the mean position in Figure  3(a)].
For weak impurity-bath coupling strengths of either sign, i.e. g BI /Ẽx = −0.2, 0.3, the impurity continuously transmits energy to its environment until the relative energy of both species eventually saturates for longer times [49,72]. This loss of impurity energy essentially causes its dissipative oscillatory behavior [cf. Figure 3(a)]. The largest energy transfer takes place in the pinning regime. Here, the impurity overcomes the bosonic medium in the left site of the double-well only once and, thereby, transfers a large amount of energy to the bosonic medium such that the impurity becomes pinned within the latter.
Appendix B: Two-body correlation in the steady bound state regime To elucidate the interplay between the bosonic medium and the impurity in the steady bound state regime in more detail we perform an analysis of the impurity-medium reduced two-body density ρ This quantity provides information about the probability of finding the impurity at position x I and one particle of the bosonic background located at x B . A snapshot of the one-body density ρ (1) σ (x σ ) with σ ∈ {B, I} and the two-body density ρ (2) B,I (x I , x B ), respectively, at t/ω −1 = 150 and for g BI /Ẽx = −0.9 is depicted in Figures 8(b) and (c). Here, we find the impurity to be localized at the left maximum of the bosonic medium corresponding to the left site of the double-well which agrees with the observations made for the single-impurity case in Figures 5(a1) and (a2). Furthermore, the two-body density indicates that the probability to find one particle of the bath at the left site of the double-well, i.e. close to the impurity, is enhanced compared to the respective probability for the right site. In particular, the diagonal of ρ  Figure 8(c)] represents the probability to capture the impurity and one particle of the environment at the same position which we will refer to as ρ   Figure 2(a1)] and designates a high probability for the impurity and one particle of the bosonic medium to be at the same location [35,37].
Appendix C: The dissipative oscillation response regime: effective mass and damping of the Bose polaron Let us also analyze the dissipative oscillation regime in the case of a single-impurity in more detail. As observed in Figure 3(a), the mean position of the impurity for weak impurity-medium couplings exhibits a damped oscillatory behavior. Therefore, in the following we compare the analytical solution of a damped harmonic oscillator with the mean position x I (t) and mean momentum p I (t) obtained within the ML-MCTDHX method (see section III). The equation of motion of a particle subjected to a damped harmonic oscillator [78] reads where γ eff denotes the effective damping constant, ω eff the effective trapping frequency and m eff refers to the effective mass of the impurity. Here, we interpret the impurity as a quasi-particle, namely a Bose polaron, which is dressed by the excitations of its surrounding and moves in an effective harmonic oscillator. The mean position for a particle obeying Eq. (C1) reads . Additionally, we assume that the particle is initially at rest, i.e. p 0 = 0, and shifted by x eff 0 = x I (0) . The corresponding mean momentum of Eq. (C1) is accordingly written as Subsequently, we fit the analytical results of Eqs. (C2) and (C3) to the mean position and mean momentum calculated from the ML-MCTDHX approach for the free parameters γ eff , ω eff and m eff [37,51].
In Figures 9(b)-(d) we present the fitted parameters for impurity-medium interactions ranging exemplary from g BI /Ẽx = −0.2 to 0.2. We find that the effective mass [28,51] of the impurity decreases for larger absolute values of g BI , whereas at the same time the damping constant increases. The increase of γ eff for increasing attractive and repulsive couplings can be explained by the corresponding growing influence of the bosonic environment. Furthermore, we find an approximately linear decrease of the effective frequency ω eff . In order to intuitively explain this behavior we show in Figure 9(a) the time-averaged effective potential for the two considered extrema of g BI (i.e. g BI /Ẽx = −0.2 and g BI /Ẽx = 0.2). As can be seen, the effective potential is deeper in the case of attractive g BI compared to the one obtained for repulsive g BI leading to a higher effective frequency ω eff in the attractive case than in the repulsive one.  purity does not solely occupy the energetically lowest eigenstate of the effective potential, but also populates energetically higher-lying states [Figures 5(c1) and (c2)]. In the main text we referred to these states as hidden excitations since they are not apparent by merely inspecting the evolution of the one-body density ρ In the following we aim at investigating these hidden excitations in more detail. To this end, we monitor the time-evolution of the density associated with the individual species functions |Ψ I j for g BI /Ẽx = 1.0 [ Figures 10(a2)-(c2)]. As it can be readily seen, the densities of the species functions remain constant once the impurity is pinned. Moreover, a careful inspection of the densities of the first three species functions reveals an ascending number of nodes which is tantamount to the existence of energetically higher-lying excitations of |Ψ I j . In order to unravel the structure of the aforementioned species functions we project them on the basis set consisting of the eigenfunctions ψ eff i of the effective potential [cf. Figure 5(c1)] and take the absolute square of the respective overlap, i.e. P i;j (t) = | Ψ I j (t)|ψ eff i | 2 , where j = 1, 2, 3. Indeed, we find that the three dominantly occupied species functions cor- respond to the three energetically lowest eigenfunctions of the effective potential, see Figures 10(a1)-(c1). In particular, |Ψ I 1 (x I , t)| 2 matches with the energetically lowest eigenstate ψ eff 1 , whereas |Ψ I 2 (x I , t)| 2 and |Ψ I 3 (x I , t)| 2 correspond to the second and third eigenstate, i.e. ψ eff 2 and ψ eff 3 .
Appendix E: Dynamical response for a triple-well trapped environment To extend our basic conclusions regarding the impurity's response, described in section II, we replace the double-well of the bosonic medium with a triple-well. However, the harmonic trap of the impurity as well as the employed quench protocol to induce the dynamics remain unchanged. The triple-well of the bosonic environment reads V B (x B ) = m B ω 2 B (x B ) 2 /2+g − (x B )+ g + (x B ), where a superimposed harmonic trap with frequency ω B /ω = 0.15 and two Gaussians g ± ( shifted by ∆ from the trap center are used. Also, the Gaussians have a width of w B /x = 0.8 and a height of h B /Ẽx = 1.8 while the displacement is ∆/x = 2.5. The system consists of N B = 20 bosons for the bosonic medium and a singleimpurity N I = 1. Additionally, we employ M = 6 species and d B = d I = 6 single-particle functions for the calculations to be presented below. Figure 11(a) shows the time-evolution of the impurity's mean position corresponding to the six identified dynamical response regimes which are labeled as I-VI in the respective crossover diagram in terms of g BI illustrated in Figure 11(b). As in the double-well case [see Figure 3(a)], we find that for strong attractive impurity-medium interactions (regime I) a localization of the impurity in the vicinity of the most left site of the triple-well occurs. This behavior is attributed to the initial large overlap of the impurity with the bath on the left site 3 . By decreasing g BI to intermediate attractive couplings, i.e. g BI /Ẽx = −0.5 corresponding to the regime II, we observe that the impurity localizes at the central site of the triple-well. Entering the weak attractive and repulsive impurity-medium coupling strengths we identify a dissipative oscillation of the impurity similar to the case in which the bosonic bath is trapped to a double-well [ Figure 3(b)]. For stronger repulsive g BI the aforementioned oscillatory character vanishes (regime IV) and the one-body density of the impurity ρ (1) I (x, t) (not shown) exhibits two humps located at the two maxima of the triple-well.
Here, x I (t) tends to zero, e.g. for g BI /Ẽx = 0.7. Note that even though the regimes II and IV show a similar behavior in terms of x I (t) we can distinguish them by evaluating the variance (s) which is in the attractive case smaller than in the repulsive one (not shown here). A further increase of the impuritybath repulsion to g BI /Ẽx = 1.2 leads to a localization of the impurity at the position between the left and central site of the triple-well [cf. regime V in Figure  11(b)]. For strong repulsive impurity-medium inter-action strengths we reach the regime VI which corresponds to the total reflection of the impurity as in the double-well case [38].
In summary, similarly to the double-well scenario we observe for a triple-well continuous transitions between the emergent dynamical response regimes of the impurity with respect to g BI . Interestingly, we find that the behavior of the impurity can be steered and controlled via the impurity-medium coupling strength, where the emerging response regimes have a character similar to the double-well case.