Entanglement-assisted tunneling dynamics of impurities in a double well immersed in a bath of lattice trapped bosons

We unravel the correlated tunneling dynamics of an impurity trapped in a double well and interacting repulsively with a majority species of lattice trapped bosons. Upon quenching the tilt of the double well it is found that the quench-induced tunneling dynamics depends crucially on the interspecies interaction strength and the presence of entanglement inherent in the system. In particular, for weak couplings the impurity performs a rather irregular tunneling process in the double well. Increasing the interspecies coupling it is possible to control the response of the impurity which undergoes a delayed tunneling while the majority species effectively acts as a material barrier. For very strong interspecies interaction strengths the impurity exhibits a self-trapping behaviour. We showcase that a similar tunneling dynamics takes place for two weakly interacting impurities and identify its underlying transport mechanisms in terms of pair and single-particle tunneling processes.

We unravel the correlated tunneling dynamics of an impurity trapped in a double well and interacting repulsively with a majority species of lattice trapped bosons. Upon quenching the tilt of the double well it is found that the quench-induced tunneling dynamics depends crucially on the interspecies interaction strength and the presence of entanglement inherent in the system. In particular, for weak couplings the impurity performs a rather irregular tunneling process in the double well. Increasing the interspecies coupling it is possible to control the response of the impurity which undergoes a delayed tunneling while the majority species effectively acts as a material barrier. For very strong interspecies interaction strengths the impurity exhibits a self-trapping behaviour. We showcase that a similar tunneling dynamics takes place for two weakly interacting impurities and identify its underlying transport mechanisms in terms of pair and single-particle tunneling processes.
In this context, especially strongly particle imbalanced mixtures have attracted a lot of interest recently. In the extreme case such systems consist of a single impurity immersed in a majority species. These setups have been studied theoretically [25][26][27][28][29][30][31][32] and experimentally [33][34][35][36], for a single impurity, serving as a simulator of polaron physics, as well as for many impurities [37][38][39][40][41][42] and are indeed a subject of ongoing research. While the ground state properties of a single impurity in a bath are to a certain extent well understood, less focus has been placed on the transport properties and the emergent collisions of the impurity through the bath [43][44][45][46]. Indeed, in these systems correlation effects, such as entanglement, are expected to be a crucial ingredient since the impurities form a few-body subsystem [47]. Moreover, the underlying trapping potential plays an important role for the behaviour of the impurity species, which has been analyzed for homogeneous systems [48][49][50], harmonic confinements [51][52][53][54][55] as well as lattice potentials [33,56,57]. The majority of the abovementioned investigations have been focusing on the case where both species are trapped in the same geometry. However, introducing different trapping potentials for each species is expected to alter sig-nificantly the observed dynamics. A setting of particular interest involves a bath of lattice trapped bosons which act as multiple material barriers for the tunneling dynamics of the impurity.
In the present work we explicitly focus on an impurity which is confined in a one-dimensional double well and interacts repulsively via contact interaction with a majority species of bosons trapped in a lattice. For single component bosons in a double well the analogue of the well-known superconducting Josephson junction can be established. The bosonic Josephson junction provides the testbed for many, also experimentally observed, intriguing phenomena, such as Josephson oscillations, macroscopic quantum self-trapping [58][59][60][61][62][63] and correlated pair tunneling [64][65][66]. Extensions of these phenomena to multicomponent setups have also been extensively studied, see for instance [67][68][69][70]. Herein, we extend these investigations by exploring the dynamics of impurities in a double well immersed in a few-body bath of lattice trapped bosons. This gives rise to an effective potential for the impurities whose shape strongly depends on the interspecies interaction strength. Depending on the latter, one can realize tunneling scenarios which are beyond the well-known regimes of Josephson oscillation and quantum self-trapping and rely on the interspecies entanglement. This can be of particular interest for future applications in atomtronics [71].
In our setup of a single impurity in a double well the dynamics is steered by the repulsive coupling to the majority species. Varying the interspecies interaction strength we unravel different dynamical response regimes of the impurity upon quenching the tilt of the double well. These regimes range from rather irregular tunneling in the double well for small interspecies interaction strengths to dynamical self-trapping in a single site for very strong couplings [6,7,55]. For intermediate coupling strengths we observe a strong impact of the density distribution of the majority species on the arXiv:1909.00998v2 [cond-mat.quant-gas] 13 Feb 2020 impurity's tunneling dynamics. The impurity initially collides with the material barrier imposed by the density of the majority species and then tunnels to the corresponding other site of the double well. This offers a controlled way of transporting the impurity within the double well. The entire tunneling process in the case of intermediate interspecies interaction strengths is accompanied by a strong entanglement between the subsystems revealing the complexity of this phenomenon. We remark that in the absence of entanglement this process does not take place. Additionally, in this case the self-trapping behaviour is altered. Surprisingly, we find that the dynamics of the impurity can be described in terms of Wannier states [41,42] which are associated with the superposition of the effective time-averaged potential induced by the density of the majority species and the double well potential. This proves to be a valuable tool that captures the dynamics of the impurity adequately, even though a strong entanglement persists throughout the dynamics [12,46]. Due to the strong correlations appearing in our system it is necessary to utilize an approach which operates beyond lowest band and mean-field approximations, such as the Bose-Hubbard model or Gross-Pitaevskii approximation. Therefore we track the emergent non-equilibrium dynamics by employing the Multi-Layer Multi-Configurational Time-Dependent Hartree Method for atomic Mixtures (ML-MCTDHX) [72][73][74] that enables us to capture all the important particle correlations.
Our work is structured as follows: In section II we present the system under investigation and the employed computational methodology. In section III we unravel the quench-induced tunneling dynamics of the impurity, revealing also the crucial role of the inter-and intraspecies correlations. Section IV is dedicated to an in-depth characterization of the microscopic effects involved in the dynamical response of the impurity. We extend our results to the case of two weakly interacting impurities in section V and conclude with a summary of our findings and a discussion of future directions in section VI.

