Many-body quantum dynamics and induced correlations of Bose polarons

We study the ground state properties and nonequilibrium dynamics of two spinor bosonic impurities immersed in a one-dimensional bosonic gas upon applying an interspecies interaction quench. For the ground state of two non-interacting impurities we reveal signatures of attractive induced interactions in both cases of attractive or repulsive interspecies interactions, while a weak impurity-impurity repulsion forces the impurities to stay apart. Turning to the quench dynamics we inspect the time-evolution of the contrast unveiling the existence, dynamical deformation and the orthogonality catastrophe of Bose polarons. We find that for an increasing postquench repulsion the impurities reside in a superposition of two distinct two-body configurations while at strong repulsions their corresponding two-body correlation patterns show a spatially delocalized behavior evincing the involvement of higher excited states. For attractive interspecies couplings, the impurities exhibit a tendency to localize at the origin and remarkably for strong attractions they experience a mutual attraction on the two-body level that is imprinted as a density hump on the bosonic bath.

The focus of the majority of the above-mentioned theoretical studies have been the stationary properties of the emergent quasiparticle states for single impurities in homogeneous systems. However, the nonequilibrium dynamics of impurities is far less explored and is expected to be dominated by correlation effects which build up in the course of the evolution [34-36, 39, 42-45]. Existing examples include the observation of self-trapping phenomena [46,47], formation of dark-bright solitons [6,42], impurity transport in optical lattices [48][49][50][51], orthogonality catastrophe events [35,52], injection of a moving impurity into a gas of Tonks-Girardeau bosons [53][54][55][56][57][58][59][60] and the relaxation dynamics of impurities [45,61,62]. Besides these investigations, which have enabled a basic description of the quasiparticle states in different interaction regimes, a number of important questions remain open and a full theoretical understanding of the dynamics specifically of Bose polarons is still far from complete.
A system of particular interest consists of two impurity atoms immersed in a Bose-Einstein condensate (BEC), where the underlying interactions between the impurities come into play. In such a system impurity-impurity correlations [10,63,64] can be induced by the BEC, even in the case where no direct interaction between the impurities is present. However, the competition between direct and induced interactions can also be expected to lead to interesting effects. It is therefore natural to investigate the dynamical response of the impurities with varying interspecies interactions (attractive or repulsive) and to identify in which regimes robustly propagating Bose polaron states exist [25,39,43]. In addition it is interesting to study the existence of bound states between the impurities [1,10], the effect of strong correlation between the impurities on the orthogonality catastrophe [35,52], phase separation between the two atomic species [65][66][67] and energy exchange processes [68,69]. Comparing the effects in systems with single and multiple impurities is an interesting task, as well as their theoretical interpretation in terms of the spin polarization (alias the contrast) which has not yet been analyzed in the case of two impurities and involves more energy channels compared to the case of a single impurity. For these reasons, we study in this work an interspecies interaction quench for two bosonic impurities overlapping with a harmonically trapped BEC. To address the correlated quantum dynamics of the bosonic multicomponent system we use the Multi-Layer Multi-Configuration Time-Dependent Hartree method for atomic mixtures (ML-MCTDHX) [70][71][72], which is a non-perturbative variational method that enables us to comprehensively capture interparticle correlations.
In this work we start by studying the ground state of two non-interacting impurities in a bosonic gas and show that for an increasing attraction or repulsion they feature attractive induced interactions, a result that persists also for small bath sizes and heavy impurities [8]. However, two weakly repulsively interacting impurities can experience a net repulsion for repulsive interspecies interactions.
When quenching the multicomponent system, we monitor the time-evolution of the contrast and its spectrum [19,25] for varying postquench interactions. We show that the polaron excitation spectrum depends strongly on the postquench interspecies interaction strength and the number of impurities while it is almost insensitive to the direct impurity-impurity interaction for the weak couplings considered herein. Additionally, a breathing motion of the impurities can be excited [73,74] for weak postquench interspecies repulsions, while for stronger ones a splitting of their single-particle density occurs. In this latter case a strong attenuation of the impurities motion results in the accumulation of their density at the edges of the bosonic gas and they mainly reside in a superposition of two distinct two-body configurations: the impurities either bunch on the same or on separate sides of the BEC, while the bath exhibits an overall breathing motion. For attractive interspecies couplings, the impurities exhibit a breathing motion characterized by a beating pattern. The latter stems from the values of the impuritie's center-of-mass and relative coordinate breathing modes, whose frequency difference originates from the presence of attractive induced interactions. Additionally, the impurities possess a tendency to localize at the trap center, a behavior that becomes more pronounced for stronger attractions [75]. Strikingly, for strong attractive interspecies interactions we show that during the dynamics the impurities experience a mutual attraction on the two-body level and the density of the bosonic bath develops a small amplitude hump at the trap center. We find that a similar dynamical response also takes place for two weakly repulsively interacting impurities but the involved time-scales are different. To interpret the observed dynamics of the impurities we invoke an effective potential picture that applies for weak couplings [35,36,73,75].
Our work is structured as follows. Section II presents our setup and introduces the correlation measures that are used to monitor the dynamics. In Sec. III we address the ground state properties of the impurities for a wide range of interspecies interaction strengths. The emergent nonequilibrium dynamics triggered by an interspecies interaction quench is analyzed in detail in Sec. IV. In particular, we present the time-evolution of the contrast and the system's spectrum [Sec. IV A-IV C] and study the full dynamics of the single-particle and two-body reduced density matrices for repulsive [Sec. IV D] and attractive [Sec. IV E] postquench interactions. We summarize and discuss future perspectives in Section V. Finally, Appendix A details our numerical simulation method and demonstrates the convergence properties.

II. THEORETICAL FRAMEWORK
A. Hamiltonian and quench protocol We consider a highly particle number imbalanced Bose-Bose mixture composed of N I = 2 bosonic impurities (I) possessing an additional pseudospin-1/2 degree of freedom [76], which are immersed in a bosonic gas of N B = 100 structureless bosons (B). Moreover, the mixture is assumed to be mass-balanced, namely m B = m I ≡ m and each species is confined in the same one-dimensional external harmonic oscillator potential of frequency ω B = ω I = ω. Such a system can be experimentally realized by considering e.g. a 87 Rb BEC where the majority species resides in the hyperfine state |F = 2, m F = 1 and the pseudospin degree of freedom of the impurities refers for instance to the internal states |↑ ≡ |F = 1, m F = 1 and |↓ ≡ |F = 1, m F = −1 [77,78]. Alternatively, it can be realized to a good approximation by a mixture of isotopes of 87 Rb for the bosonic gas and two hyperfine states of 85 Rb for the impurities. The underlying MB Hamiltonian of this system readŝ The non-interacting Hamiltonian of the bosonic gas isĤ 0 ↓} being the indices of the spin components. HereΨ σ (x) refers to the bosonic field-operator of either the bosonic gas (σ = B) or the impurity (σ = a = {↑, ↓}) atoms. Furthermore, we operate in the ultracold regime where s-wave scattering is the dominant interaction process.
Therefore both the intra-and the intercomponent interactions can be adequately modeled by contact ones. The contact intraspecies interaction of the BEC component is modeled byĤ and between the impurities viaĤ int where either a = a =↑, ↓ or a =↑, a =↓. Note also that we assume g ↑↑ = g ↓↓ = g ↑↓ ≡ g II . Most importantly, we consider that only the pseudospin-↑ component of the impurities interacts with the bosonic gas while the pseudospin-↓ is noninteracting. The resulting intercomponent interaction isĤ int In all of the above-mentioned cases, the effective onedimensional coupling strength [79] is given by and µ = m 2 is the reduced mass. The transversal length scale is a ⊥ = /µω ⊥ with ω ⊥ being the transversal confinement frequency and a s σσ denotes the threedimensional s-wave scattering length within (σ = σ ) or between (σ = σ ) the components. In a corresponding experiment, g σσ can be tuned either via a s σσ with the aid of Feshbach resonances [15,16] or by adjusting ω ⊥ using confinement-induced resonances [79]. In the following, the MB Hamiltonian of Eq. (1) is rescaled with respect to ω. As a consequence, length, time, and interaction strengths are given in units of mω , ω −1 and 3 ω m respectively. To study the quench dynamics, the above-described multicomponent system is initially prepared in its ground state configuration for fixed g BB = 0.5 and g BI = 0 and either g II = 0 or g II = 0.2. In this way, the case of two non-interacting and that of weakly interacting impurities are investigated. This initial (ground) state emulates a system prepared in the |1, −1 = |↓ 1 ⊗ |↓ 2 configuration for the spin degree of freedom i.e. where the impurity-BEC interaction is zero. Note that the spinor part of the wavefunction is expressed in the basis of the total spin i.e. |S, S z [80]. Accordingly, the spatial part |Ψ 0 BI of the ground state of the system obeys the following eigenvalue equation Ĥ −Ĥ BI |Ψ 0 BI |1, −1 = E 0 |Ψ 0 BI |1, −1 , with E 0 being the corresponding eigenenergy and H BI |Ψ 0 BI |1, −1 = 0. To trigger the dynamics we carry out an interspecies interaction quench from g BI = 0 to a finite positive or negative value of g BI at t = 0 and monitor the subsequent time-evolution. In a corresponding experiment, this quench protocol can be implemented by using a radiofrequency π/2 pulse with an exposure time much smaller than ω −1 [19]. The pulse acts upon the spin degree of freedom of the impurity, which maps the pseudospin-↓ impurities to the superposition state |ψ S i ≡ with i = 1, 2 [18]. The corresponding MB wavefunction of the system, |Ψ(t) = e −iĤt/ |Ψ 0 BI (|ψ S 1 ⊗ |ψ S 2 ) , is then given by The setup and processes addressed in our work can be experimentally realized utilizing radiofrequency spec-troscopy [9,18,22,23,43] and Ramsey interferometry [18].

