Crossover from attractive to repulsive induced interactions and bound states of two distinguishable Bose polarons

We study the impact of induced correlations and quasiparticle properties by immersing two distinguishable impurities in a harmonically trapped bosonic medium. It is found that when the impurities couple both either repulsively or attractively to their host, the latter mediates a two-body correlated behavior between them. In the reverse case, namely the impurities interact oppositely with the host, they feature anti-bunching. Monitoring the impurities relative distance and constructing an effective two-body model to be compared with the full many-body calculations, we are able to associate the induced (anti-) correlated behavior of the impurities with the presence of attractive (repulsive) induced interactions. Furthermore, we capture the formation of a bipolaron and trimer state in the strongly attractive regime. The trimer refers to the correlated behavior of two impurities and a representative atom of the bosonic medium and it is characterized by an ellipsoidal shape of the three-body correlation function. Our results open the way for controlling polaron induced correlations and creating relevant bound states.

Interestingly, it was predicted [37,39] that there is also the possibility of mediating repulsive impurity-impurity interactions when two impurities are coupled with different signs to a bosonic bath.In this sense, the underlying experimentally relevant three-component system [16,18] allows to unravel additional polaronic properties as it has been also argued by immersing impurities into a two-component pseudospinor mixture [4][5][6][63][64][65][66] in order to create, for instance, spin-wave excitations and magnetic polarons [4,5], impurities diffusive response [64] or to facilitate the detection of the dressing cloud via interferometry [4].However, quasiparticle formation in three-component systems is largely unexplored, besides the few above-mentioned recent studies.An interesting direction is to exploit the tunability of such mixtures, e.g. in terms of different intercomponent couplings, for devising the ground state quasiparticle properties such as the impurities effective mass and induced interactions.Here, it is important to understand the interplay of the latter properties and the underlying impurities' correlations.Also, the formation of relevant bound states either solely among the impurities (bipolarons) or between the impurities and the host atoms (trimers) remains elusive.To address these questions, we consider two distinguishable and non-interacting impurities that are embedded into an one-dimensional bosonic gas.The impurities' couplings with the host are individually tuned spanning the regime from attractive to repulsive interactions.Here, the effective interactions between the impurities can be only mediated in the presence of impurity-medium entanglement and bound states require the involvement of strong correlations.As such, to account for the relevant inter and intra-component correlations we employ the variational multilayer multiconfiguration time-dependent Hartree method for atomic mixtures (ML-MCTDHX) approach [67][68][69] which is well established for investigating impurity physics [36].
Inspecting the spatial two-body correlations between the two impurities we reveal that, in general, they are correlated (anti-correlated) when the two impurity-medium coupling strengths posses the same (opposite) sign.To shed more light on the impact of induced impurities' correlations we carefully monitor their relative distance [70], excluding all mean-field type contributions, for varying coupling strengths.A central result of our work is that the impurities' correlated (anti-correlated) behavior is related to a decrease (increase) of their relative distance, thus, indicating the presence of an induced attraction (repulsion) between them.This observation is additionally confirmed by constructing an effective two-body model in the weak impurity-medium coupling regime inspired from the case of indistinguishable impurities [36,52].It specifically allows to assign the impurities' induced interaction strength and sign but also other quasiparticle related properties such as their effective mass and trap frequency.
For strong impurity-medium attractions, we identify the formation of a bipolaron state involving the two distinguishable impurities.This bound quasi-particle state is characterized by the so-called bipolaron energy [32], and the size of the impurities' dimer state featuring an exponential decrease for larger attractions.Proceeding a step further, we find that for such strong attractive impurity-medium interactions the three-body correlation function features an ellipsoidal shape indicating bunching and revealing the creation of a trimer state among the two impurities and a corresponding bath atom.To further testify the existence of this trimer state we employ the Jacobi relative distances of the three distinguishable atoms [71] showing an exponentially decreasing trend for increasing impurity-medium attractions.
This work is organized as follows.In section 2, the three-component setup under consideration is introduced and in Section 3 we explain the variational method used to obtain the ground state properties of the many-body system.Section 4 elaborates on the possible ground state density configurations upon varying the impurity-medium couplings.The emergence of induced impurity-impurity correlation patterns is explicated in Section 5.The interrelation of the aforementioned induced correlations with the induced attractive and repulsive impurity in-teractions is provided in Section 6 through monitoring their relative distance and constructing an effective two-body model.Delving into the strongly attractive impurity-medium interaction regime, we demonstrate the formation of a bipolaron state among the two distinguishable impurities in Section 7 and the generation of a trimer state among the impurities and a bath atom in Section 8. We summarize our findings and discuss future perspectives in Section 9.The behavior of the logarithmic negativity in order to quantify the bipartite intercomponent entanglement is discussed in Appendix A. Appendices B and C provide supplemental information regarding the polaron characteristics and induced effective interactions.In Appendix D we comment on the impact of the impurity mass and the number of bath particles on the ground state properties of the system.Finally, in Appendix E we elaborate on the microscopic excitation processes of the system via a number state analysis.

Two distinguishable impurities in a bosonic gas
We consider a one-dimensional harmonically trapped three component mixture.It contains a bosonic medium A with N A = 15 atoms of mass m A and two distinguishable impurities B and C, i.e., N B = N C = 1, having masses m B and m C , respectively.The many-body Hamiltonian of this system reads where Ĥσ denotes the Hamiltonian of each component σ and Ĥσσ ′ represents the intercomponent interaction contribution with σ, σ ′ ∈ {A, B, C}.Specifically, Assuming that the system is at ultracold temperatures it dominantly experiences s-wave scattering processes that can be described by two-body contact interactions between particles of the same as well as of different species characterized by the generic strength g σσ ′ [14].The latter depends on the respective three-dimensional scattering lengths a 3D σσ ′ and the transversal confinement frequency ω ⊥ that are experimentally tunable via Feshbach resonances [13,14] and confinement induced resonances respectively [12].The latter would allow the tuning of interactions even in the absence of a Feshbach resonance.
For simplicity, we focus on the mass-balanced case m σ ≡ m (unless stated otherwise) and thus ω σ ≡ ω.Moreover, we rescale our Hamiltonian in harmonic oscillator units ħ hω which means that the length and interaction scales are given in ħ h/mω and ħ h 3 ω/m, respectively.Such a three-component system could be experimentally realized [16,18] e.g., by trapping three different hyperfine states of 87 Rb which can feature various Feshbach resonances.An alternative candidate may be two isotopes of Rubidium atoms with 85 Rb emulating the medium and two-hyperfine states of 87 Rb [72,73] representing the impurities.Since our main findings persist also for mass-imbalanced mixtures, see the discussion in Section D, corresponding heteronuclear mixtures of different isotopes could also be used.We also note that the experimental realization of three-component mixtures was reported in Refs.[16,18] and a proposal for a corresponding impurity system was recently made in Ref. [74].Since our aim is to understand the role of induced interactions between the impurities mediated by the medium, in the ground state of the system, it is natural to consider two non-interacting impurities setting g BC = 0, which could be realized, for instance, via magnetic Feshbach resonances [75].

