Interspecies entanglement with impurity atoms in a lattice gas

The dynamics of impurity atoms introduced into bosonic gases in an optical lattice have generated a lot of recent interest, both in theory and experiment. We investigate to what extent measurements on either the impurity species or the majority species in these systems are affected by their interspecies entanglement. This arises naturally in the dynamics and plays an important role when we measure only one species. We explore the corresponding effects in strongly interacting regimes, using a combination of few-particle analytical calculations and density matrix renormalisation group methods in one dimension. We identify how the resulting effects on impurities can be used to probe the many-body states of the majority species, and separately ask how to enter regimes where this entanglement is small, so that the impurities can be used as probes that do not significantly affect the majority species. The results are accessible in current experiments, and provide important considerations for the measurement of complex systems with using few probe atoms.

At the same time, entanglement in many-body systems [44,45] has generated a lot of interest, especially because of the information that can be extracted from entanglement in spatial modes. This is helpful in understanding a variety of phenomena, e.g., to identify topological phases [46][47][48][49] or understand the growth of local entropy during thermalisation [50,51]. Such entanglement in space has also been measured in experiments with cold atoms in optical lattices [52][53][54]. For multicomponent gases, the interplay between correlation and entanglement effects has been explored in both the Bose-Hubbard and the Hubbard model [55][56][57][58]. Recently, interspecies entanglement has been used to characterise the shift of the phase transition points in the two-component Bose-Hubbard model due to interspecies interactions [59]. Entanglement of an impurity in a few-body continuous system has also been discussed in terms of a reduced single-particle density matrix in reference [60].
The system we will investigate here is a lattice model for impurities introduced into a Bose gas loaded into the lowest Bloch band of an optical lattice [61][62][63]. This has been realised in experiments [38,64,65], and recently discussed as an important example for characterising the probing of strongly correlated systems with impurity atoms [66]. We are particularly interested in asking about the role of entanglement between the impurity atoms and the majority species, and how this impacts measurements made on a single species alone. Specifically, when momentum distributions of impurity atoms are measured in some two-species experiments [37,38] there is a notable decrease in the visibility of peaks in these momentum distributions, beyond what might be expected from an increase in the effective mass or interactions mediated by the majority species. Entanglement between two species can lead to a mixed reduced density operator for a single species, substantially affecting first-order coherence properties, including measured momentum (or quasi-momentum) distributions. These effects, as we will see below, are strongest when going away from the limit where polaron physics is studied (or indeed, where a usual polaron description is valid), i.e., it is important for strong interactions, or when there is not a large ratio between the atom numbers of the majority and minority species.
As an example, in figure 1 we depict a two-component mixture in a lattice. In figure 1(a), the two species are in a product state with no entanglement, as might be expected to occur in the ground state of a mixture with no interspecies interaction. If we trace over one species and ask what the reduced density operator is that describes the second species alone, we find it represents a pure state, with the atoms still delocalised, and with their momentum peaked at p = 0. However, in figure 1(b), we have an initial state of perfectly correlated dimers between the two species, with dimer momentum peaked at p = 0. This represents an entangled state of the two species, and when we now trace over one species, the reduced density operator for the other species contains a mixture of all possible configurations, with a completely flat momentum distribution.
Below, using numerical and analytical methods, we analyse this behaviour for different parameter regimes of both impurity and majority species particles. We find complex behaviour that exhibits particular signatures associated with the quantum phase diagram of the Bose-Hubbard model for the majority species. In this sense, understanding the impurity-majority species entanglement can be useful as an alternative route to probing the complex many-body behaviour, either by observing the impurity atoms, or by observing their effect on the state of the majority atoms. At the same time, if the impurity atoms are being used as a probe in the sense of references [8,66], it may be important to minimise direct entanglement between the species, and we analyse requirements in simple example cases to achieve this regime.
The remainder of this article is organised as follows: in section 2 we introduce the model and numerical methods we use to analyse the system. We then analyse the entanglement and momentum distributions for two particles on a lattice, one of each species, as a starting point for further investigations in section 3. In section 4, we make use of a Born-Oppenheimer approximation to help us separate, for simple examples, the effects of a change in effective mass of the impurities from effects arising from the entanglement between the impurities and the majority species. In section 5 we analyse the behaviour for systems of multiple atoms on the lattice in parameter regimes corresponding to different quantum phases of the majority species. In section 6 we then ask under which circumstances a single atom will become disentangled from the system to which it is coupled, and to this end we investigate cases where the impurity atoms and the majority species have different tunnelling rates in the lattice, and where the impurity atoms are confined to a fraction of the full length of the system. We then provide a conclusion and outlook in section 7.