B. Many-body wavefunction ansatz
To calculate the stationary properties and to track the MB nonequilibrium quantum dynamics of the multicomponent bosonic system discussed above we employ the ML-MCTDHX method [70][71][72]. This is an abinitio variational method for solving the time-dependent MB Schrödinger equation of atomic mixtures and it is based on the expansion of the total MB wavefunction with respect to a time-dependent and variationally optimized basis tailored to capture both the intra-and the interspecies correlations of a multicomponent system [35,65,81,82].
To include the interspecies correlations, the MB wavefunction (|Ψ(t) ) is first expanded in terms of D distinct species functions, |Ψ σ i (t) , for each component σ = B, I, and then expressed according to a truncated Schmidt decomposition [83] of rank D, namely Here the time-dependent expansion coefficients λ k (t) are the Schmidt weights and will be referred to in the following as the natural populations of the k-th species function. Evidently, the system is entangled [84] or interspecies correlated when at least two different λ k (t) possess a nonzero value. If this is not the case, i.e. for λ 1 (t) = 1, λ k>1 (t) = 0, the wavefunction is a direct product of two states. Therefore, in order to account for intraspecies correlations, each of the above-mentioned species functions is expressed as a linear superposition of time-dependent number-states, | n(t) σ , with time-dependent coefficients A σ i; n (t) as Each number state | n(t) σ is a permanent building upon d σ time-dependent variationally optimized single-particle functions (SPFs) |φ σ l (t) , l = 1, 2, . . . , d σ with occupation numbers n = (n 1 , . . . , n d σ ). Consecutively, the SPFs are expanded on a time-independent primitive basis. The latter refers to an M dimensional discrete variable representation (DVR) for the majority species and it is denoted by {|k }. For the impurities this corresponds to the tensor product {|k, s } of the DVR basis for the spatial degrees of freedom and the two-dimensional pseudospin-1/2 basis {|↑ , |↓ }. Accordingly, each SPF of the impurities is a spinor wavefunction of the form with B I jk↑ (t) [B I jk↓ (t)] being the time-dependent expansion coefficients of the pseudospin-↑ [↓] (see also Refs. [35,82] for a more detailed discussion).
The time-evolution of the (N B + N I )-body wavefunction |Ψ(t) governed by the Hamiltonian of Eq. (1) is obtained via solving the so-called ML-MCTDHX equations of motion [70]. The latter are determined by utilizing e.g. the Dirac-Frenkel [85,86] variational principle for the generalized ansatz introduced in Eqs. (3), (4) and (5). This procedure results in a set of D 2 linear differential equations of motion for the λ k (t) coefficients which are coupled to D( nonlinear integrodifferential equations for the species functions and d B + d I nonlinear integrodifferential equations for the SPFs.
A main aspect of the ansatz outlined above is the expansion of the system's MB wavefunction with respect to a time-dependent and variationally optimized basis. The latter allows to efficiently take into account the intraand intercomponent correlations of the system using a computationally feasible basis size. In the present case the Bose gas consists of a large number of weakly interacting particles and therefore its intracomponent correlations are suppressed. As a consequence they can be adequately captured by employing a small number of orbitals, d B < 4. Additionally, the number of impurities, N I < 3, is small giving rise to a small number of integrodifferential equations allowing us to employ many orbitals, d I , and thus account for strong impurity-impurity and impurity-BEC correlations. Therefore, the number of the resulting equations of motion that need to be solved is numerically tractable. Since our method is variational, its validity is determined upon examining its convergence. For details on the precision of our simulations see Appendix A.

C. Correlation measures
To study the quench-induced dynamics of each species at the single-particle level we calculate the one-body reduced density matrix for each species [89,90] Here,Ψ σ (x) is the σ-species bosonic field operator acting at position x and satisfying the standard bosonic commutation relations [91]. For simplicity, we will use in the following the one-body densities for each species i.e. ρ (1) σ (x; t) ≡ ρ (1) σ (x, x = x; t), which is a quantity that is experimentally accessible via averaging over a sample of single-shot images [65,92,93]. We remark that the eigenfunctions and eigenvalues of ρ (1) σ (x, x ; t) are termed natural orbitals ϕ σ i (x; t) and natural populations n σ i (t) [65,70] respectively. In this sense, each bosonic subsystem is called intraspecies correlated if more than a single natural population possess a non-zero contribution.
To unveil the role of impurity-impurity correlations following the interspecies interaction quench we calculate the time-evolution of the corresponding diagonal of the two-body reduced density matrix where a, a =↑, ↓. The two-body reduced density matrix refers to the probability of finding simultaneously one pseudospin-a boson at x 1 and a pseudospin-a boson at x 2 [65,66]. Moreover, it provides insights into the spatially resolved dynamics of the two impurities with respect to one another. Indeed, the impurities are dressed by the excitations of the bosonic gas forming quasiparticles which in turn can move independently or interact, and possibly form a bound state [8,10,42,94].
To capture the emerging effective interactions between the two bosonic impurities we monitor their relative distance [9,42] given by Here,N a with a =↑, ↓ is the number operator that measures the number of bosons in the spin-a state. Experimentally, r aa (t) can be probed via in-situ spin-resolved single-shot measurements on the spin-a state [93]. More precisely, each image gives an estimate of r aa (t) between the bosonic impurities if their position uncertainty is assured to be adequately small [93]. Subsequently, r aa (t) is obtained by averaging over several such images.

III. INDUCED INTERACTIONS IN THE GROUND STATE OF TWO BOSONIC IMPURITIES
Before investigating the nonequilibrium dynamics of the two bosonic impurities immersed in a BEC it is instructive to first analyze the ground state of two impurities interacting with the bosonic medium for varying interspecies interactions g B↑ ranging from attractive to repulsive. Note that such a configuration corresponds in our case to two impurities residing in the pseudospin-↑ state since only this state is interacting with the bath [see also Eq. (1)]. The aim of this study is to reveal the presence of induced impurity-impurity interactions mediated by the bath. As discussed in Section II A, the mass-balanced multicomponent bosonic system consists of two impurities N I = 2 immersed in a MB bath of N B = 100 atoms with g BB = 0.5 and it is externally confined in a harmonic oscillator potential of frequency ω = 1. Later on, also the mass-imbalanced and the fewbody (N B = 10) scenaria will be investigated. Below we consider either two non-interacting (g II = 0) or two weakly interacting impurities (g II = 0.2). To obtain the interacting ground state of the system as described by the Hamiltonian of Eq. (1) we employ either imaginary time propagation or improved relaxation [70,71] within ML-MCTDHX.
The relative distance [Eq. (8)] between the two impurities as well as their two-body reduced density matrix [Eq. (7)] for different values of g B↑ are shown in Fig. 1. Focusing on the case of two non-interacting impurities, g II = 0, we see that for larger attractions the relative distance between the impurities decreases (see Fig, 1 (a)) and converges towards a constant value i.e. r ↑↑ ≈ 0.1 for g B↑ < −2. The decrease in r ↑↑ for −2 < g B↑ < 0 implies that the impurities effectively experience an attraction with respect to one another. This attraction is a manifestation of the attractive induced interactions mediated by the bosonic gas since g II = 0 [8]. The impurities reside together in the vicinity of the trap center since ρ where r ↑↑ become approximately constant, the impurities come very close with respect to one another. Here, the corresponding ρ (2) ↑↑ (x 1 , x 2 ) shrinks along its anti-diagonal and its diagonal becomes elongated [see Fig. 1 (b 1 )], which is indicative of a bound state having formed between the impurities known as a bipolaron state [8,10,94].
Turning to weak interspecies repulsions 0 < g B↑ < 0.5 we find that r ↑↑ slightly increases [see Fig. 1 (a)] while the two impurities reside close to the trap center [see Fig. 1 It is important to mention that this increase in r ↑↑ does not directly imply that the impurities experience a weak repulsion mediated by the bosonic bath. Indeed, by neglecting all correlations between the impurities, i.e. by substituting ρ ↑ (x 2 )/2 into r ↑↑ we find the same tendency of r ↑↑ with even slightly larger values (see also the discussion below). Since in the limit of the non-correlated case there are no induced interactions, the fact that r ↑↑ is smaller when correlations are taken into account means that the impurities still feel an effective attractive force. Note that for the other interaction regimes presented herein such an unexpected behavior of r ↑↑ does not occur as it can also be deduced by the corresponding two-body spatial configurations building upon ρ (2) ↑↑ (x 1 , x 2 ) (see below). Furthermore, it can be seen that at g B↑ = g BB = 0.5, where the miscibility/immiscibility transition between the impurity and the BEC takes place [65,67], the behavior of r ↑↑ is suddenly altered. Indeed for g B↑ ≥ 0.5, r ↑↑ shows a decreasing tendency which indicates the presence of attractive induced interactions between the impurities. In particular, for 0.5 ≤ g B↑ < 1.1, r ↑↑ reduces and the impurities tend to bunch together at the same location. This can be confirmed by the fact that ρ (2) ↑↑ (x 1 , x 2 ) shows a populated elongated diagonal as depicted in Fig. 1 (b 4 ) for g B↑ = 0.5. Moreover for stronger repulsions g B↑ > 1.1, r ↑↑ remains almost constant. Especially so for g B↑ > 1.5, where the two impurities residing either on the left or the right edge of the Thomas-Fermi profile of the BEC. The latter can be evidenced in Fig. 1 (b 5 ) by the two strongly populated spots appearing at x 1 ≈ x 2 ≈ ±R T F with R T F denoting the Thomas-Fermi radius.
In view of the results of Ref. [35] it is tempting to interpret our above findings in terms of an effective potential, V ef f (x; g BI ). A valid candidate for such a potential can be constructed as where ρ (1) B (x; g BI = 0) refers to the equilibrium density of the BEC for g BI = 0. Equation (9) implies that ρ (1) B (x; g BI = 0) acts on the impurities just as an additional repulsive (g BI > 0) or attractive (g BI < 0) potential on top of the externally imposed parabolic trap. It is noteworthy that the simplification of the impurity problem provided by Eq. (9) neglects several phenomena that might be important for the description of the ground state of the impurity system. First, the renormalization of the impurity's mass, m I → m ef f I by the coupling with its environment is neglected and, most importantly, the possible emergence of induced interactions is not contained in Eq. (9), due to the absence of twobody terms. The latter are extremely important for the description of ρ To provide an estimate of the quantitative error obtained by the approximation of Eq. (9) we include in Fig. 1 (a), also the results for r ↑↑ within the effective potential picture. It is evident that when using V ef f (x; g BI ), r ↑↑ is always larger than the corresponding full MB result for g BI = 0 . This effect is particularly pronounced for g BI > 0.5 where r ↑↑ within Eq. (9) exhibits an increasing tendency instead of a decreasing one with g BI . Such an effect can be attributed to the vanishing off-diagonal elements of ρ (2) ↑↑ (x 1 , x 2 ) which cannot be captured within V ef f (x; g BI ), as in the latter case ρ ↑ (x 2 )/2. Indeed, the large impurity-impurity interactions within this regime render the effective potential incapable of describing the ground state of the bath impurity system within this interaction regime. Similarly, for g BI < −2, r ↑↑ using the effective potential is significantly larger than the corresponding MB result, which can be attributed to the prominent role of induced interactions in the formation of the bipolaron state [10]. 2) impurities as well as few-and many bath particles are shown (see legend) for a mass-balanced system mI = mB. r ↑↑ from the effective potential picture of Eq. (9) for two non-interacting bosonic impurities is also illustrated (see legend) with respect to g B↑ . Inset illustrates r ↑↑ of two non-interacting impurities in the case of a mass-balanced (mI = mB) and a mass-imbalanced (mI ≈ 1.53mB) system with respect to g B↑ . The corresponding two-body reduced matrix of the ground state of the two pseudospin-↑ (b1)-(b5) non-interacting and (c1)-(c5) interacting (gII = 0.2) impurities for different interspecies interactions (see legends). In (b1)-(b5) and (c1)-(c5) the mixture consists of NB = 100 bosons and NI = 2 bosonic impurities. Also, in (b4), (b5), (c4) and (c5) the dashed magenta lines indicate the location of the Thomas-Fermi radius of the bosonic gas. In all cases gBB = 0.5 and the system is trapped in a harmonic oscillator potential with ω = 1.
Considering a smaller bath consisting of N B = 10 atoms does not significantly alter the ground state properties of the two non-interacting bosonic impurities. Here, r ↑↑ [ Fig. 1 (a)] exhibits a similar behavior as for N B = 100 atoms, with the most notable difference occurring in the region of g B↑ ≈ g BB where a smoother decrease occurs when compared to the N B = 100 case. The value for which the distance becomes constant is also shifted to larger values when N B = 10. These differences can be qualitatively understood within a corresponding effective potential picture which we will discuss in Section IV D 1, see Eq. (15) and the remark [95].
A similar to the above-described overall phenomenology of the two non-interacting bosonic impurities for a varying g B↑ is also observed for the case of heavier impurities as can be seen in the inset of Fig. 1 (a). Here we consider a 87 Rb bosonic gas and two 133 Cs impurities prepared e.g. in the hyperfine states |F = 1, m F = 0 and |F = 3, m F = 2 respectively and being both confined in the same external harmonic oscillator [96,97]. Compared to the mass-balanced scenario the behavior of r ↑↑ around g B↑ ≈ g BB becomes somewhat smoother and the maximum value is also slightly shifted to larger interaction strengths. Another conclusion that can be drawn, is that heavier impurities prefer to remain closer to each other compared to the lighter ones, since r ↑↑ has smaller values in the former than in the latter case. As a consequence we can infer that heavy impurities experience stronger attractive induced interactions than light ones. These differences can also be explained in terms of the effective potential picture which will be introduced in section IV D 1, see also remark [98].
When a weak intraspecies repulsion among the impurities is introduced, g II = 0.2, see Fig. 1 (a), the ground state properties remain the same for attractive g B↑ but change fundamentally in the repulsive regime. Indeed r ↑↑ decreases for an increasing interspecies attraction, signifying an induced attraction between the impurities despite their repulsive mutual interaction, until it becomes constant for g B↑ < −2. More specifically, for −2 < g B↑ < 0 the impurities are likely to remain close to the trap center [see Fig. 1 is predominantly populated. Furthermore, for g B↑ < −2 the impurities bunch together at a fixed distance [ Fig. 1 (a)] and the two-body reduced density matrix becomes elongated along its diagonal [see Fig. 1(c 1 )], suggesting the formation of a bound state similar to the g II = 0 case. However, for g B↑ > 0, r ↑↑ exhibits an overall increasing tendency, which indicates that the two impurities are located mainly symmet- B (x; 0) with respect to g B↑ for NI = 2 and gII = 0. (c) Ground state one-body density of two non-interacting impurities as a function of g B↑ . In all cases the bath consists of NB = 100 bosons with gBB = 0.5.
rically around the trap center. This latter behavior can be directly deduced by the relatively wide distribution of the anti-diagonal of their two-body reduced density matrix [see Figs. 1 (c 3 ) and (c 4 ) for 0 < g B↑ < 1]. Moreover, and in sharp contrast to the g II = 0 case, for g B↑ > 1 the impurities acquire a large fixed distance and in particular can be found to reside one at the left and the other at the right edge of the BEC. This configuration of the impurities can be seen from the fact that solely off-diagonal elements of ρ Finally, it is worth mentioning that for two weakly repulsive impurities the induced effective attraction can never overcome their direct s-wave interaction for g B↑ > 0.
To further support the existence of attractive induced interactions between the two impurities we study the ground state energy of the system for varying g B↑ . In particular, we calculate the expected position of the polaronic resonances [9] namely ∆ N I . As it can be seen, for both, N I = 1 and N I = 2, the resonance position ∆ N I + (g B↑ ) increases for a larger g B↑ and it takes negative and positive values for attractive and repulsive interactions, respectively. Moreover, in the N I = 2 scenario ∆ N I + (g B↑ ) is found to be negatively shifted when compared to the corresponding N I = 1 case for g II = 0. This behavior indicates the presence of attractive induced interactions for both attractive and repulsive Bose polarons [8,10,63]. Focusing on g II = 0.2 and g B↑ < 0 a small decrease of ∆ N I + (g B↑ ) occurs when compared to the g II = 0 case showing that attractive induced interactions become more pronounced when direct s-wave impurity-impurity repulsions are involved. However, for repulsive polarons i.e. g B↑ > 0 the presence of s-wave impurity-impurity interactions counteracts the effect of attractive induced interactions and accordingly ∆ N I + (g B↑ ) is almost the same for N I = 2, g II = 0.2 and N I = 1, see the inset of Fig. 2 (a).
The underlying mechanism behind the abovementioned impurity-impurity induced interactions can be qualitatively understood as follows. For attractive g B↑ the presence of impurities gives rise to a small density enhancement of the BEC in the vicinity of their spatial position. This effect is captured by the deformation of the BEC density quantified by δρ B (x; 0) and shown in Fig. 2 . This density enhancement of the BEC forces the impurities to approach each other leading to the emergence of attractive impurity-impurity induced interactions. Similarly for g B↑ < 0 the impurities tend to reside in regions of lower bath density causing a density depletion of the BEC characterized by δρ The above-described density depletion of the bath gives rise to the attractive induced interactions analogously to g B↑ < 0. It is also worth commenting that for g B↑ > 0.5 ρ

IV. QUENCH INDUCED DYNAMICS
Next, we study the interspecies interaction quenched dynamics for the mass-balanced multicomponent system which is initially prepared in its ground state and characterized by g BB = 0.5 and g B↑ = 0. In this case the Thomas-Fermi radius of the BEC is R T F ≈ 4.2 and the impurities are in a superposition of their spin components described by Eq. (2). We mainly analyze the case of two non-interacting (g II = 0) impurities and briefly discuss the scenario of two weakly interacting impurity atoms in order to expose the effect of their mutual interaction in the dynamics.
To induce the nonequilibrium dynamics we perform at t = 0 a sudden change from g B↑ = 0 to either attractive [Sec. IV E] or repulsive [Sec. IV D] finite val-ues of g B↑ . To examine the emergent dynamics we first discuss the time-evolution of the spin polarization (alias contrast) and its spectrum. Consequently we discuss the dynamical response of the impurities in terms of their single-particle densities and the corresponding two-body reduced density matrix. An effective potential picture for the impurities is constructed in order to provide an intuitive understanding of the quench dynamics.

A. Interpretation of the contrast of two impurities
To examine the quench-induced dynamics of the two spinor bosonic impurities we first determine the time-evolution of the total spin polarization (contrast) which enables us to infer the dressing of the impurities during the dynamics [18].
, with σ k ab denoting the Pauli matrices. The contrast for a single impurity has been extensively studied [25,39,43,99] and it is related to the so-called Ramsey response [18] and therefore the structure factor. The time-dependent overlap between the interacting and the noninteracting states is given by where |Ψ 0 BI is the spatial part of the MB ground state wavefunction of a single impurity with energyẼ 0 when g BI = 0.Ĥ =PĤP withP being the projector operator to the spin-↑ configuration, andĤ denotes the postquench Hamiltonian [Eq. (1)]. Note also that the contrast is chosen here to take values in the interval [0, 1]. From Eq. (10) zero contrast implies that the overlap between the interacting and the non-interacting states vanishes signifying an orthogonality catastrophe phenomenon [52,99]. On the other hand, if | Ŝ (t) | 2 = 1 then the non-interacting and the interacting states coincide and no quasiparticle is formed. Therefore only in the case that 0 < | Ŝ (t) | 2 < 1 we can infer the dressing of the impurity and the formation of a quasiparticle.
When increasing the number of impurity atoms to N I > 1, | Ŝ (t) | 2 is more complex since additional spin states contribute to the MB wavefunction (see Eq. (2)). To understand the interpretation of | Ŝ (t) | 2 during the dynamics we therefore first discuss it for the case of two impurities. The contrast of two pseudospin-1/2 bosonic impurities reads where the spatial overlap between two different spin con-figurations namely |S, S z and |S , S z is defined as [80] referring to the spatial wavefunction corresponding to the spin configuration |S, S z and |Ψ 0 BI being the spatial part of the initial MB state for two impurities.
In particular in our case we consider two Recall that a quasiparticle is a free particle that is dressed by the excitations of a bosonic bath via their mutual interactions. As a consequence, Ψ 0 BI ( x B , x I ) refers to the wavefunction where no polaron quasiparticle exists since it is the ground state wavefunction of the system with g B↑ = 0. Moreover, Ψ 1,0 ( x B ; x I ) and Ψ 1,1 ( x B ; x I ) denote the wavefunctions where a single and two impurities respectively interact with the bosonic gas and therefore describe the formation of a single and two polarons, respectively. Accordingly, A(|1, 0 ; |1, −1 ) provides the overlap between the state of a single and no impurities interacting with the bath, while A(|1, 0 ; |1, 1 ) is the overlap between a single and two impurities interacting with the bath.
As a result, | Ŝ (t) | 2 = 1 means that A(|1, 0 ; |1, −1 ) = A(|1, 0 , |1, 1 ) = e iϕ where ϕ is a phase factor. The fact that |A(|1, 0 ; |1, −1 )| = 1 implies that the spatial state of a single impurity interacting with the bath is the same as the noninteracting one, except for a possible phase factor, and therefore a quasiparticle is not formed. Moreover since also |A(|1, 0 , |1, 1 )| = 1 it holds that the state of a single pseudospin-↑ interacting impurity coincides with the state of two pseudospin-↑ impurities interacting with the bath and as a consequence with a bare particle due to |A(|1, 0 ; |1, −1 )| = 1. Thus, In the former case we can deduce the occurrence of an orthogonality catastrophe phenomenon as in the single impurity case while the latter scenario is given by the destructive interference of the A(|1, 0 ; |1, −1 ) and A(|1, 0 , |1, 1 ) terms. However, for 0 < | Ŝ (t) | 2 < 1 the corresponding overlaps acquire finite values and a quasiparticle can be formed.
Notice also that in the special case of g ↑↓ = 0 and g ↓↓ = 0 (but g ↑↑ arbitrary) it can be shown that . The latter is exactly the contrast or the structure factor of a single impurity [Eq. (10)]. Indeed |Ψ 0 BI = |Ψ 0 BI ⊗ |ψ 0 I for g ↓↓ = 0 holds where |ψ 0 I is the singleparticle ground state of the impurity while |Ψ 0 BI and |Ψ 0 BI refer to the spatial part of the MB ground state wavefunction of a single (energyẼ 0 ) and two impurities (energy E 0 ) respectively. AdditionallyĤ is the postquench Hamiltonian given by Eq. (1). Consequently, the contrast in this special case acquires the simplified form Evidently, here | Ŝ (t) | depends explicitly on the structure factor S 1 (t) of a single impurity allowing for a direct interpretation of the dynamical dressing of the two impurities with respect to the single impurity case discussed in Ref. [35]. In the following, g ↑↑ = g ↓↓ = g ↑↓ ≡ g II and as a consequence g ↑↓ = 0, g ↓↓ = 0 is encountered for g II = 0 while the general case of Eq. (12) applies for the case of g II = 0.2 analyzed below.

B. Evolution of the contrast
The dynamics of the two particle contrast | Ŝ (t) | is presented in Figs. 3 (a)-(c) for both attractive and repulsive postquench interspecies interactions g B↑ . In particular, | Ŝ (t) | is shown for either two non-interacting . In all cases, six different dynamical regions with respect to g B↑ can be identified marked as R I , R II , R III , R IV , R II and R III . Focusing on the system with N B = 100 and g II = 0 these re- Fig. 3 (a)]. Specifically, within the very weakly interacting region R I the contrast is essentially unperturbed remaining unity in the course of the time-evolution and therefore there is no quasiparticle formation. For postquench interactions lying within R II or R II the contrast performs small and constant amplitude oscillations, weakly deviating from | Ŝ (t = 0) | = 1 [ Fig. 3 (f)]. This behavior indicates the generation of two long-lived coherent quasiparticles (see also section IV C). Entering the intermediate repulsive interaction region R III , | Ŝ (t) | exhibits large amplitude (0 < | Ŝ (t) | < 1) multifrequency temporal oscillations [ Fig. 3 (f)]. The latter signifies the dynamical formation of two Bose polarons which are coupled with higher-order excitations of the bosonic bath when compared to regions R II and R II as we shall expose in section IV D 1. For intermediate attractive interactions (region R III ) | Ŝ (t) | undergoes large amplitude oscillations taking values in Fig. 3 (f)]. This response of | Ŝ (t) | again signals quasiparticle formation. However, in addition to this dynamical dressing the destructive (| Ŝ (t) | = 0) and the constructive (| Ŝ (t) | ≈ 1) interference between the states of a single and two Bose polarons can be seen (see also Eq. (11) and its interpretation in section IV A).
For strong repulsive interactions lying within R IV the contrast shows a fastly decaying amplitude at short evolution times (0 < t < 2) and subsequently fluctuates around zero [ Fig. 3 (f)]. This latter behavior of | Ŝ (t) | 2 → 0 is a manifestation of an orthogonality catastrophe phenomenon of the spontaneously generated short-lived (0 < t < 2) Bose polarons. It is a consequence of the spatial phase separation between the impurity and the bosonic bath (see also Fig. 5 (h) and the discussion in section IV D 1), where the impurity prefers to reside at the edges of the BEC background, see also Fig. 2 (c). Note that this behavior is also supported by the effective potential of the impurities, see Eq. (9). Most importantly this process results in an energy transfer from the impurity to the BEC, which prohibits the revival of the dynamical state of the impurity to its initial one, implying | Ŝ (t) | 2 1. Such a mechanism has been also identified to occur for the case of a single impurity, see Ref. [35].
The emergence of the different dynamical regions in the evolution of the contrast holds equally when the size of the bath decreases to N B = 10 [ Fig. 3 (c)]. For such a few-body scenario region R II , where coherently longlived quasiparticles are formed, becomes slightly wider, i.e. 0.2 ≤ g R II B↑ < 0.6, compared to the N B = 100 case. The most notable difference between the few and the many particle bath takes place in the intermediate interaction region R III . The latter, occurs now at 0.6 ≤ g R III B↑ < 1.8, with | Ŝ (t) | performing large amplitude multifrequency oscillations implying in turn the formation of highly excited polaronic states. Note that the amplitude of the oscillations of | Ŝ (t) | here is larger than in the N B = 100 case [ Fig. 3 (a)]. Additionally, we observe that | Ŝ (t) | decreases smoothly as g B↑ increases, which is in sharp contrast to the N B = 100 case. Recall that such a smooth behavior occurring in the fewbody scenario has already been identified in our discussion of the ground state properties and in particular when inspecting the relative distance between the impurities. Also, the oscillations of | Ŝ (t) | (0 < | Ŝ (t) | < 1) for intermediate attractive interactions (region R III ) being a consequence of the destructive (| Ŝ (t) | = 0) and constructive (| Ŝ (t) | ≈ 1) interference between the states of a single and two Bose polarons are much more prevalent and regular for N B = 10 as compared to the N B = 100 case. Concluding, we can infer that the overall phenomenology of the dynamical formation of quasiparticles as imprinted in the contrast is similar for N B = 10 and N B = 100.
To test the effect of the number of impurities on the in- teraction intervals of quasiparticle formation we also consider the case of N I = 3 non-interacting, g II = 0, bosons immersed in a few-body bath of N B = 10 atoms. The dynamics of the corresponding contrast for this system following a quench from g B↑ = 0 to a finite either attractive or repulsive g B↑ is illustrated in Fig. 3 (d). As it can be seen, | Ŝ (t) | shows a similar behavior to the case of two impurities [ Fig. 3 (c)] but the regions of finite contrast become narrower. Particularly, the intermediate repulsive interaction region here occurs for 0.5 ≤ g R III B↑ < 1.5 instead of 0.6 ≤ g R III B↑ < 1.8 for N I = 2. Additionally, | Ŝ (t) | acquires lower values within the regions R III and R III for more impurities. Moreover, for N I = 3 within R III we observe a pronounced dephasing of the contrast which is absent for the N I = 2 case, see Figs.
3 (e 1 ), (e 2 ). As a consequence, we can deduce that the basic characteristics of the regions of dynamical polaron formation do not significantly change for a larger number of impurities in the regime N I N B .
Finally, we discuss | Ŝ (t) | for weakly interacting impurities. Comparing the temporal evolution of | Ŝ (t) | for g II = 0.2 [ Fig. 3 (b)] to the one for g II = 0 [ Fig. 3 (a)] we observe that the extent of the abovedescribed dynamical regions (R I , R II , R III , R IV , R II and R III ) can be tuned via g II . For instance, region R II occurs at 0.2 ≤ g R II B↑ < 0.4 for g II = 0.2 instead of 0.2 ≤ g R II B↑ < 0.5 when g II = 0, while region R III takes place at 0.4 ≤ g R III B↑ < 1.3 if g II = 0.2 and within 0.5 ≤ g R III B↑ < 1 in the non-interacting scenario. Also region R IV where the orthogonality catastrophe takes place is shifted to slightly larger interactions for g II = 0.2 compared to the g II = 0 case. Interestingly we observe that the contrast within R III and R III exhibits a decaying tendency for long evolution times t > 50 in the presence of weak impurity-impurity interactions, a behavior which is absent when g II = 0.

C. Spectrum of the contrast
To quantify the excitation spectrum of the impurity we calculate the spectrum of the contrast, namely . (14) Recall that at low impurity densities and weak interspecies interactions it has been shown that | Ŝ (t) | is proportional to the so-called spectral function of quasiparticles [18,99,100]. Figure 4 presents A(ω f ) in the case of a single and two either non-interacting (g II = 0) or weakly interacting (g II = 0.2) impurities when N B = 100 for different interspecies couplings of either sign. Evidently, for weak g B↑ belonging either to region R II with g B↑ = 0.25 [ Fig. 4 (a)] or R II with g B↑ = −0.25 [ Fig. 4 (d)] we observe a single peak in A(ω f ) located at ω f ≈ 4.27 and ω f ≈ −4.39 respectively. This single peak occurs independently of the number of impurities and their intraspecies interactions. Therefore, this peak at small g B↑ = ±0.25 corresponds to the long-time evolution of a well-defined repulsive or attractive Bose polaron respectively. Within region R III e.g. at g B↑ = 0.5 two dominant peaks occur in A(ω f ) [ Fig. 4 (b)] at frequencies ω f ≈ 8.42 and ω f ≈ 8.79 for both the N I = 1 and N I = 2 cases. Accordingly, these two peaks suggest the formation of a quasiparticle dressed, for higher frequencies, by higher-order excitations of the BEC background.
Entering the strongly interspecies repulsive region R IV a multitude of frequencies are imprinted in the impurity's excitation spectrum e.g. at g B↑ = 1.5, see Fig. 4 (c). The number of the emerging frequencies is larger for the two compared to the single impurity but does not significantly depend on g II for N I = 2. For instance, when N I = 1 mainly three predominant peaks centered at ω f ≈ 23.75, ω f ≈ 25.13, and ω f ≈ 26.26 appear in A(ω f ) whilst for N I = 2 and g II = 0 five dominantly contributing frequencies located at ω f ≈ 22.31, ω f ≈ 23.81, ω f ≈ 25.2, ω f ≈ 26.39 and ω f ≈ 27.52 occur. These frequency peaks correspond to even higher excited states of the quasiparticle than the ones within the region R III . We note that for values of g B↑ deeper in R IV a variety of low amplitude but large valued frequency peaks occur in A(ω f ). This fact indicates that the impurities tend to populate a multitude of states indicating the manifestation of the polaron orthogonality catastrophe as discussed in Refs. [35,75] (results not shown here).
Turning to intermediate attractive interactions lying within R III such as g B↑ = −0.5 a single frequency peak can be seen in A(ω f ) whose frequency is shifted towards more negative values for N I = 2 compared to N I = 1 and also for increasing g II The overall behavior of the excitation spectrum A(ω f ; g B↑ ) for N I = 2 and g II = 0 is shown in Fig.  4 (g) with varying g B↑ . Evidently, the position of the dominant quasiparticle peak in terms of ω f increases almost linearly for larger g B↑ . This behavior essentially reflects the linear increase of the energy of the initial state |Ψ(0) [Eq. (2)] directly after the quench. Moreover, comparing the position of the dominant quasiparticle peak with ∆ N I =1 + reveals that for g B↑ > 0.5, while the latter saturates, the former increases and additional peaks appear in the spectrum A(ω f ; g B↑ ). These peaks correspond to excited states of the system and already for g B↑ > 1 the ground states corresponding to ∆ N I + cease to be populated during the dynamics. In a similar fashion, such additional quasiparticle peaks occur also for attractive interactions, see Fig. 4 (g) for g B↑ < −0.5. In this case the additional quasiparticle peaks stem from the induced interactions resulting in the presence of a peak at and other ones which correspond to the occupation of higher-lying excited polaronic states with S z = 1 [Eq. (2)]. Note that such an almost linear behavior of the polaronic spectrum is reminiscent of the corresponding three-dimensional scenario but away from the Feshbach resonance regime. The latter corresponds in one-dimension to an interspecies Tonks-Girardeau interaction regime which is not addressed in the present work. We remark that in one-dimension there is no molecular bound state occuring for repulsive interactions.
Summarizing, we can infer that the quasiparticle excitation spectrum depends strongly on the value of the postquench interspecies interaction strength and also on the number of impurities outside the weakly attractive and repulsive coupling regimes [100]. However, this behavior is also slightly altered when going from two noninteracting to two weakly interacting impurities. For a relevant discussion on the lifetime of the above-described spectral features we refer the interested reader to Ref. [101]. It is also important to mention that in the weakly interacting impurity-BEC regime where the contrast is finite in the course of the evolution the spectral function A(ω f ) corresponds to the injection spectrum in the framework of the reverse rf spectroscopy [2,26].

D. Quench to repulsive interactions
Below we further analyze the dynamical response of the multicomponent system, and especially of the impurities, following an interspecies interaction quench from g B↑ = 0 to g B↑ > 0 within the above identified dynamical regions of the contrast. In particular, we explore the dynamics of the system on both the single-and the two-body level and further develop an effective potential picture to provide a more concrete interpretation of the emergent phenomena. We mainly focus on the nonequilibrium dynamics of two non-interacting impurities (g II = 0) and subsequently discuss whether possible alterations might occur for weakly interacting (g II = 0.2) impurities. Also, in the following, only the temporal-evolution of the pseudospin-↑ part of the impurities is discussed since the pseudospin-↓ component does not interact with the bosonic medium.

Density evolution and effective potential
To visualize the spatially resolved dynamics of the system on the single-particle level we first inspect the timeevolution of the σ-species single-particle density ρ (1) σ (x; t) [Eq. (6)] illustrated in Fig. 5. For weak postquench interspecies repulsions lying within the region R II e.g. g B↑ = 0.25, such that g B↑ < g BB , the impurities [see Fig.  5 (b)] exhibit a breathing motion of frequency ω I br ≈ 1.44 inside the bosonic medium [73,74]. Moreover, at initial evolution times (t < 60) the amplitude of the breathing is almost constant whilst later on (t > 60) it shows a slightly decaying tendency, see for instance the smaller height of the density peak at t = 70 compared to t = 20 in Fig. 5 (b). This decaying amplitude can be attributed to the build up of impurity-impurity correlations in the course of the evolution [42] due to the presence of induced interactions discussed later on, see also Fig. 7 (a). The breathing motion of the impurities is directly captured by the periodic contraction and expansion in the shape of the instantaneous density profiles of ρ (1) ↑ (x; t) depicted in Fig. 6 (b). On the other hand, the bosonic gas remains essentially unperturbed [ Fig. 5 (a)] throughout the dynamics, showing only tiny distortions from its original Thomas-Fermi cloud due to its interaction with the impurity.
An intuitive understanding of the observed dynamics of the impurities is provided with the aid of an effective potential picture. Indeed, the impurity-BEC interactions can be taken into account, to a very good approximation, by employing a modified external potential for the impurities. The latter corresponds to the time-averaged effective potential created by the harmonic oscillator and the density of the bosonic gas [35,51,73,75] namelȳ The averaging process aims to eliminate the emergent very weak distortions on the instantaneous density of the BEC ρ  B (x; t) with ω B br ≈ 1.82, hardly visible in Fig. 5 (a). They are canceled out in our case for T > 20. Note that ω B br < 2 is attributed to the repulsive character of the BEC background which negatively shifts its breathing frequency from the corresponding non-interacting value [102]. At g B↑ = 0.25 thisV ef f I (x) takes the form of a modified harmonic oscillator potential illustrated in Fig. 6 (a) together with the densities of its first few single-particle eigenstates. Furthermore, assuming      B (x, t) the effective trapping frequency of the impurities corresponds to ω ef f = ω 1 − g B↑ g BB . Therefore their expected effective breathing frequency would be ω ef f,I br = 2ω ef f ≈ 1.41 which is indeed in a very good agreement with the numerically obtained ω I br . The discrepancy between the prediction of the effective potential and the MB approach is attributed to the approximate character of the effective potential which does not account for possible correlation induced shifts to the breathing frequency. Moreover, in the present case the impurities which undergo a breath-ing motion withinV ef f I (x) reside predominantly in its energetically lowest-lying state E 1 , see Fig. 6 (a). It is also important to mention that this effective potential approximation is adequate only for weak interspecies interactions where the impurity-BEC entanglement is small [35,75]. Note also that the inclusion of the Thomas-Fermi approximation in the effective potential of Eq. (15) can not adequately describe the impurities dynamics when they reach the edges of the bosonic cloud, see [36] for more details. However in this case ρ (1) ↑ (x; t) lies within ρ (1) B (x; t) throughout the evolution indicating the miscible character of the dynamics for g B↑ < g BB [35,65]. Furthermore, for these weak postquench interspecies repulsions a similar to the above-described dynamics takes place also for two weakly (g II = 0.2) repulsively interacting impurities as shown in Fig. 5 (c). The impurities undergo a breathing motion within the bosonic medium in the course of the time-evolution exhibiting a slightly larger oscillation frequency than for the g II = 0 case but with the same amplitude [hardly visible by comparing Figs. 5 (b) and (c)].
For larger postquench interaction strengths g B↑ = 0.5 (region R III ), i.e. close to the intraspecies interaction of the bosonic bath g BB , the impurities show a more complex dynamics compared to the weak interspecies repulsive case [ Fig. 5 (e)]. Also, the BEC medium performs a larger amplitude breathing motion [ Fig. 5 (d)] compared to the g B↑ = 0.25 scenario but again with a frequency ω B br ≈ 1.82. Focusing on the impurities motion, we observe that at short evolution times (0 < t < 5) after the quench ρ (1) ↑ (x; t) expands and then splits into two counterpropagating density branches with finite momenta that travel towards the edges of the bosonic cloud, see Fig. 5 (e) and the profiles shown in Fig. 6 (d). The appearance of these counterpropagating density branches is a consequence of the interaction quench which imports energy into the system. Reaching the edges of ρ   ↑ (x; t) of the MB dynamics [ Fig. 6 (d)] reveals that the impurities predominantly reside in a superposition of the two lower-lying excited states (E 1 and E 2 ) ofV ef f I (x). Additionally in the case of two weakly repulsively interacting impurities, shown in Fig. 5 (f), the impurities' motion remains qualitatively the same. However, due to the inclusion of intraspecies repulsion the impurities possess a slightly larger overall oscillation frequency and the collisional patterns at the trap center appear to be modified as compared to the g II = 0 case.
Turning to strong postquench repulsions, i.e. g B↑ = 1.5 g BB which belongs to R IV , the dynamical response of the impurities is greatly altered and the bosonic gas exhibits an enhanced breathing dynamics as compared to the weak and intermediate interspecies repulsions discussed above. Initially ρ (1) ↑ (x; t = 0) consists of a density hump located at the trap center which, following the interaction quench, breaks into two density fragments, as illustrated in Fig. 5 (h), each of them exhibiting a multihump structure [see also Fig. 6 (f)]. Note that the density hump at the trap center remains the dominant contribution of ρ (1) ↑ (x; t) until it eventually fades out for t > 5, see Fig. 5 (h). This multihump structure building upon ρ (1) ↑ (x; t) is clearly captured in the instantaneous density profiles depicted in Fig. 6 (f). Remarkably, the emergent impurity density fragments that are symmetrically placed around the trap center (x = 0) perform a damped oscillatory motion in time around the edges of the Thomas-Fermi radius of the bosonic gas, see in particular Figs. 5 (g), (h).
The emergent dynamics of the impurities can also be interpreted to lowest order approximation (i.e. excluding correlation effects) by invoking the corresponding effective potential which for these strong interspecies repulsions has the form of the double-well potential shown in Fig. 6    ↑ (x; t) designate that the pseudospin-↑ impurities at initial times are in a superposition state of a multitude of highly excited states [see e.g. Fig. 6 (f) at t = 8] while for later times they reside in a superposition of lower excited states [see e.g. Fig. 6 (f) at t = 15]. We should also remark that a similar overall dynamical behavior on the single-particle level has been reported in the case of a single spinor impurity and has been also related to an enhanced energy transfer from the impurity to the bosonic bath [35,68,69,75]. Such an energy transfer process takes place also in the present case (results not shown here). Another important feature of the observed dynamical response of the impurities is the fact that they are not significantly affected by the presence of weak intraspecies interactions. This can be seen by inspecting Fig. 5 (i) which shows the time-evolution of ρ (1) ↑ (x; t) for g II = 0.2. Here, the most noticeable difference when compared to the g II = 0 scenario is that the splitting of ρ
to the time-evolution of the pseudospin-↑ impurity intraspecies two-body reduced matrix ρ (2) ↑↑ (x 1 , x 2 ; t) [Eq. (7)]. Recall that ρ (2) ↑↑ (x 1 , x 2 ; t) provides the probability of finding at time t a pseudospin-↑ boson at location x 1 and a second one at x 2 [65,66]. Most importantly, it allows us to monitor the two-body spatially resolved dynamics of the impurities and infer whether they move independently or correlate with each other [8,10,42].
↑ (x 2 , t = 0)/2 since g II = 0 and initially g B↑ = 0. As already discussed in section IV D 1 for weak interspecies postquench repulsions, namely g B↑ = 0.25 (region R II ), the impurities perform a breathing motion on the single-particle level [ Fig. 5 (b)] exhibiting a decaying amplitude for large evolution times. Accordingly, inspecting ρ (2) ↑↑ (x 1 , x 2 ; t) [Figs. 7 (a 1 )-(a 6 )] we observe that the impurities are likely to reside together close to the trap center since ρ (2) ↑↑ (−2 < x 1 < 2, −2 < x 2 < 2; t) is mainly populated throughout the evolution. In particular, at initial times ρ    Fig. 7 (a 5 )] indicating that the impurities tend to be slightly apart or at the same location respectively. This is indicative of the admittedly weak induced interactions as the breathing mode along the anti-diagonal of ρ (2) ↑↑ (x 1 , x 2 ; t) (relative coordinate breathing mode) does not possess exactly the same frequency as the breathing along the diagonal (center-ofmass breathing mode).
For larger interspecies repulsions e.g. for g B↑ = 0.5 (region R III ) the two-body dynamics of the impurities is significantly altered, see Figs. 7 (b 1 )-(b 6 ). At the initial stages of the dynamics the impurities reside together in the vicinity of the trap center as ρ is predominantly populated. However for later times two different correlation patterns appear in ρ (2) ↑↑ (x 1 , x 2 ; t) in a periodic manner. Recall that for these interactions ρ   ↑↑ (x 1 , x 2 ; t)] or they remain spatially separated with one of them residing in the left and the other in the right density branch [see the anti-diagonal elements of ρ (2) ↑↑ (x 1 , x 2 ; t)]. Moreover, during their collision at x = 0 the impurities are very close to each other as it is evident by the enhanced two-body probability in the neighborhood of x 1 = x 2 = 0 [ Fig. 7 (b 3 ), (b 5 )]. The dynamics of two weakly repulsive (g II = 0.2) impurities shows similar two-body correlation patterns to the noninteracting ones, as it can be seen by comparing Figs. 7 (b 1 )-(b 6 ) to (c 1 )-(c 6 ). This behavior complements the similarities already found at the single-particle level (see sec. IV D 1). The major difference on the two-body level between the g II = 0.2 and g II = 0 scenario is that in the former case ρ (2) ↑↑ (x 1 , x 2 ; t) is more elongated along its antidiagonal when the impurities collide at x = 0 [Figs. 7 (c 1 ), (c 3 )]. Therefore weakly interacting impurities tend to be further apart compared to the g II = 0 case, a result that reflects their direct repulsion. Other differences observed at the same time-instant in ρ (2) ↑↑ (x 1 , x 2 ; t) between the interacting and the non-interacting cases are due to the repulsive s-wave interaction that directly competes with the attractive induced interactions emanating in the system. For instance, shortly after a collision point e.g. at t = 55, shown in Figs. 7 (b 5 ) and (c 5 ), we observe that due to the repulsive s-wave interactions the attractive contribution between the impurities, see the diagonal of ρ Fig. 7 (c 5 )].
Turning to very strong repulsions, e.g. for g B↑ = 1.5 lying in region R IV , the correlation patterns of the two noninteracting impurities [Figs. 7 (d 1 )-(d 6 )] show completely different characteristics compared to the g B↑ ≤ g BB regime. Note here that in the dynamics of ρ   ↑↑ (x 1 , x 2 = x 1 ; t). Entering deeper in the evolution the aforementioned delocalization of the impurities becomes more enhanced since ρ  ↑↑ (x 1 , x 2 ; t) is inherently related to the multihump structure of ρ (1) ↑ (x; t) and suggests from a two-body perspective the involvement of several excited states during the impurity dy-namics. It is also worth mentioning that at specific time instants the diagonal of ρ