Variational wave function approach
The ground state of the three-component mixture, described by the Hamiltonian of Eq. ( 1), is determined within the ML-MCTDHX method [67][68][69]76].A central aspect of this ab-initio approach is based on the expansion of the many-body wave function on different layers using a variationally optimized time-dependent many-body basis.This leads to an efficient truncation of the underlying Hilbert space tailored to capture the relevant inter-and intracomponent correlations.Specifically, the many-body wave function is first expressed in terms of three different sets of D σ species functions as follows The time-dependent coefficients C i jk (t) bare information about the entanglement between the involved components.For instance, the bipartite entanglement between two components can be analyzed by tracing out the degrees of freedom of the third one and then apply the positive partial transpose criterion on the resulting mixed state [77] (see also Appendix A).Next, the intracomponent correlations are included into the wave function ansatz by expanding each species function as a superposition of permanents |⃗ n(t)〉 weighted by time-dependent expansion coefficients ) represents the occupation distribution of N σ particles on d σ time-dependent single-particle functions.Additionally, the single-particle functions are expanded into a time-independent discrete variable representation [78] consisting in our case of M r = 300 evenly spaced grid points.
The number of utilized species functions D σ dictates the degree of intercomponent entanglement.For instance, by providing only one species function for each component, i.e., by setting D A = D B = D C = 1, the many-body wave function reduces on its top layer to a product state, thereby, prohibiting any interspecies entanglement.Such a treatment is commonly referred to as a species mean-field ansatz (sMF) [67].For two-component mixtures the sMF ansatz is unique, however, in three-component systems there are various sMF that could be constructed.As an example, setting D σ = 1 and D σ ′ , D σ ′′ > 1, we allow for entanglement generation only between the species σ ′ and σ ′′ , whilst intercomponent correlations with species σ are suppressed.To clearly distinguish among the different possible sMF ansatzes, in the following, we abbreviate as sMFσ where σ ∈ {A, B, C} the ansatz that ignores intercomponent correlations between species σ and the remaining ones.In this sense, the sMFC is written as where only species A and B can become entangled while species C remains uncorrelated with the other species.
The ground state of the three component mixture is obtained through the imaginary time propagation method.The time-dependent coefficients of each layer, namely the species and single-particle layers, are optimally adapted to the system, e.g. by following the Dirac-Frenkel variational principle [79] in order to determine the underlying ML-MCTDHX equations of motion.The latter correspond to D A D B D C linear differential equations of motion for the C i jk (t) nonlinear integrodifferential equations for the species functions and d A +d B +d C nonlinear integrodifferential equations for the single-particle functions.This co-moving basis concept minimizes the number of required states for achieving numerical convergence.In this sense, it reduces the computational cost as compared to methods relying on time-independent basis sets, while simultaneously allows to account for all relevant correlations.The truncation of the Hilbert space is determined by the number of employed species-and single-particle functions defining the numerical configuration space Utilizing this method, it is in principle possible to describe mixtures with mesoscopic particle numbers and strong interactions.However, as the number of particles increases and correlations become enhanced a larger number of orbitals should be taken into account in order to reach numerical convergence.The latter is carefully checked by systematically increasing the numerical configuration space and ensuring that the observables of interest remain unchanged within a desired level of accuracy.As expected, this process is accompanied by a significant computational cost and in particular it is the interplay of intraand intercomponent correlations with the components atom number that limits the applicability of the method due to numerical feasibility.Elaborated discussions on the ingredients, applicability and benchmarks of this variational method to different multicomponent settings can be found in the recent reviews [36,80].
For our system, the degree of correlations in the bosonic bath, e.g. as captured by its depletion [81] 1−n A 0 with n A 0 representing the largest eigenvalue of the bath's one-body reduced density matrix is negligible within the considered interaction strength intervals.This allows us to use only a few orbitals for the medium in order to ensure convergence.On the other hand, the impurities depletion is in general larger, especially for strongly repulsive interactions, and thus we need to use more orbitals.Herewith, we have checked that employing an orbital configuration (6, 6, 6; 4, 6, 6) results in the convergence of the observables of interest, such as the species densities and intercomponent two-body correlation functions, while the amount of equations of motion are tractable.For completeness, let us note that stronger intercomponent interactions than the ones to be reported below e.g.|g AC | < 10 require a larger number of species functions and impurities orbitals which is still numerically feasible.Similarly, in order to tackle systems with stronger intracomponent bath interactions the number of the respective d A orbitals should be increased.This naturally entails more difficult convergence issues than increasing the impurities orbitals (and thus considering stronger impurity-medium interactions) since the number of the underlying equations of motion becomes larger in the former case.

One-body density configurations of the three-component mixture
To investigate the emergent spatial configurations of the three-component impurity setting arising due to different combinations of the involved interactions, we initially employ the σcomponent one-body density being normalized to unity.Namely, ρ (1) σ denotes the bosonic field operator which annihilates (creates) a σ-species atom at position x.In an experiment, the density is routinely detected through in-situ absorption imaging [82][83][84].Our understanding on the mixture spatial distributions at different interactions is also corroborated by an effective potential picture, which has been proven thus far successful in order to qualitative explicate various aspects of impurity physics in twocomponent settings [70,85,86].According to this, each σ component is subjected to an effective potential stemming from the superposition of its external harmonic trap and the density of the complementary components σ ′ weighted by the respective intercomponent interactions, (1) Figure 1: One-body σ-species density, ρ (1)  σ (x), shown together with the effective potentials [Eq.(6)] of the impurities (see legend).Two distinguishable non-interacting impurities (B, C) are considered which are individually coupled to a bosonic medium A with g AA = 0.2.The impurity-medium coupling strengths from left to right panels refer to (g AB , g AC ) = (−1.0,−0.2), (1.0, −0.2) and (1.0, 1.5).For attractive interactions the medium atoms accumulate in the vicinity of the impurities and their effective potential is attractive.Turning to repulsive couplings a tendency for impuritymedium phase-separation occurs for g Aσ > g AA with σ = B, C. i.e., Naturally, this is a sMF framework since it ignores intercomponent correlations.Moreover, it is more meaningful for the impurity subsystem since the impact of the impurity densities is suppressed for the medium.Density profiles of all three components and the impurity effective potentials are provided in Fig. 1 for characteristic impurity-medium interaction configurations, namely (g AB , g AC ) = (−1.0,−0.2), (1.0, −0.2) and (1.0, 1.5).The impurities are considered to be non-interacting among each other, i.e., g BC = 0, and the medium bosons feature throughout g AA = 0.2.As it can be seen, for an overall attractive impurity-medium coupling the bosons of the medium are placed in the vicinity of the impurities which are naturally localized at the trap center [cf. Figure 1(a)].This distribution of the medium atoms can also be understood in terms of the respective attractive impurity-medium interaction energy E int Aσ = 〈Ψ MB |H Aσ |Ψ MB 〉 for g Aσ < 0 with σ = B, C. Also, for both g AB < 0 and g AC < 0 the effective potential of each impurity corresponds to a dipped harmonic trap enforcing its localization whose degree is, of course, enhanced for stronger attractions [cf. Figure 1(a)].The value of the attractive interaction determines the degree of spatial localization, i.e., the B impurity with g AB = −1.0 is more localized than the C impurity experiencing g AC = −0.2.For sufficiently large attractive impurity-medium couplings (|g Aσ | ≫ g AA ) the impurities form a bipolaron, see for details the discussion in Section 7.
On the other hand, tuning at least one of the impurity-medium couplings towards the repulsive regime such that g Aσ > g AA is satisfied leads to the phase-separation among these components since E int Aσ > 0. In this case, the impurity forms a shell around the edges of the bath residing around the trap center [65].Such configurations can be readily observed, for instance, in Figure 1(b) where solely the B impurity is strongly repulsively coupled with the bath (g AB > g AA ) and also in Figure 1(c) where both impurities phase separate with the bath due to g AB > g AA and g AC > g AA .Notice that for strong repulsive impurity-medium couplings the underlying effective potential of the impurity has the form of a double-well potential which

