{'fraction_non_alphanumeric': 0.0392604873427421, 'fraction_numerical': 0.011507864636304108, 'mean_word_length': 4.574557315936627, 'pattern_counts': {'":': 0, '<': 9, '<?xml version=': 0, '>': 8, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Counterflow dynamics of two correlated impurities immersed in a bosonic gas

The counterflow dynamics of two correlated impurities in a double-well coupled to an one-dimensional bosonic medium is explored. We determine the ground state phase diagram of the system according to the impurity-medium entanglement and the impurities two-body correlations. Specifically, bound impurity structures reminiscent of bipolarons for strong attractive couplings as well as configurations with two clustered or separated impurities in the repulsive case are identified. The interval of existence of these phases depends strongly on the impurity-impurity interactions and external confinement of the medium. Accordingly the impurities dynamical response, triggered by suddenly ramping down the central potential barrier, is affected by the medium's trapping geometry. In particular, for a box-confined medium repulsive impurity-medium couplings lead, due to attractive induced interactions, to the localization of the impurities around the trap center. In contrast, for a harmonically trapped medium the impurities perform a periodic collision and expansion dynamics further interpreted in terms of a two-body effective model. Our findings elucidate the correlation aspects of the collisional physics of impurities which should be accessible in recent cold atom experiments.

Recently, more attention has been placed on the interplay between several impurities immersed in a quantum gas [28][29][30][31]. Among others, the coalescence of two bosonic impurities coupled to a harmonically trapped bosonic medium has been predicted [32] as well as the existence of their induced interactions [33][34][35]. In the strongly attractive interaction regime, the formation of bipolarons referring to impurity bound states was also unraveled [36,37]. Beyond these studies the non-equilibrium dynamics of quasi-particles following an interaction quench [38][39][40] has been examined unveiling, in particular, energy redistribution processes, temporal orthogonality catastrophe phenomena and the effective temperature of the impurities [39] e.g. by emulating pump-probe and Ramsey spectroscopy.
Another branch in the field of ultracold quantum gases concerns the collisional aspects of atomic ensembles. Counterflow dynamics can be triggered, e.g., by employing a magnetic field gradient separating two atomic hyperfine states [41] or releasing an ultracold quantum gas from a double-well potential into a harmonic oscillator [42]. For a single atomic species these protocols result in the oscillation of the formed dark solitons [42] or in the case of a two-component mixture in the spontaneous generation of dark-bright soliton trains [41]. Another technique to initiate atomic collisions constitutes of two counter-propagating harmonic oscillator potentials [43] which has been experimentally realized with 40 K and 87 Rb clouds utilizing two optical tweezers [44,45].
In this sense, it is intriguing to explore the counterflow correlated dynamics of impurities in combination with a superimposed superfluid background. A similar question was addressed for fermions [46,47] e.g. showing the formation of shock-waves. Thereby, of immediate interest is the influence of the background on the collisional response and the associated emergent induced interactions between the impurities [47]. The impact of the bath on the impurity dynamics is expected to depend on the confining potential of the bath, and the impurities coupling strength as well as the interaction between the impurities and the bath particles [48]. To tackle these open questions, herein we consider a minimal model of two bosonic impurities trapped in a double well and immersed in a bosonic bath. The counterflow dynamics between the impurities is induced by suddenly ramping down the potential barrier of their double well and, subsequently, let the system evolve in time for different interaction configurations.
Specifically, it is shown that already the ground state configurations depend on the chosen impuritymedium and impurity-impurity interaction strength and, importantly, on the type of the underlying trapping potential of the medium. For instance, in the case of a box-confined medium the impurities coalesce for intermediate repulsive impurity-medium couplings independently of the impurity-impurity repulsion. On the other hand, for a harmonically trapped bath the impurities separate from each other for strong impurity-impurity and impurity-medium repulsions [32]. Moreover, we find indications of bipolaron formation [36] for strong impurity-medium attractions. The dynamical response of the impurities appears to be strongly affected by the combination of the involved interaction strengths as well as by the external confinement of the bath. More precisely, in the case of a box-confined medium and intermediate repulsive impurity-medium couplings the impurities induced interactions lead to their localization at the trap center after their first collision. However, increasing the impurity-medium interaction strength the impurities experience a periodic collisional response characterized by a damped amplitude: a behavior that is argued to be governed by finite size effects determined by the size of the box potential. Employing a harmonically trapped medium the impurities localize at the trap center for attractive impurity-medium coupling strengths, while for intermediate and strong repulsions they phase separate with the medium [49,50] and thus their dressing is suppressed. Considering weakly repulsive impurity-impurity interactions a state transfer manifests from two separated to two coalesced impurities. Importantly, this process is absent in the decoupled case elucidating the role of the coupling with the bath and thus of the interspecies correlations (entanglement).
A microscopic analysis provides insights into the single-particle excitation processes and the two-body states participating in the dynamics and the aforementioned state transfer.
To describe the stationary and dynamical properties of the composite impurity-medium system we employ the multi-layer multi-configuration time-dependent Hartree method for atomic mixtures (ML-MCTDHX) [51][52][53][54]. This ab initio approach allows us to efficiently track the relevant inter-and intraspecies correlations which are anticipated to be enhanced, especially during the dynamics. This is in part due to the few-body impurity subsystem as well as the spatial inhomogeneity caused by the external potential.
This work is structured as follows. In section II and III we present the impurity model under consideration and introduce the ingredients of the variational method, respectively. We proceed in section IV with an analysis of the system's ground state and draw a phase diagram with respect to variations of the impurity-impurity and impurity-medium interaction strengths. This analysis is based on the two-body densities quantifying the correlations of the bath particles and the impurities. Next, in section V and VI the dynamical response of the system following a sudden reduction of the double-well barrier is discussed.
In particular, section V elaborates on the case of a box-confined medium in which we explicate, e.g., the localization of the impurities at intermediate repulsive impurity-medium coupling strengths and the emergence of finite size effects. The case of a harmonically trapped medium is investigated in section VI where the focus is set on the impurities excitation processes and their dependence on the impurity-medium couplings. Our results are summarized in section VII together with an outlook regarding further research directions. Appendix A elaborates on the behavior of the impurities relative distance in their ground state and in Appendix B the impact of the impurities mass on their collisions is exposed. In Appendix C we discuss the persistence of the impurities collisional features when a linear ramp is applied to the barrier height of the double-well potential.

II. IMPURITY-MEDIUM SETTINGS
The system consists of two different bosonic species B and I at ultracold temperatures. In particular, we consider N I = 2 impurities of mass m I and a bosonic bath of N B = 20 particles with mass m B . The corresponding Hamiltonian readŝ is the interaction Hamiltonian of species σ ∈ {B, I}. Each component is subject to a different external potential V σ (x σ i ), a scenario that can be achieved via species selective optical potentials [55,56]. It is also restricted to one spatial dimension [57] that can be realized experimentally, e.g., by freezing out the transverse degrees of freedom using a strong harmonic confinement [58,59].
Since we are operating in the ultracold regime it is sufficient to take into account only s-wave scattering processes and thus the interaction between two particles of the same species is modeled with a contact interaction potential [9] determined by the one-dimensional effective coupling strength parameter g σσ . Analogously, the coupling between the impurities and the bath is described through a contact inter- , where g BI denotes the impurity-medium interaction strength. Due to the fact that g σσ with σ, σ ∈ {B, I} depends, among others, on the three-dimensional swave scattering length it can be experimentally adjusted, e.g., via Feshbach resonances utilizing either magnetic or optical fields [1,10,60,61]. Below, we consider a bosonic medium of 87 Rb atoms and 133 Cs impurities. Thus, the mass ratio is m I /m B = 133/87 [5,55,62,63].
At t = 0 the system is prepared in its ground state with a specific combination of interaction strengths (g σσ ). The impurities are initially confined in a double well V dw ) which is the superposition of a harmonic oscillator potential Figure 1. Schematic representation of the considered setup and the quench protocol. Two bosons (red circles) are coupled to a bosonic medium (blue shaded area) which is either confined in (a) a box potential or (b) a harmonic oscillator. First, the two impurities are loaded into the double well and the equilibrium state of the system is obtained for a specific set of inter-and intraspecies interaction parameters (gBB, gII , gBI ). Subsequently, a counterflow dynamics between the two impurities is induced by suddenly ramping down the barrier of the double-well potential.
with frequency ω I and a Gaussian of width w I and height h I [64,65]. For the external potential of the bosonic medium we consider two cases: a box potential of size L B = 1 with V B (x) = 0 for −L B /2 < x < L B /2 and V B (x) = ∞ elsewhere, and a harmonic os- In the former scenario we considerx box = L B /10 andẼ box = 2 m Bx 2 box as characteristic length and energy scales, respectively [66]. Thus, the time and interaction strength are expressed in units oft box = m Bx 2 box andg box = 2 m Bxbox , respectively. However, in the case of a harmonically trapped medium it is more convenient to express the energy in units of the medium's harmonic oscillatorẼ ho = ω ho where ω B /ω ho = 1. It follows that the length, interaction strength and timescales are given in terms of , respectively. To construct the impurities' double well we employ ω I /ω box,ho = 0.6, h I /Ẽ box,hox −1 box,ho = 3.0 and w I /x box,ho = 0.7.
The ground state of the composite system is acquired for a set of values of interaction strengths (g BB , g II , g BI ). In the case of a box-confined medium the impurities always exhibit a finite spatial overlap with the bath particles for the considered interaction strengths. However, for a harmonically trapped bath the impurity-medium overlap vanishes as long as g BI > g BB (see also the discussion below). Subsequently, the dynamics is triggered by suddenly reducing (at t = 0) the barrier of the impurities' double-well potential [see Figure 1]. Consequently, in the course of the time evolution the impurities collide and experience a harmonic oscillator potential V ho I (x) with frequency ω I . As we shall argue below, the emerging collisional correlated dynamics depends strongly on the initial phase of the system determined by the interaction parameters (g BB , g II , g BI ).

III. VARIATIONAL APPROACH AND WAVE FUNCTION ANSATZ
To determine the time-dependent solution of the problem described by the Hamiltonian of Eq. (1) we invoke the ML-MCTDHX method [51][52][53][54]. This approach is an ab-initio one and optimizes a chosen basis e.g. in terms of the Dirac-Frenkel variational principle [67]. In particular, the basis set which underlies the many-body wave function |Ψ MB (t) is characterized by a time-dependent and multi-layered structure with individual truncation parameters [68]. Firstly, the many-body wave function is expanded into distinct sets of species functions {|Ψ σ i (t) } Dσ i=1 with D σ denoting the number of the latter for species σ ∈ {B, I}. Since here we consider a two-component mixture, |Ψ MB (t) is firstly expressed into two such basis sets and, thus, can be written in the form a truncated Schmidt decomposition [69][70][71] where D = D B = D I and |Ψ σ i (t) are the so-called natural species functions [53]. The time-dependent Schmidt coefficients λ i (t) determine the population of the i-th natural species function and provide information about the interspecies entanglement [71,72]. For instance, in the case that only a single Schmidt coefficient λ i (t) is nonzero, the system is described by a direct product ansatz of species functions indicating the absence of entanglement. On the other hand, the two species are considered to be entangled when more than one Schmidt coefficients are nonzero.
In the next step of the many-body wave function |Ψ MB (t) truncation, each species function is expanded into time-dependent permanents Here, each permanent represents one of the Nσ+dσ−1 Nσ possible configurations to distribute N σ particles on d σ single-particle functions |ϕ σ j (t) . A further imposed condition is that the number of occupied singleparticle functions for each permanent has to be equal to N σ (indicated by n|N σ ). This expansion enables us to account for intraspecies correlations. Finally, the time-dependent single-particle functions |ϕ σ j (t) are expanded into a time-independent discrete variable representation [73], which we choose here to consist of 300 grid points in an interval {−5, 5} in units of x box,ho . Additionally, in this work we employ D = 6 species functions and d A = 4, d B = 6 single-particle functions for the accurate calculation of the considered systems.
Especially, the multi-layered architecture and the time-dependent basis of the many-body wave function mainly contribute to the high degree of flexibility of the method which enables |Ψ MB (t) to approach the accurate solution for each time instant with a high fidelity even for systems containing a mesoscopic particle number. In this way, the ML-MCTDHX method keeps the number of required wave function coefficients within a computational feasible limit and, at the same time, accounts for the relevant inter-and intraspecies correlations.

IV. CHARACTERIZATION OF THE GROUND STATE
In the following we provide an overview of the ground state characteristics of two impurities trapped in a double-well and coupled to a bosonic bath confined either in a box potential or a harmonic oscillator. These ground states will subsequently serve as a starting point for examining the counterflow impurity dynamics immersed in a medium which will be discussed below in Section V. Unless stated otherwise, the interaction strength between the bath particles is kept fixed to g BB /g box = g BB /g ho = 0.5.

A. Main observables of interest
Let us first introduce the quantities that will be employed for the identification of the ground state phases and the quench dynamics of the two interacting impurities in the cases of a box-confined and a harmonically trapped medium. The distinction between the emergent ground state configurations is performed with respect to the two-body density distributions of the impurities and the bath particles at t = 0. The reduced two-body density of two particles of the same species is given by whereΨ ( †) σ (x σ 1 ) denotes the bosonic field operator which annihilates (creates) a particle of species σ ∈ {B, I} at position x σ 1 . In fact, ρ σσ (x σ 1 , x σ 2 ) is the probability of finding one particle at x σ 1 and, simultaneously, another particle of the same species at x σ 2 . In the following, we will drop the time parameter since for the ground state t = 0.
In both considered external confinements of the medium, an increase of the repulsive impuritymedium coupling strength leads to the development of interspecies correlations (entanglement) which eventually impact the ground state configurations [74]. A common measure for quantifying entanglement in a bipartite system is the von Neumann entropy [71,75], defined as Recall that D denotes the number of the employed species functions and λ i are the Schmidt coefficients [cf. Eq. (2)]. For a maximally entangled mixture the von Neumann entropy obtains its maximum value, S vN max = ln D and λ i = 1/D. This value corresponds to S vN max = 1.79 in our case. In contrast, a vanishing von Neumann entropy indicates a decoupled (i.e. nonentangled) mixture such that the total many-body wave function can be written as a direct product state of the two individual species wave functions.
Moreover, in order to judge the degree of miscibility among the impurity and medium clouds we calculate the interspecies spatial overlap [76,77] which is quantified through Here, ρ |Ψ MB being the one-body density of σ ∈ {B, I} species [78].
B. Decoupled case gBI = 0 Before addressing the ground state properties of the coupled mixture, we focus on the simpler scenario where the bath and the impurities are decoupled from each other (g BI = 0) and thus they can be treated individually. Accordingly, the impurity-bath entanglement is vanishing, i.e. S vN = 0. Then, the system reduces to two interacting bosons in a double well [79][80][81] with the bath being homogeneous or harmonically trapped [cf. Figure 2(a) and Figure 3(a) for g BI = 0]. Here, we distinguish between weakly and strongly coupled impurities. In the former case, the two impurities are delocalized over the two sites of the double well, see the dominant population of the off-diagonal compared to the diagonal elements of ρ  Figure 2 On the other hand, for larger impurity-impurity repulsion g II the density maxima along the diagonal vanish and only density peaks at the off-diagonal remain [cf. inset of Figure 2  II (x I 1 , x I 2 ) is described by the conditional probability of finding one impurity at the left and one impurity at the right site of the double well or vice versa. In this sense, the impurities tend to separate from each other and are anti-correlated within the same site of the double well. Thus, they reside in a Mott-type state. Notice that for a decoupled mixture the above described impurity configurations corresponding to regime (II) and (III) occur independently of the particular trapping geometry of the bath [cf. g BI = 0 in Figures 2(a) and 3(a)]. Thereby, the medium extends almost homogeneously over the box potential as also reflected by the shape of its two-body density [see inset of Figure 2 Having analyzed the spatial configurations of the interacting impurities for a suppressed impurity-bath coupling we then discuss the ground state properties of the composite system when the bath is confined in a box potential and the interspecies interaction strength g BI becomes finite. The respective ground state phase-diagram is mapped out and presented in Figure 2(a) based on the underlying two-body configurations identified in ρ . All crossovers among the involved phases are smooth. For each inset the range of the color map is optimally chosen and maximally extends from 0 to 2 [0.2] for ρ (2) In panels (b1)-(d2) the ground state two-body densities of the bath and the impurities are presented in terms of a weakly interacting bath (gBB/g box = 0.1) and strongly interacting impurities (gII /g box = 2.0). The interspecies interaction strength from top to bottom row is gBI /g box = 5.0, 0.2, −2.0, respectively.
Overall, we find that upon varying g BI and g II , the system deforms smoothly across the different phases, which are analyzed in detail in the following, see also [82].
In the case of strong repulsive g BI , the mixture enters regime (I) [see Figure 2(a)] [83]. Here, the impurities two-body density exhibits two peaks along its diagonal meaning that the two impurities occupy simultaneously a single site of the double-well. Such a behavior is referred to as the coalescence of the impurities and has been observed also for the case where the impurities and bath particles are both harmonically confined [32]. Intuitively, we explain this behavior as follows. One impurity lying at a specific site of the double well repels the bath particles and, thereby, creates an effective hole which attracts the other impurity [84,85]. On the other hand, the impurities impact accordingly the bath. This backaction manifests, for instance, in the off-diagonal parts of the medium's two-body density which exhibits strongly suppressed spatial regions at the location of the impurities [see inset in regime (I) of Figure 2(a)]. Indeed, the probability to find two bath particles at positions corresponding to opposite double-well sites is vanishing. This is due to the fact that the impurities lie both either at the left or at the right double-well site as it becomes apparent from their reduced two-body density. However, a configuration where two bath particles reside simultaneously at the same double-well site is still conceivable, assumed that the impurities are at the opposite site, thereby, avoiding the bath particles [note the non-vanishing density at the diagonal of ρ in the inset of Figure 2(a)]. Increasing the impurity-medium interaction strength within regime (I) for a fixed g II /g box ∈ [0, 2] we observe two prominent features appearing in terms of ρ . Firstly, the two-body density holes at the off-diagonal of ρ 2 ) become more pronounced for increasing g BI and, secondly, for a g BI /g box 2.5 two bath particles are correlated at the most right and most left or at opposite sites of the bath cloud [see the outermost density peaks at the diagonal and off-diagonal elements of ρ Figure 2(a)]. From this latter behavior we can conclude that the bath particles exhibit two-body long distance correlations. Moreover, we note that in case of strong impurity-impurity repulsions, e.g. for g II /g box = 2.0, the two-body state of the impurities begins to fermionize and the diagonal peaks of ρ (2) II (x I 1 , x I 2 ) broaden and, eventually, fragment [86] [see corresponding inset of Figure 2 Similarly to the coalescence of the impurities in the repulsive case [regime (I)] also in the attractive scenario the two impurities simultaneously occupy either the left or right site of the double well, as it can be deduced from their diagonal and highly localized twobody density configuration ρ [see the inset of Figure 2(a)]. In both the repulsive and the attractive cases, the bath mediates an induced attractive interaction between the impurities such that the latter coalesce and tend to occupy the same double-well site. For a more detailed discussion regarding the presence of the attractive induced interactions between the impurities via their relative distance, see Appendix A. Furthermore, due to the attractive interactions the bath particles localize in the vicinity of the impurities such that also the two-body density of the medium exhibits two dominant peaks along the diagonal [ Figure  2(a)]. Also, we remark that regime (II) in Figure 2(a) includes ground states corresponding to delocalized impurities, i.e., where the diagonal and off-diagonal elements of ρ are simultaneously populated. However, with varying g II and g BI the particular density peaks are deformed compared to the depicted insets of regime (II) in Figure 2(a). For instance, in the case of g II /g box = 2.0 and g BI /g box = 1.0 corresponding to regime (II) in Figure 2(a) the two-body density of the impurities exhibits peaks at its off-diagonal elements [similar to regime (III)] and fragmented density humps at its diagonal [as observed in regime (I)].
Moreover, increasing (decreasing) the interspecies interaction to large repulsive (attractive) values leads to a noticeable growth of the von Neumann entropy. Namely, the impurities become highly entangled with the bath [cf. Figure 3(b)]. At the same time the impurities and the bath share a finite spatial overlap with each other for all the considered values of g BI and g II [ Figure 3(c)]. Therefore, since these two features constitute a basic requirement for the formation of quasi-particles, e.g. as discussed in Refs. [24,39,50,87], in principle, the impurities can be dressed by the excitations of the bath and, thus, form Bose polarons.
Additionally, we investigated the ground state of the system for a weakly interacting bath, e.g. for g BB /g box = 0.1. As we will argue, the enhanced compressibility of the bath can alter the ground state configurations and this is evident, among other observables, in terms of the two-body density. In Figures 2(b1)-(d2) we exemplary present the two-body densities of the impurities and the bath particles for g BB /g box = 0.1 and g II /g box = 2.0. Here, the impurities two-body densities exhibit a qualitatively similar structure with the ones corresponding to a moderately interacting bath [ Figure 2(a)]. For instance, a localization at the sites of the double-well is observed for strongly attractive g BI [ Figure 2(d1)] as well as an anti-correlated behavior for weak impuritymedium couplings [ Figure 2(c1)] and the coalescence of the impurities for strong g BI [ Figure 2(b1)]. Similarly, the two-body density of the weakly interacting medium resembles the one of a moderately interacting bath in the cases of weak attractive and repulsive as well as strong attractive g BI [Figures 2(c2) and (d2)]. However, for strong repulsive g BI the medium's twobody density is modified for weak g BB , i.e., the offdiagonal parts of ρ are depopulated and only the diagonal ones are occupied [ Figure 2(b2)]. We attribute this behavior to the increased compressibility of the bath which suppresses correlations between two bath particles residing at longer distances e.g. the opposite edges of the cloud as observed in Figure 2(a). Additionally, this behavior is accompanied by strong anti-correlations between the impurities and the bath particles, i.e., the respective two-body density ρ exhibits only peaks at its off-diagonal (not shown).
Furthermore, we have found that in the corresponding ground state phase diagram of a weakly interacting medium (g BB /g box = 0.1) with N B = 20, the regimes (II) and (III) shrink as compared to the g BB /g box = 0.5 case shown in Figure 2(a). Specifically, their phase boundaries are shifted towards the line corresponding to g BI = 0. A similar, but less pronounced, shift of the phase boundaries is observed when increasing the number of bath particles to N B = 30 but keeping g BB fixed. Summarizing, both decreasing g BB or increasing N B while considering fixed all other parameters leads to an enhancement of the magnitude of the attractive induced interactions between the impurities when g BI is switched on towards finite attractive or repulsive values.

D. Harmonically trapped medium
We then proceed to analyze in more detail the system consisting of a harmonically trapped medium. This change of the external confinement reduces the mobility of the bath particles which are then naturally bounded by the harmonic oscillator around the trap center. The respective phases presented in Figure 3(a) feature smooth crossovers among them besides the one between the regimes (I) and (III) which is abrupt. To testify the "smoothness" of the underlying crossover regions we track, as in the box-confined scenario, the behavior of the impurities two-body densities (cf. [82]).
An increasing impurity-medium repulsion such that g BI > g BB leads to a phase separation between the impurities and the bath particles as it is captured by the diminishing spatial overlap depicted in Figure  3(b). In this case the impurities are no longer dressed by the excitations of the bath and, thus, the quasiparticle notion is essentially lost [39,49,50]. Thereby, we distinguish between two cases according to the impurity-impurity interaction strength. In the case of weak g II and strong g BI corresponding to regime (I) in Figure 3(a) the impurities coalesce in a similar manner as described above (see also Ref. [32]). However, as the impurity-impurity interaction strength becomes large enough or the impurity-medium repulsion sufficiently small, regime (III) is accessed in which the impurities spatially separate. This is identified by the exclusive population of the off-diagonals of ρ (2) Figure 3(a)]. Turning to strongly attractive impurity-medium interaction strengths assigned as regime (IV) in Figure  3(a) a localization of the impurities in the barrier of the double-well is observed, see in particular the elongation of ρ (2) II (x I 1 , x I 2 ) along its diagonal. This property is related to the formation of a bipolaron, referring to a dimer bound state consisting of two polarons [35,36,39]. We base our argument of bipolaron formation on the following observations which have also been used in Refs . [36] to expose the existence of such states in three-dimensions. The continuous decrease of the so-called bipolaron energy E BP = E 2 − E 1 + E 0 has been verified for increasing impurity-medium attraction, where E i denotes the total energy of the bosonic gas containing i = 0, 1, 2 impurities. In the same manner, also the size of the dimer state quantified in our case by 1/ r 2 II increases for larger g BI (not shown here).
Concluding, let us mention in passing that similarly to the case of a box-confined medium, the main requirements for the formation of Bose polarons are also fulfilled in the presence of a harmonic trap. The only exception consists of the region of phase-separation among the impurities and the medium at g BI > g BB . Thus, the impurities response, to be presented below, can be interpreted as the counterflow correlated dynamics of two quasi-particles, here Bose polarons.

V. COLLISIONAL MANY-BODY DYNAMICS FOR A BOX CONFINED MEDIUM
Next, we investigate the time evolution of the composite system upon suddenly ramping down the central barrier of the impurities' double-well potential such that they are henceforth externally confined in a harmonic oscillator and, thus, their counterflow dynamics is triggered. In a decoupled mixture (g BI = 0) this quench results in an undamped periodic impurities motion where they collide and subsequently expand repeatedly. Turning to finite impurity-medium couplings the response is substantially altered and depends strongly on the trapping potential of the medium, as we will argue below.
A. Response through the time-evolution of the density We monitor the counterflow dynamics of two impurities coupled via g BI to a bosonic medium trapped in a box potential. After ramping down the potential barrier of the double well, the impurities are left The insets depict the two-body densities at the same time instants and scales, but for two strongly interacting impurities, i.e. gII /g box = 2.0. At those time instants, the one-body densities for two weakly and strongly interacting impurities reveal similar features, thus, allowing the comparison on the two-body density level.
to evolve in the resulting harmonic trap. As a first step, we categorize the emergent dynamical response regimes by inspecting the one-body densities ρ I (x, t) and ρ Figures 4(a1)-(c1) and (a2)-(c2), respectively. The impurity-impurity coupling is kept fixed g II /g box = 0.2 and only the impuritymedium interaction is varied. We are able to identify four distinct dynamical response regimes taking place at strong attractive (g BI /g box < −0.8), weak attractive and repulsive (−0.8 g BI /g box < 0.8), intermediate repulsive (0.8 g BI /g box 2.0) and strongly repulsive (2 < g BI /g box ) values of g BI . These dynamical regimes are, of course, inherently related to the corresponding phases unraveled in the ground state of the system, see Figure 2(a). Note that the behavior of the one-body densities in the respective regions do not qualitatively alter for varying g II from small to large repulsive values at least in the range of 0 ≤ g II /g box ≤ 2.0 considered herein. Only by inspecting higher-body observables, such as the reduced two-body density, reveals significant alterations of the impurities response regarding variations of g II .

B. Dynamics for weakly attractive and repulsive impurity-medium couplings
In the case of either weakly attractive or repulsive g BI the impurities one-body densities feature a periodic motion consisting of a collision and an expansion of their cloud [cf. Figure 4(a1) for g BI /g box = −0.2]. As a consequence, the bath is only weakly perturbed from its initial homogeneous profile showing small amplitude distortions at the instantaneous location of the impurities [see Figure 4(a2)]. This response emerges when considering initial configurations corresponding to the interaction regimes (II) and (III) discussed in Figure 2(a). Interestingly, the timeevolution of the one-body density does not depend strongly on variations of the impurity-impurity interaction strength [88]. Therefore, one has to rely on two-body observables, such as the two-body density ρ (2) II (x I 1 , x I 2 ), e.g. presented in Figure 4(a3) for g BI /g box = −0.2 and g II /g box = 0.2 at t/t box = 10.5 and in the respective inset for two strongly interacting impurities with g II /g box = 2.0. In the former case, the two weakly interacting impurities are initially and throughout the evolution delocalized over both sites of the double well since both the diagonal and the offdiagonal elements of ρ are non-vanishing. However, for strongly interacting impurities we find that in the course of the evolution a pronounced correlation hole occurs [39], i.e., solely the off-diagonal of ρ (2) Figure 4(a3)]. In other words, due to the strong repulsion the impurities reside at spatially opposite positions and avoid each other during the dynamical evolution. The respective dynamics of the impurities is characterized by the periodic expansion and contraction of their cloud around the trap center while avoiding to reside at the same location, see the correlation hole of ρ (2) II (x I 1 , x I 2 ). However, independently of g II the impurities remain within the medium, thus, forming a polaron due to the finite g BI .

C. Time-evolution for intermediate interspecies repulsions
Increasing the impurity-medium repulsion to intermediate values (g BI /g 2.0) and, thereby, entering regime (I) in Figure 2(a) a comparatively altered response is realized. Indeed, once the two impurities collide at the trap center, they remain localized [89], see Figure 4(b1) for g BI /g box = 1.5. Consequently, because of the repulsive character of the employed g BI , the bath is pushed towards the edges of the box with each density branch undergoing weak amplitude oscillations due to its reflection from the walls of the confining box [ Figure 4(b2)]. Simultaneously, the bosonic medium becomes highly depleted, meaning that higher-lying natural orbitals (being the eigenvalues of the reduced one-body density matrix ρ (1) B (x, x , t)) are macroscopically populated. Accordingly, the bosonic gas is correlated and deviates from a perfect BEC. This behavior is in contrast to the case of g BI /g box < 0.8 where the first orbital is dominantly occupied.
However, the corresponding two-body density ρ (2) II (x I 1 , x I 2 ) when the impurities collide at the trap center at t/t box = 4.5 is shown in Figure 4(b3). The elongated shape of ρ (2) II (x I 1 , x I 2 ) along the diagonal indicates the presence of an attractive interaction between the impurities induced by the coupling with the bath [32,39,50]. Apparently the strength of this induced attraction is larger than for g BI /g box = 0.2. Moreover, by comparing this case to the one of strongly interacting impurities (g II /g box = 2.0), we find once again significant differences only on the behavior of their two-body density. At t = 0, the shape of ρ (2) II (x I 1 , x I 2 ) corresponds to two coalesced impurities whose density peaks lay on its diagonal and are fragmented [cf. inset of Figure 2(a)]. In the course of the dynamics the impurities collide at the center where they remain in the course of the evolution [see inset of Figure 4(b3)]. Thereby, the strong repulsion between the impurities hinders a population at the exact diagonal of the two-body density. Notice their small spatial overlap with the bath hinting towards their suppressed dressing [90].

D. Dynamical response for strong impurity-medium repulsions
For even stronger impurity-medium repulsions (g BI /g > 2.0), the impurities once they collide around x = 0 they drift apart from each other and collide again [see Figure 4(c1) for g BI /g box = 5.0]. This behavior is repeated in the course of time with a damped collision amplitude. Thereby, the population on the diagonal of ρ [see the inset of Figure 4(c3)] with the two respective density fragments exhibiting each four faint maxima which persist until an evolution time of t/t box = 8. During the evolution the one-body density of the bath allocates at the edges of the box when the impurities collide and reoccupies the trap center when the latter drift apart, see also below for a more detailed discussion.
We attribute the impurities' density splitting for large g BI , even though the induced attractive interaction is higher in this case than for intermediate repulsive couplings, to the finite size of the considered box potential. In order to elucidate the underlying mechanism we show in Figure 5(a) the time-dependent spatial variance of the medium X 2 B (t) which serves as a measure for the instantaneous spatial extension of the medium's cloud [19,91]. At t = 0 and for sufficiently strong repulsive g BI the impurities reside both either at the right or at the left site of the double well while the bath particles avoid these pairs [cf. regime (I) in   Subsequently, after ramping down the barrier of the double well, the impurities collide around x = 0 enforcing the medium to depopulate the trap center, a process that results in the increase of X 2 B (t). To facilitate further our discussion, we provide specific profiles of the impurities' one-body density and their effective potential [24,39,92]. The effective potential is constructed from the superposition of the impurities' (post-quench) harmonic oscillator and the one-body density of the bath weighted by the impurity-medium coupling strength. It reads, B (x I , t). (7) A maximum of X 2 B (t) is reached, i.e. the spatial extend of the medium is largest, when the impurities allocate at the trap center [see Figures 5(b) and (c)]. Thereby, the impurities transfer energy to the medium leading to an increased interaction energy between the particles of the latter, i.e. Figure 5(a)]. After reaching a maximum of Ĥ BB the bath reoccupies the trap center as indicated by the reduction of X 2 B (t) and the impurities' density splits again, a behavior that is repeated in the course of time. As argued below, this dynamical response can be attributed to finite size effects stemming from the size of the medium's box potential.

E. Impact of the barrier height, atom number and box size on the impurity dynamics
In order to check the robustness of the observed dynamical response regimes of the impurities against the system parameters we have additionally varied the height of the double well h I , the size of the box potential L B and the number of bath particles N B . For small g BI and fixed g II /g box = 0.2, increasing the height of the double well from h I /Ẽ boxx −1 box = 2 to 7 leads to a crossover of the impurities, i.e. from a superfluid to a Mott-insulating phase. The former phase corresponds to a two-body density ρ where in all quarters prominent density peaks are present. In the latter case only the off-diagonal of ρ is populated [compare with the ground state configurations in regime (III) of Figure 2(a)]. Further increasing the impurity-medium coupling to intermediate repulsive values the impurities coalesce again independently of the barrier height of the double-well. Thereby, the impurities' ground state corresponds to the one of regime (I) depicted in Figure 2(a). Regarding the dynamical response of the above-described system, we did not find a qualitatively different behavior when h I is varied but rather a shifting of the identified regimes towards a larger value of g BI .
In Section IV C it was mentioned that the ground state phase diagram is altered with respect to variations of g BB and N B . As such, also the dynamical response of the system is affected upon tuning these parameters by means that the specific dynamical features are realized for smaller values of g BI . This is attributed to the fact that the impurities feature an enhanced magnitude of attractive induced interactions for either a decreased g BB (e.g. g BB /g box = 0.1 and especially in the interval |g BI | < g BB ) or an increased N B (for instance at N B = 30). Consequently, they localize around the trap center after their first collision for a smaller g BI than for g BB /g box = 0.5 or a larger N B such as N B = 30. This effect is arguably more prominent in the former scenario and during the dynamics manifests by the enhanced impurities localization at the trap center.
An opposite behavior is observed for varying only the size of the medium's box potential L B . In this case, an increasing L B broadens the medium such that the impurities are less affected by the presence of bath particles which, eventually, for intermediate repulsive g BI increases the relative distance between the impurities for a larger L B (not shown here). However, when the box size is of the order of the distance of the double well minima, the bath particles localize at the trap center (between the two double-well sites) for intermediate repulsive impurity-medium interaction strengths. Additionally, for strong repulsive g BI an increasing box size leads to the localization of the impurities around x = 0 and a dynamical response similar to the one observed in Figure 4(c1) is absent. This holds also when we simultaneously increase the number of bath particles and the size of the box potential while keeping the ratio N B /L B fixed [93]. B (x, t) Figure 6. Spatiotemporal evolution of the one-body densities for two weakly interacting impurities (gII /g ho = 0.2) and a weakly interacting bath (gBB/g ho = 0.5) trapped in a harmonic oscillator. The dynamics is induced by ramping down the barrier of the double well in which the impurities initially reside. In each column a different impuritymedium interaction strength is considered, which is from left to right gBI /g ho = −0.2, 0.2, 0.6.

VI. IMPURITIES DYNAMICS FOR THE HARMONICALLY TRAPPED BATH
Next, we examine the counterflow dynamics (induced by the same quench protocol) of the two impurities coupled to a harmonically confined bath. Since in this case the medium tends to localize at the trap center a phase separation between the two species is facilitated for intermediate to strong repulsive impuritymedium interaction strengths which eventually prohibits a subsequent dynamical mixing of the species. Indeed, we find that for values larger than g BI /g ho = 0.6 > g BB /g ho , corresponding to a vanishing spatial overlap at t = 0 [cf. Figure 3(b)], the initial phase separation between the impurities and the medium persists also in the course of the propagation. However, for g BI /g < 0.6 an intriguing response is observed.
A. Collision features in terms of the one-body density Figure 6 illustrates the time evolution of ρ (1) B (x, t) and ρ (1) I (x, t) for weak impurity-impurity couplings, i.e. g II /g ho = 0.2, and for varying impurity-medium interaction strength. As it can be exemplary inferred from Figures 6(a1) and (b1) for weak attractive or repulsive g BI the impurities perform a periodic motion within the harmonic trap and induce only small deformations to the bath density associated with soundwave emission of the latter [Figures 6(a2) and (b2)]. In the former case, a persisting breathing dynamics of the initially localized impurities takes place being somewhat similar to the one which has been previously discussed for a medium confined in a box potential [ Fig. 4(a1)]. In the latter scenario corresponding to Figure 6(b2), the originally spatially separated impurities collide around x = 0 and then split in a periodic fashion (see also the discussion below).   Moreover, we show the respective density evolution in the phase-separation regime, i.e. for g BI /g ho = 0.6 > g BB /g ho [Figures 6(c1) and (c2)]. It can be readily deduced that here the impurities are already initially phase separated with the bath and remain in this state also in the course of the evolution while performing small amplitude oscillations due to their collisions with the bath edges [94]. In particular, for g BI /g ho > 0.6 the impurities remain in the coalescence regime if the initial chosen values for g BI and g II coincide with regime (I) [ Figure 3(a)] and are spatially separated if the values for g BI and g II correspond to regime (III) [ Figure 3(a)]. Furthermore, for very strong attractive g BI the impurities and the medium localize together at the trap center where they remain throughout the time-evolution (not shown), see also Ref. [95] for a similar dynamics. We finally remark that as in the case of a box-confined medium increasing the number of bath atoms, e.g. to N B = 30, does not lead to significant alterations of the observed dynamical response.

B. Two-body density evolution for weak impurity-medium repulsions
Let us now focus on the dynamical properties of two impurities which are weakly repulsively coupled to the medium, e.g., via g BI /g ho = 0.2 and interacting among each other with g II /g ho = 0.2 [ Figure  6(b1)]. Specifically, we are interested in the dynamical evolution of the impurities' reduced two-body density ρ (2) II (x I 1 , x I 2 ) depicted in Figure 7. At t = 0, corresponding to the ground state in which the two impurities are confined in a double well, their two-body density exhibits two dominant density peaks across its off-diagonal and two suppressed peaks at its diagonal elements [see Figure 7(a)]. Considering the two diago- nal peaks as sufficiently small, this two-body configuration can be interpreted as the probability for the impurities to occupy opposite double-well sites. Following the quench, the impurities collide at the trap center and their cloud starts to contract and expand with a frequency corresponding to the periodic motion of the one-body density [ Figure 6 Figures 7(f)]. Notice that this dynamical response is inherently related to a process which is hidden on the one-body level [cf. Figure 6(b1)].

C. Single-particle dynamical excitation processes
To obtain insights into the underlying microscopic processes in the course of the impurity dynamics we project the many-body wave function onto basis functions consisting of the generalized Wannier functions φ I i (x I ) of the initially considered (pre-quenched) double-well potential [96,97]. In this way we can retrieve the probabilities for the impurities to occupy certain localized states of this basis and distinguish between the left and right double-well sites. The Wannier functions are constructed as a superposition from the six energetically lowest eigenfunctions of the onebody HamiltonianĤ (1),dw = − 2 2m I ∂ 2 (∂x I ) 2 + V dw I (x I ) and are provided in Figure 8(a) together with their associated eigenenergies i . In particular, the Wannier state corresponding to i = 1 (i = 2) is the energetically lowest one at the left (right) site. Analogously, i = 3, 5 (i = 4, 6) signify the first and second excited Wannier states at the left (right) site. Note that, even though for the analysis a basis of a doublewell potential is utilized, the impurities' dynamics still takes place within a harmonic oscillator. The respective two-body probabilities for the impurities to simultaneously occupy the i-th and j-th Wannier state are given by Here 1 B is the unit operator defined in the subspace of the bath and |φ I i φ I i | are the one-body projectors of the i-th Wannier state acting on a single impurity. The quality of the basis is tested by summing up all probabilities P (2) i,j (t) for each time instant and verifying that i,j P (2) i,j (t) > 0.97 holds until t/t ho = 70. The above-described two-body probabilities for g BI /g ho = 0.2 and g II /g ho = 0.2 are presented in Figure 8(b). All probabilities show an oscillatory behavior stemming from the periodic collision and expansion of the impurities [cf. Figure 6(b)]. Beyond this rapid motion, a decay and revival of P 2,2 (t)) Wannier state is maximized as demonstrated in Figure 6(b). This observation implies that a single-particle intraband excitation process takes place. Moreover, also energetically higher-lying Wannier states contribute to the ground state configuration of the impurities as well as to their dynamical response. For instance, the second and third (first and fourth) Wannier states contribute with P 2,3 (0) = P 1,4 (0) = 19 % to the ground state configuration. Therefore, the initial state is a superposition of different single-particle states. This is to be contrasted with the discussion below [see Figure 9] where a two-body basis is employed accounting in a more natural way for effects stemming from impurityimpurity interactions and the coupling to the bath.
For comparison the case of two weakly interacting impurities g II /g ho = 0.2 which are decoupled from the bath (g BI = 0) is showcased in Figure 8(c). It can be readily seen that the impurities do not perform a state transfer similar to the one depicted in Figure 8(b) but rather retain their delocalized configuration. In particular, at time-instants corresponding to an expansion of the impurities cloud the twobody probabilities associated with the energetically lowest Wannier states lying at opposite and the same sites significantly contribute to the impurities manybody wave function, thus confirming the former statement. Otherwise, the impurities response is characterized by excitations to energetically higher-lying states. Next, we inspect the case of strongly interacting impurities (g II /g ho = 2.0) which are weakly coupled to the bath (g BI /g = 0.2) [ Figure 8(d)]. The dynamics begins with initially separated impurities, viz. P 1,2 (t) obtains a maximum at t = 0, and continues with the collision of the impurities at the trap center where they both dominantly populate the same energetically lowest left or right Wannier state [cf. P 2,2 (t)]. Subsequently, the strong impurityimpurity repulsion enforces the impurities to occupy again opposite double-well sites. This scheme repeats itself during the evolution and, in particular, lasts until T /t ho = 200. Again a state transfer process occurs as it can be seen from the competition of P 1,2 (t) and P 1,1 (t). However, in this case the transfer is less transparent and not as dominant as for g II /g ho = 0.2, implying that an increasing g II /g ho results in the suppression of this process.
Therefore, the intraband excitation process observed for g BI /g ho = 0.2 and g II /g ho = 0.2 proves to be sensitive to the impurity-impurity interaction strength and, most importantly, requires a finite impurity-medium coupling. In this manner, we can conclude that this state transfer of the impurities is induced by the presence of the bath. Moreover, we have verified the absence of this mechanism for a species mean-field ansatz (D = 1 in Eq. (2)), i.e. when the entanglement is not taken into account. Thus, we can deduce that many-body effects and, in particular, the impurity-medium entanglement play a crucial role for the realization of such processes.

D. Effective two-body impurity mechanisms
To further understand the participating excitation processes we consider a projection of the many-body wave function onto a two-body basis set. In particular, we choose for this investigation the ground state and the first four energetically lowest excited states of an effective Hamiltonian. This effective Hamilto-nianĤ (2),eff describes two weakly interacting impurities (g II /g ho = 0.2) trapped in the effective potential defined in Eq. (7) with g BI /g ho = 0.2. In this manner, we take the backaction induced by the medium into account. In Figure 9(c)-(f) we present the impurities two-body density for the ground state as well as the first three excited states |Φ I i ofĤ (2),eff . As such, we associate the ground state (i = 0) with the twobody state (|LR + |RL )/ √ 2 where |L (|R ) represents a single-particle state corresponding to the left (right) site of the double well. Analogously, we refer to the first (i = 1) and second (i = 2) excited states as the configurations (|LL − |RR )/ √ 2 and (|LL + |RR )/ √ 2, respectively. The corresponding eigenenergies are shown in Figure 9(b). Note that the first (i = 1) and second (i = 2) excited eigenstates are II (x I 1 , x I 2 ) Figure 9. (a) Time evolution of the fidelity between the many-body wave function for gBI /g ho = 0.2 and gII /g ho = 0.2 and the excited states of the effective two-body Hamil-tonianĤ (2),eff (see also Eq. (7)). (b) Sum of the impuritymedium interaction energy Ĥ BI and the energy of the impurities Ĥ I in the many-body approach for the ground state as a function of the impurity-medium coupling (blue line) depicted together with the eigenenergies ofĤ (2),eff (grey horizontal dashed lines). The two-body density of the ground state and the first three excited eigenstates are provided in (c)-(f).
approximately degenerate. In order to support the validity of this two-body approach for two impurities coupled to a larger medium, we additionally provide the sum of the impurity energy and the interaction energy ( Ĥ I + Ĥ BI ) at t = 0 as predicted within the many-body approach, namely when the impurities are still trapped in a double-well potential. Since this energy matches at g BI /g ho = 0.2 the ground state energy of the effective approach [cf. intersection of Ĥ I + Ĥ BI with i = 0 in Figure 9(b)] we conclude that the effective potential adequately accounts for the presence of the medium at t = 0.
As a next step, we calculate the fidelity of the two-body eigenstates |Φ I i with the time-dependent species functions of the impurities being coupled to the medium. In this way, the probabilities of the contributing two-body configurations are revealed. Therefore, we estimate the absolute square of the projection of D j=1 |Ψ B j (t) ⊗|Φ I i on the many-body wave function |Ψ MB (t) defined in Eq. (2), which reads as The dynamics of the fidelity with respect to the ground state and the first four excited eigenstates of H (2),eff is provided in Figure 9(a). Analogously with the analysis regarding the Wannier states, we observe besides a fast periodic motion a slower decay and re-vival of the ground state |Φ I 0 associated with two separated impurities. This behavior is accompanied by a complementary increase of the second excited state |Φ I 2 = (|LL + |RR )/ √ 2 associated with the coalescence of the impurities. Since the first excited state, corresponding to the anti-symmetric configuration (|LL −|RR )/ √ 2, is strongly suppressed we conclude that the impurities undergo the two-body state transfer from (|LR +|RL )/ √ 2 to (|LL +|RR )/ √ 2. Moreover, we note that during the impurities' collision an appreciable amount of higher excited states need to be taken into account as indicated, for instance, by the non-negligible occupation of the third exited state |Φ I 3 .
Concluding we have explicated the microscopic mechanisms on both the one-and the two-body level taking place during the collision dynamics of interacting impurities coupled to a harmonically confined bath. For instance, the single-particle intraband excitation process appears to be sensitive with respect to the impurity-impurity interaction strength and requires a finite coupling to the bath. Moreover, an analysis with respect to a two-body basis deciphered the transitions among particular two-body configurations.

VII. CONCLUSIONS AND OUTLOOK
We have investigated the ground state and correlated dynamics of two interacting bosonic impurities confined in a double well and immersed in a bosonic medium. The latter either experiences a box potential or it is confined in a harmonic trap. We establish the phase diagram of the ground state for varying impurity-impurity and impurity-medium coupling strengths. Thereby, the emergent ground states have been characterized with the aid of the two-body densities and impurity-medium entanglement. An analysis of the impact of different trapping geometries on the formation of these phases has been performed. For instance, we explicate that the coalescence of the impurities at strong (repulsive or attractive) impuritymedium interaction strengths is preserved for different impurity-impurity repulsions when the bath is in a box. However, in the case of a harmonically trapped bath the impurities separate from each other for strong impurity-impurity repulsion residing in a Mott-type configuration. Moreover, in the latter scenario, we observe at strong impurity-medium attractions indications for the formation of a bipolaron.
Focusing on a specific interaction-dependent ground state configuration we trigger the dynamics by suddenly ramping down the potential barrier of the impurities' double well. Firstly, the dynamical response regimes of the impurities coupled to a box confined medium are unraveled with respect to their associated one-and two-body densities. In particular, for intermediate impurity-medium repulsions a localization of the impurities at the trap center after the original collision is realized. The impurities' two-body density features an elongated shape along the diagonal for weak impurity-impurity repulsion which suggests the presence of attractive induced interactions mediated by the bath. This induced localization of the impurities persists also when the coupling strength between the impurities is further increased. This observation together with the existence of a spatial overlap for finite impurity-medium interaction strengths support the formation of quasi-particles, i.e. two interacting polarons. We have attested the robustness of the above phenomena with respect to variations of the number of bath particles and the size of the box potential. However, for strong repulsions finite size effects of the medium's box potential come into play and govern the dynamical response of the system. Essentially, after the impurities' collision at the trap center they drift apart and then the medium reoccupies the center.
By considering a harmonically confined bath the impurities' response is qualitatively altered. Due to the spatial localization of the bath at the trap center the impurities and the medium undergo a phase separation already for intermediate impurity-medium repulsions as it was also observed on the ground state level. The response becomes especially intriguing for weak impurity-medium couplings where the impurities are able to perform a breathing motion within the bath. Specifically, for weak impurity-medium repulsions we observe a state transfer of the impurities starting with two spatially separated ones located at different double-well sites and evolving into a coalesced configuration, i.e. the impurities cluster. Interestingly, this state transfer process does not emerge for strongly interacting impurities and, most importantly requires a finite impurity-medium interaction strength, viz. it is induced by the coupling to the bath. Moreover, it is shown that this mechanism can be well understood in terms of a single-particle Wannier basis of the double-well. Additionally, we reveal the participating two-body states in this process using an analysis in terms of a two-body basis which consists of the low-lying excited states of a corresponding effective two-body Hamiltonian.
There are several possible extensions of our results. An immediate one will be to investigate the collision features of the impurities immersed in a spatially extended bosonic gas with the aim to unveil their possible damping mechanisms and appreciate the corresponding drag force. In another context, it would be worth including a spin degree of freedom for the impurities. Here, the dynamics of the emergent spin-spin correlations is of interest especially when the impurities localize around the trap center. Moreover, it would be intriguing to consider two impurities with different masses, e.g. a light and heavy one, coupled to a background. This way it would be feasible to investigate the influence of the mass on the emergent collisional aspects of the impurities and their induced interactions as well as trigger specific population transfer channels by considering a Rabi-coupling term.   Figure 10. Ground state relative distance between the impurities rII (in units ofx ho ) with respect to the impuritymedium interaction strength gBI for different impurityimpurity couplings (see legend). The relative difference ∆ rII between the impurity-impurity distance predicted within the many-body approach ( rII ) and the effective Hamiltonian of Eq. (7) with V dw I (x I ) = 0 ( r eff II ) is also provided (see main text). In the following, we examine the impurities relative distance r II which can serve as an indicator for the existence of their induced interactions mediated by the bath [38,39,98]. This quantity, which is accessible through in-situ spin-resolved single-shot measurements on the impurities state [99], reads The impurities relative distance is presented in Figure 10 for the ground state of the system (t = 0) for varying impurity-medium interaction strength and different impurity-impurity couplings. In the case of two weakly interacting impurities as well as for two strongly interacting ones we observe that with an increasing absolute value of g BI the relative distance between the impurities reduces. This behavior implies an induced attraction mediated by the coupling to the bath. Furthermore, it is evident that for larger impurity-impurity interactions r II is enhanced when compared to the one of weakly interacting impurities. This behavior is caused by the increased intraspecies impurities direct repulsion compensating their induced attraction.
To further justify the presence of induced interactions, we compare the resulting impurities relative distance as obtained using the many-body approach ( r II ) with the one ( r eff II ) predicted within the effective two-body HamiltonianH (2),eff . The latter was introduced in section VI for a harmonically confined medium. Notice that here this effective Hamiltonian describes the interplay between two interacting particles confined in an effective potential constructed by the one-body density of a box-trapped medium [see Eq. (7)]. As such, the entanglement between the impurities and the medium is neglected while effects stemming from the backaction to the bath are taken into account. By comparing r II between these two methods we can determine whether the decrease of r II for increasing |g BI | [see Figure 10] originates from an entanglement-assisted induced interaction or it is due to an alteration of the effective potential.
Inspecting ∆ r II = ( r eff II − r II )/ r eff II depicted in Figure 10 a large deviation among r II and r eff II becomes evident for increasing |g BI |. This confirms the presence of attractive induced interactions between the impurities. The fact that ∆ r II is finite can be traced back to the shape of the impurities two-body densities as obtained within the aforementioned approaches for large g BI [corresponding to regime I in Figure 2(a)]. While in the many-body scenario the impurities coalesce (namely only the diagonal of ρ (2) II (x I 1 , x I 2 ) is occupied), in the effective potential case both the diagonal and the off-diagonal elements of ρ (2) II (x I 1 , x I 2 ) are equally populated. This naturally leads to a larger relative distance. Concluding, the comparison with the effective model underlies the importance of considering correlations in the system and reveals the presence of attractive induced interactions between the impurities, see more details in Refs. [39,100].

Appendix B: Collisional properties of two heavy impurities
Throughout this work we have considered a system consisting of 87 Rb bath particles and 133 Cs impurities corresponding to a mass ratio of m I = 133/87m B . Below, we employ a mixture characterized by a mass ratio m I = 5m B in order to elaborate on the dynamical response of two heavy impurities [38] immersed in a bath confined in a box potential [101]. The time evolution of the corresponding one-and two-body densities of the impurities and the bath particles are depicted in Figures 11(a1)-(c1) and (a2)-(c2) for interaction parameters as the ones employed in the main text [ Figure 4]. In particular, we assume two weakly interacting impurities (g II /g box = 0.2) and vary the impurity-medium interaction strength g BI . Regarding the one-body density evolution we do observe a qualitatively similar behavior as compared to the case B (x, t) Figure 11. Time evolution of the one-body density of (a1)-(c1) two heavy impurities and (a2)-(c2) a box-confined bosonic bath. Each column represents the dynamics for a fixed impurity-medium interaction strength which is from left to right gBI /g box = −0.2, 1.5, 5.0. The impurityimpurity coupling remains constant being gII /g box = 0.2.
(a3)-(c3) Snapshots of the impurities two-body density. The insets show the two-body densities at the same time instants, but for two strongly interacting impurities, i.e. gII /g box = 2.0.
of lighter impurities. Only in the case of weak attractive g BI [ Figure 11(a1)] the heavy impurities perform a more pronounced breathing oscillation with a larger oscillation period. Inspecting a two-body density snapshot reveals that the impurities are spatially separated from each other and oscillate along the offdiagonal of ρ II (x I 1 , x I 2 ) [cf. Figure 11(a3)]. Therefore, they are not delocalized as their lighter counterparts [ Figure 4(a3)], and this behavior persists for two strongly interacting impurities [cf. inset of Figure  11(a3)].
Additionally, for intermediate impurity-medium repulsions, i.e. g BI /g box = 1.5, the one-body densities of the impurities and the medium [Figures 11(b1) and (b2)] as well as the impurities two-body density [ Figure 11(b3) and its inset] do not reveal a qualitatively different response with respect to the case of lighter impurities [cf. Figures 4(b1)-(b3)]. Indeed the impurities remain very close throughout the timeevolution exhibiting a more pronounced localization trend around the trap center as compared to lighter ones. The same holds also for the case of strong impurity-medium interactions [see Figures 11(c1)-(c3) for g BI /g box = 5.0] where the impurities feature multiple collisions with a dissipative amplitude. Only, the fragmentation in terms of the diagonal of the twobody density for g II /g box = 2.0 is more prominent in the case of heavy impurities than for lighter ones [cf. inset of Figure 11(c3)] and becomes visible even on the one-body density level (not shown here). The impurities are weakly interacting with gII /g box = 0.2 and are coupled to a box-confined medium with an interspecies coupling gBI /g box = 1.5. An increasing ramp time maintains an almost constant magnitude of entanglement until the impurities collision. (b) The corresponding onebody density evolution of the impurities for τ /t box = 30. As it can be seen the linear ramp delays the impurities first collision event when compared to Fig. 4 (b1) referring to the corresponding quench dynamics.
Appendix C: Impurity dynamics after a linear ramp of their double-well potential In the following we address the robustness of the impurities dynamical response when ramping down their barrier height h I in a time-dependent manner and not suddenly as in the main text. Specifically, we apply the following linear protocolh I (t) = h I − h I t/τ if 0 ≤ t ≤ τ whilsth I (t) = 0 as long as τ < t. The ramp time τ is defined as the one at which the barrier height vanishes, i.e. whenh I (τ ) = 0. For our purposes, we restrict our study to finite values of τ which deviate from the instantaneous quench but also do not refer to an adiabatic linear ramping.
To visualize the impact of the time-dependent protocol on the impurities collision process we present in Figure 12 the time-evolution of the von Neumann entanglement entropy [Eq. (5)] and the impurities one body density. We follow a linear ramp of the impurities barrier height characterized by a finite ramp time τ /t box = 30. Notice that the latter appreciably deviates from the quench scenario (τ /t box = 0). Moreover, we exemplarily invoke the system where the medium is confined in a box potential while the relevant interactions are g BB /g box = 0.5, g II /g box = 0.2 and g BI /g box = 1.5. Recall that in the main text it has been shown that in this regime the impurities localize at the trap center after their first collision event, see also Figure 4(b1). As it can be readily seen [ Fig-ure 12(b)], the impurities dynamical response in terms of ρ (1) I (x, t) remains qualitatively unchanged as compared to the quench [ Figure 4(b1)]. The most prominent difference regards the timescale at which the impurities collide and subsequently localize at the trap center. As expected the initial collision can be delayed for an increasing ramp time. Turning to the evolution of the respective von Neumann entropies for varying ramp times we observe an interesting behavior. At short timescales the impurities and the medium are highly entangled [ Figure 12(a)]. The magnitude of the entanglement is maintained in the course of the evolution until the impurities collide and then localize at the trap center where it suddenly decreases. This drop of the entanglement depends strongly on the ramp time and in particular it takes place faster for smaller τ since in this case the collision event is accelerated.
Concluding, we remark that the linear protocol affects the remaining response regimes, occurring for other interspecies interaction strengths, in a similar vein. Namely the main features as described in Section V do not substantially alter but rather the underlying timescales change. For instance, considering weak attractive or repulsive g BI referring to the impurities breathing motion [ Figure 4(a1)], we find a decreasing tendency of the breathing amplitude and frequency for larger τ (not shown). This is attributed to the fact that for increasing τ the collision of the impurities is less violent thus producing a less pronounced breathing. Analogous effects are observed for a harmonically trapped medium where, for instance, also in this case the amplitude and frequency of the underlying breathing motion depend strongly on the ramp time τ .