Using entanglement to discern phases in the disordered one-dimensional Bose-Hubbard model

We perform a matrix product state based density matrix renormalisation group analysis of the phases for the disordered one-dimensional Bose-Hubbard model. For particle densities N/L = 1, 1/2 and 2 we show that it is possible to obtain a full phase diagram using only the entanglement properties, which come"for free"when performing an update. We confirm the presence of Mott insulating, superfluid and Bose glass phases when N/L = 1 and 1/2 (without the Mott insulator) as found in previous studies. For the N/L = 2 system we find a double lobed superfluid phase with possible reentrance.

Introduction. -The study of bosons in one dimension has been of great interest in both theoretical and experimental physics for many years due in part to the existence of a quantum phase transition from a superfluid to insulator at zero temperature [1]. The introduction of disorder causes a further phase transition into a localised Bose glass phase, which is insulating but remains compressible [2]. The experimental study of phase transitions in bosonic systems is possible using Helium in porous media [3,4], Josephson junction arrays [5], thin films [6,7] and, more recently, optical lattices [1,8]. It is now possible to introduce disorder in a controlled manner to optical lattices using speckle potentials [9,10] to study these transitions directly [11,12].
Analytical results even for clean systems are limited. There is an approximate Bethe-ansatz solution [13], where the maximum number of bosons per site is set to two. For disordered systems Giamarchi and Schulz used renormalization group (RG) techniques to determine the weak disorder physics given the Luttinger parameter K [14,15]. There are further real-space RG results for the case of strong disorder [16,17]. Numerical approaches provide some of the most effective means of garnering information. Quantum Monte Carlo has been employed in 1, 2 and 3 dimensions [18][19][20][21][22][23], but these methods become difficult in the limit of zero temperature. An ideal method for analysing one dimensional systems is the density-matrix RG (DMRG) [24]. It has been applied with great success to a number of physical systems from quantum chemistry [25] to quantum information [26], including the disordered Bose-Hubbard model [27,28]. The phase diagrams obtained using these methods for the one dimensional case [19,28] remain inconclusive; whilst qualitatively similar, they are quantitatively very different. This difference could be down to the choice of different observables and the difficulties each has with finite size effects.
In recent years the use of entanglement properties as a means of deciphering phase has become commonplace [29][30][31][32][33][34]. Entanglement is a measurement of a wavefunction's non-locality and as such it is an ideal means of analysing various phases. Modern numerical techniques such as tensor networks and DMRG obtain entanglement information as part of the update algorithms, so large amounts of information about the phase is gathered automatically in the course of the RG iterations [35]. In this paper we perform a DMRG simulation of the disordered Bose-Hubbard model in the form of a variational update of a matrix product state (MPS) [35,36] implemented in the ITensor libraries [37]. The disordered Bose-Hubbard model is made up of bosonic creation, b † i , and annihilation operators, b i , on sites of a linear lattice. The Hamiltonian is [28] p-1 where n i = b † i b i is the local occupation or number operator that gives the number of bosons on site i. The potential disorder is modelled via i.i.d. µ i ∈ [−∆µ/2, ∆µ/2] values. For ease of comparison, we have adopted prior conventions [28] for t and U and throughout the rest of the analysis we keep t = 1.
Observables. -The Mott insulator can be differentiated from the Bose glass phase by the existence of the Mott gap, E g , between the ground and first excited state. While DMRG ordinarily finds the ground state of the system, low lying excited states have to be constructed iteratively by orthogonalising with respect to the lower lying states [35]. For the Bose-Hubbard chain it is numerically more convenient to use the fact that the energy of the excited state is equal to the difference in energy between the chemical potential for particle, µ p = E N +1 − E N , and hole, µ h = E N − E N −1 , excitations [28]. This means that E g can be found by calculating the energies E N +1 , E N and E N −1 of the N + 1, N and N − 1 particle sectors, respectively, as Hence the determination of E g requires to run a DMRG for three different particle numbers and each set of parameters.
The superfluid phase is determined by a non-zero superfluid fraction ρ s . This is defined as the difference between the ground state energies of a chain with periodic boundaries and anti-periodic boundaries where L is the chain length and N the number of bosons [28]. For Mott insulator and Bose glass phases, we have ρ s = 0, so a finite ρ s indicates superfluidity in the phase diagram. From a computational point of view, ρ s is not an easy quantity to determine as it requires the use of periodic boundaries, which are well-known to converge slower and be less accurate than for open systems when using DMRG [35]. Furthermore, as ρ s is the difference between two energies, two such periodic DMRG calculations have to be performed for each set of parameters. The two-point correlation function b † i b j -also referred to as the one-particle density matrix [28] or bosonic Green's function [38] -provides information regarding the localisation of the wavefunction. For the Bose glass and Mott insulating phases the correlation function decays exponentially [38] b where ξ is the correlation length and . . . denotes the expectation value when averaged over all pairs of sites separated by |i − j| and all disorder realisations. Extended phases like the superfluid are not localised so ξ diverges in the thermodynamic limit. In the absence of disorder the superfluid phase will be described by Luttinger liquid theory [39], hence the correlation function will admit a power law decay where K is the Luttinger parameter. K takes the value 2 for a Kosterlitz-Thouless (KT) transition from superfluid to Mott insulator [40,41]. By utilizing an RG approach, Giamarchi and Schulz [14] showed that disorder scales to zero in the weak disorder regime when K > 3/2, giving a superfluid phase. On the other hand, disorder grows for K < 3/2 signifying a Bose glass. This was later extended [15] to the medium disorder case (U ∼ ∆µ). Instead of the infinite-system size result (5) we use the conformal field theory (CFT) expression [38] for an open chain of size L, The expression has to be averaged over all i and j with separation |i−j| and, of course, also averaged over disorder realisations. Calculating correlation functions does not require multiple DMRG runs, but requires the calculation of an expectation value for each combination of i and j, of which there are L(L − 1)/2. Furthermore, the accuracy of locating the KT transition from correlation functions for the Bose-Hubbard model has previously been questioned [38,42].
In each DMRG run, bipartitioning the chain into system and environment blocks is done routinely to compute singular-value decompositions [35]. These singular values, s a , can themselves be used to obtain quantitative measures of phase boundaries [31][32][33] without the need for multiple runs, thus saving substantial numerical costs. The most common such measure is the entanglement entropy or Von Neumann entropy defined as which gives the entanglement between regions A and B [35]. The reduced density matrix, ρ A , for region A is the density matrix obtained by tracing over degrees of freedom from region B. Its eigenvalues are given as squares of the s a 's. Hence S A|B is a measure of the spread of the s a values. If there is one non-zero singular value then the regions are in a product state of the two regions. The other extreme is if all singular values are equal, in which case the subsystems are maximally entangled. In the subsequent analysis we shall average the entanglement entropy over all possible bipartitions along the chain. This averaged entanglement entropy is very good at distinguishing between phases with high and low entanglement, for example the superfluid and Mott insulating phases. p-2 Entanglement and the disordered 1D Bose-Hubbard model Deng et. al. [33] used the entanglement spectral parameter, ζ, to obtain the phase diagram for an extended Bose-Hubbard model. The ζ parameter is defined as the sum of the difference between the first and second, and third and fourth, respectively, eigenvalues s 2 a of ρ A when averaged over all bipartition positions such that , a = 1, 2, 3, 4, the bipartition-averaged a-th eigenvalue. In fig. 1 we show the typical behaviour of the four lowest s a values in the superfluid, fig. 1(a), Mott insulator, fig. 1(b), and Bose glass, fig. 1(c), regimes. We see that the entanglement spectra of the superfluid phase is somewhat noisy, but without too much variation along L A . One of the striking features of the entanglement spectrum for the Mott insulator regime is that effectively corresponds to a product state for all bipartitions, even in the presence of disorder. This means that the average entanglement entropy will be nearly zero with a negligible deviation. Last, entanglement spectra of the Bose-glass show pronounced localised peaks, resulting in a much larger variation of ζ (and S A|B ) than in the superfluid phase. These findings suggest that the spatial variations of S A|B and ζ might also be used to distinguish the phases of the disordered Bose-Hubbard model.
Results. -We collected most data for a modest system size of L = 50 in order to allow disorder-averaging over 100 samples using as our DMRG implementation the ITensor libraries [37]. We increase L up to 200 for highprecision estimates of phase boundaries at certain (U, ∆µ) points. For disordered systems, getting stuck in local minima is particularly problematic, so we use a relatively large bond dimension χ = 200 and perform 20 DMRG sweeps of the chain for each sample. Our truncation error is less than 10 −10 . We also introduce a small noise term for the first few sweeps; this perturbs a perhaps bad initial wavefunction, letting the system converge faster into the desired target states.
Bosons do not obey the Pauli exclusion principle and hence can condense in principle onto a single site. In order to capture such a behaviour, the one-site basis dimension should to be as large as the number of particles in the system. This is numerically infeasible and it is necessary to introduce a finite maximum number of bosons that can occupy each site. We set the maximum number to five 1 . This is consistent with Ref. [41] who find that a higher particle number does not effect the results appreciably for U > 0 [27,28,38].