G
(2) Another interesting phenomenon reflecting the richness of three-component systems arises upon considering distinct interactions between each impurity and the bath.Indeed, varying the impurity-medium coupling for a specific impurity affects the shape of the bath accordingly and, in turn, impacts the distribution of the other impurity.This is visualized in Figures 1(a) and (b) where g AC is the same while g AB is modified from attractive to repulsive values ultimately altering the spatial localization of impurity C, see in particular the peak of ρ Therefore, it is possible to implicitly manipulate the distribution of one impurity by adjusting the coupling of the other impurity with the bath and importantly in the absence of direct impurity-impurity interaction.This property, as it will be discussed below, can be proved crucial for controlling impurity-impurity induced interactions.

Intercomponent (induced) correlations and entanglement
Next, we shed light on the associated intercomponent correlation patterns with a particular emphasis on the existence of induced correlations between the impurities mediated by the bosonic gas.The intercomponent two-body spatial correlations, or two-body coherence, can be quantified through [84], Here, we subtract the probability of independently detecting a σ and a σ ′ atom at positions x σ 1 and x σ ′ 2 from the probability to simultaneously measure one at x σ 1 and the other at The latter is provided by the reduced two-body density which is normalized to unity.In this sense, the two particles are correlated or bunched (anticorrelated or antibunched) if G (2) ) is positive (negative); otherwise, they are referred to as two-body un-correlated [84,87].Experimentally the two-body correlation function is accessible through analyzing the respective single-shot images, see e.g.Refs.[88][89][90][91][92].

Characteristic correlation patterns
First, we study the emergent two-body correlation patterns between the B impurity and the medium for different intercomponent interactions [Figures 2(a1)-(c1)].For attractive g AB < 0 and g AC < 0 the B impurity is correlated with a bath atom at the same position, see the diagonal of G (2) 2 ) > 0, while these two particles are anti-correlated when symmetrically placed with respect to the trap center as it is shown from the anti-diagonal of G (2) In this sense, the B impurity prefers to occupy the same spatial region with the bath.Turning to repulsive g AB > 0 and independently of g AC ≶ 0, the above-discussed two-body correlation distributions are inverted and the B impurity features an anti-bunched (bunched) behavior at the same (different) location with a bath particle as can be deduced by the diagonal (anti-diagonal) of G (2) Let us now discuss the induced correlations among the non-interacting impurities.When both impurities are attractively coupled to the bath they exhibit a bunching tendency which is, of course, mediated by the bosonic gas, see the diagonal of G (2) Otherwise, the impurities are anti-bunched when residing at different locations with respect to x = 0.This two-body configuration of the impurities manifests the presence of their attractive induced interactions regulated by the impurity-medium attractive interactions as we will discuss in Section 6.Note also that a further increase of the impurity-bath attraction can result in the formation of a bipolaron state which we analyze in detail within Section 7. A similar two-body impurities correlation pattern occurs when they both repulsively couple with their bath [Figure 2(c2)].However, in this case the impurities cluster either at the left or the right side of the bath, while the probability to reside at opposite sides is suppressed . This trend which is inherently related to the impurity-medium phase-separation has also been observed for two indistinguishable impurities and it is known as their coalescence [42].In sharp contrast, if one impurity couples repulsively and the other attractively to the bath the reverse to the above-described correlation behavior is observed.Namely, the impurities anti-bunch (bunch) at the same (different) location in terms of the trap center, see Figure 2(b2).This scenario manifests the flexibility offered by three component mixtures and it is connected to the emergence of repulsive impurity-impurity induced interactions, a phenomenon that can not occur in two-component systems and we analyze in Section 6.

Emergent correlation regimes
To provide an overview of the two-body correlation behavior stemming from the interplay of the distinct impurity-medium couplings, we inspect the spatially integrated over [−∞, 0 ] (due to symmetry) correlation function It quantifies the amount of intercomponent correlations or anti-correlations by means that it is positive (negative) when the particles prefer (avoid) to occupy the same region with respect to the trap center. 1 The phase diagrams of the impurity-medium C AB and impurity-impurity C BC integrated correlations as a function of g AB and g AC are depicted in Figure 3(a) and (b) respectively.Recall that since g BC = 0 all emerging impurity correlations are induced by their coupling to the bath.An anti-correlated (correlated) behavior between the B impurity and the bath occurs for g AB > 0 (g AB < 0) and varying g AC , see also Figures 2(a1)-(c1).Notice also the un-correlated tendency for strongly attractive g AC and repulsive g AB [Figures 3(a), (b)].Indeed, due to the large g AC < 0 both the bath A and the C impurity localize at the trap center minimizing their spatial overlap with the B impurity since g AB > 0 and thus C AB is suppressed.Naturally, a less attractive g AC enhances the overlap between impurity B and the bath leading to an anticorrelated behavior.The largest degree of anti-correlation as captured by C AB is reached when g AB > g AA and g AC > g AA where both impurities form a shell around the bath and coalesce [cf.corresponding region in Figure 3(a)].
Turning to the impurities' correlations, we observe that as long as they both couple either repulsively or attractively to the bath it holds that C BC > 0, implying that they are correlated [see also Figures 2(a1) and (c1)].However, when the couplings g AB and g AC have opposite signs, with one lying in the weak and the other in the strong interaction regime, then mostly C BC < 0, i.e., the impurities are anti-correlated [cf. Figure 2(b1)].A notable exception takes place if one of the impurities couples strongly repulsively to the bath (e.g.g AB > g AA ) and the other strongly attractively (e.g.|g AC | > g AA ).This leads to a suppressed spatial overlap among the bath and the repulsively interacting impurity2 and thus the bath is only correlated with the attractively coupled impurity, see also the discussion above.Together with the fact that the impurities are spatially separated in this interaction region, if mediated impurity correlations occur they have to be nonlocal.This is indeed the case since the impurities are found to be anti-correlated, C BC < 0, see the two parameter regimes in Figure 3(b) enclosed by the dashed lines.

Quantification of impurities induced interactions
Below, we examine how the mediated correlations among the distinguishable impurities alter their relative distance and, subsequently, relate the induced impurity-impurity correlation patterns with an effective induced interaction strength.The latter as it will be argued can be either attractive or repulsive due to the genuine three-component nature of the system and it is further quantified via an effective two-body model.

