Entanglement structure of the two-component Bose-Hubbard model as a quantum simulator of a Heisenberg chain

We consider a quantum simulator of the Heisenberg chain with ferromagnetic interactions based on the two-component 1D Bose-Hubbard model at filling equal to two in the strong coupling regime. The entanglement properties of the ground state of the two-component Bose-Hubbard model are compared to those of the effective spin model as the interspecies interaction approaches the intraspecies one. A numerical study of the entanglement properties of the two-component Bose-Hubbard model is supplemented with analytical expressions derived from the effective spin Hamiltonian. When the pure ferromagnetic Heisenberg chain is considered, the entanglement properties of the effective Hamiltonian are not properly predicted by the quantum simulator.

The Bose-Hubbard model is almost ubiquitous nowadays in the interpretation of ultracold atomic gases experiments with optical lattices 1 . It provides the prime ingredient that allows ultracold atomic setups to mimic well-known many-body problems 1,2 . In particular, it makes these systems extremely competitive for building quantum simulators of a wide range of notably difficult physical problems 3,4 . A particularly relevant example is the use of a two-component Bose-Hubbard (TCBH) model as a quantum simulator of spin models [4][5][6][7] . As pointed out in refs 5,6 . different spins, e.g. 1/2, 1, etc, can be simulated depending on the filling factor of the two species in the chain. In the present article we concentrate on the specific configuration of filling one for both species, i.e. equal number of atoms of both species in the chain, in which case the TCBH can be mapped into a ferromagnetic Heisenberg spin −1 model 5 .
In general terms, a quantum simulator is defined as an experimentally feasible and versatile setup which is able to mimic a target Hamiltonian in different parameter regimes. In this way, for instance, the low energy physics of both the target Hamiltonian and the quantum simulator should be very similar. One could, however, be interested in other properties of the target system such as simulating the quantum correlations present in the ground state. The entanglement spectrum, characterizing the pairwise entanglement present in the system, provides a powerful witness of the presence of quantum correlations 8 . This seems a relevant goal for near-future applications of quantum entanglement to a variety of quantum technologies, see for instance 9 .
In this article we consider the question: To what extent does the quantum simulator exhibit similar entanglement properties to the simulated Hamiltonian? In particular, we focus on critical regimes where specific entanglement properties universally characterize the phase of the system. The analysis will be performed in the strongly interacting regime, where the interaction strength of both species is equal and much larger than the tunneling rate. We will study the entanglement properties of the system as the interspecies interaction is increased towards the point where all interactions are equal. In this way, the effective spin model goes from an anisotropic Heisenberg model into the isotropic Heisenberg one. Analytical results using perturbation theory will be complemented with numerical calculations using DMRG (density matrix renormalization group). In this way we can compare the entanglement present in the TCBH with that of the spin model, paying particular attention to the critical phases which appear in the latter.

Model
We consider two bosonic species with contact-like interactions in a 1D optical lattice at zero temperature. We assume the system to be described by the TCBH Hamiltonian, are the annihilation (creation) bosonic operators at site i = 1, …, L for species α = A, B, and α n i are their corresponding number operators. We have assumed equal tunneling strength, t > 0, and equal repulsive intra-interaction strength, U > 0, for both components. For the rest of the work we set the energy scale to t = 1. The ground state (GS) of Eq. (1) in the strong-coupling regime (U ≫ U AB , t) is a Mott insulator (MI) with a total filling ν = N A /L + N B /L ≡ ν A + ν B . In this work we fix ν A = ν B = 1.
We define the entanglement properties through the reduced density matrix obtained tracing out the right half of the system ρ L/2 = Tr R |ψ〉〈ψ|, where |ψ〉 is the ground state of the TCBH. The amount of entanglement is quantified with the von Neumann entropy S E = −Trρ L/2 logρ L/2 . Finally, the entanglement spectrum (ES) 10 is defined in terms of ξ i = −logλ i , where λ i are the eigenvalues of the reduced density matrix.

