Doping a lattice-trapped bosonic species with impurities: From ground state properties to correlated tunneling dynamics

We investigate the ground state properties and the nonequilibrium dynamics of a lattice trapped bosonic mixture consisting of an impurity species and a finite-sized medium. For the case of one as well as two impurities we observe that, depending on the lattice depth and the interspecies interaction strength, a transition from a strongly delocalized to a localized impurity distribution occurs. In the latter regime the two species phase separate, thereby forming a particle-hole pair. For two impurities we find that below a critical lattice depth they are delocalized among two neighboring outer lattice wells and are two-body correlated. This transition is characterized by a crossover from strong to a suppressed interspecies entanglement for increasing impurity-medium repulsion. Turning to the dynamical response of the mixture, upon quenching the interspecies repulsion to smaller values, we reveal that the predominant tunneling process for a single impurity corresponds to that of a particle-hole pair, whose dynamical stability depends strongly on the quench amplitude. During the time-evolution a significant increase of the interspecies entanglement is observed, caused by the build-up of a superposition of states and thus possesses a many-body nature. In the case of two bosonic impurities the particle-hole pair process becomes unstable in the course of the dynamics with the impurities aggregating in adjacent lattice sites while being strongly correlated.

We investigate the ground state properties and the nonequilibrium dynamics of a lattice trapped bosonic mixture consisting of an impurity species and a finite-sized medium. For the case of one as well as two impurities we observe that, depending on the lattice depth and the interspecies interaction strength, a transition from a strongly delocalized to a localized impurity distribution occurs. In the latter regime the two species phase separate, thereby forming a particle-hole pair. For two impurities we find that below a critical lattice depth they are delocalized among two neighboring outer lattice wells and are two-body correlated. This transition is characterized by a crossover from strong to a suppressed interspecies entanglement for increasing impurity-medium repulsion. Turning to the dynamical response of the mixture, upon quenching the interspecies repulsion to smaller values, we reveal that the predominant tunneling process for a single impurity corresponds to that of a particle-hole pair, whose dynamical stability depends strongly on the quench amplitude. During the time-evolution a significant increase of the interspecies entanglement is observed, caused by the build-up of a superposition of states and thus possesses a many-body nature. In the case of two bosonic impurities the particle-hole pair process becomes unstable in the course of the dynamics with the impurities aggregating in adjacent lattice sites while being strongly correlated.

I. INTRODUCTION
Ultracold atomic physics offers an excellent testbed for probing the static properties and in particular the nonequilibrium quantum dynamics of multicomponent systems for both fermions and bosons [1,2]. It provides an exquisite level of control of several system parameters including, for instance, the intra-and intercomponent scattering lengths via Feshbach resonances [1,3], the shape of the external trapping potential [4,5] as well as the particle number with remarkable experimental achievements especially in one spatial dimension [6,7].
However, the majority of these lattice trapped impurity investigations have been mainly focusing on the case that only the impurities experience the lattice potential and the host resides in a homogeneous environment. Moreover, they have been predominantly restricted to the single impurity case [53] and operated within the lowestband approximation [54,55]. Thus, the situation where both the impurities and their medium are trapped in the same lattice potential remains an open question. In such a setting the impurities act as defects possessing a particle character and it would be intriguing to study the different phases that arise in the ground state of this composite system for variable impurity-medium interactions and unveil their underlying correlation properties. Recall that lattice trapped particle-balanced bosonic mixtures exhibit quantum phases [56][57][58][59][60] being absent in their single-component counterpart. This is, partly, caused by the non-negligible presence of interspecies correlations [61]. For instance, modifications of the Mott-insulator (MI) to the superfluid (SF) phase transition [56][57][58][59] have been reported due to the existence of a second component leading to the so-called paired and counterflow superfluid states [62,63], quantum emulsion states [57,64] as well as losses of the intracomponent coherence [65]. In this sense, it is natural to investigate the existence and interplay of the different phases in the particle imbalanced scenario with respect to the interspecies interaction strength. As a prototypical example, henceforth we consider one or two bosonic impurities immersed in a majority species of bosons with both components being lattice trapped in one-dimension.
Having established the ground state properties of this setup another fruitful prospect is to inspect its corresponding nonequilibrium dynamics by quenching the system between the different emergent phases. Here the analysis and consequent control of the tunneling dynamics of the impurities is of particular importance since it might give rise to a variety of complex transport phenomena, self-trapping events and formation of (repulsively) bound pairs [48-51, 66, 67]. Furthermore, the identification of the correlated many-body nature of the different tunneling processes will allow us to infer their microscopic origin which is certainly of interest. To track the static properties and the quench dynamics of the particle imbalanced Bose-Bose mixture we utilize the variational Multi-Layer Multi-Configuration Time-Dependent Hartree method for atomic mixtures (ML-MCTDHX) [68][69][70] which enables us to capture all the relevant interparticle correlations of this multicomponent setup.
Regarding the ground state of a single and two bosonic impurities immersed in a majority bosonic species we find a transition from a SF to a MI phase of the composite system (doped insulator) for a specific lattice depth and increasing interspecies repulsion. This transition takes place for weaker impurity-medium interactions for deeper lattices, a result which is more pronounced in the two impurity case. Within the SF phase the impurity and the majority species show a delocalized behavior with the interspecies entanglement being enhanced and the medium being characterized by strong two-body correlations. However entering the MI state of the mixture, the species phase separate forming a particle-hole pair and their entanglement is suppressed [48,71]. In this latter case the many-body state of the system exhibits a twofold degeneracy [73]. Moreover, for two impurities we observe that below a critical lattice depth the impurities are delocalized among two neighboring outer wells of the lattice and are two-body correlated.
Turning to the dynamical response of the mixture, upon quenching the interspecies repulsion from a MI to a SF phase, we reveal that the predominant tunneling process for a single impurity corresponds to that of a particle-hole pair [72], whose dynamical stability depends strongly on the quench amplitude. More specifically, the initially localized impurity becomes spatially delocalized in the course of the evolution while it gradually tunnels from one side of the lattice to the other. On the other hand, the majority species particles tend to avoid the impurity in the course of the tunneling. During the time-evolution a significant increase of the interspecies entanglement is observed, which is due to the build-up of superposition of states and thus possesses a many-body nature. Additionally, strong correlations occur between the particles of the majority species. In the case of two bosonic impurities the particle-hole pair process becomes unstable during the evolution with the impurities aggre-gating in adjacent lattice sites while being strongly correlated.
This work is structured as follows. In sec. II we introduce our setup and discuss the variational many-body approach. Sec. III presents the ground state properties in a finite lattice for a single and two impurities immersed in a strongly interacting majority bosonic species with filling smaller than unity. The nonequilibrium dynamics of the impurities by quenching the interspecies interaction strength from the doped insulator to the SF phase is analyzed in sec. IV. We summarize our results and provide an outlook in sec. V. Appendix A elaborates on the lattice trapped ground state phase diagram of an impurity in a unit filling majority species.