Effect of the induced impurity-impurity correlations on their relative distance
A reliable measure for this purpose, that has also been utilized in two-component settings [70,87] and can be experimentally monitored via in-situ spin-resolved single-shot measurements [93], is the relative distance between the impurities Specifically, in order to extract the contribution stemming from genuine impurity-medium correlations we estimate the modified relative distance at different correlation levels as dictated by the respective truncation of the many-body (MB) wave function (see also Section 3), namely Here, sMF stands for the general species mean-field case where all intercomponent correlations are neglected, while sMFB (sMFC) refers to the case at which only intercomponent correlations between the B (C) impurity and the medium are ignored [36,65].Excluding the sMF contribution as well as the ones corresponding to the entanglement between the bath and impurity C or B [cf. last four terms of Eq. ( 11)] from the relative distance where all correlations are included, i.e., 〈r MB BC 〉, we are able to distill the effects originating from the mutual correlation among the impurities and the bosonic gas by tracking ∆〈r BC 〉.As such, ∆〈r BC 〉 captures the genuine effects of the induced correlations as described by C BC [Figure 3(b)].We interpret a value of ∆〈r BC 〉 which is positive (negative) as the signal of emergent repulsive (attractive) impurities' induced interactions.
The modified relative distance, ∆〈r BC 〉, is presented in Figure 4(a) with respect to the g AC coupling and for characteristic fixed g AB values.In general, we find an induced attraction between the impurities when they both couple either attractively or repulsively to the medium, while they feature a mediated repulsion if one of them couples attractively and the other repulsively to the bosonic gas.Since ∆〈r BC 〉 is closely related to C BC , an induced correlation (anti-correlation) between the impurities can be associated to their attractive (repulsive) BC predicted within the many-body approach (see main text).(c) Fidelity F BC of the impurities wave function as found in the many-body method and the effective two-body model with respect to the impurity-medium couplings g AB and g AC .We consider two non-interacting but distinguishable impurities immersed in a bosonic gas of N A = 15 atoms with g AA = 0.2.induced interaction and vice versa [cf.Figures 3(b) and 4(a)].For instance, considering repulsive g AB and tuning g AC to weak attractions, ∆〈r BC 〉 becomes positive denoting an induced repulsion between the impurities.However, for stronger repulsive g AC ∆〈r BC 〉 is negative and thus attractive induced interactions occur maximizing in the coalescence regime where g AB and g AC are both strongly repulsive, see also the inset of Figure 3(a).Furthermore, in the case of suppressed mediated correlations between the impurities (C BC ≈ 0), i.e., in the trivial case where g AB = 0 or for strong attractive g AC and repulsive g AB [cf. Figure 3(b)], also ∆〈r BC 〉 vanishes (see Figure 4(a) for strong attractive g AC and g AB = 0.2, 1.0).In the last scenario, the gradually increasing g AC attraction leads to a reduction (enhancement) of the correlation between the bath and the B (C) impurity whose interplay impedes the development of mediated impurity correlations and therefore induced interactions.
In the case of an attractively coupled impurity B, e.g.g AB = −1.0,−0.2, ∆〈r BC 〉 decreases when g AB is tuned to strong attractive values, a phenomenon also occurring for Here, increasing the attraction between impurity C and the bath enhances their correlation, while at sufficiently strong attractive g AC the correlation between the bath and the impurity B begins to slightly decrease for constant attractive g AB (cf. Figure 3).This competition between the different impurity-medium correlations suggests an interesting interplay between the individual intercomponent correlations and could in principle hinder the bath to mediate correlations between the impurities leading eventually to the observed reduction of the induced impurity-impurity correlation/interaction.Such an interplay of intercomponent correlations is indicative of a more intricate and generic correlation transfer process among the species [36], that is an exciting future perspective but lies beyond the focus of our study.However, note that for decreasing g AB = g AC results in a saturation of the impurity-impurity correlation, a fact that will also become important later in the discussion regarding the bipolaron formation in Section 7.
Finally, notice that a similar qualitative behavior of the intercomponent correlations and thus also of ∆〈r BC 〉 takes place for either increasing the number of atoms of the bosonic medium or the bare mass of one of the impurities, see Appendix D. In fact, both scenarios lead for repulsive g AB and g AC to an amplified impurities entanglement and to a stronger attractive induced interaction.

Effective two-body model
To determine the strength of induced impurity-impurity interactions, we reduce the threecomponent many-body system to an effective two-body model consisting of two interacting quasi-particles.This is a common approach to identify polaron properties from many-body simulations and has been successfully applied to two indistinguishable impurities [52] but not to distinguishable ones.Here, the effective two-body model employs the effective potential V eff σ (x σ ) [defined in Eq. ( 6)] for each impurity and thus neglects impurity-medium correlations.Also, the underlying impurities induced interactions are represented by a contact potential of strength g eff BC (a treatment with finite range interactions leads to similar results as it is demonstrated in Appendix C).Specifically, the corresponding effective two-body Hamiltonian reads The effective potential accounts for the effective mass and frequency of each impurity [94].
These effective parameters originate from the polaron picture where the impurity becomes dressed by the excitations of the bath, see Appendix B for a more detailed discussion.
In order to deduce the effective interaction strength g eff BC , we minimize ∆G (2) BC and G (2),eff BC are the impurities' two-body correlation functions calculated from the many-body three-component mixture and the effective two-body model, respectively. 3By estimating the value of g eff BC which minimizes ∆G (2) BC , we are able to associate the emergent induced correlation pattern between the impurities described in Fig. 3(b) with a corresponding induced interaction strength g eff BC .The resultant behavior of g eff BC provided in Figure 4(b) for fixed g AB and varying g AC agrees qualitatively with the observations made for ∆〈r BC 〉.The impurities experience an induced attraction when they both couple either attractively or repulsively to the bath, corresponding to an induced correlation, otherwise they feature an induced repulsion related to their anti-correlated tendency. 4To testify the validity range of the effective two-body model [Eq.(12)] for describing the impurities, we calculate the fidelity F BC of their ground state wave function as obtained from H (2),eff (|Φ BC eff 〉) and the full three-component mixture (| ΨBC i 〉). 5 The fidelity is provided in Figure 4(c) as a function of g AC and for different fixed values of g AB .It becomes apparent that H (2),eff is not valid for g AA < g Aσ where the respective impurity phase separates with the bath.We further note that especially in the regime where the impurities are anti-correlated and share no significant spatial overlap, an effective treatment considering a contact interaction potential fails to describe the full many-body calculations.Instead, in this interaction regime, due to the presence of non-local correlations, a more appropriate choice to model effective impurity-impurity interactions would be a long-range interaction potential, such as the one used in Appendix C. 3 We find ∆G (2) BC ≲ 10 −5 for all considered interaction strengths g AB and g AC . 4Note that g eff BC = 0 if one of the impurities does not interact with the bath which further confirms the validity of the effective model predictions since in this case no correlations are mediated. 5For this reason we use the Schmidt decomposition where the λ i correspond to the Schmidt coefficients [95,96].As such the fidelity is expressed as Still, within this effective two-body model different observables for the impurities such as their residue and correlation functions can be extracted and shown to exhibit a qualitative correct behavior.Deviations from the full many-body results are mostly traced back to the absence of intracomponent correlations of the bath and impurity-medium ones.