Model
We consider an ensemble of bosonic atoms loaded into the lowest band of an optical lattice. We denote the majority species as species 1, with N 1 atoms, and the impurity species as 2, with N 2 N 1 atoms. For sufficiently low temperatures and where interactions are smaller than the energy separation between Bloch bands, this situation is generally well described by a multi-species Bose-Hubbard model ( ≡ 1) [36], where b † σ,i b σ,l and n σ,l are the creation(annihilation) operator and number operator for species σ ∈ {1, 2} on the lth lattice site. Each species has nearest-neighbour tunnelling rate J σ and intra-species onsite interaction U σ . The on-site interspecies interaction energy shift is then denoted U 12 . We will generally take J 1 = J 2 ≡ J unless otherwise specified, and we usually take the same 1D lattice length M. In section 6 we will consider the impurity particles to be confined to a lattice length M 2 M, and denote M 1 as the full length for the majority species.
Throughout this work we will mainly restrict our calculations to one dimension, in order to simplify the computations. However the basic principles we discuss here and the qualitative behaviour of the entanglement in different parameter regimes is expected to transfer directly to higher dimensions.
In what follows, we will use analytical methods to obtain exact results for a few atoms, and exact diagonalisation methods for small lattice sizes, especially to obtain values for the von Neumann entropy of entanglement. If we compute the reduced density matrix ρ σ for either of the two species, where |Ψ denotes the state of the total system, and σ here denotes the opposite species to σ, then the von Neumann entropy of entanglement can be computed as Note that if the total state of the system is pure, then S vN is independent of the choice of σ. The entropy of the reduced density matrix for one species in this case entirely represents the entanglement between the two species. If S vN = 0 the reduced density matrix is a pure state. This occurs when the entanglement is zero, and the total state is a product state of the two species. As a guide to larger system behaviour, we employ mean-field methods based on the bosonic Gutzwiller ansatz [67], where the ground state of the two species Bose-Hubbard Hamiltonian in (1) on an M-site chain is written as, Here f (l) n 1 ,n 2 is the amplitude associated with different number states for each particle on the lth site. To provide additional information on the many-body physics beyond this, we employ density matrix renormalisation group (DMRG) methods based around matrix product states [68,69] to determine the ground state. In each case, we ensure that the results are properly converged in the bond dimension of the matrix product state, D.