Two-body dynamics within the effective potential picture
To further expose the necessity of taking into account the intra-and the interspecies correlations of the system in order to accurately describe the MB dynamics of the impurities we next solve the time-dependent Schrödinger equation that governs the system's dynamics relying on the previously introduced effective potential picture [Eq. (15)] via exact diagonalization [103]. Thus our main aim here is to test the validity ofV ef f I (x) at least to qualitatively capture the basic features of the emergent nonequilibrium dynamics of the two impurities. We emphasize again thatV ef f I does not include any interspecies correlation effects that arise in the course of the temporal-evolution of the impurities. Within this approximation the effective Hamiltonian that captures the quench-induced dynamics of the impurities reads whereΨ ↑ (x) is the bosonic field-operator of the pseudospin-↑ impurity and g ↑↑ denotes the intraspecies interactions between the two pseudospin-↑ impurity atoms. Recall that the intercomponent contact interaction of strength g B↑ and the intraspecies interaction between the bath atoms are inherently embedded intō V ef f I [Eq. (15)]. In particular, withinV ef f I we account for the correlated Thomas-Fermi profile of the BEC since ρ (1) B (x; t) is determined from the MB approach. Below, we exemplarily study the dynamics of two noninteracting impurities and therefore we set g ↑↑ = 0 in Eq. (16). Moreover, in order to trigger the nonequilibrium dynamics we consider an interspecies interaction quench from g B↑ = 0 (t = 0) to a finite repulsive value of g B↑ . Such a sudden change is essentially taken into account via a deformation ofV ef f I [Eq. (15)]. The corresponding instantaneous two-body reduced density matrix of the impurities within H ef f is depicted in Fig. 8 for distinct values of g B↑ . Focusing on weak postquench interactions, e.g. g B↑ = 0.25, we observe that at the initial times the two-body dynamics of the impurities is adequately described within ↑↑ (x1, x2; t), of the two pseudospin-↑ non-interacting (gII = 0) bosonic impurities within the effective potential picture when considering an interspecies interaction quench to (a1)-(a6) g B↑ = 0.25, (b1)-(b6) g B↑ = 0.5 and (c1)-(c6) g B↑ = 1.5. The harmonically trapped bosonic mixture consists of NB = 100 atoms with gBB = 0.5 and NI = 2 impurities and it is prepared in its corresponding ground state configuration.
correct shape of ρ (2) ↑↑ (x 1 , x 2 ; t) and more precisely its deformations occuring along its diagonal or anti-diagonal [see Figs. 8 (a 4 )-(a 6 )] which stem from the build up of higher-order correlations during the dynamics.
Increasing the repulsion such that g B↑ = 0.5, deviations between the effective potential approximation and the correlated approach become more severe. For instance, at the initial times the sharp two-body probability peak of ρ (2) ↑↑ (x 1 , x 2 ; t) in the vicinity of x 1 = x 2 = 0 arising in the MB dynamics [ Fig. 7 (b 1 )] becomes smoother within H ef f [ Fig. 8 (b 1 )] although the overall shape of ρ (2) ↑↑ (x 1 , x 2 ; t) remains qualitatively similar. Moreover, the observed elongations along the diagonal of ρ ↑↑ (x 1 , x 2 ; t) of two different two-body configurations occurring at specific time-instants is also predicted at least qualitatively via H ef f , see Figs. 8 (b 2 ), (b 4 ) and (b 6 ). We remark that the differences in the patterns of ρ (2) ↑↑ (x 1 , x 2 ; t) between H ef f and the correlated approach are even more pronounced when g II = 0.2 (results not shown).
Strikingly for strongly repulsive interactions, g B↑ = 1.5, H ef f completely fails to capture the two-body dynamics of the impurities. This fact can be directly inferred by comparing ρ Even at the initial stages of the dynamics the effective potential cannot adequately reproduce the correct shape of ρ (2) ↑↑ (x 1 , x 2 ; t), compare Fig. 8 (c 1 ) with Fig. 7 (d 1 ). Note, for instance, the absence of the central two-body probability peak in the region −2 < x 1 , x 2 < 2 within H ef f which demonstrates the correlated character of the dynamics. More precisely, ρ  ↑↑ (x 1 , x 2 ; t) in the MB dynamics [see e.g. Fig. 7 (d 4 )-(d 6 )] is a pure correlation effect and a consequence of the participation of a multitude of excited states in the impurity dynamics which is never captured within H ef f .