Bipolaron formation
Strong attractive induced interactions between two dressed impurities, commonly occurring for strong attractive impurity-medium direct interactions, eventually lead to the formation of a bound dimer quasi-particle state, the so-called bipolaron [32,45].In order to probe the presence of such a dimer impurity bound state in our setup, we study the bipolaron energy, Here, E(g AB , g AC ) denotes the total energy of the system including the two distinguishable impurities, E 0 is the energy of the bosonic gas in the absence of impurities and is the energy of one impurity coupled to the bath.The bipolaron energy is presented in Figure 5(a) covering a wide range of attractive and repulsive impurity-medium interactions, g AB and g AC .It features a rapid decrease when both impurities couple attractively to the medium, thereby, evincing the formation of a bound state. 6 A complementary observable used for the identification of the bipolaron is the spatial size of this dimer state.This is naturally captured by σ ∼ 〈r 2 BC 〉, where 〈r 2 BC 〉 is the squared relative distance [cf.Eq. (10)] between the impurities B and C [32].Specifically, in the following, we track σ/σ 0 with σ 0 being the distance in the uncoupled scenario, i.e., at g AB = g AC = 0, such that we explicitly estimate the impact of the impurity-medium interactions on the dimer size.This is depicted in Figure 5(a) as contour dashed lines along which σ/σ 0 is constant in the g AB -g AC plane on top of the bipolaron energy.It can be readily seen that for increasing magnitude of the attractive impurity-medium couplings, i.e., g AB and g AC , the size of the dimer state shrinks further, see in particular the dashed lines in Figure 5 which from bottom left to top right correspond to σ/σ 0 ≈ 0.18, 0.29, 0.65.
The bipolaron dimer state refers to the bunching behavior of the impurities at the same spatial region which manifests in the elongated shape of their two-body density ρ (2) 2 ) along the diagonal.In the non-interacting case, i.e., g AB = g AC = 0, ρ 2 plane and becomes gradually elongated for larger attractions due to the mediated attraction between the impurities, see e.g.Figures 5(b) and (c) for the cases (g AB , g AC ) = (−0.5,−0.5) and (−1.5, −1.5), respectively, also marked as gray dots in Figure 5(a).To quantify the degree of the aforementioned elongation, we fit the half maximum of the impurities' two-body density, 7 i.e. ρ (2) BC (0, 0)/2 to a rotated ellipse [see white dotted lines in Figures 5(b) and (c)] and determine the corresponding eccentricity e = 1 − b 2 /a 2 where a (b) denotes the semi-major (semi-minor) axis marked by the black lines of the ellipse. 8Apparently for e = 0, ρ 2 ) is circularly symmetric while in the case of e < 1 it is elongated having the shape of an ellipse.
The eccentricity of the impurities' two-body density is depicted in Figure 5(d) for g AB = g AC .By tuning the impurity-medium coupling from the non-interacting limit towards strong attrac- 6 The bipolaron energy decreases exponentially if both impurity-medium couplings (g AB , g AC ) are equally varied from the non-interacting limit to the strongly attractive regime, i.e., along the diagonal in Figure 5(a). 7We remark that choosing ρ (2) BC (0, 0)/2 for the fitting is employed for convenience.Indeed, also other density values were used, e.g.ρ (2) BC (0, 0)/4, verifying the same behavior of the eccentricity. 8For the fitting we use the general ellipse equation αx 2  1 + β x 1 x 2 + γx 2 2 + δx 1 + εx 2 + φ = 0, which in the frame of the ellipse reduces to x1 2 /a 2 + x2 2 /b 2 = 1.(2) 2 ) for (g AB , g AC ) = (−0.5,−0.5) and (−1.5, −1.5), respectively, in units of mω/ħ h [see also corresponding gray dots in panel (a)].The region where ρ (2) σσ ′ (0, 0)/2 is fitted to an ellipse (white dotted line) and shown together with the semi-minor and semi-major axis (black lines).The corresponding eccentricity is depicted in panel (d) assuming g AB = g AC .The transition to a bipolaron state where the eccentricity saturates for increasing impurity-medium attractions and the size of the dimer state is σ/σ 0 ≈ 0.29 occurs at g AC = −1.5 (gray dashed line).We consider two non-interacting but distinguishable impurities immersed in a bosonic gas of N A = 15 atoms with g AA = 0.2.
tions, e increases from e ≈ 0 at g AB = g AC = 0 to finite positive values until it saturates at around g AB ≈ −1.5.A larger attraction leads only to an additional shrinking of the dimer size, see in particular the exponential decrease of σ/σ 0 in Figure 5(d), leaving the shape of ρ 2 ) almost unchanged.In this sense, we deduce that the bipolaron state is formed at g AB = g AC ≈ −1.5 corresponding to σ/σ 0 ≈ 0.29 [vertical gray dashed line in Figure 5(d)].This observation allows us to generalize our conclusions for the bipolaron formation also in the case of g AB ̸ = g AC from the critical size of the dimer state being σ/σ 0 ≲ 0.29, which corresponds to the central contour dashed line in Figure 5(a).
We remark that the above-described behavior of both E bip (g AB , g AC ) and σ/σ 0 is in accordance with previously studied two-component systems containing two indistinguishable bosonic impurities that form a bipolaron 9 in the strongly attractive coupling regime [32].However, our results generalize these findings demonstrating the existence of a bipolaron in the case of two distinguishable impurities and suggesting that this bound state is robust to individual variations of g AB or g AC as indicated by the contour lines in Figure 5. Another aspect that we have addressed is that increasing the mass of one impurity, e.g.considering m B = 2, leads to a faster reduction of the dimer state size as well as the bipolaron energy for decreasing g AB = g AC while the eccentricity saturates at smaller impurity-medium attractions as compared to the mass-balanced case.This suggests, as expected, that a heavier impurity facilitates bipolaron formation.
3 -planes.The spatial coordinates x σ are expressed in units of ħ h/mω, whereas ρ ABC are presented in units of (mω/ħ h) 3/2 .For visualization purposes we only show the data whose correlation or density value is larger than 0.2 of the respective maximum value.The region corresponding to ρ ABC (0, 0, 0)/2 is fitted to an ellipsoid rotated in space (part of the fitted ellipsoid is marked by the white dashed lines).The three semi-axis are denoted by the green lines in panel (d).(e) Eccentricities calculated from the semiaxis (see main text) for attractive g AB = g AC .(f) Jacobi relative distances 〈r  16)] for g AB = g AC .We mark the transition to a trimer state at g AC = −1.5 [gray dashed line in panels (e) and (f)].For the three-component setup two non-interacting but distinguishable impurities immersed in a bosonic gas of N A = 15 atoms with g AA = 0.2 are considered.

Three-body correlations and trimer state
In the following, we aim to shed light on the existence of three-body correlations appearing in the ground state of the two distinguishable impurities embedded into the bosonic gas.For this purpose, we construct as a first step the normalized reduced three-body density which represents the spatially resolved probability of finding at the same time a representative atom of the medium at position x A 1 and the impurities B and C at positions x B 2 and x C 3 [97,98].Experimentally, the three-body density could be obtained by detecting simultaneously the positions of the three particles of interest, here, the two impurities and one bath atom, and then average over a sample of experimental absorption images [99].Having defined the three-body density, we construct the spatially resolved three-body correlation function as a straightforward extension of the two-body one defined in Eq. ( 7), i.e., According to this measure, the three participating particles are correlated (anti-correlated) if G (3) 3 ) = 0 implies that they are uncorrelated.Note, that this measure still contains two-body correlation effects since only the product of one-body densities has been subtracted from the three-body density.
The three-body correlation function is depicted in Figures 6(a) and (b) for the case of strong repulsions between impurity B and the bath (g AB = 1) and either weak attractive or repulsive couplings between the bath and the C impurity, namely g AC = −0.2 and 0.2, respectively.Moreover, for visualization and completeness issues, we additionally showcase within the 3 planes the underlying two-body correlation functions G (2) 3 ) and G (2) 3 ), respectively. 10Focusing on g AC = −0.2, it becomes evident that G (3) 3 ) fragments into two correlated and two anti-correlated parts.The correlated segments indicate that it is likely for one bath atom and the C impurity to reside at the same side with respect to the trap center while the repulsively coupled impurity B favors to be on the opposite side.On the other hand, the anti-correlated fragments suggest that a configuration where the impurities and a bath atom are at the same location is not favorable.The spatial arrangement of these fragments is altered in the three-dimensional space if the sign of g AC is inverted, in a sense that the correlated and anti-correlated regions are rotated by roughly 90 • around the x B 2 direction.In such a configuration the impurities are located at the same side in terms of the trap center and a bath atom lies on the opposite side.The corresponding two-body correlation functions G (2) 2 ) preserves its pattern, see the contours in Figures 6(a) and (b).Subsequently, we turn to strongly attractive impurity-medium interactions with g AB = g AC .Here, the three-body density ρ 3 ) becomes elongated exhibiting an ellipsoidal shape, see e.g. Figure 6(d) for (g AB , g AC ) = (−1.5,−1.5).Thereby, the three-body density is stretched along the (x A 1 , x B 2 , x C 3 )-direction, i.e., the diagonal of the coordinate system, demonstrating a bunching behavior of the two impurities and a representative atom of the bath species.In particular, the corresponding three-body correlation function, presented in Figure 6(c), features a correlated pattern along the diagonal around which a shell-like structure consisting of anti-correlated fragments is formed.
To quantify the deformation of the three-body density, we fit its half maximum, i.e., ρ ABC (0, 0, 0)/2, to a rotated ellipsoid (see white dashed lines in Figure 6(d) corresponding to a profile of the ellipsoid).Specifically, we fit the ellipsoid equation x2 , where xi refers to the coordinate system of the ellipsoid spanned by its semi-axis with lengths a, b and c [green lines in Figure 6(d)].From the semi-axis we determine three eccentricities, namely e a b = 1 − b 2 /a 2 , e ac = 1 − c 2 /a 2 and e bc = 1 − c 2 /b 2 with a ≥ b ≥ c.These eccentricities are depicted in Figure 6(e) together with the relative deviation, er r, from the ellipsoid function for varying g AB and assuming g AB = g AC .In the non-interacting case, i.e., g AB = g AC = 0, the eccentricities are already finite indicating a deviation from a spherical shape, which is in contrast to the bipolaron [cf. Figure 5(d)].This is attributed to the presence of finite intraspecies interactions among the bath particles causing the observed spatial deformation.Importantly, the eccentricities show an increasing tendency for stronger attractive values of g AB = g AC , meaning that the elongation of the ellipsoid is enhanced until it saturates at around g AB = g AC ≈ −1.5. 10 As an example, notice that the contours in the x A 1 -x B 2 and x B 2 -x C 3 planes of Figure 6(c) correspond to the