Two particles on a 1D lattice
We can obtain an intuition for the behaviour we expect by considering the entanglement of two distinguishable particles in an optical lattice. Using an exact solution for the ground state of equation (1) [70][71][72], we can quantify how the entanglement and momentum distribution for each particle change as a function of the interaction strength between the bosons.
As a useful starting point, we consider the limiting cases. Because the ground state of a single particle on the lattice is a state with quasi-momentum p = 0, for two particles the non-interacting ground state when U 12 = 0 is given by When we compute the reduced density matrix, we then obtain which is a pure state, for which the resulting quasi-momentum distribution n(p) is peaked at p = 0. If, however, the particles are interacting such that |U 12 | J, then for attractive interaction, we obtain instead If we now compute the reduced density matrix, we then obtain which is a mixed state with S vN = log 2 M. This mixed state with particles localised on each site arises in a sense because the state of the second species contains information on the locations of the state of the first species. The resulting momentum distribution is completely flat, despite the fact that the doublon momentum distribution is peaked at p = 0.
To analyse the behaviour for arbitrary interaction strengths we look at the general solution for the two-particle wavefunction with complex coefficients ψ(x, y). Taking periodic boundary conditions, we can separate the centre of mass R = (x + y)/2 and relative r = x − y coordinates, and determine an analytical solution [70][71][72], for which we provide more details in appendix A. We show the analytical calculation for this relative wavefunction in figure 2(a). We see that the peak of the bound pair solutions becomes sharper as the interaction strength is increased. The single particle momentum distributions are shown for comparison in figure 2(b), and show clearly the effect of interactions. As expected for U < 0, with increasing interaction strength the two particles become more tightly bound and this leads to broadening of the single-particle momentum distribution. In this general setting, the entanglement between the two particles causes mixedness in the single species reduced density matrix, implying a reduction in the first-order coherence of this species when measured alone. Therefore, the momentum distribution, although being a single particle observable, is affected by the entanglement between species. In the rest of the article we will use the height of the momentum peak or visibility for each species σ, V σ , as an indicator of the corresponding changes in momentum distributions.
In figure 2(c) we then look directly at the von Neumann entropy of entanglement. For attractive interactions (U 12 < 0), the entanglement grows very sharply as a consequence of the direct pairing of the particles in position space that creates the bound state and reaches the saturation value for small U 12 /J (which is log 2 M as noted above). The entanglement in position space is generated by repulsion (U 12 > 0) between the particles, which makes it energetically unfavourable for them to be present on the same lattice site. This does not provide as strong entanglement as in the attractive case, but still increased with increasing U 12 , towards an asymptotic value. The visibility of the momentum distribution peak V 1 is shown in figure 2(d) and directly mirrors the behaviour of the entanglement, falling very sharply on the attractive side U 12 < 0 as the particles become highly entangled and the momentum distribution for a single particle tends rapidly to a flat distribution. Similarly, on the repulsive side (U 12 > 0), the slight drop of the visibility profile followed by a steady value reflects the corresponding entanglement of the particles.