E. Quench to attractive interactions
Next we discuss the dynamical behavior of both the BEC medium and the bosonic impurities on both the oneand the two-body level after an interspecies interaction quench from g B↑ = 0 to the attractive regime of g B↑ < 0. To explain basic characteristics of the dynamics of the impurities an effective potential picture is also employed. As in the previous section we first examine the emergent time-evolution of two non-interacting impurities (g II = 0) and then compare our findings to that of two weakly interacting (g II = 0.2) ones.

Single-particle dynamics and effective potential
To investigate the spatially resolved dynamics of the multicomponent system after an interaction quench from g B↑ = 0 to g B↑ < 0, we first analyze the spatiotemporal evolution of the σ-species single-particle density ρ (1) σ (x; t). The dynamical response of ρ (1) σ (x; t) triggered by the quench is presented in Fig. 9    B (x; t) characterized by a predominant frequency ω I br ≈ 2.76 and a secondary one ω I br ≈ 2.88 thus producing a beating pattern. These two distinct frequencies stem from the center-of-mass and relative coordinate breathing modes of the impurities, whose existence originates from the presence of attractive induced interactions in the system. We remark that the breathing frequency of the center-of-mass can be estimated in terms of the corresponding effective potential of the impurities, see also Eq. (17). In particular for g B↑ = −0.5, ω I br = 2 √ 2.06 ≈ 2.87 (see also the comment in Ref. [104]) which is in very good agreement with ω I br . The relevant contraction of ρ (1) ↑ (x; t) can be inferred by its increasing amplitude that takes place from the very early stages of the nonequilibrium dynamics [ Fig. 10 (b)]. The beating pattern can be readily identified e.g. by comparing the maximum height of ρ (1) ↑ (x; t) during its contraction at initial and later stages of the dynamics, see e.g. ρ (1) ↑ (x; t) at t = 10 and t = 40 in Fig. 9 (b). Moreover, as a consequence of the motion of the impurity and the relatively weak interspecies attraction, i.e. g B↑ = −0.5, the Thomas-Fermi cloud of the bosonic gas becomes slightly distorted. In particular, a low amplitude density hump is imprinted on ρ  ↑ (x; t) as shown by the white colored region in Fig. 9 (a) in the vicinity of x = 0 [75]. An almost similar effect to the above-mentioned breathing dynamics is present also for the case of two weakly interacting impurities [ Fig. 9 (c)]. Here, the secondary frequency manifests itself at later evolution times resulting in turn in a slower beating of ρ (1) ↑ (x; t) compared to the g II = 0 scenario [hardly visible in Fig. 9 (c)]. This delayed occurrence is attributed to the presence of intraspecies repulsion which competes with the attractive induced interactions.
For a larger negatively valued interspecies coupling, e.g. for g BI = −1 within region R III , ρ (1) ↑ (x; t) becomes more spatially localized and again performs a decaying amplitude breathing motion, the so-called beating identified above, but with a larger major frequency, ω I br ≈ 3.2, compared to the g B↑ = −0.5 case [ Fig. 9 (e)]. Notice that the observed beating motion of the impurities persists while being more dramatic for this stronger attraction [compare Figs. 9 (b) and (e)]. This enhanced attenuation of the breathing amplitude together with the strong localization of the impurities is a direct effect of the dominant presence of interspecies attractions between the impurity and the bath, see also Refs. [75]. Also, due to the stronger g B↑ and the increased spatial localization of ρ ↑ (x; t) exhibits a sech-like form tending to be more localized for a larger interspecies attractions g B↑ , see e.g. ρ  . We should remark that for large negative g B↑ the system becomes strongly correlated and the BEC is highly excited. The latter is manifested by the development of an overall weak amplitude breathing motion of the bosonic gas, see The above-mentioned dynamics can also be qualitatively explained in terms of a corresponding effective potential approximation [35,73,75]. Yet again, the effective potential experienced by the impurities consists of the external harmonic oscillator V (x) and the single-particle density of the BEC background. Importantly, since ρ (1) B (x; t) is greatly distorted from its original Thomas-Fermi profile due to the motion of the impurities, we invoke a time-averaged effective potential. Consequently, the effective potential of the impurity reads where T = 100 denotes the corresponding total propagation time. We remark that for the considered negative values of g B↑ the shape ofV ef f I (x) does not significantly change after averaging over T = 60. A schematic illustration ofV ef f I (x) and the densities of its first few singleparticle eigenstates at g B↑ = −1 is presented in Fig. 10  (a), see also remark [104]. The observed localization tendency of ρ  ↑ (x; t) for a larger |g BI | and thus its increasing localization tendency.