Conclusions and perspectives
We have studied the correlation properties in the ground state of two non-interacting distinguishable impurities immersed in a bosonic bath with the entire three-component system being harmonically trapped.The impurities become dressed by the excitations of the bosonic gas generating quasiparticle states, herein Bose polarons, having characteristic properties such as effective mass and featuring induced correlations.In order to appreciate the impact of interand intracomponent correlations we rely on the variational ML-MCTDHX method whose flexible wave function truncation ansatz allows to operate at different correlation orders.An emphasis is placed on the high tunability of the three-component setting unveiling rich density and correlation patterns, the manipulation of both the sign and the strength of impurities induced interactions as well as the formation of bound impurity states.Specifically, we demonstrate that upon varying the involved impurity-medium couplings, both impurities can either localize at the trap center (attractive intercomponent interactions), form a shell around the bosonic gas (repulsive interactions), i.e., phase-separate, or one of them localize and the other phase-separate (alternating signs of impurity-medium couplings).These density configurations can be understood at least qualitatively in terms of an effective potential picture for the impurities which refers to a dipped harmonic oscillator (double-well) for attractive (repulsive) intercomponent interactions.
A detailed characterization of the induced correlations is provided in a wide range of impurity-medium interactions aiming to expose their intricate role.Inspecting the two-body intercomponent correlation functions we find that the bosonic gas mediates anti-correlations among the impurities if one of them couples repulsively and the other attractively to it.In contrast, induced two-body correlations occur as long as both impurities couple either attractively or repulsively to their medium.The origin of the aforementioned correlation patterns is traced back to the spatial configurations of each component.This means that if the impurities have a finite spatial overlap with the bath the latter mediates two-body correlations between them.Interestingly, there is also the possibility that the impurities are not overlapping but can be still correlated implying that non-local correlations are in play.To quantify the strength and sign of the induced interactions we employ the relative two-body distance among the impurities extracting all contributions stemming from mean-field effects.In this sense, it is demonstrated that induced two-body correlations (anti-correlations) are related to mediated attractive (repulsive) impurity interactions.These findings are further supported by an effective two-body model containing the impurities effective trapping potential and their induced interactions.Importantly, this approach allows to determine the strength and sign of the effective interactions mediated between the impurities through a comparison with the full many-body results.Moreover, by constructing an effective one-body Hamiltonian enables us to estimate the effective mass and trapping frequency of each distinguishable impurity (polaron), see Appendix B.
Evidences regarding bipolaron formation are provided, when both impurities are strongly attractively coupled to the bosonic gas, by means that the bipolaron energy and the size of the underlying dimer state rapidly decrease for stronger attractions.Interestingly, we determine the intercomponent three-body correlation function according to which overall weak threebody correlations exist and become enhanced for strongly attractive impurity-medium interactions signaling the formation of trimers among the impurities and an atom of the medium.
In this investigation we have restricted ourselves to the ground state of the threecomponent mixture.Further understanding on the character of the impurities induced interactions and in particular their nonlocal character and their dependence on the statistics of the medium are interesting perspectives.In this context, a systematic finite size scaling analysis with respect to the number of bath particles in order to infer the persistence of our findings e.g. in terms of the crossover of the impurities induced interactions (see also Appendix D) and in general the build-up of intercomponent correlations would be desirable as well.Also, the emulation of spectroscopic schemes that will allow the identification of the ensuing polaron states and excitations [24,102] constitutes an intriguing direction.Furthermore, studying the behavior of impurities induced interactions and bound states in different external trapping potentials is also an interesting direction.Here, a setup of immediate interest would be to load the bath atoms in a ring potential and investigate the formation of impurities bound states in both the attractive and the repulsive impurities-medium interaction regimes.Another straightforward extension would be to explore the nonequilibrium impurities dynamics in order to understand the build-up of induced correlations.An additional fruitful research direction is to understand the Bose polaron formation when indistinguishable impurities are immersed in an attractive two-component gas forming a droplet.Certainly, studying correlation effects in particle-balanced three component settings with an emphasis on the few-to many-body crossover and in particular close to the pair immiscibility threshold is worth to be pursued.

A Behavior of the bipartite entanglement
A standard measure to estimate the bipartite entanglement of mixed states that exist in a multi-component system 11 is encapsulated in the logarithmic negativity [56][57][58][59]62,104,105].
It is based on the partial transpose of the two-body species reduced density matrix, 12 which, e.g.referring to species σ and σ ′ , is obtained by integrating out the degrees of freedom of species σ ′′ leading to ρ 55,77].
Its partial transpose T σ with respect to species σ is calculated by exchanging the indices i and l associated with species σ, i.e., ρ T σ and in particular summing up its negative eigenvalues µ i yields the so-called negativity, N σσ ′ = i |µ i |.Subsequently, the logarithmic negativity reads This measure exploits the fact that for a separable mixture, e.g.ρ (2),spec σσ ′ = i p i ρ(1),spec , the partial transpose does not alter the spectrum of ρ (2),spec σσ ′ and, hence, all eigenvalues remain positive.In this sense, the presence of negative eigenvalues guarantees the existence of entanglement.However, this statement can not be inverted, i.e., even if the logarithmic negativity is zero the species σ and σ ′ can still be entangled [103].
The logarithmic negativity between the bath and the B impurity, E AB , as well as among the impurities, E BC , is illustrated in Figures 7(a) and (b) respectively within the g AB -g AC plane.As expected it overall captures the main features of the integrated correlation functions shown in Figures 3(a) and (b).For instance, E AB vanishes for strongly attractive g AC and strongly repulsive g AB [Figure 7(a)], while the parameter region referring to the impurities coalescence is in a similar way pronounced in E BC as it has been observed for C BC , compare Figures 3(a) and (b) for repulsive g AB and g AC .Recall that while E σσ ′ provides only a quantitative diagnostic for the bipartite entanglement and does not describe the correlated or anti-correlated behavior as C σσ ′ it still gives insight into the entanglement content of the many-body system.As such, for large g AB < 0 the logarithmic negativity uncovers that the bath and the B impurity are strongly MB denotes the one-body distribution of the full threecomponent many-body system, whereas ρ entangled especially so in the repulsive g AC > 0 region, while varying g AB towards the weakly attractive regime and for |g AC | > 1 entanglement is reduced [Figure 7(a)].This is attributed to the simultaneous increase of E AC , 13 unveiling a competition between the intercomponent entanglement of individual impurities with the medium.Finally, in line with the predictions of C BC , E BC demonstrates that entanglement is finite when both impurities are either weakly attractive or strongly repulsively coupled to the medium, see Figure 3(d).