Comparison of the effects of entanglement with mediated interactions and increased effective mass
In some sense, this decrease in the visibility of a momentum distribution for impurity atoms, when interacting with a majority species, might be expected. After all, the majority species will tend to mediate additional interactions between the minority atoms, and also potentially to increase the effective mass. This has been discussed in the past, e.g., for polarons in a lattice system immersed in a continuous reservoir gas [61,73]. Here we are generally interested in limits where we go away from the usual regimes of polarons, but for certain limiting cases it is possible to directly extract approximate mediated interactions and effective masses for one species of atom, and comparing the visibility profile resulting from this effective model, from the full calculation including the reduction in first-order coherence due to entanglement.
This can be approached in the limit where the species of interest (i.e., impurities) are much heavier than the second species, so that they have much lower tunnelling rate (J 2 J 1 ). In this limit, we can treat the problem in a Born-Oppenheimer (B-O) approximation, which allows us to extract effective interactions as a function of distance. By considering a single impurity confined to two lattice sites, we can also determine effective tunnelling rates that reflect any increased effective mass, allowing us to generate an effective model for the behaviour of the impurities. We will do this first for a single particle of the second species per impurity, and then repeat the calculation for a case where the second species becomes the majority species.
In the B-O approximation, the full Hamiltonian (1) is divided into two parts. The first part is the tunnelling of the impurities, H T = − l,j J 2 b † 2,l b 2,j , which is considered to be very small on the timescale relevant for the dynamics of the second species of atoms. The second part of the Hamitonian then governs the resulting configuration where the impurities are essentially taken to be motionless. With the impurities separated by a distance R (in terms of lattice parameter), the corresponding Hamiltonian can be written as which is independent of the overall position on the lattice, for periodic boundary conditions. An eigenstate of the full Hamiltonian, |ψ , with eigenenergy E, can be written in terms of the eigenfunctions of Applying the full Hamiltonian on |ψ and using orthonormalisation of the set of |φ k,R lead to where E k,R are the eigenenergies of H R . The B-O approximation is valid when these eigenenergies are well separated. In this case, the off-diagonal coupling terms φ k,R |H T |φ k ,R can be neglected. Now, E k,R in the uncoupled equation (12) can be interpreted to be the long-range interaction potential between the two impurities, mediated by the majority atoms tunnelling faster. We can now find the ground state of the full Hamiltonian by diagonalising equation (12) with lowest energy states of H R .
We first numerically investigate the case where two atoms from each species are considered in long chains (in order to reduce finite size effects) for different values of R, J 2 /J 1 , U 1 , and U 12 . The findings are shown in figure 3 where we take a chain with M = 50 sites and fix J 2 = 0.01J and U 1 = 2J. The validity of the B-O approximation is checked and as shown in figure 3(a), the spectrum is found to be properly gapped only for sufficiently attractive U 12 . In the attractive case, as |U 12 | grows larger than J 1 , the ground state is a bound state where the two majority particles are localized at the positions of the static impurities, with an energy gap opens to a state where only one majority particle is localized, and the other one is delocalized around that site. For repulsive U 12 the spectrum is gapless as transition to excited states is caused by any tunnelling of the lighter majority particles, therefore, we only consider attractive interactions in the following.
The potential energy surfaces E R for this system are then computed. The attractive long-range effective interaction U eff 2 is shown in figure 3(b), as a function of the distance R. Here, the U eff 2 values are the ground state energies of H R , and the zero of the energy is chosen at the largest R. Next, we calculate the effective tunnelling J eff 2 of an impurity in the presence of delocalised majority species bosons. This can be done by confining the impurity in a double-well in the middle of the chain with majority species bosons with small enough U 1 . The effective tunnelling is then given by half of the energy difference between the lowest even and odd states. This is shown in figure 3(c) as a function of U 12 . With these parameters computed, the resulting visibility profile for two impurities is shown as the red line in figure 3(d). This reflects the mediated interactions and increased effective mass captured by the B-O approximation. The blue line in figure 3(d) displays the total visibility in the actual ground state of the full four particle system, which shows an additional decrease in visibility from the B-O model.
We now carry out similar calculation for two impurities interacting with relatively larger number of majority bosons, in the context of the systems we consider in this work. In figure 4 we show the results for a 12-site chain with 12 majority species bosons and 2 impurity bosons. For reasons stated before, we only consider attractive U 12 . We also ensure that the coupling between the ground state and the excited states in the Born-Oppenheimer treatment stays small in the parameter ranges we are interested in. This is shown in figure 4(a) where the maximum, taken over all R, of the ratio of tunnelling coupling between ground state and a few excited states and the energy difference between them, is plotted, for J 2 = 0.1J. We have taken U 1 = 2J to ensure that the majority species is delocalised. In order to not disrupt the ground state of the majority species too much we use |U 12 | < U 1 . With these choices the effective parameters J eff 2 and U eff 2 for the two impurities are calculated as before. The resulting visibility profile is shown in figure 4(b), as the blue line. To show the corresponding visibility for this system due to the combined effects, with entanglement, we compute the many-body ground state, using DMRG methods. The visibility of the impurity species is shown as the red line in figure 4 With this understanding we now focus on the interspecies entanglement and its effects on larger systems beyond the two particle system. Below we will see analogous behaviour of that presented in section 3 in the many-body case, made somewhat more complicated by the dynamics of interacting particles in the majority species.

Many-body case for impurities on a lattice
In this section, we now investigate the interspecies entanglement in regimes where many-body dynamics play a key role. We identify the corresponding effects of interactions on the visibility of a p = 0 peak in the momentum distribution, and use this to understand signatures of the many-body phase diagram of the majority species in the dynamics of the impurity atoms.