A. Setup and Hamiltonian
Our setup consists of two different species of bosons A and B, also referred to as the majority species and the impurity, respectively, which interact repulsively via a contact potential of strength g AB . Each species is confined in a different onedimensional optical potential at zero temperature. Experimentally this can be realized by preparing e.g. 87 Rb atoms in two different hyperfine states, i.e. |F = 2, m F = −2 and |F = 1, m F = −1 , thereby obtaining a two-species bosonic mix-ture. Utilizing the so-called 'tune-out' wavelength [75,76] it is possible to create species-dependent potentials, such that the two species experience different optical potentials [3]. The majority A species, composed of bosons of mass m A and interacting repulsively via a contact interaction of strength g AA , is trapped in a six-well lattice potential. The minority B species on the other hand, consisting of N B impurities of mass m B interacting repulsively via a contact interaction of strength g BB , resides in an initially tilted double well potential. The resulting many-body Hamiltonian of the system reads Here, the lattice potential V A = V 0 cos 2 (k 0 x A i ) of the majority species is characterized by its depth V 0 and wave vector k 0 = π/l where l denotes the distance between two successive minima of the potential. The double well of the impurities is constructed by the combination of a harmonic oscillator potential with frequency ω B and a Gaussian potential characterized by a width w and a height h. Additionally, we superimpose a linear tilting potential V tilt = αx B i to the double well leading to an asymmetry between the two wells, whose degree can be controlled by the parameter α. Assuming zero temperature we can model the inter and intraspecies interaction potential between the atoms via a bare delta potential with effective cou- where σ and σ refer to the corresponding species A and B [4]. Here, M AB = m A m B m A +m B represents the reduced mass and a ⊥ = h M AB ω ⊥ the transversal length scale which is steered by the frequency of the transversal confinement ω ⊥ perpendicular to the one-dimensional Bose gas. Apart from varying ω ⊥ , it is possible to control the coupling strength g σσ through the free space, threedimensional scattering length a σσ 0 which can be tuned via Feshbach resonances in magnetic or op-tical fields [5,[77][78][79][80].
Throughout this work we consider a fixed number of bosons for the majority species N A = 8 and set m = m A = m B . As mentioned above, our setup can be experimentally realized by considering two hyperfine states of 87 Rb. Note that we have also simulated the corresponding dynamics of a mass imbalanced system consisting e.g. of a 87 Rb bosonic ensemble and a 133 Cs impurity. For this latter case we confirmed that an overall similar phenomenology compared to the mass balanced case occurs but the emerging tunneling regimes to be presented below take place at smaller interspecies interaction strengths. The energy scales for the Hamiltonian in Eq. (1) are given in units of the recoil energy E r =h 2 k 2 0 /(2m), whereas the length and time scales are expressed in units of k −1 0 and ω −1 For the lattice potential of the majority species we use a depth of V 0 /E r = 8. The harmonic part of the double well potential of the impurities has a harmonic oscillator frequency of ω/ω r = 0.1 · √ 2 and the barrier height and width are h/E r k −1 0 = 2 and w/k −1 0 = 1, respectively. Furthermore, the intraspecies interaction strength among the bosons of the majority species is kept fixed to the value g AA /E r k −1 0 = 1. Hard-wall boundary conditions are imposed at x/k −1 0 = ±3π. In the following, we present the quench protocol which induces the tunneling dynamics. A sketch of the employed procedure is depicted in Figure  1. First, we obtain the many-body ground state of our system, assuming the above mentioned parameters. Here, the tilting strength of the double well is set to α/E r k −1 0 = 0.1 (the effect of a smaller tilting strength is analyzed in the Appendix), such that the impurities localize in the left well of the asymmetric double well potential. To trigger the tunneling dynamics of the impurities the system is quenched to a geometry, constituting a symmetric double well, i.e. the tilting strength is set to α = 0. Varying the interspecies interaction strength g AB , we explore the dependence of the system dynamics on g AB .

B. Approach to the correlated many-body dynamics
To unravel the dynamics of the system we employ ML-MCTDHX [72][73][74]. As explicated below, this ab initio method gains its efficiency from the time-dependent and with the system co-moving basis set. In the first step, the total many-body wave function |Ψ MB (t) is expanded with respect to M different species functions |Ψ σ (t) for each of the species σ and expressed according to the following Schmidt decomposition [81] Here, the Schmidt coefficients λ i (t), in decreasing order, provide information about the degree of population of the i-th species function and thereby signify the degree of entanglement between the two species. In the case that only one Schmidt coefficient is non-zero, the species A and B are not entangled with each other and the system can be described by a species mean-field ansatz (M = 1). However, in general it is necessary to provide several species functions for the expansion of the total many-body wave function, since entanglement might prove crucial for the adequate description of the systems dynamics.
Furthermore, the species wave functions |Ψ σ (t) describing an ensemble of N σ bosons are expanded in a set of permanents, namely Such an expansion allows us to take intraspecies correlations of the σ-species into account. Moreover, in this expression the vector n σ = (n σ 1 , n σ 2 , ...) describes the occupations of the time-dependent single-particle functions (SPF) of the species σ, which are further expanded in terms of a timeindependent discrete variable representation [82]. The notation n σ |N σ indicates that for each n σ i the particle number conservation condition i n σ i = N σ has to be fulfilled. For the time propagation of the many-body wave function we employ the Dirac-Frenkel variational principle δΨ MB | (i∂ t − H) |Ψ MB [83][84][85] with the variation δΨ MB and obtain the corresponding equations of motion [74,86].
In conclusion, the ML-MCTDHX method takes all inter-and intraspecies correlations into account and gives us access to the complete manybody wave function. In contrast to standard approaches, where the wave function for solving the time-dependent Schrödinger equation is built upon time-independent Fock states with time-dependent coefficients, the ML-MCTDHX method takes a comoving time-dependent basis into account, where the Fock states, spanned by the SPFs, as well as the coefficients are time-dependent. This concept of a time-dependent basis reduces not only the required number of basis states and, hence, improves the computational effort, but it also provides at the same time an accurate description of the system's many-body state. We note here that in order to ensure the convergence of our many-body simulations, to be presented below, we have employed M = 6 (M = 10) species and d A = 6, d B = 6 (d A = 6, d B = 6) single-particle functions respectively for the case of a single (two) impurity atom(s). In this context we define the orbital configuration C = (M, d A , d B ) which determines the size of the truncated Hilbert space.

III. CORRELATED TUNNELING DYNAMICS OF A SINGLE IMPURITY
In the following we consider a mass-balanced bosonic mixture described by the Hamiltonian of equation 1 where the majority species consists of N A = 8 and the impurity species of N B = 1 particles. We initially prepare our system in its ground state with a tilting strength of α/E R k −1 0 = 0.1 for different interspecies interaction strengths g AB . Due to the initial tilt the impurity is found to be well localized in a single site of the double well potential. Moreover, we find that the impurity species exhibits a rather large spatial overlap with the majority species for small g AB [see Figure 1 (b)] which of course reduces with increasing repulsive g AB [see Figure 1 (c)]. In particular, for small g AB the majority species occupies all sites of the lattice potential, such that the impurity strongly overlaps with it [cf. Figure 1 (b)]. However, for strong repulsive interactions the majority species depopulates the well of the lattice potential in which the impurity tends to localize, resulting in a weak spatial overlap of the two species [cf Figure 1 (c)]. Upon quenching the tilting strength to α/E R k −1 0 = 0 towards a symmetric double well we initiate the tunneling dynamics, thus favoring the tunneling of the impurity to the right well as the corresponding energy offset between the two wells vanishes, see also Fig. 1 (a). As a consequence the impurity becomes mobile, thereby colliding with the lattice trapped majority species which in general acts as a material barrier for the impurity dynamics. Varying the interspecies interaction strength we find four different regimes for the dynamical response of the impurity (see below).
As a first step, we quantify these regimes by monitoring the time evolution of the one-body density ρ (1) σ (x, t) = Ψ MB (t)|Ψ † σ (x)Ψ σ (x)|Ψ MB (t) of the corresponding subsystems σ. The spectral decomposition of the σ-species one-body density is given by where n σj (t) are the so-called natural populations and Φ σj (x, t) the corresponding natural orbitals. The dynamics of ρ (1) σ (x, t) is presented in Figure 2, for different interspecies interaction strengths g AB . As it can be seen ρ (1) σ (x, t) exhibits four distinct dynamical response regimes. For small interspecies interaction strengths, in our case g AB /E R k −1 0 = 0.2, the impurity undergoes a rather complex tunneling dynamics to the other site of the double well [ Figure 2 (a)]. This is a single-particle effect caused by the strong initial tilt and is therefore also present for g AB = 0. For short evolution times, i.e. 0 < t/ω −1 r < 50, the impurity performs oscillations in the initial well and then tunnels [see ellipse in Figure 2 (a)] to the other well. Here, the oscillations within each of the two wells, which still persist even for g AB = 0 (not shown here), are caused by the rather strong initial tilt of the double well and are not present for smaller tilts 1 [see Figure  11 (a)]. In this sense, the majority species barely affects the tunneling dynamics of the impurity and exhibits weak amplitude modulations from its initial profile due to the finite g AB [ Figure 2 However, for larger coupling strengths the impurity is strongly influenced by the density distribution of the majority species, e.g. see Figure 2 (b),(f). The majority species distributes over the lattice such that ρ (1) A (x, t) is accumulated close to the minima of the lattice potential. Due to the repulsive interspecies interaction the impurity has to overcome on top of the double well barrier these additional material barriers imposed by the accumulation of the density of the majority species. 1 We note that this tunneling behaviour differs from that of a single particle in a double well potential in the case of smaller tilts, yielding a single frequency Rabi tunneling. This leads to an oscillation of the impurity through the neighboring density maximum of the A species [see white rectangle in Figure 2 (b)]. This tunneling through the material barrier imposed by the majority species we will refer to as material barrier tunneling in the following. Throughout this enduring oscillation process the impurity performs a transport [87,88]   and as a result we enter the self-trapping regime. We remark that this self-trapping behavior is caused by the presence of the majority species, in sharp contrast to the wellknown case of interacting bosons in a double well. Note also that the impurity undergoes dipole-like oscillations within the left site of the double well. Also, we have checked that this self-trapping behavior [cf. Figure 2 (d)] of the impurity persists up to t/ω −1 r = 400 evolution times (not shown here). Considering the behavior of the majority species A, we observe the development of excitations of ρ (1) A (x, t) as a back-action of the tunneling pro-cess of the impurity [12]. In particular, ρ Predominantly, this is the case for the inner four wells. This behaviour of the majority species is caused by the repulsive interspecies interaction which leads in the course of the impurity tunneling to a shift of the density of species A, thereby reducing the overlap between the species. In the extreme case [cf. g AB /E R k −1 0 = 4.0] where the impurity remains localized in one site of the double well, the majority species redistributes such that a density hole is formed in one lattice site [cf. Figure 2 (h)], in order to avoid the impurity. Here, the overall density of the majority species barely changes in time due to the absence of the impurity's tunneling.
In order to quantify the dynamical response of the system even further it is convenient to analyze how strongly the time-dependent many-body wave function deviates from the initial state |Ψ 0 at t = 0 in the course of time. This is well captured by the fidelity F (t) = | Ψ MB (t)|Ψ 0 | 2 which is defined as the overlap between the time-dependent and the initial wave function. Figure 3 shows the fidelity F (t) for various interspecies interaction strengths corresponding to the four different tunneling regimes identified in the time evolution of the one-body densities in Figure 2. We clearly observe that the behavior of the fidelity is qualitatively driven by the one-body density distribution of the impurity over time. For the cases in which the impurity tunnels to the other site of the double well [ Figure 2 (a),(b)], the fidelity deviates significantly from unity, i.e. |Ψ MB (t) deviates from the ground state |Ψ 0 . However, in the regimes where the tunneling of the impurity is suppressed the fidelity remains close to unity, e.g. see F (t) for g AB = 2.0, 4.0. In this sense,   . Nevertheless, using solely the fidelity it is not possible to distinguish between the different tunneling mechanisms. In order to discern between the above-mentioned four possible regimes of the impurity's dynamical response it is useful to consider the integrated one-body density of the impurity where L is the size of the system. This quantity provides the probability of finding the impurity in one half of the initially populated well of the double well. Indeed, the integrated density allows to distinguish between the emergent tunneling dynamics since it incorporates the effect of the material barrier. In Figure 4 (a) we show the temporal evolution of this quantity for the correlated many-body approach 2 for different g AB . We clearly observe four distinct regimes for the response of the impurity which correspond to the one-body densities in Figure 2 (a)-(d). Indeed, in regime I an irregular oscillatory pattern of the integrated density is found. Regime II exhibits regular oscillations whose intensity decays in time. This corresponds to the material barrier tunneling with a final transfer of the impurity to the other site of the double well [see 2 (b)]. In regime III the oscillations of the integrated density remain stable in time which is due to the material barrier tunneling of the impurity in the initially populated well without a transfer to the other site. Finally, in regime IV we find a higher-frequency oscillatory behavior with a finite amplitude throughout the evolution. This behavior corresponds to the selftrapping regime [see Figure 2 (d)]. Regarding the dependence of the tunneling behavior on the different system parameters see Appendix B.
However, so far we did not get insight into the degree of the system's correlation throughout the dynamics. To unravel the degree of correlations which accompanies the tunneling dynamics of the impurity we distinguish between inter-and intraspecies correlations. The former are described by the Schmidt coefficients (Eq. 3), which provide a measure for the degree of entanglement between the subsystems, whereas the latter can be inferred from the natural populations (cf. Eq. 5). Since the B species consists of a single particle, the natural populations of the B species coincide with the Schmidt coefficients. Therefore, in the following we invoke the deviation 1 − n B1 (t) as a measure of entanglement between the subsystems. Accordingly, 1 − n A1 (t) indicates the degree of intraspecies correlations of the majority species. The temporal evolution of the depletion 1 − n σ1 (t) of the most populated natural orbital of the A and the B species is illustrated in Figure 5 for different g AB upon quenching the tilting strength. We observe that for small interspecies interaction strengths, i.e. g AB = 0.2, the subsystems are mainly disentangled throughout the dynamics, since 1 − n B1 (t) ≈ 0. Increasing the interspecies interaction strength to g AB /E R k −1 0 = 1.0 the subsystems become strongly entangled in the course of time, i.e. 1 − n B1 (t) > 0. This can be associated with tunneling of the impurity to the other site of the double well and the involved increasing interspecies interaction between the subsystems. Naturally, the motion of the impurity through the majority species has an impact on the natural populations of the A species which is connected to the intrinsic tunneling processes of the A species in the lattice potential [cf. Figure 2 (e)-(f)]. Indeed, 1 − n A1 > 0 independently of g AB and it is maximized in the above-described third tunneling region. Interestingly, the rather strong degree of entanglement remains in the self-trapping regime for g AB /E R k −1 0 = 4.0, even though the impurity barely overlaps with the majority species.
In order to emphasize the importance of the entanglement between the subsystems for the tunneling behavior of the impurity, we additionally perform calculations assuming only a single product state |Ψ MB = |Ψ A ⊗ |Ψ B in Eq. (3), thereby neglecting all interspecies correlations. The dynamics of the σ-species one-body densities employing a species mean-field ansatz, i.e. assuming a single product state between the species but still including intraspecies correlations, are shown in Figure 6. For g AB /E R k −1 0 = 0.2 we find no visible differences between the full many-body approach and the species mean-field calculations. This is  We consider a minority species consisting of NB = 1 particle and a majority species of NA = 8 particles which interact repulsively with gAA/Erk −1 an expected result, as the degree of entanglement is rather small for these interactions (cf. Figure  5). However, as soon as entanglement becomes important, we find strong deviations in the corresponding one-body densities. In particular, for Figure 6 (b)] we do not observe the previously predicted tunneling to the other site of the double well [cf. Figure 2 (b)]. Furthermore, the onebody density of the impurity for g AB /E R k −1 0 = 2.0 exhibits additional oscillation frequencies in the species-mean field scenario [cf. Figure 6 (c)] compared to the full many-body case [2 (c)]. In the selftrapping regime, g AB /E R k −1 0 = 4.0, the species mean-field calculations seem to capture the dynamics quite well at first glance. However, on a closer inspection of the one-body density it turns out that the spatial position of the impurity differs compared to the complete many-body approach, while the temporal oscillations of the density are also of different amplitude and frequency, see Figures 2 (d) and 6 (d). This general difference in the tunneling dynamics is well captured by the integrated density shown in Figure 4. Indeed, the species mean-field ansatz is not able to recover regime II in Figure 4 (a), while the self-trapping regime III in Figure 4 (b) is strongly altered compared to that one in Figure 4 (a) corresponding to regime IV in the many-body treatment. Even regime III in Figure 4 (a) is quantitatively changed when using the species mean-field ansatz [cf. Figure 4 (b) regime II]. In this sense, entanglement between the impurity and the majority species plays a crucial role, in order to describe the dynamics correctly.

IV. CHARACTERIZATION OF THE IMPURITY DYNAMICS
To analyze the tunneling behavior of the impurity and the accompanying correlations due to the presence of the majority species (cf. Figure 2) we next develop an effective potential model for the impurity. This effective potential is obtained by superimposing the time-averaged density of the majority species to the external double well potential. To adequately describe the dynamical response of the majority species we employ the associated Wannier functions. In particular, we project the complete many-body wave function obtained via ML-MCTDHX onto these Wannier functions in order to analyze the behavior of the impurity in a fixed basis set.
A. Construction of the effective potential Initially, we prepare our system such that it is given by the ground state of the Hamiltonian (eq. 1) with an underlying asymmetric double well. With respect to the quenched Hamiltonian (t > 0) our system and in particular the impurity is energetically excited due to the tilting. This enables the impurity to tunnel through the potential barrier of the double well into the right well. However, as mentioned in section III, for specific interspecies interaction strengths this residual energy appears to be not large enough to overcome the potential barrier. Therefore, the impurity B rather performs a tunneling in the initial site of the double well through the material barrier imposed by the onebody density of the majority species [cf. Figure 2 (c)].
In the following, we aim at understanding this tunneling behavior using an effective potential for the impurity. We remark that this effective potential serves only as a tool for an analysis of the underlying tunneling processes. Integrating out the majority species we arrive at the following effective potential for the impurity This effective potential is composed by the double well potential V B and the one-body density of the majority species ρ A being weighted by the number of particles N A and the interspecies interaction strength g AB . Note that ρ (1) A (x B , t) is calculated within the correlated many-body approach and thereby includes all necessary correlations. ρ (1) A (x B , t) cannot be recovered within a mean-field approach. To proceed, we average this effective potential over time and obtain a Time-Averaged Effective Potential (TAEP) where T denotes the total propagation time. We can justify this time-average by the small one-body density deformations of the majority species over time. Furthermore, we remark that Eq. 6 is a species mean-field effective potential and, therefore, only assumes a single product state. Even though we have seen in Figure 5 that the entanglement between the subsystems plays a crucial role this ansatz turns out to be a powerful tool to analyze the basic aspects of the tunneling behavior of the impurity and gives an intuitive understanding [12,46,91]. The time-averaged effective potentials are depicted in Figure 7 (a) for the four interspecies interaction strengths g AB /E r k −1 0 = 0.2, 1.0, 2.0, 4.0 corresponding to the four tunneling regimes already discussed in section III. For weak interspecies interaction strengths, e.g. g AB /E r k −1 0 = 0.2, the TAEP resembles the shape of the double well potential V B , and only small deviations are visible raised from the superimposed one-body density of the majority species. Moreover, the height of the central potential barrier of the TAEP (≈ 0.89E r ) is approximately the same as for the double well V B with a value of h/( √ 2πw) ≈ 0.80 at x B = 0. Therefore, one can assume that the tunneling dynamics of the impurity in the TAEP would differ only marginally compared to the behavior of a single particle confined in the double well V B , i.e. for g AB = 0. Since in the TAEP the one-body density ρ (1) A of the majority species is weighted by the interspecies interaction strength g AB , the spatial distribution of ρ (1) A becomes more pronounced for increasing g AB . Consequently, for a larger g AB we observe the appearance of six maxima on top of the double well structure, which stem from the majority species trapped in the minima of the lattice potential V A . We find that for an interspecies interaction strength of g AB /E r k −1 0 = 1.0 the density of the majority species distributes such that we obtain a nearly parity symmetric TAEP with respect to x = 0 [cf. Figure 7 (d)].
In the following, we refer to the first, second, etc. maximum of the TAEP ordered from left to right excluding the case of g AB /E r k −1 0 = 0.2 due to the small deviations of the TAEP from the double well structure. Increasing the interspecies interaction strength eventually breaks the spatial symmetry of the TAEP w.r.t. x = 0, which can be readily seen e.g. at g AB /E r k −1 0 = 2.0 and g AB /E r k −1 Figure 7 (e), (f)]. In particular, the TAEP for g AB /E r k −1 0 = 4.0 exhibits a distinct asymmetry. Here, the second maximum of the TAEP is strongly suppressed compared to the other maxima. This can be attributed to the fact that the one-body density of the majority species ρ A is strongly depopulated compared to the other wells of V A [cf. Figure 2 (h)] which leads to the suppression of the second maximum of the TAEP. Focusing on the maxima of the TAEP especially, on the third and fourth maximum, we can interpret these maxima as the potential barriers that the impurity has to overcome in order to tunnel from the left to the right side of the TAEP. We observe that the corresponding maxima heights increase with increasing g AB , which indicates that the effective potential barrier for the impurity also increases with g AB .
In the following, we investigate this effective potential barrier that separates the left from the right side of the TAEP. For this purpose, we calculate the relative difference ∆ max,i = (Λ eff i − Λ eff 2 )/Λ eff 2 between the second maximum height Λ eff 2 and third and fourth maximum height, Λ eff 3 and Λ eff 4 , of the TAEP, where i = 3, 4. Figure 7 (b) shows the relative difference ∆ max,i , which serves in the following as a measure for the effective potential barrier, in dependence on the interspecies interaction strength g AB . For values below g AB /E r k −1 0 = 1.0 the effective potential barrier ∆ max,i ≈ 0 for i = 3, 4 meaning that the maxima barely deviate. For values above g AB /E r k −1 0 = 1.0 the third and fourth maximum height of the TAEP become larger than the second one, breaking in this manner the symmetry of the TAEP. For large interspecies interaction strengths, e.g. g AB /E r k −1 0 = 4.0, the effective potential barrier ∆ max,i abruptly increases which is also due to the absence of ρ (1) A in the second well of the lattice potential and, subsequently, the lowering of the second maximum height of the TAEP. The abrupt increase of the effective potential barrier ∆ max,i intuitively leads to the assumption that for large interspecies interaction strengths, e.g. g AB /E r k −1 0 = 4.0, a tunneling of the impurity from the left to the right TAEP should be strongly suppressed, as already seen in the one-body density [cf. Figure 2 (d)] obtained within the complete many-body approach.
Let us also describe the tunneling behavior of the impurity in terms of states that are highly localized in the minima of the TAEP. Since, for sufficiently strong interspecies interaction strengths g AB the TAEP resembles a lattice with five sites [cf. Figure 7 (a)] we calculate five of those functions. For this purpose, we construct an effective HamiltonianĤ eff (x B ) using the TAEP. For the localized functions we use the notion of generalized Wannier functions [89,90] which have the advantage that they can be also obtained in the presence of a non-periodic potential. To obtain five generalized Wannier functions w   Figure 7 (c)] we find, compared to the results for larger g AB , the largest overlap between the Wannier functions. This indicates that those Wannier functions are rather ill-defined. The reason can be found in the TAEP which resembles for small interspecies interaction strengths, e.g. g AB /E r k −1 0 = 0.2, more a double well than a lattice potential. Increasing g AB [cf. Figure 7 (d)-(f)], the Wannier functions become more localized in the minima of the TAEP and, therefore, are suited for the further analysis of the many-body wave function.
In summary, we have developed an effective one-body Hamiltonian using a TAEP for the purpose of constructing generalized Wannier functions from the eigenfunctions of this effective one-body Hamiltonian. This procedure resulted in functions which are highly localized in the wells of the TAEP for sufficiently large interspecies interaction strengths g AB . Projecting these Wannier functions onto the full many-body wave function, obtained in the course of our numerical simulations, we are in the following able to get a deeper insight into the tunneling dynamics of the impurity.