Results
Perturbative regime. In the strong-coupling regime (U ≫ t), the ES can be obtained perturbatively following 11 . In order to organize the ES we introduce the quantum numbers δN α = N α,L/2 − L/2, with α = A, B. They measure the excess (δN α > 0) or absence (δN α < 0) of bosons with respect to the MI with filling ν A = ν B = 1 on the left subsystem (of size L/2). In Fig. 1 we report the obtained entanglement spectrum as a function of the interspecies interaction, U AB , for a fixed, large, value of U = 50. For U AB = 0 the ES of the single-component Bose-Hubbard model is recovered, with clearly separated levels corresponding to different orders in perturbation theory on 1/U 11 . For non-zero values, U AB > 0, some entanglement values exhibit an explicit dependence on this interaction. The entanglement values associated with δN A = ±1; δN B = 0 and δN A = 0; δN B = ±1 are given by, and do not show an explicit dependence on U AB at the order studied. Furthermore, these ones are completely analogous to the first ones obtained for the single-component Bose-Hubbard. The lowest entanglement value associated with δN A = δN B = 0 gets a contribution ξ = U 8/ 0 (2) 2 due to the renormalization of the wavefunction. Genuine second-order contributions are of two different kind: (+), corresponding to δN A = δN B , which favor the movement of two different bosons through the boundary in the same direction and, (−), δN A = −δN B , which favor the hopping of two different bosons through the boundary in opposite directions. Unlike ξ 1 (2) , these ones are absent in the single-component Bose-Hubbard model as they are directly related to the presence of two different components. Specifically, configurations with δN A = −δN B are associated to the phase separation of the two components through the boundary. An analytic formula can also be obtained, The two different branches, (+) and (−), have a very different behavior as U AB is increased. The (−) one is seen to decrease as U AB /U → 1. This is expected as, in absence of tunneling, the system becomes highly degenerate at www.nature.com/scientificreports www.nature.com/scientificreports/ U AB = U. We predict a closing of the entanglement gap -difference between the two lowest entanglement values 10 -at U AB = U − 2/U, see Fig. 1 (left panel). The fact that the entanglement gap closes is reflected in an increase of the von Neumann entropy S E , see right panel in Fig. 1. The analytical predictions are in very nice agreement with the DMRG calculations. Note also that the structure of the ES changes dramatically as we approach this point, with higher order processes becoming comparable to the lowest entanglement value. These higher order processes also show a logarithmic dependence as found for ξ − (2) with a slope that indicates the order of the perturbation theory at which they are found. At this point we expect the system to enter in a critical regime. In this regime the von Neumann entropy S E increases as we commented previously, but it also starts to depend on the system size, see Fig. 1.
the tCBH as a quantum simulator. Interesting physics appears as (U − U AB ) ~ 1/U. As discussed above, the system enters in a critical regime which cannot be described by the simple perturbation theory. Instead one has to consider a degenerate perturbation theory. For the general case of n atoms per site, the low-energy Hilbert space is described by an effective spin S ≡ n/2 where A and B are taken as a pseudo-spin 1/2. The effective Hamiltonian describing this low-energy spin subspace is given by superexchange processes at second-order in the hopping 5,7,12 In the limit of zero hopping, t = 0, and isotropic interactions, U = U AB , the spectrum of the Hamiltonian (1) has different degenerated subspaces separated by an energy scale of order U. The introduction of a small hopping t/U ≪ 1 and a small anisotropy |U − U AB | ≪ U breaks the degeneracy of these subspaces. Considering degenerated second-order perturbation theory in the hopping, one can write an effective Hamiltonian (4) which describes the degeneracy breaking inside the lowest energy subspace. In this way, the Hamiltonian (1) is mapped into the Hamiltonian (4). This is what allows one to term the TCBH a quantum simulator of the Heisenberg model. But what happens with observables? We deal with this question by using degenerated perturbation theory, see e.g. 19 .
We split the TCBH Hamiltonian in two pieces, H = H int + H t , corresponding to the interaction part, H int , and the tunnelling piece, H t , respectively. The lowest energy degenerated subspace of H int is the one defining the subspace P, with energy E p 0 . The remaining states, which are separated by a gap proportional to U, define the Q subspace. In our case, the effective spin Hamiltonian (4) acts on the P space, while the TCBH acts on the complete Hilbert space, P ⊕ Q. Given the complete wave function |ψ〉, one can obtain the projection on the P subspace, ψ ψ | 〉 = | 〉 P P , where the projector reads, α α = ∑ | 〉〈 | α∈ P P , and |α〉 are eigenstates of H int . The inverse problem is formally written as ψ ψ | 〉 = Ω| 〉 P , where Ω is called wave operator. One can write expressions for Ω at a specific order in perturbation theory. Since the mapping between H and H eff has been obtained at second-order in the hopping, it is sufficient to express Ω at first-order in the hopping 19 , where β E 0 are the eigenenergies of H int corresponding to the eigenstate |β〉. In terms of the wave operator one can write the explicit expression = ΩˆĤ PH P eff which ensures that if H|φ〉 = E|φ〉 then H eff |φ P 〉 = E|φ P 〉, with φ φ | 〉 = | 〉 P P . Finally, the reduced density matrix ρˆL /2 associated to the groundstate of the complete Hamiltonian (1) |ψ GS 〉 can be expressed in terms of the groundstate of the effective one (4) ψ | 〉 P GS as, Thus, the question of how well the entanglement properties are reproduced in a quantum simulator is rewritten as: can Ω introduce some extra structure which affects the universal entanglement properties? entanglement in the critical regime. The scaling of the von Neumann entropy can be used to characterize the different phases of the system. From a CFT description this is a well-known result 20,21 and the magnitude ΔS = S E (L) − S E (L/2) captures the scaling behavior properly 22 . Following the known behavior of the effective model, Eq. (4), one expects to go from ΔS → 0 in the large −D phase, to ΔS = (c/6)log2 in the critical XY phase with c = 1. This is exactly what is seen in Fig. 2, where we observe the crossover between the two regimes in the TCBH as we vary U(U − U AB ) ~ D/J, with a very nice agreement with the results obtained for the effective spin model 23 . Furthermore, these results are mostly independent of U, for sufficiently large U. Therefore, we can conclude that the transition in the spin picture from a large −D to a critical XY ferromagnetic phase is captured by the transition observed in the TCBH. On the other hand, as U − U AB → 0 a dependence on U starts to appear. , the ground state of the system is ψ , which is a superposition of all spin configurations in the chain satisfying that the total magnetization is S T z . Therefore, considering a bipartition of the system A of length l the ES is organized by eigenstates with well defined magnetization = − S Sl m A z in the subsystem A with eigenvalues m with m = 0, …, 2Sl, which is a natural extension of the results presented in 24,26,27 . From Eq. (7) an asymptotic expression for the von Neumann entropy S E can be obtained considering l = L/2 and = S 0 Notice that in the thermodynamic limit any small anisotropy D > 0 will restore the conformal symmetry. For finite systems a smooth crossover between the CFT and the scale invariant prediction (8) is expected 28 , see Fig. 2. In this region is where a non-universal behavior of the TCBH model is observed and we obtain different scalings of S E for different values of the interaction U. Furthermore, we observe that in the limit (U − U AB ) → 0 the scaling mostly depends on the value of the interaction U and does not coincide with the value predicted by the spin model, Eq. (8).
In order to understand the dependence of the scaling of the entanglement entropy on the interaction U at U = U AB , we examine the ES of the TCBH model (1) and compare it with the analytical prediction for the spin model (7). The ES represented in Fig. 3 displays a parabolic dependence as a function of δN A − δN B (which is analogous to δ = − S S S z T z A z in the effective spin model). This parabolic dependence is also expected from the spin picture, see Eq. (7), but the curvature is considerably different. This curvature is directly related with the entanglement gap and one can observe that in both situations it depends linearly on the inverse of the system size L, see Fig. 3. But this linear dependence is different in the two models. Specifically, from the spin picture we obtain that δ → 4/L, so it closes in the thermodynamic limit L → ∞. Conversely, in the TCBH model the gap does not close in the thermodynamic limit for finite values of the interaction. Furthermore, the ES predicted by the spin model (7) has a well defined magnetization δS z , meaning that for each value of the magnetization there is a unique entanglement value ξ δS z . On the other hand, the ES of the TCBH model shows a richer structure with different parabolic envelopes for the same magnetization. Focusing on this extra structure we observe that the second parabolic envelope has associated a half-integer magnetization δS z , unlike the first one which has integer magnetization. This can be understood expanding the wave operator at first-order, see Eq. (5), where H t is the hopping term of the Hamiltonian (1). The second envelope is obtained by the application of ψ | 〉 H t P GS over the frontier which defines the bipartition of the system used to compute the ES. Therefore, these entanglement eigenstates correspond to having an extra particle or hole δN = ±1 for any of the two species which explains the half-integer nature of δS z . Notice that this component of the ground state wavefunction is reminiscent of the first entanglement eigenstates with eigenvalue ξ 1 (2) , Eq. (2). But now, due to the non-trivial entangled structure of the ground state, for each value of the subsystem magnetization we have this particle-hole excitation over the frontier which gives a large number of states, of order L. We have checked that the gap between the first two parabolic envelopes closes like 2logU and does not show an explicit dependence on the system length L.
The effect of including ψ | 〉 H t P GS in the wavefunction has large effects on the von Neumann entropy. The main reason is that the number of entanglement states given by ψ | 〉 H t P GS is of order L which is the same than the number of entanglement states in ψ | 〉 P GS . Therefore, the contribution of both parts to the von Neumann entropy is logL and  www.nature.com/scientificreports www.nature.com/scientificreports/ we can estimate the total contribution as S E ∝ (1/2 − A/U 2 )logL, with A some constant value. In order to verify that, we define the slope c eff (L) = 6(S E (L) − S E (L 0 ))/(log(L/L 0 )) with a reference size L 0 = 50, for which finite size effects will be reduced 29 . In Fig. 4 we see that there is always a logarithmic behavior and the slope c eff (L) shows a clear dependence on 1/U 2 which confirms our predictions.

Conclusions
The extent to which a quantum simulator of a well-known spin system captures the entanglement properties of the ground state of the effective Hamiltonian has been scrutinized. We have considered the entanglement properties of the ground state of the two-component 1D Bose-Hubbard model in the strong-coupling regime for total filling ν A = ν B = 1. This model acts as a quantum simulator of the spin 1 Heisenberg model with ferromagnetic interactions. In the regime in which the spin system is in a critical XY phase (U − U AB ~ t 2 /U) the two-component Bose-Hubbard model shows a universal (independent of the interaction U) scaling of the von Neumann entropy, which matches the CFT prediction expected for the effective spin system. On the other hand, we observe that this universality is lost as we approach the isotropic point U = U AB where the effective spin model loses the conformal invariance. By comparing the ES of the quantum simulator with the effective spin model, which has been analytically obtained, we observe large discrepancies between the two of them for large values of the interaction U. In particular, magnitudes which should display a universal behavior (like the slope of the scaling in the entanglement entropy) strongly depend on the interaction U. This dependence has been analytically predicted constructing the wavefunction of the two-component Bose-Hubbard model using the wave operator.