Two-body correlation dynamics and comparison to the effective potential approximation
Having described the time-evolution of the impurities on the single-particle level, we next analyze the dynamical response of the pseudospin-↑ component by invoking the corresponding two-body reduced density matrix ρ (2) ↑↑ (x 1 , x 2 ; t) [see also Eq. (7)]. The time-evolution of ρ (2) ↑↑ (x 1 , x 2 ; t) is depicted in Figs. 11 (a 1 )-(a 5 ) for two non-interacting (g II = 0) impurities following an interspecies interaction quench from g B↑ = 0 to g B↑ = −0.5 (region R III ). Before the quench the impurities lie together in the vicinity of the trap center since ρ (2) ↑↑ (x 1 = 0, x 2 = 0; t = 0) shows a high probability peak [ Fig. 11 (a 1 )]. However as time evolves the two bosons start to occupy a relatively smaller spatial region as can be deduced by the shrinking of the central twobody probability peak across the diagonal at t = 10 in Fig. 11 (a 2 ). Then they move either opposite to each other [see the elongated anti-diagonal in Figs. 11 (a 3 ), (a 5 )] or tend to bunch together at the same location [see the pronounced diagonal of ρ (2) ↑↑ (x 1 , x 2 = x 1 ; t = 60) in Fig. 11 (a 4 )]. This latter behavior of the impurities is the two-body analogue of their wavepacket periodic expansion and contraction (relative coordinate breathing motion) discussed previously on the single-particle level [ Fig. 9 (b)].
The dynamics of two weakly repulsively interacting (g II = 0.2) impurities [Figs. 11 (b 1 )-(b 5 )] shows similar characteristics to the above-described non-interacting scenario. Indeed, initially [ Fig. 11 (b 1 )] and at short times [ Fig. 11 (b 2 )] the impurities reside close to the trap center while later on they repel [see e.g. Fig. 11 (b 3 )] or attract [ Fig. 11 (b 4 )] each other as a result of their breathing dynamics [see also Fig. 9 (c)]. The major difference between the weakly interacting and the non-interacting impurities is that their distance which is   given by the anti-diagonal distribution of their two-body reduced density matrix is slightly different, see Fig. 11 (e 1 )-(e 5 ). For instance at t = 40 the non-interacting impurities are further apart from each other as compared to the case of interacting impurities, while this situation is reversed at t = 90. The aforementioned difference owes its existence to the distinct relative coordinate breathing frequencies. This can be directly inferred from the fact that ρ (2) ↑↑ (x 1 , x 2 ; t) possesses a larger spatial distribution when g II = 0.2 and it is attributed to their underlying mutual repulsion. For instance, even initially ρ (2) ↑↑ (x 1 , x 2 ; t = 0) for g II = 0.2 [ Fig. 11 (b 1 )] is slightly deformed towards its anti-diagonal compared to the g II = 0 case [ Fig. 11 (a 1 )]. This behavior persists also during the evolution independently of the expansion or the contraction of the impurity cloud, as can be seen by comparing Figs. 11 (b 4 ) to (a 4 ) and Figs. 11 (b 5 ) to (a 5 ).
To reveal the importance of both intra-and interspecies correlations for the impurity dynamics we then utilize the effective potential,V ef f I (x), introduced in Eq. (17) and solve numerically the time-dependent Schrödinger equation of the impurities via exact diagonalization. We remark once more thatV ef f I neglects the interspecies correlations of the multicomponent system but includes the density profile of the BEC determined by the MB approach. In particular, we construct the effective Hamiltonian H ef f of Eq. (16) but using thē V ef f I (x) of Eq. (17). For brevity we focus on the case of g ↑↑ = 0.2 and analyze the dynamics after an interspecies interaction quench from g B↑ = 0 (t = 0) to g B↑ = −0.5. As explained in Sec. IV D 3 within the effective potential picture this quench scenario accounts for the deformation ofV ef f I . Snapshots of ρ  11 (c 1 )-(c 5 )] significant deviations occur between the two methods. Indeed, during the time-evolution the correlation patterns visible in ρ (2) ↑↑ (x 1 , x 2 ; t) calculated via H ef f exhibit similar overall characteristics to the ones taking place in the correlated approach but at completely different time-scales. In fact, ρ (2) ↑↑ (x 1 , x 2 ; t) shows elongated shapes along its diagonal [ Fig. 11 (c 3 )] or anti-diagonal [ Fig. 11 (c 4 )] implying that the impurities tend to be relatively close or apart from one another respectively. The latter is again a manifestation of the breathing motion of the impurities at the two-body level. However H ef f fails in general to adequately capture the correct spatial shape of ρ Finally, turning to strong postquench attractions within R III , e.g. for g B↑ = −1 presented in Figs. 11 (d 1 )-(d 5 ), we observe that the two-body dynamics of the impurities is drastically altered with respect to the weakly attractive case g B↑ = −0.5 described above. Initially, at t = 0, the two bosons bunch together in the vicinity of the trap center since ρ Fig. 11 (d 1 )]. Subsequently the two-body distribution of ρ (2) ↑↑ (x 1 , x 2 ; t) spatially shrinks exhibiting a highly intense peaked structure around −0.2 < x 1 , x 2 < 0.2 as shown in Figs. 11 (d 2 ), (d 3 ). For longer evolution times ρ (2) ↑↑ (x 1 , x 2 ; t) deforms possessing an elongated shape across its diagonal [see Figs. 11 (d 4 ), (d 5 )] which indicates that the impurities experience a mutual attraction. This latter behavior suggests the appearance of attractive induced interactions between the impurities mediated by the bosonic gas.