Effects on the majority species
For each of the cases in this section, we compute the ground states of the Hamiltonian equation (1), which we compute using ED (exact diagonalisation), DMRG, and Gutzwiller mean-field methods as discussed in section 2. As discussed above, we refer to the height of the peak of the quasi-momentum distribution per particle, denoted by V σ for the species σ, as the visibility of the momentum distribution. Similarly, the von Neumann entropy S vN shows the effects of entanglement between the two species. For the ED calculations (with periodic boundary condition) the lattice consists of 6 sites (M) and for the mean field and DMRG calculations we have used M = 32 and M = 16 respectively. In all cases the number-dominant reservoir species is at unit filling (N 1 = M) and the impurity species is at quarter filling (N 2 = M/4), except for M = 6 where we have taken N 2 = M/3.
In general, as the two species become more entangled with increasing interactions, the visibility of the momentum distribution clearly decreases. However, there are a number of important many-body phenomena that are visible at specific points in the visibility profiles as a function of interspecies interaction. The passage from a delocalised phase to a localised one for species 1, as the intra-particle interaction increases, can cause a rise in the visibility, while entanglement keeps increasing. Alternatively, there are points in the parameter space where the ground state goes through abrupt structural changes, causing a sudden increase in the von Neumann entropy, which cannot be captured by observing the momentum distribution alone. In the following we analyse these features in more detail.
In figure 5 we show the entanglement of the species 1 in terms of the von Neumann entropy S vN as a function of interspecies interaction from the ED results in figure 5(a), and the visibility profiles from the DMRG calculations in figures 5(b) and (c) and mean field calculations in figure 5(d). We notice an increase in entanglement and decrease in visibility as the interspecies interaction U 12 is increased from zero, as was seen in the previous section for the system of two bosons on an optical lattice. A change in U 2 does not have significant effect on the general entanglement or the visibility profiles, so we fix the value to be U 2 = 32J.
When the majority species 1 particles are non-interacting (U 1 = 0) they are delocalised at U 12 = 0 and we see a high visibility of the momentum peak at p = 0 in this case. When U 12 is increased, we see the decrease in the visibility, in a form that is largely familiar from the two-particle case in the previous section.
For repulsive interactions between majority atoms, we first look at U 1 = 2J where the particles of species 1 are still largely delocalised at U 12 = 0 as they still are in the superfluid regime of the single-species Bose-Hubbard model. The finite U 1 value, however, results in a slight decrease of the visibility V 1 . As U 12 is increased from zero we see similar behaviour for both the visibility and von Neumann entropy to that seen in the U 1 = 0 case but now taking U 12 > 0 further reduces the delocalisation of species 1 atoms. This leads to a small further decrease of visibility compared with the non-interacting case, but for U 12 < 0 the decrease in visibility is correspondingly less than the U 1 = 0 case.
We see strong features of the many-body physics of the majority species entering the dynamics as we further increase U 1 , so that U 1 /J is larger than the critical value for the superfluid-Mott insulator transition. In 1D this critical value has been reported as, (U/J) c ≈ 3.3 [74]. For U 1 = 8J we see that the visibility has a very different shape, which is characteristic now of the behaviour when the majority species is in a Mott insulator regime for U 12 = 0. We see that the visibility V 1 here has a minimum at zero interspecies interaction. This can be understood as being due to the fact, in the Mott insulator the particles are exponentially localised at each lattice site. This results in a broadening of their momentum distributions and causes the dip in the height of the momentum peak. For low interspecies interaction the particles from the species 2 do not have sufficiently strong interactions to excite species 1 atoms out of the Mott insulator entirely, but they do lead to some delocalisation through virtual amplitudes to create such excitations. This can also be seen in an increase of the von Neumann entropy. If we look at the U 1 = 8J line in figure 5(b) the subsequent local maxima on the both sides of the minimum occurring at U 12 = 0 happen due the increase in U 12 where the effect of the presence of a second species becomes stronger. As U 12 becomes comparable to U 1 , the energy input due to the presence of a species 2 particle disrupts the localized phase as the energy penalty for having a double occupation of species 1 is comparable to the energy required to put two particles from the different species on a single site. Thus the species 1 particles begin to delocalise and the visibility increases substantially. However further increasing U 12 imposes a restriction on this delocalisation process which causes a drop in the visibility again.
The interplay between the different interaction parameters gives rise to some particularly interesting individual features at particular points in parameter space, most notably a surprisingly large increase in the value of the von Neumann entropy that occurs for U 1 = 2J at around U 12 = −9J, as can be seen in figure 5(a) (green line). This happens due to drastic changes in the nature of the ground states, and shows how sensitive this measure can be to such structural changes, in a regime where this cannot be detected via momentum distribution changes. Around U 12 = −10J it is energetically favourable to have all the species 1 and species 2 particles at one single site. In figure 5(a), this can happen in 6 possible ways as we look at a 6 site system with periodic boundary condition. On the other side of the peak-like structure, around U 12 = −8J, it is energetically favourable to have both the species 2 particles on adjacent sites, and this configuration can also achieved in 6 different ways. The von Neumann entropy is therefore indeed log 2 6 on both sides of the peak. Now around the peak, which is near U 12 = −9J all the 12 configurations become important and the von Neumann entropy becomes log 2 12. Carrying out a Schmidt decomposition between the species reveals that the ground state is very close to a maximally entangled states with 6 almost equal singular values for U 12 = −10J and U 12 = −8J, and 12 almost equal singular values for U 12 = −9J. The other singular values are suppressed by at least three orders of magnitude. Looking at the energy levels of the composite system we can also see that the lowest six levels are very close to each other at U 12 = −10J and U 12 = −8J whereas there is an avoided crossing with second lowest six levels at around U 12 = −9J. For a general M-site system with two impurities with large intra-species repulsion, this jump of S vN from log 2 M to log 2 2M would occur at particular value of attractive U 12 , determined by the system parameters. A simple estimation of energies in the two distinct configurations shows that this happens around |U 12 /J| ≈ U 1 N 1 /2 + U 2 /N 1 .
For U 1 = 32J the particles in species 1 are in the deep Mott insulator regime and in the range of U 12 that we are looking at here the energy input in the system by the presence of the particles of species 2 cannot affect the Mott insulator as U 12 is always much smaller than U 1 . Since varying the interaction strength does not entangle species 1 with the other the von Neumann entropy stays at zero. The visibility of these highly site-localised species 1 particles also stays constant at a very low value which is even much smaller for the mean field treatment as we can see in figure 5(d). This is because in mean field treatment the spatial correlations in a Mott insulator are exactly zero and in a numerically exact treatment they fall exponentially with the distance in space.
In figure 6 we look more closely at the visibility profile which changes from going through a maximum at zero interspecies interaction (for example, the U 1 = 1.7J line in figure 6(a)) to a minimum (for example, the U 1 = 2.9J line in figure 6(a)) as a function of the U 1 value. Here we notice the transition like feature which occurs in the regime where particles become more localised, but note that it occurs before the superfluid to Mott insulator phase transition in 1D, as it occurs at around U 1 = 2J when computed using DMRG.