B Effective mass and trap frequency of a single impurity
In the following, we approach the three-component impurity setting as a polaron problem since each individual impurity via its coupling to the bosonic gas is dressed by the excitations of the latter.In this sense, we aim to capture the effective behavior of the B and C impurity with the effective one-body model [94], where m eff σ and ω eff σ denote the polaron effective mass and trapping frequency with σ ∈{B, C}. 14  To identify the values of the effective mass and frequency, we minimize the cost function L σ = ∆ρ (1)  σ In this expression, the first term refers to ∆ρ (1)  σ = dx σ ρ (1),MB for the characteristic interaction configurations (g AB , g AC ) = (−0.2,0.1) and (0.2, 0.1), respectively.For comparison we additionally provide the one-body density ρ As it can be readily seen, the one-body densities predicted by the two effective one-body models are in excellent agreement with the one corresponding to the full three-component many-body system.Deviations start to become evident for strong repulsive impurity-medium couplings (not shown) where the impurity and the medium phase separate [36,94].Recall that the effective model is by definition valid for weak intercomponent repulsions where the impurity does not probe the edges of the bosonic cloud.
The effective masses and frequencies of the B and C impurities after minimization of the cost function given by Eq. (B.2) are represented in Figures 8(c) and (d) with respect to the impurity-medium couplings.It is important to point out that both the effective mass and frequency of a specific impurity, e.g. the B one, primarily depend on its coupling with the bath g AB .The interaction strength of the other impurity (C) with the bath, e.g.g AC , has almost no impact on the effective parameters of impurity B. For instance, this conclusion can be drawn from the nearly constant behavior of m eff B and ω eff B for varying g AC shown in Figure 8(c), or the fact that m eff C and ω eff C remain almost intact for fixed g AC and different g AB , see Figure 8(d).For an attractively coupled impurity with the bosonic gas, the effective mass and frequency become larger than their bare values [gray dashed lines in Figures 8(c  in Figure 8(a)].On the other hand, in the case of a repulsively coupled impurity the effective trapping frequency is still tighter than the original value, but the effective mass becomes smaller than its bare value [cf.m eff B , ω eff B for g AB = 0.2 in Figure 8(c) as well as m eff C , ω eff C for g AC > 0 in Figure 8(d)].In particular, the effective mass is small enough to compensate the increased effective frequency meaning that the underlying harmonic trap is eventually broadened [cf.V ho−eff B in Figure 8(b)].Additionally, the comparatively smaller effective mass is related to 14 Recall that within the effective two-body model described by Eq. ( 12) we implicitly account for the effective mass and frequency via the effective potential V eff B,C Eq. ( 6).Indeed, beyond mean-field corrections imprinted on ρ

BC
, which correspond to the overlap of the impurities two-body wave function obtained from the full many-body approach and the effective model of Eq. ( 12) containing either an exponential or a contact-type interaction potential, respectively.(b) Difference between ∆G (2),exp BC and ∆G (2),contact BC , referring to the variance of the two-body correlation function calculated within the effective two-body model using either the contact or the exponential interaction potential with respect to the full three-component system.For both quantities the relative deviations are minor, testifying the validity of both effective interaction potentials.a spatial delocalization of the impurity cloud. 15In this way, the effective one-body model captures the effects imprinted on the impurity in the three-component system.

C Modelling the effective impurity interactions with an exponential potential
To verify the validity of the contact interaction potential for describing the induced interactions between the impurities [Eq.( 12)], we next exemplify that our results do not change if one instead uses an exponential potential.The latter has been derived in Refs.[38,39] and holds in the homogeneous case and for immobile impurities residing at distances satisfying A (0) ≈ 0.6 being the healing length of the bath.In particular, we replace the interaction term in Eq. (12) with A (0) . 16 As discussed in Section 6.2, we judge the quality of the effective twobody model by estimating the fidelity, F BC , between the impurities two-body wave function as extracted from the full many-body system and the effective two-body model containing either a contact or an exponential interaction potential.Subsequently, we determine the difference F exp BC − F contact BC which as shown in Figure 9(a) testifies deviations at most of the order 10 −4 .Proceeding one step further, we determine the overlap between the respective two-body correlation functions of the impurities determined within the full three-component system and the effective two-body model.Namely, we track ∆G where G (2),exp BC denotes the two-body correlation function obtained within the effective twobody model (see also Section 6.2) with an exponential interaction.To infer the deviations among the exponential and contact effective interactions at the two-body correlation level, we calculate the difference ∆G Therefore, the contact and exponential effective interaction potentials lead essentially to the same description regarding the impurities properties.This outcome was not a-priori expected since the exponential potential is originally derived in the homogeneous case.

D Impact of mass-imbalanced impurities and the atom number of the bosonic gas
Let us demonstrate the generalization of our results in the main text when the impurities are mass-imbalanced or the bosonic medium contains a larger number of particles.For this purpose, we focus on the behavior of the intercomponent correlations which can be quantified through the integrated correlation function [Eq.( 9)] presented in Figure 10 for different system parameters.
In general, increasing the mass of an impurity disturbs the cloud of the bosonic gas to a larger degree which should eventually lead to an enhanced impurity-medium correlation.This is indeed evident in Figure 10(a) where the integrated correlation function, C AB , is increased as compared to the mass-balanced case, thus testifying an overall larger degree of entanglement.Furthermore, since the correlation between the C impurity and the bath is not affected by the change of m B [Figure 10(b)], the larger C AB leads to a stronger mediated correlation between the impurities, see e.g.C BC in Figure 10(c).The latter naturally leads to an amplified impurities' induced interaction for increasing m B .In particular, for g AB = 0.2 and strong repulsive g AC , where C BC features the largest increase.
Next, we concentrate on the mass-balanced system but consider a larger number of bath particles and in particular N A = 30, while maintaining the same mean-field interaction, i.e., N A g AA = const.As it can be seen, the impurity-medium correlations, as captured by C AB and C AC , are reduced compared to the reference case N A = 15, g AA = 0.2 [Figure 10(a), (b)].This is attributed to the smaller intra-species coupling strength g AA = 0.1 resulting in a decrease of the respective intra-species correlations among the bath particles.However, the mediated correlations among the impurities B and C are clearly enhanced when g AB and g AC are both repulsive, see Figure 10(c).In this sense, a larger number of bath particles featuring a decreasing intraspecies interaction is associated to a reduction of intraspecies correlations of the bath and impurity-medium ones but enhances to a certain degree the mediated correlation between the impurities.This behavior hints towards a complicated correlation transfer mechanism to the impurity-impurity subsystem which deserves further future investigations.Nevertheless, a systematic finite size scaling analysis in terms of the atom number in the bath is required in order to deduce the robustness of our findings.However, we expect that the main features of the impurities dressing, e.g. the crossover from a correlated to an anti-correlated behavior (associated to attractive and repulsive induced interactions as discussed in Sections 5 and 6), and the existence of the impurities bound states for attractive interactions are retained for larger number of bath atoms.

E Estimating the importance of correlations on the many-body wave function
To expose the impact of intercomponent correlations at different interaction regimes on the level of the many-body wave function we analyze the fidelity 〈Ψ sMF |Ψ MB 〉 2 , see Figure 11(a).
Here, |Ψ MB 〉 denotes the full many-body wave function where all emergent inter-and intracomponent correlations are taken into account, while |Ψ sMF 〉 refers to the species mean-field wave function which ignores all intercomponent correlations.Naturally, the fidelity is unity when the species are non-interacting, i.e., g AB = g AC = 0, since in this scenario intercomponent correlations are a-priori prohibited.However, the fidelity decays for increasing impurity-medium coupling strengths as intercomponent correlations are triggered in this case.The largest deviation between the many-body and species mean-field wave functions occurs in the parameter region corresponding to the coalescence of the impurities, i.e., for strongly repulsive g AB and g AC .Further understanding of the respective correlation mechanisms can be delivered by identifying the participating microscopic configurations.For this reason we construct the species function eigenbasis |ψ A i 〉|ψ B j 〉|ψ C k 〉 obtained by calculating the eigenfunctions of an effective k 〉, constructed from the eigenstates of an effective species Hamiltonian (see main text), with the many-body wave function |Ψ MB 〉.Apparently, energetically higher-lying excited states possess substantial contribution.Probability amplitudes which remain below 0.02 within the interaction range −2.0 ≤ g AC ≤ 2.0 are shown as gray lines.The harmonically trapped three component system consists of two non-interacting but distinguishable impurities immersed in a bosonic gas of N A = 15 atoms with g AA = 0.2.
species Hamiltonian [cf.Eq. ( 2)] characterized by the effective potential defined in Eq. ( 6). 17s basis for the bath we take the ground and the energetically two lowest excited states of the effective potential into account, while for the two impurities we consider the corresponding energetically lowest six eigenstates leading to a total number of 108 three-component basis states |ψ A i 〉|ψ B j 〉|ψ C k 〉. g AB = 1.0).In general, it is observed that finite interactions yield a non-negligible population of energetically higher-lying excited states.Importantly, this behavior becomes enhanced in the coalescence regime, i.e., for strong repulsive g AB and g AC .This means that there are several macroscopically occupied basis states reflecting the significant intercompoment entanglement (cf.Figures 2 and 7).

) −2 0 2 gFigure 3 :
Figure 3: (a)-(b) Phase diagram of the intercomponent (see legends) spatially integrated correlation functions C σσ ′ [Eq.(9)] in the parametric plane of the impuritymedium interaction strengths (g AB , g AC ).A value of C σσ ′ < 0 (C σσ ′ > 0) indicates an anti-correlated (correlated) behavior between the atoms of species σ and σ ′ , while C σσ ′ = 0 denotes the absence of two-body correlations (see also main text).The gray circles correspond to the interaction combinations (g AB , g AC ) depicted in Figures1 and 2. The regions enclosed by the dashed lines in panel (b) indicate the interaction regions where the impurities do not overlap but are still two-body anti-correlated.The harmonically trapped three component system consists of two non-interacting but distinguishable impurities immersed in a bosonic gas of N A = 15 atoms with g AA = 0.2.