V. SUMMARY AND CONCLUSIONS
We have investigated the ground state properties and the interspecies interaction quench quantum dynamics of two spinor bosonic impurities immersed in a harmonically trapped bosonic gas from zero to finite repulsive and attractive couplings. For two non-interacting impurities, we have shown that for an increasing attraction or repulsion their overall distance decreases indicating the presence of attractive induced interactions. Moreover, at strong attractions or repulsions the impurities acquire a fixed distance and bunch together either at the trap center or at the edge of the Thomas-Fermi profile of the bosonic gas respectively. For two weakly repulsive impurities we find that their ground state properties remain qualitatively the same for attractive couplings, but for repulsive interactions they move apart being located symmetrically with respect to the trap center. A similar to the above-described overall phenomenology takes place for smaller system sizes and heavier impurities.
Regarding the quench dynamics of the multicomponent system we have analyzed the time-evolution of the contrast and its spectrum. We have revealed the emergence of six different dynamical response regions for varying postquench interaction strength which signify the existence, dynamical deformation and the orthogonality catastrophe of Bose polarons. We have also shown that the extent of these regions can be tuned via the intraspecies repulsion between the impurities, the impurity concentration and the size of the bath. Moreover, we have found that the polaron excitation spectrum depends strongly on the postquench interspecies interaction strength and the number of impurities but it is almost insensitive on the impurity-impurity interaction for the weak couplings.
Focusing on weak postquench interspecies repulsions the non-interacting impurities perform a breathing motion manifested as a periodic expansion and contraction of their density on both the one-and two-body level. For an increasing repulsion the impurities single-particle density splits into two counterpropagating density branches that travel to the edges of the BEC medium where they are reflected back towards the trap center and subsequently collide, repeating this motion in a periodic manner. Here the impurities mainly reside in a superposition of two distinct two-body configurations, namely they either reside together or each one lies at a specific density branch, while during their collision they tend to remain very close to each other. In the strong repulsive regime we have observed that the density of the impurities shortly after the quench breaks into two fragments which are symmetric with respect to the origin and which exhibit a multihump structure and perform a damped oscillatory motion close to the Thomas-Fermi radius of the bosonic gas. This multihump structure leads to a spatially delocalized behavior of the corresponding twobody correlation patterns and suggests the involvement of higher excited states. In all cases the bosonic gas exhibits a breathing motion whose amplitude becomes more pronounced for an increasing repulsion.
Turning to attractive interspecies couplings, the impurities show a beating breathing motion and experience a spatial localization tendency at the trap center on both the one-and two-body level, a behavior that becomes more pronounced for larger attractions. Strikingly, for strong attractive interactions we unveil that gradually the impurities experience a mutual attraction on the two-body level. This effect demonstrates the pronounced presence of induced interactions for attractive interspecies ones. As a result of the impurities motion the density of the bosonic bath deforms, developing a low amplitude density hump located at the origin. The occurrence of this hump is a direct consequence of the presence of induced interactions.
In all cases investigated in the present work, an in-tuitive understanding of the dynamics of the impurities is provided via an effective potential picture which is shown to be an adequate approximation for weak couplings where correlations are negligible. However, for increasing interaction strengths this effective model largely fails to adequately describe the dynamics on both the one-and two-body level due to the presence of both induced attraction and higher-order correlations. Finally, in all of the above-mentioned cases we showcase that a similar dynamical response takes place for two weakly repulsive impurities but the corresponding time-scales are slightly altered due to the competition between their mutual repulsion and the developed attractive induced interactions.
There is a multitude of fruitful possible extensions of the present effort that can be addressed in future works. A intriguing aspect would be to examine whether thermalization of the impurities dynamics takes place for strong repulsions in the framework of the eigenstate thermalization hypothesis [105]. An imperative prospect is to study the robustness of the emergent quasiparticle picture in the current setting in the presence of temperature effects [106,107]. Moreover, the study of induced interactions of two bosonic impurities immersed in a Fermi sea would be an interesting prospect especially in order to expose their dependence on the different statistics of the medium. Additionally, the generalization of the present results to higher-dimensional settings would be highly desirable. Another interesting direction would be to investigate the collisional dynamics of subsonically or supersonically moving impurities in a lattice trapped bosonic gas. Here, one could unravel the properties of the emergent quasiparticles, such as their lifetime, residue, effective mass and induced interactions with respect to the interspecies interaction strength.
Appendix A: Remarks on the many-body simulations To solve the underlying time-dependent MB Schrödinger equation of the considered multicomponent system we invoke the Multi-Layer Multi-Configurational Time-Dependent Hartree Method for Atomic Mixtures (ML-MCTDHX) [70,71]. As discussed in Section II B it constitutes a variational approach for calculating the stationary and most importantly the nonequilibrium quantum dynamics of bosonic and fermionic multicomponent mixtures [35,36,65] including spin degrees of freedom [9,35,82]. A key advantage of the method is that it assumes the expansion of the total MB wavefunction in terms of a time-dependent and variationally optimized basis. Such a treatment enables us to capture both the intra-and intercomponent correlation effects by employing a computationally feasible basis size. The latter flexibility allows to span the relevant subspace of the Hilbert space efficiently for each time-instant which is in contrast to numerical methods relying on a time-independent basis.  (5)]. Additionally, within our implementation a sine discrete variable representation (sine-DVR) is utilized as the primitive basis for the spatial part of the SPFs with M = 600 grid points. The latter intrinsically introduces hard-wall boundary conditions at both edges of the numerical grid imposed herein at x ± = ±50. We have ensured that the position of the hard-walls does not affect the presented results by assuring that no appreciable density occurs beyond x ± = ±20. The eigenstates of the composite MB system are obtained by means of the so-called improved relaxation method [70,71] implemented in ML-MCTDHX. In order to simulate the nonequilibrium dynamics we propagate in time the wavefunction [Eq. (3)] utilizing the appropriate Hamiltonian within the ML-MCTDHX equations of motion.
To infer the convergence of our MB simulations we ensure that all observables of interest, e.g. | Ŝ (t) |, ρ ↑ (x; t), become to a certain degree insensitive upon varying the employed orbital configuration space chosen herein to be C = (D; d B ; d I ) = (12; 3; 10). Below, we exemplarily showcase the convergence behavior of the contrast during evolution for a system composed of N B = 100 bosons with g BB = 0.5 and N I = 2 non-interacting (g II = 0) impurities. More precisely, we investigate its absolute deviation between the C = (10; 3; 10) and other orbital configurations C = (D; d B ; d I ) during the nonequilibrium dynamics, namely The time-evolution of ∆ |S(t)| C,C is illustrated in Fig.  12 after an interspecies interaction quench from g B↑ = 0 to intermediate repulsions e.g. g B↑ = 1 [ Fig. 12 (a)] and strong ones such as g B↑ = 4 [ Fig. 12 (b)]. As it can be readily seen by inspecting ∆ |S(t)| C,C , a systematic convergence of | Ŝ (t) | can be achieved in both cases. At intermediate postquench repulsions, e.g. g B↑ = 1, ∆ |S(t)| C,C e.g. between the C = (12; 3; 10) and C = (10; 3; 8) [C = (8; 3; 8)] orbital configurations acquires a maximum value of the order of 3% [7%] at large propagation times as shown in Fig. 12 (a). As expected, an increasing g B↑ yields a larger relative error [ Fig. 12 (b)] but still remaining at an adequately small degree. Indeed, turning to strong repulsions such as g B↑ = 4 we observe that the deviation ∆ |S(t)| C,C with C = (12; 3; 10) and C = (12; 3; 8) [C = (10; 3; 8)] lies below 5% [9%] throughout the evolution, see Fig. 12 (b). Finally, we should mention that a similar analysis has been performed for all other interspecies interaction strengths and observables discussed in the main text and found to be adequately converged (results not shown here for brevity).