Effects on the impurity species
In figure 7 we further investigate the effects on the impurity particles (species 2). As noted above, the properties of these impurity particles depend strongly on the many-body state of the species 1 particles and therefore should be affected by the choice of U 1 values, allowing them to be used as probes for the physics of the species 1 particles. The visibility of species 2 as a function of U 12 in general has a peak around U 12 = 0 that decreases on each side. This peak-like structure starts broadening as we keep increasing U 1 starting from U 1 = 0. The visibility profiles are quite similar to those seen in the previous section for the system of two bosons on an optical lattice in terms of the mechanisms that create a maximum at U 12 = 0 and a slower decrease for repulsive U 12 . The value of the maximum visibility also follows a similar trend as a function of U 2 and falls sharply for attractive U 2 whereas it falls much more slowly on the repulsive side. This particular behaviour is seen in figure 5(c). For very large and positive U 1 the majority species particles are localised and for small interspecies interaction the entanglement remains very small. This phase however is disrupted at sufficiently strong U 12 , causing steep rise in entanglement and subsequent decline in impurity visibility profiles. These interesting features are reported in the following paragraphs and figures 7(a), (b) and (d). In figure 7(c) we show the visibility V 1 of the species 1 for identical parameters, computed with DMRG calculations.
For U 1 = 32J and unit filling, species 1 particles are in the deep Mott insulator regime. For U 2 = 0 we expect the impurities to behave similarly to free particles on a modified lattice and as a result there is little entanglement with species 1, causing a strong visibility around U 12 = 0. As U 12 is increased, it eventually becomes energetically favourable for the impurity particles to be found together at one single lattice site and push out the species 1 particle, creating a hole in the Mott insulator. This causes an abrupt structural change in the localised phase of species 1, in the sense that the impurity atoms now participate in a delocalised form in the Mott insulating phase, which will reflect the z-antiferromagnet phase of a general two-species Bose-Hubbard model [75,76] (which will generally exhibit phase separation). The impurity particles become localised in position space by the species 1 particles through this process and therefore we see a sharp rise in entanglement that results in a sudden decrease in the visibility of species 2. The value of U 12 at which this happens depends on the number of impurity particles and expectedly we notice this value to be U 1 /N 2 in figure 7 as that would be the interspecies energy required to have one species 1 particle and all the species 2 particles at the same site that would match the excitation energy of the Mott insulator. Now on the attractive side of U 12 around the same magnitude (U 1 /N 2 ) it also becomes energetically favourable to create a hole in the Mott insulator and to have all the impurity particles on that neighbouring site of the hole where the species 1 particle has tunnelled to. Due to this similar localisation the two species become highly entangled and the visibility V 2 again falls drastically. For an infinite lattice system (which is the case when one treats the problem in mean field theory) this localisation process causes the visibility to completely vanish, as shown in figure 7(d). For a finite system ( figure 7(b)) the visibility falls down and takes a constant value that decreases as we increase the system size. Now for U 2 = 2J the visibility at U 12 = 0 will be smaller than in the previous case as the impurity particles repulsively interact among themselves. This decrease in height of the visibility persists as we increase the U 2 value but not arbitrarily as we have discussed before. For small values of U 12 the visibility again remains unchanged and we also see the same localisation effect causing a drop in visibility at sufficiently high U 12 as before. The magnitude of U 12 at which the drop happens increases with increase in U 2 as the repulsion between the impurity particles also needs to be overcome. We see this for U 2 = 2J, 8J, and 32J in figure 7(b).
For attractive U 2 we expect a similar plateau-like visibility profile but with much smaller values to start with (at and around U 12 = 0) and we see that in figure 7(b) for U 2 = −2J, although in this case the plateau-like structure is hardly visible due to such small value of the peak. The value of visibility (for all the profiles) after the drop tends to go to 1/M which is the lower limit for the peak of a single particle momentum distribution on a lattice with M sites (one can think about two particles on a lattice with very strong inter-particle interactions where the single particle momentum distribution is completely flat in the ground state).

