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

Our computation approach is the ab-initio multi-layer multi-configuration time-dependent Hartree method for bosonic (fermionic) mixtures (ML-MCTDHX) [69–71], which accounts for all the relevant correlations of the atomic mixture [50, 72, 75–78]. 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 [79]

$$\begin{equation}\vert {{\Psi}}_{\text{MB}}\left(t\right)\rangle =\sum _{i=1}^{M}\sqrt{{\lambda }_{i}\left(t\right)}\vert {{\Psi}}_{i}^{A}\left(t\right)\rangle \otimes \vert {{\Psi}}_{i}^{B}\left(t\right)\rangle .\end{equation} \tag{ 1 }$$

Here, the Schmidt coefficients $\sqrt{{\lambda }_{i}}$, in decreasing order, provide information about the degree of population of the ith species function and thereby determine the degree of entanglement between the impurities and the majority species. In case that λ₁ = 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

$$\begin{equation}\vert {{\Psi}}_{i}^{\sigma }\left(t\right)\rangle =\sum _{{\to {n}}^{\sigma }\vert {N}_{\sigma }}{C}_{\sigma \to {n}}\left(t\right)\vert {\to {n}}^{\sigma };t\rangle ,\end{equation} \tag{ 2 }$$

where the vector ${\to {n}}^{\sigma }=\left({n}_{1}^{\sigma },\enspace {n}_{2}^{\sigma },\dots \right)$ denotes the occupations of the time-dependent single-particle functions of the species σ. The notation ${\to {n}}^{\sigma }\vert {N}_{\sigma }$ indicates that for each $\vert {\to {n}}^{\sigma };t\rangle$ we require the condition ${\sum }_{i}{n}_{i}^{\sigma }={N}_{\sigma }$. The time propagation of the many-body wave function is achieved by employing the Dirac–Frenkel variation principle $\langle \delta {{\Psi}}_{\text{MB}}\vert \left(\mathrm{i}{\partial }_{t}-\mathcal{H}\right)\vert {{\Psi}}_{\text{MB}}\rangle$ [80–82] 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 [83, 84]. Besides investigating the quantum dynamics it allows us to determine the ground (or excited) states by using either imaginary time propagation or improved relaxation [85], thereby being able to uncover also possible degeneracies of the many-body states. We remark that in commonly used approaches for solving the time-dependent 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 [86].

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 single-particle functions spanning the time-dependent number states $\vert {\to {n}}^{\sigma };t\rangle$ (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.

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 (figure 1). The impurity species with N_B = 1, 2 particles is denoted as species B and the majority species containing N_A = 4 (main text) or N_A = 5 (appendix A) particles is referred to as species A. This setup lies within reach of current experimental techniques [87, 88]. 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}, reads

$$\begin{align}\hfill {\hat{H}}_{\sigma }& ={\int }_{-L/2}^{L/2}\mathrm{d}\text{x}\enspace {\hat{{\Psi}}}_{\sigma }^{{\dagger}}\left(\text{x}\right)\left[-\frac{{\hslash }^{2}}{2{m}_{\sigma }}\frac{{\mathrm{d}}^{2}}{{\mathrm{d}\text{x}}^{2}}+{V}_{0}\enspace {\mathrm{sin}}^{2}\left(\frac{\pi k\text{x}}{L}\right)+{g}_{\sigma \sigma }\enspace {\hat{{\Psi}}}_{\sigma }^{{\dagger}}\left(\text{x}\right){\hat{{\Psi}}}_{\sigma }\left(\text{x}\right)\right]{\hat{{\Psi}}}_{\sigma }\left(\text{x}\right),\hfill \end{align} \tag{ 3 }$$

where ${\hat{{\Psi}}}_{\sigma }^{{\dagger}}$ is the field operator of species σ, m_σ their mass and V₀ 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. ⁸⁷Rb atoms in two different hyperfine states [89]. The interaction between the species A and B is given by

$$\begin{equation}{\hat{H}}_{AB}={g}_{AB}{\int }_{-L/2}^{L/2}\mathrm{d}\text{x}\enspace {\hat{{\Psi}}}_{A}^{{\dagger}}\left(\text{x}\right){\hat{{\Psi}}}_{A}\left(\text{x}\right){\hat{{\Psi}}}_{B}^{{\dagger}}\left(\text{x}\right){\hat{{\Psi}}}_{B}\left(\text{x}\right),\end{equation} \tag{ 4 }$$

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}_{\alpha }^{3D}$. By assuming the above-mentioned strong transversal confinement with the same trapping frequencies ${\omega }_{\perp }^{\sigma }={\omega }_{\perp }$ for both species σ ∈ {A, B} it is possible to integrate out frozen degrees of freedom, leading to a quasi one-dimensional model with ${g}_{\alpha }=2\hslash {\omega }_{\perp }{a}_{\alpha }^{3D}$.