B. Dynamical response in terms of Wannier states
In the following discussion, we analyze the results of the correlated many-body calculations utilizing the Wannier functions derived in the previous section more specifically. Therefore, we project the i-th Wannier function w (B) i onto the manybody wave function and thereby receive the timedependent probability P 1,B (w ) that the impurity occupies the i-th well of the TAEP. This probability is defined as Before discussing the results, let us comment on the Wannier functions as a basis set for the species wave function of the impurity. By summing up the five occupation probabilities of the Wannier functions, we obtain for each time instant of the evolution a measure for the accuracy of this basis representation. For our results we find that ) > 0.95 so that we can consider the Wannier functions as an adequate basis set for describing the tunneling dynamics of the impurity. The time evolution of the probability P 1,B (w  Figure 8. In the following, we aim at understanding the respective one-body densities ρ right side of the TAEP is more uniform than at g AB /E r k −1 0 = 0.2. Initially, we observe an exchange of probability between the first and second Wannier states [cf. Figure 8 (b1)], which represents the material barrier tunneling process in the initial site of the double well. Eventually, probability is transferred from the first and second to the fourth and fifth Wannier state [cf. Figure 8 (b2)], reflecting the controlled tunneling behavior observed in Figure 2 (b). Partially, this can be understood in terms of the TAEP which exhibits a lattice structure on top of the double well, while being still spatially symmetric with respect to x = 0. For stronger interspecies interaction strengths, e.g. g AB /E r k −1 0 = 2.0 and g AB /E r k −1 0 = 4.0, the first two Wannier functions are predominantly occupied. As shown in Figure 8 (c), for times directly after the quench the first Wannier function is the most occupied one, whereas for larger times the occupation probability for the second Wannier function becomes the dominant one. We can understand this intuitively by inspecting the corresponding TAEP depicted in Figure 7 (e). Here, the TAEP exhibits a global minimum in the second well which makes it energetically favorable for the impurity to reside there. Furthermore, the first and second occupation probability exhibit a strong counterwise oscillation behavior. Comparing this with the oscillation of the one-body density of the impurity around a potential barrier imposed by the majority species [see Figure 2 (c)], we find very good agreement. Moreover, we find that the occupation probability for the other Wannier functions, i.e. the third fourth and fifth, are strongly suppressed, but seem to continuously increase over time. Therefore one might assume that for sufficiently long evolution times the impurity eventually tunnels to the right side of the TAEP.
Turning to the occupation probabilities at g AB /E r k −1 0 = 4.0 [ Figure 8 (d)], we observe that the occupation probabilities of the third, fourth and fifth Wannier function are close to zero during the entire time propagation. We can understand this in terms of an effective potential barrier that separates the left side of the TAEP from the right side w.r.t. x = 0. Here, the value for the relative difference ∆ max,i for g AB /E r k −1 0 = 4.0 is much larger than for the other considered interspecies interaction strength g AB [see Figure 6 (b)] and, therefore, it is very unlikely that the impurity will tunnel to the right side of the TAEP even for later times. Additionally, the occupation probabilities of the first and second Wannier function perform a counterwise oscillation, likewise to the probabilities at g AB /E r k −1 0 = 2.0 [ Figure 8 (c)]. In contrast to the aforementioned oscillation of the occupation probabilities, we observe for g AB /E r k −1 0 = 4.0 almost twice the number of oscillation periods during the dynamics as well as a reduction of the amplitudes. A hint for understanding this gives again the corresponding TAEP [see Figure 6 (f)]. Here, the second maximum height is strongly suppressed and the first and second Wannier functions have a large overlap. From this we can infer that the species wave function of the impurity consists mainly of a superposition of the first and second Wannier function whose contributions oscillate over time.
In conclusion, we generated localized Wannier functions associated with the TAEP and projected them onto the time-dependent many-body wave function. Having this at hand we were able to describe the many-body impurity dynamics in terms of the evolution of the occupation probabilities of the respective Wannier functions. As a natural next step, we shall investigate the two-body correlations between the impurity and the majority species. More precisely, we determine the discrete correlation function associated with the impurity occupying a specific well of the TAEP and a single particle of the majority species occupying a certain well of the lattice potential V A .
First we generate Wannier functions associated with the single-particle Hamiltonian for the majority speciesĤ , following the procedure explained above 3 . Furthermore, we determine the one-body probability for a particle of the majority species to occupy the jth well of the lattice potential V A by projecting the many-body wave function onto the Wannier function w (A) j , thereby constructing the probability P 1,A (w ij |Ψ MB of finding a single particle of the majority species in the i-th well of V A and at the same time the impurity in the j-th well of the TAEP is defined as the expectation value of the following operator with respect to the many-body wave function |Ψ MB . The summation runs over the number of particles of the subsystem A. The discrete twobody correlation function is then given by 3 Note that the Wannier functions are sorted from the left to right regarding the sites of the lattice potential.
This function provides information about the correlation between a single particle of the majority species localized at the i-th well of V A and the impurity species localized at the j-th well of the TAEP. In case the discrete correlation function g 2,AB (w j ) equals unity the particle of the majority species and the impurity are termed uncorrelated since the conditional probability equals the unconditional one. However, if g 2,AB (w j ) is larger (smaller) than unity the impurity and the particle of the majority species are said to be correlated (anti-correlated) [15,73].
The time evolution of the discrete correlation function for the impurity occupying the second Wannier function w (B) 2 of the TAEP and the particle of the majority species occupying one of the six Wannier functions of the lattice potential V A are shown in Figure 9. We choose to present only results for g 2,AB where the impurity is occupying the second Wannier state w  Figure 9 (a)], we observe that initially the system is rather uncorrelated and develops stronger correlations for larger time. Here, the particle of the majority species in the fourth/fifth well of V A and the impurity species in the second well of V B eff show an increasing correlation amplitude over time, whereas an anti-correlation between the particle of the majority species in the second/third well of V A and the impurity species in the second well of V B eff occurs. The particle of the majority species being in the first and sixth well exhibits a similar correlation behavior with the impurity in the second well. We note that the strongest correlation (anticorrelation) occurs within the time intervals where the impurity tunnels to the right side of the TAEP [cf. Figure 2 (b)].
Turning to the discrete correlation functions at g AB /E r k −1 0 = 2.0, 4.0, depicted in Figure 9 (b1), (b2) and (c1), (c2) we find that the system is in both cases slightly two-body correlated. Predominantly, this is the case for g AB /E r k −1 0 = 4.0 where the impurity exhibits a self-trapping behavior, showing only a weak distinct correlated or anti-correlated behavior between the impurity in the second well of the TAEP and one particle of the majority species in a specific well of the lattice potential [cf. Figure 9 (c1), (c2)]. However, in Figure 9 (b1) we observe an oscillating correlation behavior between the majority species in the third well of V A and the impurity species in the second well of V B eff . Here, g 2,AB (w 2 ) oscillates between the anti-correlated and uncorrelated case which might be associated with the material barrier tunneling of the impurity in the initial site of the double well [cf. Figure 2 (c)]. Concluding, we observe an overall decrease of the discrete correlation with increasing interspecies interaction strength, which appears to be related to the manifestation of the self-trapping of the impurity.

V. CORRELATED TUNNELING DYNAMICS OF TWO IMPURITIES
So far we have investigated the tunneling dynamics of a single impurity coupled to a majority species and found that the tunneling behavior can be steered by varying the interspecies interaction strength g AB . In the following we unravel whether a similar controlled tunneling process can be realized employing two impurities. For this purpose, we consider two weakly interacting impurities, i.e. they interact repulsively via a contact potential of strength g BB /E r k −1 0 = 0.2, embedded in the same potential setup as shown in Figure 1. The parameters for the majority species remain unchanged compared to the single impurity case, i.e. N A = 8 and g AA /E r k −1 0 = 1.0. The tunneling process is induced by performing the same quench protocol as in the case of a single impurity. Due to a tilting potential V tilt the minority species is initially prepared in the left side of the double well V B . Setting V tilt suddenly to zero we monitor the respective tunneling dynamics. Figure 10 (a)-(d) presents the one-body densities ρ (1) B (x, t) of the two impurities for varying interspecies interaction strengths. As it can be seen, the dynamics of ρ  . For weak interspecies interaction strengths, i.e. g AB /E r k −1 0 = 0.2, we again observe a rather irregular tunneling dynamics of the impurity species [cf. Figure 10 (a)]. Increasing g AB , also for N B = 2 the impurity species performs a material barrier tunneling within the initially populated well and finally tunnels to the other well [cf. Figure  10 (b)]. A further increase of g AB finally yields a self-trapping behaviour of the impurity species due to the strong repulsion [cf. Figure 10 (d)]. As a result also a system with two impurities exhibits the aforementioned four tunneling regimes. However, introducing an additional impurity to the system leads to a shift of the tunneling regimes [cf. Figure 10 (b), (c)] with respect to the interspecies interaction strength (g AB /E r k −1 0 = 0.2, 0.8, 1.3, 4.0). Therefore, we can conclude that also for two impurities we are able to control the quench induced tunneling process via the coupling to the majority species.
Having identified the existence of the four tunneling regimes of the impurity species the question that arises is whether the impurities tunnel pairwise through the potential landscape, which we would refer to as pair tunneling, or whether they tunnel individually, which we would call single particle tunneling. To expose the underlying mechanisms we present in Figure 10  . Thus, we focus on the two cases where a material barrier tunneling of the impurity species takes place or where the latter process is accompanied by a subsequent tunneling to the other site of the double well. In Figure 10 (e) we observe a decreasing probability P 2,BB (w ) to find both impurities in the second Wannier state, whereas the probability to detect both impurities in the fourth Wannier state increases. This probability transfer indicates a pair tunneling from the left to the right side of the TAEP w.r.t. x = 0. Additionally, the impurities perform the material barrier tunneling as a pair, which can be inferred from the alternating increase and decrease of P 2,BB (w 2 ). A schematic representation of this process is illustrated in Figure 10 (i) assuming the TAEP at g AB /E r k −1 0 = 0.8. For the investigation of the single particle tunneling we show in Figure 10 (g), (h) the conditional probability to find one impurity in the second Wannier state and the other impurity in another Wan- increases. This suggests a tunneling of one impurity from the first well of the TAEP to the fourth well, whereas the other impurity remains in the second well. This process is depicted in Figure 10 (j). We note that this is one of many single particle tunneling processes that can take place. In this sense, the tunneling process of the impurity species is rather complex, consisting of single particle and pair tunneling processes.
Concluding we have realized the four tunneling regimes which we previously identified in Figure 2 (a)-(d) also for two weakly interacting impurities coupled to a majority species. This implies that it is also possible to control the tunneling process of two impurities via the interspecies interaction strength. Eventually, we have characterized the tunneling processes underlying the dynamical response of the impurity species in terms of single particle and pair tunneling processes [55].

VI. CONCLUSIONS AND OUTLOOK
We have investigated the correlated tunneling dynamics of impurities trapped in a double well potential and immersed in a lattice trapped majority species. The tunneling dynamics was initiated by implementing an initial tilt of the double well, thereby localizing the impurity species in one of the wells, and quenching this to a symmetric potential configuration. In case of a single impurity we have identified four different tunneling regimes w.r.t. the interspecies interaction strength. For very weak interspecies interaction strengths the tunneling of the impurity can be characterized as rather complex, exhibiting no regular or repetitive structure. However, increasing the coupling to the majority species leads to a regular tunneling behavior of the impurity, which consists of an initial material barrier tunneling due to the presence of the majority species and is followed by a transfer of the impurity to the other site of the double well. Additionally, this effect is accompanied by a strong entanglement between the subsystems. A further increase of the interspecies interaction strength leads to a sole material barrier tunneling in the initial site of the double well for long time intervals and finally for very large couplings forces the impurity to localize in the initially populated well and being self-trapped.
In order to gain insight into the underlying microscopic processes of the emergent correlated tunneling dynamics, we have constructed a timeaveraged effective potential (TAEP) based on the one-body density of the majority species. Depending on the interspecies interaction strength, this effective potential exhibits an additional structure in each site of the double well, thus explaining the material barrier tunneling. Increasing the coupling to the majority species, the TAEP is predominantly formed by the one-body density of the majority species and the presence of the double well is of minor consequence, resulting in the observed self-trapping of the impurity. Moreover, the generalized Wannier states associated with this potential allowed for a characterization of the impurity's dynamical response as well as the involved correlations. We concluded our study with an investigation of two weakly repulsively interacting impurities which we prepared analogously to the case of a single impurity. We were able to identify the previous four tunneling regimes for smaller interspecies interaction strengths, being shifted to g AB /E r k −1 0 = 0.2, 0.8, 1.3, 4.0 respectively, compared to the scenario of a single impurity. Employing again the TAEP we have developed an understanding of the tunneling dynamics, which consists of a superposition of pair tunneling as well as single particle tunneling processes.
There are various interesting research directions that prove to be promising for future investigations relying on the findings of the current work. A direct extension involves the inclusion of spin degrees of freedom between the impurities. Here, the possible formation of an analogue of a Cooperpair in the course of the tunneling dynamics would be of immediate interest. Another straightforward direction would be to consider quench protocols which also include a variation of the interspecies interaction strengths. For example, one might think of a subsequent interaction quench after a transfer of the impurity in order to prevent tunneling to the initially populated site. Also, dynamically driving the corresponding parameters of the system might be useful for transferring the impurity species in a more controlled and systematic manner. In the following we demonstrate that a certain minimal tilting strength α is necessary for observing the tunneling dynamics as in Figure 2 where α/E R k −1 0 = 0.1. Figure 11 shows the temporal evolution of the one-body densities of the impurity [ Figure 11 (a)-(d)] and the majority species [ Figure 11 (e)-(h)] using a tilting strength α/E R k −1 0 = 0.01, within the full many-body approach. Analogously to the previous discussion in section III, we induce the dynamics by initially tilting the double well V B of the impurity with a tilting strength α and let the system evolve in time for α = 0. However, in the present case lowering the initial tilting strength to α/E R k −1 0 = 0.01 leads to a smaller initial energy offset between the sites of V B .
For weak g AB , i.e. g AB /E r k −1 0 = 0.2, 1.0, we find a rather regular tunneling of the impurity from the left to the right side of the double V B [cf. Figure 11 (a), (b)]. Comparing this with the dynamical response of an impurity for an initial tilting strength of α/E R k −1 Figure 2 (b)] we find no material barrier tunneling triggered by the density of the majority species. The difference between the two initial tilting strengths is also evident for larger interspecies interaction strengths, e.g. g AB /E r k −1 0 = 2.0, 4.0. Here, the impurity essentially remains localized throughout the dynamics and does not perform any oscillations [cf. Figure 11 (c), (d)]. Furthermore, the one-body density of the majority species behaves accordingly and does not exhibit a distinctive dynamics compared to the α/E R k −1 0 = 0.1. Namely, ρ A remains well localized at the sites of the lattice potential during the propagation [cf. Figure 11 (e)-(h)]. These observations lead to the conclusion that indeed a sufficiently high initial tilting strength α is needed in order to observe a material barrier tunneling with a subsequent controlled transfer of the impurity to the other side of the double well.
Appendix B: Dependence of the tunneling process on the system parameters We have analyzed the tunneling behavior of the impurity species for a specific choice of the intraspecies interaction strength g AA of the majority species and the barrier height h. In Figure 12 we show that the qualitative behavior of the tunneling dynamics discussed in the main text can be recovered for significantly varying g AA and h. As a measure for the characteristic dynamical response of the impurity we again investigate the temporal evolution of the integrated one-body density of the impurity  B (x, t)dx. The latter enables us to distinguish between the different tunneling regimes, for a fixed interspecies interaction strength of g AB /E R k −1 0 = 1.0. This value lies in regime II in Figure 4 (a), where we observe a material barrier tunneling within the initially populated well with a final transport of the impurity to the other site of the double well. In Figure 12 (a) we observe that an increase of g AA leads to a faster revival of the material barrier tunneling in the initially populated well. Decreasing g AA leads to a temporal prolongation of the material barrier tunneling and thereby a delayed transfer of the impurity to the other site of the double well.
A similar process can be observed when increasing the height of the double well barrier. However, at a certain height of the barrier, e.g. h/E r k −1 0 = 5.0, the impurity is barely able to tunnel to the other site of the double well within the considered time interval and solely performs the material barrier tunneling in the initial well. In contrast, a sufficiently small barrier height, e.g. h/E r k −1 0 = 1.0, leads to a less dominant material barrier tunneling within the initial well since the impurity can be easily transferred to the other well.
Finally, we have investigated the dependence of the impurity's dynamical response on the initial tilt α [cf. Figure 12 (c)]. For small tilts we find almost no oscillations in the initially populated well and the impurity is directly transferred to the other site of the double well [see also Figure 11 (b)]. For increasing tilts α the oscillations due to the mate-rial barrier tunneling become more prominent and thereby the tunneling to the other well is delayed. In this sense, we find that the material barrier tunneling can be recovered for various parameters of the system.