Effect of the reservoir size
Up to now we have looked at strong effects of entanglement between the two species. For certain applications we would also like to identify regimes where this entanglement can be made small. This includes proposed experiments in which the impurities can be used to measure correlation functions of the majority species [8,66]. We note that it is important that we do not reduce the coupling strength to values so small that impurities cannot be used as probes on reasonable experimental timescales. Instead, if we keep the interaction strength constant, but increase the size of the reservoir, the overall effect of the impurity on the reservoir becomes small, and this is reflected in the interspecies entanglement and corresponding observables associated with the reservoir. For example, when an impurity in an optical lattice is immersed in a Bose-Einstein condensate much larger in size compared to that of the lattice, then in the limit of weak coupling there is no visible effect of the interspecies entanglement. We note that this this limit is also a necessary (but not sufficient) requirement for making a Markov approximation when the impurity atom is coupled to a reservoir in an open quantum system [77,78]. It would be an interesting future study to determine the degree of entanglement between impurities and the majority species as a function of time, and how this changes with the coupling strength between markovian and non-markovian regimes, e.g., in the setting of reference [66], in which impurity atoms are used to probe a majority species of bosons in an optical lattice.
To estimate the effect of the size of the reservoir with non-interacting species, we can consider the two-particle problem where the impurity species 2 is confined to a small part of the lattice, with M 2 lattice sites, as opposed to the full length, M 1 , and we also allow for different tunnelling rates for the two species.
In figure 8, we examine the solution to this problem, and show the characteristic differences in visibilities as a function of the ratio of the lattice sizes, and for different tunnelling rates. As can be seen in figure 8(a), V 2 increases steadily for smaller values of M 1 for a given M 2 and after M 1 /M 2 = 3 reaches the value obtained in the absence of interactions.
In figure 8(b), we explore the effect of having a heavier impurity particle in the lattice system so that its tunnelling rate J 2 is smaller than J 1 . In particular, we would like to find out how small J 2 should be compared to J 1 , so that we can treat the tunnelling term for the impurity particles in perturbation theory. We can solve the two-particle problem both perturbatively, in the limit of small J 2 /J 1 , and also numerically. Figure 8(b) shows the corresponding visibilities, V 2 , as solid lines found with exact diagonalisation, against M 1 /M 2 for descending values of J 2 /J 1 . For J 2 /J 1 = 10 −3 we observe good agreement with results from perturbation theory.
We see that, consistent with our expectations, increasing the overall size of the majority species lattice while restricting to U 12 > 0 rapidly achieves a regime where the entanglement and the effects on the visibility both become small. This effect is interestingly enhanced for equal interactions, and heavy impurities require a larger lattice for the majority species to reach the same regime.