Zoom In Zoom Out Reset image size

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πℏ)²/2m_Aλ² is the recoil energy and λ = 2L/k the optical lattice wavelength.

In the following, we explore the ground state of the system containing a single impurity, i.e. N_B = 1, for varying V₀ 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

$$\begin{equation}{\rho }_{\sigma }^{\left(1\right)}\left(x\right)=\langle {{\Psi}}_{\text{MB}}\vert {\hat{{\Psi}}}_{\sigma }^{{\dagger}}\left(x\right){\hat{{\Psi}}}_{\sigma }\left(x\right)\vert {{\Psi}}_{\text{MB}}\rangle .\end{equation} \tag{ 5 }$$

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

$$\begin{equation}{S}_{AB}=-\sum _{i}{\lambda }_{i}\enspace \mathrm{ln}\left({\lambda }_{i}\right)\end{equation} \tag{ 6 }$$

as a measure for the entanglement between the subsystems A and B, where λ_i are the Schmidt coefficients defined 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₀. 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₀/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.

Zoom In Zoom Out Reset image size

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 figures 2(c) and (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 figures 2(c) and (d) and its parity-symmetric (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₀/E_R = 10.5 we observe a minor population of the impurity in the outer wells. However, this spatial species distribution in ${\rho }_{B}^{\left(1\right)}\left(x\right)$ is more pronounced for smaller values of V₀. This means that the corresponding one-body density is increased in the outer wells. Indeed, for V₀/E_R = 7.3 and g_AB/E_Rλ = 0.047 (as compared to V₀/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 [90–93] 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 [73], 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

$$\begin{equation}{\rho }_{\sigma {\sigma }^{\prime }}^{\left(2\right)}\left({x}_{1}^{\sigma },{x}_{2}^{{\sigma }^{\prime }}\right)=\langle {{\Psi}}_{\text{MB}}\vert {\hat{{\Psi}}}_{\sigma }^{{\dagger}}\left({x}_{1}^{\sigma }\right){\hat{{\Psi}}}_{{\sigma }^{\prime }}^{{\dagger}}\left({x}_{2}^{{\sigma }^{\prime }}\right){\hat{{\Psi}}}_{\sigma }\left({x}_{1}^{\sigma }\right){\hat{{\Psi}}}_{{\sigma }^{\prime }}\left({x}_{2}^{{\sigma }^{\prime }}\right)\vert {{\Psi}}_{\text{MB}}\rangle ,\end{equation} \tag{ 7 }$$

where σ, σ' ∈ {A, B}. This measure corresponds to the probability of finding a particle of species σ at the position ${x}_{1}^{\sigma }$ and another particle of species σ' at the position ${x}_{2}^{{\sigma }^{\prime }}$.³ Figure 3(a₁) 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₀/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 two-body density ${\rho }_{AB}^{\left(2\right)}$ (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. ${\rho }_{AB}^{\left(2\right)}\left(x,x\right)\approx 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₁) we find that the probability of detecting two particles of species A at the same position is approximately zero, i.e. ${\rho }_{AA}^{\left(2\right)}\left(x,x\right)\approx 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

$$\begin{equation}\vert {{\Psi}}_{\text{MB}}\rangle \approx \vert 1,1,1,1,0{\rangle }_{A}\otimes \vert 0,0,0,0,1{\rangle }_{B}.\end{equation} \tag{ 8 }$$

Here, the number state $\vert {\to {n}}^{\sigma }\rangle =\vert {n}_{1}^{\sigma },{n}_{2}^{\sigma },{n}_{3}^{\sigma },{n}_{4}^{\sigma },{n}_{5}^{\sigma }{\rangle }_{\sigma }$ is constructed with a generalized Wannier basis of the lowest band [94, 95]. This means that e.g. ${n}_{1}^{\sigma }$ 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.

Zoom In Zoom Out Reset image size

In order to unravel the role of correlations for the ground state distribution in the strong interaction regime, we calculate the noise correlation [64] between particles of species σ and σ' which is defined as

$$\begin{equation}{g}_{\sigma {\sigma }^{\prime }}^{\left(2\right)}\left({x}_{1}^{\sigma },{x}_{2}^{{\sigma }^{\prime }}\right)={\rho }_{\sigma {\sigma }^{\prime }}^{\left(2\right)}\left({x}_{1}^{\sigma },{x}_{2}^{{\sigma }^{\prime }}\right)-{\rho }_{\sigma }^{\left(1\right)}\left({x}_{1}^{\sigma }\right){\rho }_{\sigma }^{\left(1\right)}\left({x}_{2}^{{\sigma }^{\prime }}\right).\end{equation} \tag{ 9 }$$

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}_{\sigma {\sigma }^{\prime }}^{\left(2\right)}=0$ and in the latter case ${g}_{\sigma {\sigma }^{\prime }}^{\left(2\right)}\ne 0$. In figure 3(a₂) we show ${g}_{AA}^{\left(2\right)}=0$ corresponding to the previously discussed two-body density for two particles of the A species for V₀/E_R = 10.5 and g_AB/E_Rλ = 0.12. Let us remark that in case of ${g}_{\sigma {\sigma }^{\prime }}^{\left(2\right)}{< }1{0}^{-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 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}_{AA}^{\left(2\right)}$, 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.

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₀, 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 figures 4(e) and (f)), similar to the case of a single impurity. However, for smaller values of the lattice depth (e.g. V₀/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₀/E_R = 5.7 and V₀/E_R = 8.9 and large g_AB (cf figures 4(c)–(f)) are two-fold degenerate with a parity-symmetric counterpart (with respect to x = 0).

Zoom In Zoom Out Reset image size

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

$$\begin{equation}\langle \vert {x}_{1}-{x}_{2}\vert \rangle ={\int }_{-L/2}^{L/2}{\int }_{-L/2}^{L/2}\mathrm{d}{x}_{1}\mathrm{d}{x}_{2}\vert {x}_{1}-{x}_{2}\vert {\rho }_{BB}^{\left(2\right)}\left({x}_{1},{x}_{2}\right),\end{equation} \tag{ 10 }$$

with ${\rho }_{BB}^{\left(2\right)}\left({x}_{1},{x}_{2}\right)$ being the previously introduced two-body 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₀ the distance of the impurities decreases and saturates to a specific value for even larger g_AB. Also, for smaller V₀ 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₁ − x₂|⟩ 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₀. 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–20, 74]. 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₁)) 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 ${\rho }_{BB}^{\left(2\right)}$ of finding two impurities at specific positions x₁ and x₂ accumulates at a single outer site of the lattice geometry (see figure 5(e₁)]. 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 ${\rho }_{AB}^{\left(2\right)}$ (cf figure 5(f₁)]. 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₁)). The latter fact can also be observed in the two-body density ${\rho }_{AB}^{\left(2\right)}$ for two particles of different species (see figures 5(c₁) and (f₁)). 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(a)).