A. Treatment of Many-Body Correlations and Dynamics
Our computation approach is the ab-initio Multi-Layer Multi-Configuration Time-Dependent Hartree method for bosonic (fermionic) Mixtures (ML-MCTDHX) [68][69][70], which accounts for all the relevant correlations of the atomic mixture [50,71,[74][75][76][77]. As a first step, the total many-body wave function |Ψ MB (t) is expanded in M species functions |Ψ σ (t) of species σ and written as a Schmidt decomposition [78] |Ψ Here, the Schmidt coefficients √ λ i , in decreasing order, provide information about the degree of population of the i−th species function and thereby determine the degree of entanglement between the impurities and the majority species. In case that λ 1 = 1 the species A and B are not entangled and the system can be described with a species mean-field ansatz corresponding to a single product state (M = 1).
Furthermore, the species wave functions |Ψ σ (t) describing an ensemble of N σ bosons are expanded in a set of permanents where the vector n σ = (n σ 1 , n σ 2 , ...) denotes the occupations of the time-dependent single-particle functions of the species σ. The notation n σ |N σ indicates that for each | n σ ; t we require the condition i n σ i = N σ . The time propagation of the many-body wave function is achieved by employing the Dirac-Frenkel variation principle δΨ MB | (i∂ t − H) |Ψ MB [79][80][81] with the variation δΨ MB . ML-MCTDHX provides access to the complete many-body wave function which allows us consequently to derive all relevant characteristics of the underlying system. In particular, this means that we are able to characterize the system by projecting onto number states with respect to an appropriate single-particle basis [82,83]. Besides investigating the quantum dynamics it allows us to determine the ground (or excited) states by using either imaginary time propagation or improved relaxation [84], thereby being able to uncover also possible degeneracies of the many-body states. We remark that in commonly used approaches for solving the timedependent Schrödinger equation, one typically constructs the wave function as a superposition of time-independent Fock states with time-dependent coefficients. Instead, it is important to note that the ML-MCTDHX approach considers a co-moving time-dependent basis on different layers, meaning that in addition to time-dependent coefficients the single particle functions spanning the number states are also time-dependent. This leads to a significantly smaller number of basis states and configurations that are needed to obtain an accurate description of the system under consideration and thus reduces the computation time [85].
The degree of truncation of the underlying Hilbert space is given by the orbital configuration C = (M, d A , d B ). Here, M refers to the number of species functions in the Schmidt decomposition (cf. equation 1), while d σ with σ ∈ {A, B} denote the number of singleparticle functions spanning the time-dependent number states | n σ ; t (cf. equation 2). The orbital configuration C = (7, 7, 7) has been employed for all many-body calculations presented in the main text, yielding a converged behavior of our observables.

B. Lattice Trapped Bosonic Mixture
Our system consists of a mixture of two bosonic species which are trapped in a one-dimensional lattice with hard wall boundary conditions at its endpoints. The impurity species with N B = 1, 2 particles is denoted as species B and the majority species containing N A = 4, 5 particles is referred to as species A. This setup lies within reach of current experimental techniques [86,87]. Furthermore, we introduce a coupling HamiltonianĤ AB between the two species. Both subsystems are confined along the longitudinal spatial direction, accounting for the one-dimensional character of our setup. Excitations in the corresponding transversal direction are energetically suppressed in the scenario under investigation and can therefore be neglected in our setup. This results in a Hamiltonian of the formĤ =Ĥ A +Ĥ B +Ĥ AB . The Hamiltonian of the species σ, with σ ∈ {A, B}, readŝ whereΨ † σ is the field operator of species σ, m σ their mass and V 0 the lattice depth. Also, g σσ refers to the intraspecies interaction strength of the two-body contact interaction among the σ atoms, k is the number of lattice wells and L is the length of the system, while x ∈ [−L/2, L/2]. Moreover, we assume equal masses for the species m A = m B . Experimentally this can be achieved by preparing e.g. 87 Rb atoms in two different hyperfine states [88]. The interaction between the species A and B is given bŷ where g AB is the effective one-dimensional interspecies interaction strength. The interaction strengths g α (α ∈ {A, B, AB}) can be expressed in terms of the three dimensional s-wave scattering lengths a 3D α . By assuming the above-mentioned strong transversal confinement with the same trapping frequencies ω σ ⊥ = ω ⊥ for both species σ ∈ {A, B} it is possible to integrate out frozen degrees of freedom, leading to a quasi one-dimensional model with g α = 2hω ⊥ a 3D α . Throughout this work we consider a k = 5 well lattice, while the interaction among the majority atoms is set to a value where the particles distribute in a Mott-like state for large lattice depths, namely g AA /E R λ = 0.04. Here, E R = (2πh) 2 /2m A λ 2 is the recoil energy and λ = 2L/k the optical lattice wavelength.

III. GROUND STATE PROPERTIES
Let us analyze the ground state properties of the Bose-Bose mixture with respect to the lattice depth V 0 and the interspecies coupling strength g AB for N B = 1 and N B = 2 impurities. We calculate the many-body ground state of the Bose-Bose mixture using ML-MCTDHX, which enables us to obtain the resulting full many-body wave function. In order to be able to interpret the wave function, we shall analyze reduced quantities such as the   von Neumann entropy and the one-and two-body densities of each species on the basis of the numerically obtained many-body wave function. As a result, we are able to gain an in-depth insight into the spatial distribution of the two species in the lattice potential and the accompanying intra-and inter-species correlations.

A. Single Impurity
In the following, we explore the ground state of the system containing a single impurity, i.e. N B = 1, for varying V 0 and g AB , while keeping fixed the intraspecies interaction strength to g AA /E R λ = 0.04. As a first step, we analyze the spatial distribution of the two species in terms of the one-body density of the ground state |Ψ MB of the species σ, which is defined as Additionally, in order to deepen our understanding of the ground state of the binary mixture, we investigate the degree of correlations and, in particular the entanglement between the impurity species and the majority species. For this purpose, we introduce the von Neumann entropy as a measure for the entanglement between the subsystems A and B, where λ i are the Schmidt coefficients de-fined in equation 1. Recall that in the case of a single contributing product state in equation 1, the subsystems are disentangled and the von Neumann entropy is S AB = 0, whereas any deviation from this value indicates entanglement between the A and the B species. Figure 2 (a) shows the von Neumann entropy of the ground state as a function of the interspecies interaction strength g AB and the lattice depth V 0 . For small g AB we observe that S AB ≈ 0, which indicates that our system is well described by a single product state. Increasing g AB leads to a growth of the von Neumann entropy. For sufficiently large lattice depths S AB is maximized for a specific value of g AB (cf. V 0 /E R = 7.3 and g AB /E R = 0.047), while any further increase of the latter leads to a sudden reduction of the entropy becoming subsequently close to zero which again corresponds to a single product state representation of the many-body wave function.
In order to understand the relationship between the particle distribution of each species and the von Neumann entropy it is useful to investigate the one-body density of the σ-species as a function of the interspecies interaction strength [see Figure 2 (c), (d)]. As it can be seen, up to a specific value g AB the majority species A is distributed over the whole lattice geometry, with a slight decrease (increase) of the density in the central (outer) well(s) for a larger g AB . From this critical value onward the majority species forms a hole in one of the outer wells and it is now only distributed over the four remaining wells. Similarly, we observe for the impurity that up to this critical value of g AB it is distributed over the central three sites, showing an increasing density in the outer wells for larger g AB . However, as soon as a hole is formed in the majority species the impurity localizes in a single outer well which is unoccupied by the majority species. The latter is accompanied by the formation of a two-fold degeneracy in the ground state. In this sense, the ground state is given by the density distribution for large g AB as depicted in figure 2 (c), (d) and its paritysymmetric (with respect to x = 0) counterpart. Consequently, the densities shown for large g AB correspond to only one of the two energetically degenerate ground states. Focusing on the above-described critical value of g AB , e.g. for V 0 /E R = 10.5 we observe a minor population of the impurity in the outer wells. However, this spatial species distribution in ρ (1) B (x) is more pronounced for smaller values of V 0 . This means that the corresponding one-body density is increased in the outer wells. Indeed, for V 0 /E R = 7.3 and g AB /E R = 0.047 [as compared to V 0 /E R = 10.5 and g AB /E R = 0.02 in Figure 2 (d)] the impurity is largely distributed over the lattice geometry and exhibits an increased density in the outer wells [see Figure 2 (b)]. Correspondingly, the density of the majority species is smaller in the outer wells compared to the central ones. This large overlap between the two species on the level of the one-body densities is responsible for the maximized von Neumann entropy at the critical value of g AB . In turn, this explains the sudden decrease of S AB which is associated with a phase-separation of the two species [89][90][91][92] and the formation of a doped insulator for the composite system. Concluding, an increase of the von Neumann entropy is associated with a strong delocalization of the impurity, thereby leading to a large overlap of the two species. From a critical value of g AB onward we observe a phase-separation of the two species which is accompanied by a decrease of the von Neumann entropy to S AB ≈ 0. This can also be viewed as the formation of a particle-hole pair [72], where the majority species forms the hole in one of the outer wells.
To obtain a deeper understanding of the particle distribution, as a next step, we inspect the two-body reduced density which reads where σ, σ ∈ {A, B}. This measure corresponds to the probability of finding a particle of species σ at the position x σ 1 and another particle of species σ at the position x σ 2 [95]. Figure 3 (a 1 ) illustrates this quantity for two particles of the majority species and Figure 3 (b) for one particle of the majority species and one impurity for V 0 /E R = 10.5 and g AB /E R λ = 0.12. This parameter set corresponds to the case where the two species form a particle hole-pair and thereby phase-separate. This fact can also be observed in the behavior of the interspecies twobody density ρ (2) AB [see Figure 3 (b)]. Indeed, the probability of finding a particle of species A and a particle of species B at the same position is approximately zero, i.e. ρ (2) AB (x, x) ≈ 0. Instead, measuring the impurity in one of the outermost wells, we may find the majority species localized in any other well with approximately the same probability. This indicates that the majority species is equally distributed over four out of the five wells. In this context, the question arises how the majority species particles distribute among each other. As shown, in Figure  3 (a 1 ) we find that the probability of detecting two particles of species A at the same position is approximately zero, i.e. ρ (2) AA (x, x) ≈ 0. This means that two particles of the majority species tend to avoid each other and do not occupy the same well. From this we can conclude, that they form a Mott insulator-like state on the four populated wells. In this sense, the complete wave function of the system in this regime can be well described as follows Here, the number state | n σ = |n σ 1 , n σ 2 , n σ 3 , n σ 4 , n σ 5 σ is constructed with a generalized Wannier basis of the lowest band [93,94]. This means that e.g. n σ 1 describes the number of σ atoms in the left localized Wannier state of the lowest band. This state essentially describes the doped insulator configuration of the composite system.
In order to unravel the role of correlations for the ground state distribution in the strong interaction regime, we calculate the noise correlation [63] between particles of species σ and σ which is defined as The noise correlation is a measure for the deviation of the conditional probability of finding two particles at specific positions from the unconditional one given by the product of two single-particle events. In this sense, it gives insight into whether two particles can be viewed as independent from each other or not and therefore suggest the occurrence of beyond single-particle processes in the system. In the former case g (2) σσ = 0 and in the latter case g (2) σσ = 0. In Figure 3 (a 2 ) we show g (2) AA = 0 corresponding to the previously discussed two-body density for two particles of the A species for V 0 /E R = 10.5 and g AB /E R λ = 0.12. Let us remark that in case of g (2) σσ < 10 −2 , two-body processes are negligible and therefore we set the color of the colormap to white. Moreover, we do not show the interspecies noise correlation, since it does not exhibit any significant structures due to the ground state being well described by the wave function of equation 8. The measurement of one particle of species A does not depend on the previous measurement of the particle in species B and vice versa. However, as one might expect the noise correlation between two particles of species A shows a more involved structure. In particular the conditional probability of finding two particles of the majority species at the same position deviates strongly from the unconditional one. This is a Figure 3. The two-body density for (a1) two particles of the majority species ρ (2) AA and (b) for one particle of the majority species and one being an impurity ρ (2) AB . (a2) The noise-correlation for two particles of the majority species ρ (2) AA . We consider V0/ER = 10.5 and gAB/ERλ = 0.12. The particle number of the respective species is NA = 4 and NB = 1.
clear signature of the Mott insulator-like state formed by the A species. Furthermore, we find a difference between the two probabilities in the off-diagonal elements of g (2) AA , which again emphasizes the correlated nature of this state.
In summary, we have found that our system containing a single impurity interacting repulsively with a majority species undergoes a transition regarding the ground state of the system in dependence of the interspecies interaction strength and the lattice depth. This transition manifests itself in an increase of interspecies entanglement with increasing g AB , which is accompanied by a delocalization of the impurity followed by a sudden decrease of the entanglement. The latter is due to the phase-separation (or particle-hole formation) between the species constituting a doped insulator, which takes place for sufficiently large g AB and exhibits an energetic degeneracy of the ground state.

B. Two Bosonic Impurities
In the following, we investigate the ground state of the above-discussed system, but with an additional impurity, i.e. N B = 2. In order to focus on the effect of the interspecies interaction we set the intraspecies interaction strength among the impurity particles to zero, i.e. g BB = 0. Analogously to the analysis above we first examine the interspecies entanglement and the corresponding one-body densities of each species. Evidently, in Figure 4 (a) we observe an alteration of the crossover diagram, given by the von Neumann entropy as a function of g AB and V 0 , compared to the case of a single impurity [see also Figure 2 (a)]. More specifically, for large lattice depths an increase of entanglement between the species takes place up to a critical value of the interspecies interaction strength followed by a sudden reduction of S AB to zero for a further increase of g AB . This increase of entanglement is again related to the delocalization of the impurities in the lattice, thereby enhancing the overlap between the species on the one-body density level [see Figure 4 (f)]. The sudden reduction of S AB in turn is a result of the phase-separation of the two species [see Figure 4 (e),(f)], similar to the case of a single impurity. However, for smaller values of the lattice depth (e.g. V 0 /E R = 5.7) we observe a slightly different behavior of the entanglement and the related one-body densities, as compared to the case of a single impurity. We still find a critical g AB which is associated with the delocalization of the impurity species, while for larger values of g AB the entanglement does not drop to zero. Instead, S AB saturates towards a finite value, which means that the two species remain entangled and the ground state cannot be described by a single product state. The reason for this can be traced back to the distribution of the one-body densities of the two species. For large g AB the impurity species distributes over two of the outer wells with an increased density in the well which is closer to the center [see Figure 4 (d)]. The majority species in turn exhibits a residual density in this very well [see Figure 4 (c)], which leads to a finite overlap of the two species and thereby to a finite entanglement. Let us again remark here that the above-discussed two ground states for V 0 /E R = 5.7 and V 0 /E R = 8.9 and large g AB [cf. Figure 4 (c)-(f)] are two-fold degenerate with a parity-symmetric counterpart (with respect to x = 0).
In the next step, we want to understand how the two impurities are distributed. As a first measure, we investigate the relative distance [30,32] of the two impurities  BB (x 1 , x 2 ) being the previously introduced twobody density of two impurities. Figure 4 (b) shows the impurity distance as a function of the interspecies interaction strength for various lattice depths. We observe that for increasing g AB and fixed V 0 the distance of the impurities decreases and saturates to a specific value for even larger g AB . Also, for smaller V 0 we identify a larger impurity distance. Apparently, in spite of the delocalization of the impurity species in the one-body density up to a critical value, the distance |x 1 − x 2 | decreases, which   suggests that the nature of the delocalization is not completely intuitive and needs to be inspected in detail (see below). The saturation of the impurity distance towards a finite value for large g AB can be associated with the localization of the impurities in a single outer well for very deep lattices or in two adjacent sites for smaller V 0 . The fact that the distance for the latter case should be larger as compared to the former one, becomes evident in the clear separation of the corresponding lines for large g AB [cf. black diamonds and green dashed line in Figure  4 (b)]. In general, note that the decay of the relative impurity distance is indicative of an induced attractive interaction between the impurities [18][19][20]73]. However, this quantity does not allow for a detailed insight into the actual spatial distribution of the individual impurities. This can be gained by investigating the spatially resolved two-body density of the impurities (see equation 7).
Let us in the following focus on the two cases of large g AB , where the impurity species either localizes on a single site [cf. Figure 4 (f)] or on two adjacent sites [cf. Figure 4 (d)], a result that depends on the lattice depth. For the case where the impurities localize on a single site we find that the majority species exhibits mostly a Mott insulator like structure [cf. Figure 5 (d 1 )] as in the single impurity scenario. This implies that the particles of the majority species mostly tend to avoid each other and each one occupies a single distinct lattice site. Moreover, in this context the probability ρ BB of finding two impurities at specific positions x 1 and x 2 accumulates at a single outer site of the lattice geometry [see Figure 5 (e 1 )]. This means that the two impurities tend to occupy the same site. The previously observed phase-separation of the two species is then also reflected in the corresponding two-body density ρ (2) AB [cf. Figure 5 (f 1 )]. In general, we can conclude that two impurities behave similarly to a single impurity for large lattice depths and large interspecies interaction strengths. However, the situation is not that intuitive when the impurity species occupies two adjacent sites for smaller lattice depths. Here, the majority species does not exhibit a Mott insulator-like structure, where the particles avoid each other. Instead, two particles of the majority species can be mostly found either at the same site or on two different ones, while the occupied sites are dominantly those which are not populated by the impurity species [see Figure 5 (a 1 )]. The latter fact can also be observed in the two-body density ρ (2) AB for two particles of different species [see Figure 5 (c 1 ) ,(f 1 )].Indeed, to a large extent the two species avoid each other, which indicates a phase separation where the impurity species occupies two adjacent outer wells and the majority species populates the other unoccupied wells. Nevertheless, this phase-separation is not complete and there is a finite, but small, overlap of the two species which we already discussed in the framework of the von Neumann entropy [see Figure 4  Coming back to the question of how the individual impurities distribute, in Figure 5 (b 1 ) it becomes clear that the impurity species occupying two adjacent sites may either have two particles at the same site or on two different sites. In this sense, the impurities are delocalized over these two sites. We also note that we have a slightly increased probability of finding the impurities at the same site which is adjacent to the sites occupied by  BB and (c1), (f1) for one particle of the majority species and one impurity ρ   the majority species.
In order to understand whether the findings in the two-body densities are dominated by two-body processes we subsequently investigate the noise-correlation for two particles of the same species (see equation 9) for V 0 /E R = 5.7, V 0 /E R = 8.9 and g AB = 0.12. For V 0 /E R = 8.9 the structure of the noise-correlation for two particles of the A species [ Figure 5 (d 2 )] is similar to the one for a single impurity in Figure 3 (a 2 ) which is to be expected since the particles of the majority species are in a Mott insulator-like state. Interestingly, in this scenario the noise-correlation for two impurities only shows a very weak structure of the order of 10 −2 [see Figure 5 (e 2 )]. This means that the accumulation of the two impurities in a single site is well described by the one-body densities such that the measurement of one impurity is independent of the previous one of the other impurity. However, this is not the case for V 0 /E R = 5.7 and g AB = 0.12, where the impurities accumulate in adjacent sites. Here, the conditional probability of finding two particles of the impurity species at the same position deviates strongly from the unconditional one [see Figure 5 (b 2 )]. This is also the case for finding the impurities at different sites. Also, for two particles of the majority species we observe that g (2) AA = 0 in the relevant occupied sites [see Figure 5 (a 2 )], with a structure similar to Figure 5 (d 2 ). Therefore, we can conclude that for V 0 /E R = 5.7 not only the particles of the A species exhibit two-body correlation effects but also the two impurities as compared to the case of V 0 /E R = 8.9.
In summary, we have found that for large lattice depths and strong interspecies interaction strengths two impurities accumulate in a single outer lattice well, a behavior that is similar to the case of a single impurity. As a result the majority species occupies the other four lattice sites in a Mott insulator-like state and the composite system forms a doped insulator configuration. However, for smaller lattice depths the impurity species distributes over two adjacent lattice sites which in turn leads to a larger overlap of the two species. Also, the majority species now dominantly populates the three remaining unoccupied wells. This change in the distribution is accompanied by the presence of non-negligible correlations not only among the majority species atoms, as in the previous case, but also between the impurity atoms.

IV. CORRELATED TUNNELING DYNAMICS
Having analyzed in detail the ground state properties of a lattice trapped bosonic mixture, we subsequently   study its dynamical response upon quenching the interspecies interaction strength. To this end, we prepare our system, for one as well as two impurities, in its ground state for large lattice depths and strong interspecies interaction strengths such that the impurity species occupies one of the outer wells. In this regime the two species phase separate and form a particle-hole pair, as discussed in the previous sections [96]. By suddenly lowering the interspecies interaction strength we aim at initiating a tunneling process of the impurity species through the lattice geometry. We start by examining the tunneling properties in the case of a single impurity and then extend it to two impurities.

A. Transporting a Single Impurity
In the following we prepare our system in its ground state with a lattice depth V 0 /E R = 10.5 and an interspecies interaction strength g AB /E R λ = 0.04. This choice leads to a doped insulator density distribution of the two species as depicted in Figure 2 (c), (d). In this sense the two species phase-separate and thereby form a particle-hole pair where the impurity plays the role of the particle. Moreover, the majority species atoms occupy each well separately. To trigger the dynamics, we quench the interspecies interaction strength to a smaller value such that we cross the transition from S AB ≈ 0 to large values of S AB regarding the ground state entangle-ment crossover diagram shown in Figure 2 (a).
As a representative example of the emergent tunneling dynamics of each species in Figure 6 (a), (b) we present the temporal evolution of the corresponding one-body densities following a quench to g AB /E R λ = 0.025, while keeping fixed V 0 /E R = 10.5. In this case the impurity tunnels through the lattice geometry and ends up well localized in the opposite outer well [see Figure 6 (a)]. On the other hand the hole of the majority species in the initial well vanishes, which means that particles of the majority species travel towards this well which was initially solely populated by the impurity. Finally, for t/(h/E R ) ≈ 600 a hole can again be found at the opposite outer well where now the impurity resides [see Figure 6 (b)].
In order to appreciate the involved tunneling processes we next examine the probability of finding an impurity particle in a specific Wannier state. To this end, we construct the operatorÔ where |w B l w B l | projects the particle of the B species onto the l−th Wannier state of the lowest band. Evaluating this operator with respect to the complete many-body wave function yields the probability P l = Ψ MB |Ô (1) l |Ψ MB of detecting an impurity particle in the l−th Wannier state. Note that the complete many-body wave function is obtained via ML-MCTDHX and subse-quently we analyze this high-dimensional object by eval-uatingÔ (1) l with respect to the wave function. In the following the Wannier states are ordered from left to right, i.e. |w 5 describes the Wannier state which is associated with the initially (t = 0) populated well. The Wannier states prove to be a suitable basis set, since in all cases analyzed in the following, we find that l P l ≈ 99.9%. Figure 7 (b), shows the probability P l of finding the impurity in the l−th Wannier state upon quenching to g AB /E R λ = 0.025. The energetically ordered associated Wannier states are depicted in Figure 7 (a). The probability of finding the impurity in the initially occupied well decreases in the course of time to zero, while the other wells are populated such that at the end of the process a maximum probability in the left Wannier state occurs. This behavior clearly leaves an imprint on the relative energy of the subsystems σ Ĥ rel The reason for the latter will be discussed below.
First, let us turn back to the superposition of the three central Wannier states in the course of time. Interestingly, the propagation of the impurity species cannot be understood as a subsequent tunneling from one well to the adjacent one. The impurity rather delocalizes over the lattice geometry with a maximum delocalization indicated by the white dashed line in Figure 6 (a), (b). The one-body density for this scenario is depicted in Figure  6 (d). Here, the impurity strongly localizes in the central well, with around half of that population in the two directly adjacent sites and a minor density in the outer wells. Apparently, in order to avoid overlap, the majority species exhibits a density minimum in the central well with increasing density towards the outer wells. This can also be observed in the population of the corresponding Wannier states [see Figure 7 (b)]. It is also worth noticing, that the one-body densities of both species clearly deviate from the ones discussed in the context of the ground states in Figure 2 (b) where the roles of the two species are reversed with respect to the distribution of the one-body densities. However, this increased overlap not only increases the interspecies energy [see Figure 7 (c)] but additionally leads to a drastic increase of the von Neumann entropy as in the static case for the ground states [see Figure 6 (e) dashed blue line]. Initially, we start with a von Neumann entropy of S AB ≈ 0 which then increases to a maximum and decreases again when the impurity species resides in the outermost well. This decrease is due to the fact that the two species again phase-separate at later evolution times. Let us remark here that S AB does not drop back to zero since a minor residual density of the impurity remains in the second well.
This tunneling process also has an impact on the correlations among the particles of the majority species. To unravel the role of these correlations, we define a measure for the correlations which are present in each subsystem itself. The spectral decomposition of the one-body density of species σ reads where n σj (t) in decreasing order, obeying j n σj = 1, are the so-called natural populations and Φ jσ (x, t) the corresponding natural orbitals. In this sense, the natural orbitals are the eigenstates, while the natural populations are the corresponding eigenvalues [68,69,84], which are determined by diagonalizing the one-body density matrix. Similar to the Schmidt coefficients the natural populations serve as a measure for the correlations in a subsystem. In this spirit, we define the fragmentation [49,76,77] in the subsystem σ as S σ (t) = − j n σj (t) ln(n σj (t)).
Here, the case of S σ = 0 means that the subsystem σ is not depleted, implying that all particles occupy the same single particle state, i.e. n σ1 = 1. As just discussed, the delocalization of the impurity species is accompanied by a delocalization of the majority species. We observe that this is also reflected in the fragmentation S A of the majority species [see Figure 6 (f)]. At t = 0 the majority species already exhibits a non-negligible fragmentation which is due to the fact that this species resides in a Mott insulator-like state. During the transport process, the majority species delocalizes which leads to an increase of S A , thereby maximizing the latter. The subsequent decrease of the fragmentation can be associated with the transition back to a Mott insulator-like state. Since this transition is not complete due to residual densities in the wells, the fragmentation does not drop back to its initial value.
In order to clarify that this tunneling process is not unique to a specific post-quench interspecies interaction strength, we show in Figure 6 (c) the temporal evolution of the center of mass of the impurity which is defined as The initial value of the impurity's center of mass is given by X B /λ ≈ 1, while the final position when the impurity is completely transported to the opposite outer well is X B /λ ≈ −1. For various post-quench g AB we find that the impurity can be transported to the opposite outer well, while the time it takes for the impurity to reach that well decreases with smaller g AB . We also note that the impurity tunnels back to the initially occupied site afterwards, which shall not be the focus of further discussions. Nevertheless, it becomes additionally clear that it needs a certain minimal post-quench interspecies interaction strength in order to observe the impurity tunneling within the depicted time interval. E.g. for post-quench g AB /E R λ = 0.03 the impurity rather resides in the initially occupied well, instead of tunneling through the lattice geometry, whereas for smaller g AB the tunneling process takes place. As in the case of g AB /E R λ = 0.025, which we discussed in detail, also for smaller g AB the von Neumann entropy as well as the fragmentation of the majority species are affected by the tunneling of the two species. However, the well separated initial increase and the subsequent decrease in both measures become less clear or not even evident anymore, which indicates that also the correlated tunneling takes place in a less structured manner with decreasing post-quench g AB .
In order to understand the above-described difference with varying post-quench g AB , we first analyze the temporal evolution of the two-body densities for an interspecies interaction strength to g AB /E R λ = 0.025. These observables provide a more comprehensive insight into the involved correlated processes indicated by the von Neumann entropy and the fragmentation of the majority species. In Figure 8 we present temporal snapshots of ρ (2) AB and ρ (2) AA . The first column corresponds to t = 0, while the second one refers to the case where the impurity species strongly delocalizes over the lattice geometry, thereby maximizing S AB . The third column represents the time instance when the impurity resides in the opposite outer well with X B /λ ≈ −1. Focusing now on the majority species we clearly identify that throughout the tunneling process two particles of this species avoid to occupy the same site. Even in the case of maximum delocalization of the majority species two particles do not reside in the same site. We can understand this process as a superposition of all number states where four particles are distributed separately among the five Wannier states, e.g. |1, 1, 1, 1, 0 A , |1, 1, 1, 0, 1 A etc. . In this sense, the hole in the majority species is delocalized over the lattice geometry. Moreover, we can infer that measuring the impurity at a specific site we will most probably not find a majority species atom at the same site. Even in the case of maximum delocalization these two particles will avoid each other. In this manner, the increase of S AB in the context of the delocalization is not a trivial single-particle effect. It is rather due to the strong avoidance of the two species and thereby the involved superposition of many contributing states, e.g. |1, 1, 1, 1, 0 A ⊗ |0, 0, 0, 0, 1 B , |1, 1, 1, 0, 1 A ⊗ |0, 0, 0, 1, 0 B etc., which leads to the increase of the von Neumann entropy. This in turn leads to the delocalization of the two species on the level of the one-body densities when integrating out the respective degrees of freedom (see equation 5). Based on these findings we have strong indications that the tunneling of the two species manifests itself in the effective transport of a particle hole-pair from one of the outer wells to the opposite one.
To further quantify this tunneling process we determine the probability of finding a particle hole pair in any lattice site in the course of time. We define this as the with respect to the system's wave function, i.e. Ô (2) ph . Here, |w B l w B l | projects the particle of the B species onto the l−th Wannier state and |w iA l w iA l | projects the i−th particle of the A species onto the l−th Wannier state. 1 σ are the unity operators of the respective subsystems σ. Figure 9 (a) presents the temporal evolution of the probability of finding a particle hole pair in any lattice site for various post-quench interspecies interaction strengths. Since our system is prepared in its ground state for all post-quench g AB , initial particle-hole pair probability of Ô (2) ph (t = 0) ≈ 98.8% holds. This means that we can safely assume that our system exhibits a particle-hole pair.
Focusing on the case of g AB /E R λ = 0.025, which corresponds to the two-body densities depicted in Figure 8, the quench leads to a slight reduction of Ô (2) ph which recovers again for larger times. Closely inspecting the particle hole-pair probability [see Figure 9 (b)], we can deduce that the strongest decrease happens at the time interval in which the two species strongly delocalize over the lattice geometry. For this reason, we can assume that our system is not solely described by the superposition of number states exhibiting a particle-hole pair, e.g. |1, 1, 1, 1, 0 A ⊗|0, 0, 0, 0, 1 B , |1, 1, 1, 0, 1 A ⊗|0, 0, 0, 1, 0 B etc. . Nevertheless, the admixture of other states is rather small since in the case with the largest deviations we still find a particle hole-pair with a probability of Ô (2) ph ≈ 93.2%. Moreover, we observe that Ô (2) ph increases again in the course of time, if not completely to its initial value, which is associated with the impurity then occupying dominantly the opposite outer well. In contrast, when quenching to smaller values of g AB , the stability of the particle-hole pair cannot be guaranteed anymore. For a quench to g AB /E R λ = 0.01 we lose up to 20% of the particle-hole pair which becomes even more for a quench close to zero, i.e g AB /E R λ = 0.002. In the latter case Ô (2) ph drastically decreases to 65% in the first few time steps. In this sense, in order to maintain a stable particle-hole pair during the transport of the impurity to the opposite outer well it is necessary not to quench too strongly. On the other hand a quench which is too weak allows for a stable particle-hole pair, which resides in the well initially occupied by the impurity, e.g. g AB /E R λ = 0.03 [see Figure 9 (a) red line].
In conclusion, we have found that quenching the binary mixture starting in a phase separated state leads to a very controlled correlated tunneling dynamics. On the one hand the majority species tunnels such that its particles do not occupy the same lattice site. On the other hand we observe the effective transport of a particle-hole pair from one of the outer wells to the opposite one whose stability strongly depends on the presence of a finite interspecies interaction strength. Moreover, the tunneling dynamics is accompanied by a strong entanglement between the species as well as correlations among the majority species atoms.

B. Transport Properties of Two Bosonic Impurities
Having understood the correlated tunneling dynamics of a single impurity interacting repulsively with a majority species we now turn to the case of two non-interacting bosonic impurities. We prepare our system in its ground state characterized by a lattice depth V 0 /E R = 8.9 and interspecies interaction strength g AB /E R λ = 0.05, while g AA /E R λ = 0.04. This leads to a density distribution of the two species as depicted in Figure 4 (e), (f), such that the two species phase-separate. Moreover, the two impurities accumulate at the same single lattice site, while the majority species atoms strongly avoid each other, occupying the lattice sites each one separately. In this sense, our initial state is similar to the one with a single impurity, apart from the additionally added impurity.
Quenching the interspecies interaction strength to g AB /E R λ = 0.002 we observe a tunneling of the impurity species to the opposite outer well [see Figure 10 (b)]. However, compared to the case of a single impurity a larger portion of density remains in the wells which are traversed during the dynamics. This means that compared to the single impurity case the transport portion of two impurities is less complete due to some density remaining in the traversed lattice sites. In turn the majority species tunnels in a counterflow into the opposite direction as compared to the direction of the impurities. In this context, one might again ask whether this process can be described by an effective transport of a particle hole-pair. Figure 10 (a) shows the corresponding particle hole-pair probability Ô (2) ph during the evolution.
Initially, the particle hole-pair process is the dominant tunneling channel with a probability of Ô (2) ph ≈ 98.6% which later on drastically decreases to Ô (2) ph ≈ 63.8%. The subsequent increase of this probability, which is associated with the impurities residing in the opposite outer well, does not revive to its initial value, indicating the loss of the particle hole-pair. This effect is similar to the one discussed for a single impurity when quenching to very small values of g AB , where the stability of the particle hole-pair cannot be guaranteed either. Apart from that, the system still exhibits a rather pronounced increase of correlations in the course of time. This involves for example the increase of the von Neumann entropy which can be traced back to the residual impurity density in the remaining wells during their transport to the other side of the lattice geometry. As a result the spatial overlap between the species is increased and thereby also S AB [see Figure 10  Furthermore, it is of interest to analyze in which manner the two bosonic impurities propagate in relation to each, i.e. whether they move as a pair or they delocalize in the course of time. In order to answer this question we inspect the two-body density of the impurities ρ (2) BB upon quenching to g AB /E R λ = 0.002 [see Figure 10 (g)-(i)]. The initially accumulated (in an outer well) impurities delocalize over next-neighbour sites, which means that it is possible for the two impurities to either reside at the same site or in adjacent ones [see Figure 10    A (x) of the species A upon quenching the interspecies interaction strength from gAB/ERλ = 0.05 to gAB/ERλ = 0.002. Temporal evolution of (d) the von Neumann entropy SAB, (e) the fragmentation SA and (f) the fragmentation SB upon quenching to various gAB. (g)-(i) Snapshots of the two-body density of the impurities ρ (2) BB upon quenching to gAB/ERλ = 0.002. The particle number of the respective species is NA = 4 and NB = 2 and we set the lattice depth to V0/ER = 8.9.
stances at which the probability of finding the impurities at the same site is more pronounced as compared to detecting them in adjacent ones (not shown here). From this we can conclude that in general the impurities delocalize during the propagation until they reach the opposite outer well where they eventually strongly localize [see Figure 10 (i)].
Let us finally remark that in comparison to the single impurity case, for two impurities it is necessary to quench to lower interspecies interaction strengths in order to achieve a significant transport of the impurity species to the opposite outer well. For smaller quench amplitudes we either find no tunneling of the two species [cf. Figure 10 (d)-(f) green crosses] or the impurities delocalize over the lattice geometry without accumulating in the opposite outer well. Therefore, a transport of the impurity species to the opposite outer well as in the case of a single impurity is only possible for large quench amplitudes. However, this leads to a loss of the initial particle-hole pair and a less structured development of correlations.

V. CONCLUSIONS AND OUTLOOK
We have investigated the ground state properties of a binary bosonic mixture trapped in a one-dimensional lattice geometry consisting of a majority species which is doped by one or two bosonic impurities. We demonstrate the existence of a crossover diagram of the interspecies entanglement as a function of the lattice depth and the interspecies interaction strength. The transition from strong interspecies entanglement to S AB = 0 is accompanied by a crossover from a spatially delocalized impurity species to its strong localization in one of the outer lattice wells. For large lattice depths we find this transition for the case of a single as well as two bosonic impurities. Moreover, analyzing the corresponding two-body densities we can conclude that the majority species occupies a Mott insulator-like state, while the two species phaseseparate, thereby forming a particle hole-pair. For two impurities we additionally observe in the case of smaller lattice depths a localization of the impurity species in two adjacent sites. This phenomenon also manifests itself in the impurity-impurity correlations, which we do not find when the impurities localize in a single site.
Having understood the ground state distributions, in the next step we aimed at transporting the impurity species through the lattice by performing an interspecies interaction quench from a Mott to the superfluid phase of the composite system. For a single impurity the initial particle hole-pair tunnels to the opposite outer lattice well, while we are able to control the stability of the particle hole-pair by varying the post-quench interspecies interaction strength. Here, it is important to perform Figure 11. (a) The von Neumann entropy SAB as a function of the interspecies interaction strength gAB and the lattice depth V0. The two-body density for (b) two particles of the majority species ρ (2) AA and (c) for one particle of the majority species and one impurity ρ (2) AB . We consider V0/ER = 10.5 and gAB/ERλ = 0.12. The particle number of the respective species is NA = 5 and NB = 1.
weak amplitude quenches in order to maintain the initial particle hole-pair. However, for two impurities we cannot guarantee a stable transport of this pair. This is due to the fact that only large quench amplitudes can initiate a significant transport of the impurities to the opposite outer well at all. As in the case of a single impurity such large quench amplitudes lead to a strong loss of the particle hole-pair probability.
The understanding of the crossover between a spatially delocalized to a localized bosonic impurity species serves as a perfect starting point for even more complex setups, e.g. by introducing more lattice atoms and sites. Indeed, there are several possible directions of future investigations. For example, it would be of interest to allow for a spin degree of freedom in both species, such that the particle hole-pair may carry an additional effective spin. In this spirit, a quench of the interspecies interaction strength might lead to a redistribution of the spins in the system or an effective spin transport. Another fruitful perspective is to investigate the ground state phase diagram of the setting considered herein but including dipolar interactions within and between the species. In this way, it would be possible to generate more phases due to the long-range character of the interactions.
Appendix A: Doping a unit-filling bosonic majority species To explicitly showcase the generalization of our findings regarding the ground state properties of the considered lattice trapped bosonic mixture, we consider a unit filling of the majority species, doping it with a single impurity, i.e. N A = 5 and N B = 1. We explore the ground state of the system for varying V 0 and g AB , while fixing the intraspecies interaction strength to g AA /E R λ = 0.04. Figure 11 (a) shows the von Neumann entropy of the ground state as a function of the interspecies interaction strength g AB and the lattice depth V 0 . As in the case of N A = 4 and N B = 1 [ Figure 2 (a)], we find a region of increased interspecies entanglement, which decreases to zero for a further increase of g AB . This is again associated with a strong delocalization of both species for large S AB and a phase separation in case of S AB ≈ 0 occurring for large g AB . However, there is a slight change compared to N B = 4 in the corresponding two-body densities for V 0 /E R = 10.5 and g AB /E R λ = 0.12, where an interspecies phase-separation takes place. Since the majority species exhibits an additional particle, the latter distributes over the four lattice sites which are occupied by the majority species in case of an interspecies phaseseparation. For this reason a non-negligible, but small, probability of finding two particles of the A species at the same site occurs. In ρ (2) AB the phase separation is still reflected with an increased probability of finding the majority species in the two central occupied sites, which is due to the additional majority species particle. In this sense, also for N A = 5 we find the formation of a particle hole-pair for large g AB , while the additional particle in the majority species delocalizes over the sites occupied by the latter.