Appendix C: Convergence of the many-body simulations
In the following we briefly discuss the convergence behavior of our results. As discussed in section II the size of the Hilbert space is given in terms of the orbital configuration C = (M, d A , d B ).  Here, M describes the number of species functions in the Schmidt representation (see Eq. 3), while d σ with σ ∈ {A, B} refer to the number of single-particle functions building the time dependent number states | n σ (t) (see Eq. 4). In the process of increasing the number of species functions and single-particle functions it is possible to recover the solution of the many-body wave function with an increasing accuracy. Due to the exponentially increasing size of the Hilbert space it is computationally prohibitive to use too many species and single-particle functions. How-ever, we are able to obtain numerical solutions which incorporate all the necessary correlations and go beyond mean-field approximations utilizing ML-MCTDHX. We determine the effect of the truncation of the Hilbert space by investigating as a representative example the integrated one-body density of the impurity −L/4 −L/2 ρ (1) B (x, t)dx upon varying the orbital configuration C. In Figure 13 we show the latter for an interspecies interaction strength of g AB /E R k −1 0 = 1.0 [cf. Figure 2 (b)]. Note that g AB /E R k −1 0 = 1.0 lies in the interval where the degree of correlations is maximized (cf. Figure 5). As it can be seen, increasing the size of the Hilbert space systematically it is possible to achieve convergence. Based on these findings the orbital configuration C = (6, 6, 6) has been employed in all many-body calculations presented in the main text. F. T. and K. K. contributed equally to this work.