Zoom In Zoom Out Reset image size

Coming back to the question of how the individual impurities distribute, in figure 5(b₁) 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 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₀/E_R = 5.7, V₀/E_R = 8.9 and g_AB/E_Rλ = 0.12. For V₀/E_R = 8.9 the structure of the noise-correlation for two particles of the A species (figure 5(d₂)) is similar to the one for a single impurity in figure 3(a₂) 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 4} (see figure 5(e₂)). 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₀/E_R = 5.7 and g_AB/E_Rλ = 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₂)). This is also the case for finding the impurities at different sites. Also, for two particles of the majority species we observe that ${g}_{AA}^{\left(2\right)}\ne 0$ in the relevant occupied sites (see figure 5(a₂)), with a structure similar to figure 5(d₂). Therefore, we can conclude that for V₀/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₀/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.

In the following we prepare our system in its ground state with a lattice depth V₀/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 figures 2(c) and (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 entanglement crossover diagram shown in figure 2(a).

As a representative example of the emergent tunneling dynamics of each species in figures 6(a) and (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₀/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/(ℏ/E_R) ≈ 600 a hole can again be found at the opposite outer well where now the impurity resides (see figure 6(b)).

Zoom In Zoom Out Reset image size

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

$$\begin{equation}{\hat{O}}_{l}^{\left(1\right)}=\vert {w}_{l}^{B}\rangle \langle {w}_{l}^{B}\vert ,\end{equation} \tag{ 11 }$$

where $\vert {w}_{l}^{B}\rangle \langle {w}_{l}^{B}\vert$ projects the particle of the B species onto the lth Wannier state of the lowest band. Evaluating this operator with respect to the complete many-body wave function yields the probability ${P}_{l}=\langle {{\Psi}}_{\text{MB}}\vert {\hat{O}}_{l}^{\left(1\right)}\vert {{\Psi}}_{\text{MB}}\rangle$ of detecting an impurity particle in the lth Wannier state. Note that the complete many-body wave function is obtained via ML-MCTDHX and subsequently we analyze this high-dimensional object by evaluating ${\hat{O}}_{l}^{\left(1\right)}$ with respect to the wave function. In the following the Wannier states are ordered from left to right, i.e. |w₅⟩ 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 ∑_lP_l ≈ 99.9%. Figure 7(b), shows the probability P_l of finding the impurity in the lth 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 σ

$$\begin{equation}\langle {\hat{H}}_{\sigma }^{\mathrm{r}\mathrm{e}\mathrm{l}}\rangle \left(t\right)=\langle {\hat{H}}_{\sigma }\rangle \left(t\right)-\langle {\hat{H}}_{\sigma }\rangle \left(t=0\right),\end{equation} \tag{ 12 }$$

which we define as the deviation of the expectation values of the individual Hamiltonians Ĥ_σ (with σ ∈ {A, B, AB}) deviating from their initial value at t = 0. The energy of the impurity first decreases during the transport to the opposite outer well and then increases again (see figure 7(c)). This can be easily understood in terms of the population of the corresponding Wannier states. The energy of the three central Wannier states is smaller as compared to the two outer ones (see figure 7(a)). Since the former are strongly populated in a superposition in the course of time, we observe a decrease of the impurity energy. The revival of the energy is in turn associated with the population of the other outer well which has again an increased energy (due to the hard-wall boundary conditions). The reduction of the impurity energy is accompanied by an increase of the energy of the majority species and the interspecies energy. The reason for the latter will be discussed below.

Zoom In Zoom Out Reset image size

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 figures 6(a) and (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

$$\begin{equation}{\rho }_{\sigma }^{\left(1\right)}\left(x,t\right)=\sum _{j}{n}_{\sigma j}\left(t\right){{\Phi}}_{\sigma j}^{{\ast}}\left(x,t\right){{\Phi}}_{\sigma j}\left(x,t\right),\end{equation} \tag{ 13 }$$

where n_σj(t) in decreasing order, obeying ∑_jn_σ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 [69, 70, 85], 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, 77, 78] in the subsystem σ as

$$\begin{equation}{S}_{\sigma }\left(t\right)=-\sum _{j}{n}_{\sigma j}\left(t\right)\enspace \mathrm{ln}\left({n}_{\sigma j}\left(t\right)\right).\end{equation} \tag{ 14 }$$

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

$$\begin{equation}\langle {\hat{X}}_{B}\rangle ={\int }_{-L/2}^{L/2}\mathrm{d}x{\rho }_{B}^{\left(1\right)}\left(x\right)x.\end{equation} \tag{ 15 }$$

The initial value of the impurity's center of mass is given by $\langle {\hat{X}}_{B}\rangle /\lambda \approx 1$, while the final position when the impurity is completely transported to the opposite outer well is $\langle {\hat{X}}_{B}\rangle /\lambda \approx -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 ${\rho }_{AB}^{\left(2\right)}$ and ${\rho }_{AA}^{\left(2\right)}$. 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 $\langle {\hat{X}}_{B}\rangle /\lambda \approx -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.

Zoom In Zoom Out Reset image size

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 expectation value of the operator

$$\begin{equation}{\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}={\mathbb{1}}_{A}\otimes {\mathbb{1}}_{B}-\frac{1}{{N}_{A}{N}_{B}}\sum _{l=1}^{5}\vert {w}_{l}^{B}\rangle \langle {w}_{l}^{B}\vert \otimes \sum _{i}^{{N}_{A}}\vert {w}_{l}^{iA}\rangle \langle {w}_{l}^{iA}\vert ,\end{equation} \tag{ 16 }$$

with respect to the system's wave function, i.e. $\langle {\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}\rangle$. Here, $\vert {w}_{l}^{B}\rangle \langle {w}_{l}^{B}\vert$ projects the particle of the B species onto the lth Wannier state and $\vert {w}_{l}^{iA}\rangle \langle {w}_{l}^{iA}\vert$ projects the ith particle of the A species onto the lth Wannier state. ${\mathbb{1}}_{\sigma }$ 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 $\langle {\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}\rangle \left(t=0\right)\approx 98.8\%$ holds. This means that we can safely assume that our system exhibits a particle–hole pair.

Zoom In Zoom Out Reset image size

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 $\langle {\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}\rangle$ 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 $\langle {\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}\rangle \approx 93.2\%$. Moreover, we observe that $\langle {\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}\rangle$ 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 $\langle {\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}\rangle$ 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.

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₀/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 figures 4(e) and (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 $\langle {\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}\rangle$ during the evolution. Initially, the particle hole-pair process is the dominant tunneling channel with a probability of $\langle {\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}\rangle \approx 98.6\%$ which later on drastically decreases to $\langle {\hat{O}}_{\mathrm{p}\mathrm{h}}^{\left(2\right)}\rangle \approx 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(d)). Furthermore the initially uncorrelated impurity pair, S_B(t = 0) ≈ 0, develops correlations during the propagation to the opposite outer well which can be attributed to the incomplete transfer of the impurities to this very site (see figure 10(e)). Similarly, the increase of S_A is attributed to the incomplete transfer of the effective hole into the opposite direction (as compared to the propagation direction of the impurities) (see figure 10(f)).

Zoom In Zoom Out Reset image size

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 ${\rho }_{BB}^{\left(2\right)}$ upon quenching to g_AB/E_Rλ = 0.002 (see figures 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(h)). However, during the propagation there are time instances 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 figures 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.