Density = 1.
For particle density N/L = 1, in the clean case, the system is in a superfluid phase for small U but transitions into a Mott insulating phase at a critical U c . Introducing disorder enables the existence of a localized Bose glass phase [2]. The possibility of a direct transition from superfluid to Mott insulator has been discussed extensively (see references in [43]). In one dimension it was shown [44] that the transition necessarily goes via the Bose glass phase. This is now also the accepted picture for any dimension [43].
We show our results based on ζ and S A|B for L = 50 in fig. 2, (a) and (b), respectively. The superfluid, small U 1.5, and the Mott insulator, U 2, are clearly distinguishable in both panels. The boundary of the superfluid to the Bose glass is less well defined and it is not even clear that there is a Bose glass region between the Mott insulator and the superfluid. Thus, the in principle very different wavefunctions give similar average entanglement entropy. Following in from our discussion of fig. 1 in the last section, we also plot in fig. 2(c) the standard error 2 of ζ, ∆ζ, and, similarly, (d) ∆S A|B . In these plots the phases become clear and their boundaries are consistent with earlier work [28]. In particular, a Bose glass phase can be easily identified between Mott insulator and superfluid. Furthermore, we see that the contours for ζ and S A|B in fig. 1 are qualitatively similar, just as those for ∆ζ and ∆S A|B . We emphasize that for the entanglement-based measures present here, it is in fact possible to discern all of the phases with just a single DMRG run for each (U, ∆µ, disorder realisation) data point. This is a clear advantage in terms of numerical costs when compared to calculations based on E g , ρ s or K.
In order to augment the finite-size phases identified in fig. 2, we now perform runs with larger L and perform finite-size scaling. In order to find the superfluid-Bose glass transition in the thermodynamic limit we calculate K for various points along the boundary for system sizes L = 30, 50, 100 and 150. The transition is of KT type at K = 3/2. The corresponding points in (U, ∆µ) which are shown as filled circles in fig. 2. For reference, we also plot the points where K = 3/2 for L = 50 (stars). Similarly, the superfluid-Mott insulator transition point U c is the point on the zero disorder axis where K = 2. We estimate it value as U c = 1.634 ± 0.002. We also calculate E g for the same system sizes and perform finite size scaling to find the Mott insulator-Bose glass boundary. These are indicated as stars in fig. 2.
The RG analysis of Refs. [14] and [15] suggests that there may be a further Anderson glass phase in the low U < ∆µ region of the diagram. This would imply a critical point along the superfluid boundary at which point K at the transition becomes disorder dependent in a similar manner to the strongly disordered quantum rotor model [16,17]. Our entanglement analysis along with E g and K shows no sign of a transition. It does however become very difficult to observe a clear KT transition in the finite size scaling of K when U ≪ ∆µ and the truncation of the basis becomes more problematic in this region. Density = 1/2. The clean case for N/L = 1/2 remains a superfluid for all values of U [2]. When ∆µ is increased, our entanglement measures indicate the eventual emergence of a Bose glass phase as shown in fig. 3(a+b). Still, the superfluid phase for L = 50 seems to extend up to ∆µ 1 for U 5 as shown by all four entanglement measures. The Giamarchi-Schulz criterion [14,28] implies that the Bose-Hubbard model should be in a Bose glass phase for K < 3/2. In fig. 3(a+b) we show that the resulting boundaries indicate that the superfluid phase extends as far as U K=3/2 = 3.5 ± 0.1, i.e. it ends somewhat earlier for low ∆µ than suggested by our entanglement measures. In order to explore this region further, we have also calculated ρ s for fixed ∆µ = 0.5 and sizes L = 50, 100, 150, and 200 as shown in fig. 4(a). The results for ρ s have been computed for increased bond dimension χ = 400 with 40 DMRG sweeps and 20 disorder configuration to offset the reduction in precision due to periodic boundaries. The figure shows that for U 3, ρ s decreases when increasing L as expected in the Bose glass phase. However, the decrease is very slow and, for the system sizes attainable by us, even seems to saturate at non-zero values. These results suggest that for finite systems, the K = 3/2 criterion significantly underestimates the extent of the superfluid phase, while our four entanglement measures and ρ s pre- dict a much larger region. Performing a finite-size scaling analysis for K as shown in fig. 4(a) we find the U values, for which K = 3/2 in the limit L → ∞, converge towards a limiting value of U c = 3.09 ± 0.01 (see also fig. 3). This again indicates that in an infinite system, we expect the superfluid to Bose glass transition to take place at much lower values of U c than observed for L = 50. The relevance of this result is of course that experimental realizations of the Bose-Hubbard model are typically in cold atom systems, which are limited to finite system sizes, currently a typical lattice dimension is ∼ 50 − 100 [1]. For values of ∆µ 1, the situation is less severe and we see in fig. 3(a+b) that our entanglement-based measures again qualitatively agree with the Giamarchi-Schulz criterion, both for L = 50 and estimated via finite-size scaling at L → ∞. Density = 2. To the best of our knowledge, the phases for N/L = 2 have not been shown before in the literature. Due to our numerical restriction of five bosons per site, this regime is close to the limit of what can be studied reliably, particularly for small U where the occupancy per site should be large. For large U , one might expect that we will have a Mott insulator of boson pairs, while a superfluid of boson pairs emerges for small U and small ∆µ. Similarly, we envisage a disordered Bose glass phase for large ∆µ. With more particles per site than in the N/L = 1 case, we could furthermore expect that onset of the Mott transition at ∆µ = 0 is at larger values of U , since there is a larger energy penalty to pay for a doubly occupied site [2]. Similarly, as the cost for two boson pairs to go onto the same site is 2U , we expect the 2U = ∆µ/2 line to characterize the superfluid phase as in the N/L = 1 case. In addition, one might conjecture to see a remnant of the U = ∆µ/2 condition.
In fig. 3(c+d), we show that our expectations are largely validated. In particular, a double lobe shape for the superfluid phase emerges and allows a possible re-entrant  behaviour given a suitable cut across parameter space. The gradient of the Mott insulating phase boundary is shallower (∼ 4/3) when compared to N/L = 1. Furthermore, both the ζ and S A|B based entanglement measures, as well as their errors, ∆ζ and ∆S A|B , capture the phases equally well and agree with the K and E g estimates. Note that for N/L = 2, the KT superfluid-to-Mott transition at ∆µ = 0 corresponds to K = 2 and we finite-size scale the Luttinger parameter to find U c = 2.75 ± 0.03 as shown in fig. 5.
We emphasize that the points for small U , see top left of the phase diagram in fig. 3(c+d), should be viewed with caution as the basis truncation will affect the results. Conclusion. -We have analysed the phase diagrams of the disordered Bose-Hubbard model for fillings N/L = 1/2, 1 and 2 using the entanglement-based measures ζ, S A|B , ∆ζ and ∆S A|B . We found that these measures are an excellent means of quickly identifying the different phases of the system while reducing the number of multiple DMRG runs per measurement or removing the need for special boundary conditions. However, ζ and S A|B alone did not always faithfully reproduce the phase diagrams and were sometimes misleading regarding the positions of the phase boundaries. The error-based measures, ∆ζ and ∆S A|B , provide a much clearer picture -the distributions of the values contain more information regarding the nature of the phase than the mean values alone. For N/L = 1 the phase diagram is as calculated earlier using K, E g and ρ s [19,28], providing strong support for the numerically simpler approach based on using ζ, S A|B , ∆ζ and ∆S A|B . For N/L = 1/2 the diagram shows strong finite size effects and the critical U defined by the Giamarchi-Schulz criterion is not apparent for these finite systems. For N/L = 2 the superfluid phase has a double-lobed appearance giving rise to re-entrance phenomena.