Summary and outlook
In this article we have investigated the many-body entanglement between two bosonic species in an optical lattice. We quantify how the interspecies entanglement introduces additional effects, beyond mediated (a) The visibility of the momentum distribution peak for a single impurity particle V 2 for U 12 = 8J, normalised to the value for U 12 = 0, and plotted for different M 2 values while increasing M 1 . We see that for U 12 = 8J, V 2 increases up to the non-interacting value when M 1 /M 2 is around 3. (b) The visibility of the impurity particle V 2 for U 12 = 8J 1 shown for different M 1 /M 2 values while also changing the relative tunnelling J 2 /J 1 . We see that for J 2 /J 1 = 10 −3 the calculation matches well with analytical calculations using the tunnelling term for the impurity particle as a perturbation (black circles). We have taken M 2 = 5 for this calculation. In both (a) and (b) we have used periodic boundary conditions on both lattice systems for exact diagonalisation.
interactions or a changed effective mass, also in limits where it is possible to distinguish those quantities. In particular, the change of the momentum distribution of each, as a function of interaction strength, could be used as a direct probe of the majority species, either by observation of the impurity atoms or by observation of their effects on the majority species. We have also identified regimes where the entanglement is small, which would be useful in more complex probe experiments [8,66].
The interspecies entanglement could also be directly measured using interference techniques with multiple copies of the state in a quantum gas microscope, alternately performing the scheme of references [50,[52][53][54] for one or both atomic species. This also opens the path towards future studies of impurity atoms and non-Markovian dynamics when they are immersed in a strongly interacting reservoir gas [66,78].