gFigure 4 :
Figure 4: (a) and its inset: Modified relative distance [Eq.(11)] reflecting the effects on 〈r MB BC 〉 which are exclusively caused by the induced impurities correlation as a function of g AC and for different fixed g AB .(b) Induced interaction strength between the two Bose polarons estimated by maximizing the overlap between the two-body correlation functions G (2),eff BC obtained from the effective two-body model

Figure 5 :
Figure 5: (a) Bipolaron energy, E bip , as a function of the intercomponent coupling strengths g AB and g AC .The dashed lines represent contours along which the size of the dimer state σ remains fixed and in particular from bottom left to top right correspond to σ/σ 0 ≈ 0.18, 0.29, 0.65.(b), (c) Reduced two-body impurities' density ρ

Figure 7 :
Figure 7: (a)-(b) Diagram of the intercomponent (see legends) logarithmic negativity E σσ ′ [Eq.(A.1)] as a function of the impurity-medium couplings (g AB , g AC ).The harmonically trapped three component system consists of two non-interacting but distinguishable impurities immersed in a bosonic gas of N = 15 atoms with g AA = 0.2.

( 1 )
,ho−eff B are calculated using the effective one-body Hamiltonians composed of either an effective harmonic oscillator with an effective mass and frequency [cf.Eq. (B.1)] or the effective potential defined in Eq. (6), respectively.Effective mass and trapping frequency of the dressed (b) B and (c) C impurity, respectively, as deduced from the effective polaron model defined of Eq. (B.1).

σ
(x σ ) − ρ(1),ho−eff σ (x σ )2 with ρ(1),MB σ and ρ(1),ho−eff σ being the one-body density as predicted from the full three-component system and the effective one-body model, respectively.The second contribution of the righthand side in Eq. (B.2) designates the energy difference∆E σ = E MB σ − E ho−eff σ 2, whereE MB σ = 〈Ψ MB | Ĥσ |Ψ MB 〉 is the σ impurity energy and E ho−eff σ = 〈φ| Ĥ(1),ho−eff σ |φ〉 = 12 ω eff σ is the energy of the effective one-body model and |φ〉 the corresponding ground state.Note that in order to uniquely estimate m eff σ and ω eff σ one needs to adequately describe both the density and the energy of the impurity.Figures8(a) and (b) showcase the one-body densities ρ (1),MB B and ρ (1)ho−eff B ) and (d)], see in particular m eff B , ω eff B when g AB = −0.2 in Figure 8(c) and m eff C , ω eff C for g AC < 0 in Figure 8(d).As such, the emergent Bose polaron experiences a narrower trapping potential, thereby, reflecting the localization of the impurity at the trap center [cf.ρ (1),ho−eff B and V ho−eff

( 1 )Figure 9 :
Figure 9: (a) Relative deviation between the fidelities F exp BC and contact BC, see Figure9(b).Also here, only small deviations of the order 10 −5 are identified.

2 Figure 10 :
Figure 10: Integrated two-body correlation function [Eq.(9)] among (a) the B purity and the medium, (b) the C impurity and the medium and (c) between the impurities as a function of the intercomponent interaction strength g AC .In all panels, we consider fixed g AB = −0.2,0.2 as well different masses of the B impurity (simultaneously setting ω B = m A /m B ) and atom numbers of the medium (see legend), while keeping constant the mean-field interaction N A g AA .The gray dashed line in panel (c) marks ∆ 〈r BC 〉 = 0.

Figure 11 :
Figure 11: (a)  Fidelity between the many-body wave function, |Ψ MB 〉 (including all emerging intra-and intercomponent correlations) and the species mean-field wave function |Ψ sMF 〉 where intercomponent correlations are neglected.The reduction of the overlap from unity for finite interactions evinces the participation of intercomponent correlations.(b) Probability amplitude P i jk denoting the overlap of a threecomponent time-independent basis |ψ A i 〉|ψ B j 〉|ψ C k 〉, constructed from the eigenstates of an effective species Hamiltonian (see main text), with the many-body wave function |Ψ MB 〉.Apparently, energetically higher-lying excited states possess substantial contribution.Probability amplitudes which remain below 0.02 within the interaction range −2.0 ≤ g AC ≤ 2.0 are shown as gray lines.The harmonically trapped three component system consists of two non-interacting but distinguishable impurities immersed in a bosonic gas of N A = 15 atoms with g AA = 0.2.

2 ,
A i |〈ψ B j |〈ψ C k | Ψ MB 〉 with |Ψ MB 〉 being the full many-body wave function, are presented in Figure 11(b) for g AB = 1.0 and varying g AC .Notice that the state |ψ A 0 〉|ψ B 0 〉|ψ C 0 〉, denoting the case in which each species occupies the ground state of the effective species Hamiltonian, represents the three-body ground state obtained with a sMF ansatz.Consequently, P 000 = 〈Ψ sMF |Ψ MB 〉 2 (cf.Figures 11(a) and (b) for