Superfluid density and quasi-long-range order in the one-dimensional disordered Bose-Hubbard model

We study the equilibrium properties of the one-dimensional disordered Bose-Hubbard model by means of a gauge-adaptive tree tensor network variational method suitable for systems with periodic boundary conditions. We compute the superfluid stiffness and superfluid correlations close to the superfluid to glass transition line, obtaining accurate locations of the critical points. By studying the statistics of the exponent of the power-law decay of the correlation, we determine the boundary between the superfluid region and the Bose glass phase in the regime of strong disorder and in the weakly interacting region, not explored numerically before. In the former case our simulations are in agreement with previous Monte Carlo calculations.


Introduction
The Bose-Hubbard model is a widespread paradigm in the field of strongly-correlated many-body systems [1]. Over the years it has found applications in the physics of superconducting thin-films [2], Josephson junction arrays [3], optical cavity arrays [4], and cold atoms in optical lattices [5]. The model describes bosons hopping on a tight-binding lattice and interacting through an on-site repulsion. The non-trivial properties of the model stem from the competition between the hopping (that delocalizes the particles) and the local repulsion (that tends to suppress local particle number fluctuations). In the absence of disorder such a competition dictates a quantum phase transition at integer filling: if the hopping dominates, the ground state of the system is superfluid; in the opposite regime, there is a gap in the excitation spectrum and the system is in a Mott insulating state. At a critical value of the ratio between the hopping and the repulsion parameters, an insulator to superfluid quantum phase transition takes place at zero temperature. In optical lattices this transition has been experimentally demonstrated in the seminal work by Greiner et al. [6]. The phase diagram of the Bose-Hubbard model has been studied by a number of analytical and numerical methods, a recent account of the properties of this model together with a review of the experimental advances with cold atoms can be found in the reviews [7,8].
In one dimension (1D), which is the case of interest for the present work, the quantum phase transition between the Mott insulating (MI) phase and the superfluid (SF) phase at fixed commensurate lattice filling is known to be of the Berezinkii-Kosterlitz-Thouless (BKT) type [1]. The critical properties of the quantum model map onto that of a classical two-dimensional XY model. The hallmark of a BKT transition [9][10][11] is a universal jump [12] in the superfluid density (or stiffness or helicity modulus [13]). Extensive studies of the stiffness in the classical 2D XY model are reported for example in [13,14]. Monte Carlo simulations of the one-dimensional quantum model (the mapping onto the XY model becomes accurate at large fillings) were performed originally in [15], more recent results are comprehensively discussed in [16]. On the other hand, the Density Matrix Renormalization Group alias Matrix Product States (MPS) [17] has been applied to the model under consideration [18][19][20][21], although the superfluid stiffness is difficult to extract with high precision because these methods are not ideally suited to deal with periodic boundary conditions (PBC) [22].
The presence of disorder enriches the phase diagram with a new phase, the Bose glass (BG) [1], surrounding a superfluid lobe in the interaction vs. disorder plane of the phase diagram. In absence of interactions, indeed, the disorder induces Anderson localization of single particle states and therefore an insulating phase. This is then lifted by a moderate amount of interactions that tend to reestablish coherence and consequently superfluidity. At large interactions and large disorder, finally, strong correlations should lead back to an insulating phase, distinct from the above mentioned Mott insulator in the clean setup. Although this picture is qualitatively well-understood, the full quantitative characterization of the disordered Bose-Hubbard model remains challenging: considerable efforts have been devoted to study its phase diagram, in particular the location and the nature of the glassy phase [23][24][25][26][27][28][29][30][31][32][33][34]. Developments on the issue, with a special emphasis on numerical simulations, have been reviewed in [35]. Recently, an experiment with optical lattices has made an important step forward in mapping the phase diagram from weak to strong interactions [36].
In the present work we apply a Tree Tensor Network (TTN) method [37][38][39][40][41] (specifically our recently introduced gauge-adaptive TTN technique [42]), a natural ansatz for the simulation of one-dimensional quantum many-body systems with PBC, to the disordered Bose-Hubbard model. We compute the superfluid stiffness and correlation functions, avoiding the detrimental boundary effects arising in simulations with open boundary conditions. Furthermore, by analyzing the statistics of the decay of superfluid correlators we are able to locate the superfluid to Bose glass transition line, to extend it into the previously numerically unexplored regime of small on-site interactions, and to confirm the Monte Carlo results of Ref. [26] for strong interactions.
The paper is organized as follows: In the next section we introduce the model and give a brief overview of the used numerical technique (a detailed discussion can be found in [42]) and we anticipate the main result of the paper, the phase diagram of the disordered BH model at unit filling. In section 3 we analyze the superfluid stiffness. As a benchmark we first consider the case without disorder where we are able to perform a very accurate finite-size scaling to locate the BKT transition. We compare the obtained transition point with the results from other well-known methods for determining the MI-SF transition, such as numerically extracting the Luttinger parameter from the hopping correlation functions [18], finding very good agreement. Moreover, we present our results for the average stiffness in the disordered case. In order to get a precise location of the transition points we study the hopping correlation decay, as presented in section 4. A detailed analysis of the statistics of the exponent of the power-law decay allows us to map out the complete SF-BG transition line in the phase diagram. Finally, we summarize the main results of the work in section 5.

The disordered Bose-Hubbard model
The one-dimensional Bose-Hubbard model is a system of spinless bosons moving in a single-band tight-binding lattice, described by the Hamiltonian The first term accounts for the hopping between neighboring sites (labeled by the index j), the second for the on-site repulsion and the last for the effect of a random inhomogeneous chemical potential. The couplings t, U and V j are the hopping amplitude, the on-site repulsion, and the random potential, respectively; b j (b † j ) are bosonic annihilation (creation) operators obeying [b j , b † j′ ] = δ jj ′ , with n j = b † j b j being the corresponding number operators. The random on-site disorder V j is independent on different sites and uniformly distributed in the interval [−∆, ∆]. In order to compute the stiffness we introduce a twist φ in the boundary conditions. This amounts to considering periodic boundary conditions and performing a Peierls shift in every hopping matrix element where N s is the number of sites of the ring, as schematically depicted in Fig. 1 a) where the twist is explicitly interpreted as a magnetic flux φ piercing the ring and acting on electrically charged particles. Our TTN ansatz is sketched in Fig. 1 b), showing the natural capability to represent states on 1D lattices with PBC, despite its loop-free geometry. The absence of loops in the network is crucial for the applicability of efficient energy minimization techniques needed for the ground state search. For algorithm details we refer the reader to Ref. [42], however it is important to emphasize that with our TTN approach we will be able to compute the superfluid density for large system sizes (up to N s = 256 in the absence of disorder) and therefore perform a reliable finite-size scaling analysis. As we will show, at equilibrium in one dimension this approach provides matching results to Monte Carlo calculations.
The TTN algorithm we use has abelian global symmetries built-in [43], allowing us to work in a canonical ensemble (at zero temperature) by targeting exclusively N-particles states, in accordance with the U(1) particle conservation symmetry [H, j n j ] = 0 of the Hamiltonian H. In the rest of the paper we will in particular consider the case of unit filling, i.e. N = N s . We simulate system sizes of up to  Figure 2. Phase diagram of the BH model with disorder in 1D at unit filling: Bose glass (white region), superfluid lobe (red region), Mott insulator (gray region). The boundaries of the SF lobe are determined via three different criteria: superfluid jump (green squares), Luttinger parameter K = 2/3 (red circles) and algebraic versus exponential decay of correlation functions (blue circles). The blue dashed line is a fit through the blue circles according to µ U ν , with fit parameters µ = 2.2 ± 0.1 and ν = 0.4 ± 0.1 The red dashed line follows the conjectured behavior ∆ c ∼ U 3/4 of the phase boundary for U ≪ 1 [35]. The gray shaded MI region is data from Ref. [26]. N s = 256 sites, with a truncated local boson basis up to a dimension of d = 20 in the low U regime where high local populations are possible. In the disordered scenario we average all observables over a number of n r = 100 different disorder realizations.
In Fig. 2 we anticipate the main result of this work, the phase diagram of the disordered BH model. The technical details and the discussion of the results are presented in the following sections, here we highlight that for strong interactions (the right edge of the lobe) our results obtained with three different criteria are in agreement with state-of-the-art Monte Carlo results. More interestingly, we are able to extend our numerical analysis into the low-interaction regime (left edge of the lobe): two independent criteria strongly agree down to values of the interactions of the order of U ∼ 0.4, where the numerical results nicely match the theoretical conjecture for very small U [35].

Superfluid stiffness
The superfluid stiffness ρ s is defined through the sensitivity of the ground state energy E 0 (here we consider the zero-temperature case) to a twist in the boundary conditions and we set the energy scale by choosing t = 1. Therefore, ρ s = 1 for vanishing on-site repulsion U and vanishing disorder ∆, while ρ s = 0 in the Mott-insulator and Bose glass phases.  As a benchmark of our approach, we first consider the clean case ∆ = 0. At the BKT critical point the superfluid density jumps from a finite value ρ c s to zero. The mapping onto the classical XY model yields the critical value ρ c s = (2/π 2 ) U c in resemblance with the relation P c s = (2/π) T c , known to hold for the classical XY model in 2D [12]. We determine the transition point U c by extracting the superfluid density as a function of the on-site repulsion U for different ring sizes N s . We then perform a Weber-Minnhagen fit according to the theoretically predicted functional dependence of the superfluid density at the BKT transition [13]: with b the only fit parameter. When fitting ρ s (N s ) = 2U/π 2 (1 + [2 log N s + b] −1 ) for various U, the root-mean-square (RMS) error of the fit displays a pronounced minimum at the transition point [13]. Figs. 3 a)-c) show the result of this procedure, resulting in U c = 3.394 ± 0.03 (or, equivalently, t c /U = 0.295 ± 0.003). We now introduce disorder and repeat the study for the averaged superfluid density ρ s . Since the SF-BG transition is also expected to be of the BKT type [23], at least on the right side of the SF lobe (i.e. at large U), ρ s should display a jump at the phase transition in the thermodynamic limit. An example of our numerical results is shown in Fig. 3 d). In qualitative terms, for both sides of the lobe a drop of the superfluid density for increasing system sizes can be observed when entering the BG phase. Like in the clean system, we localize the critical point by means of a Weber-Minnhagen fitting procedure, this time with ρ c s as an additional fit parameter because the height of the jump for N s → ∞ is not known for arbitrary disorder. While on the right boundary of the SF lobe we are able to successfully apply this technique employing system sizes and disorder realization numbers accessible to our numerics, this is not the case on the left boundary. Here we do not achieve sufficient numerical precision to estimate a reliable minimum in the RMS error when fitting Eq. (3), which ultimately makes this method of determining the phase transition unfeasible for the small U regime. This might be as well an indication of a different nature of the SF-BG phase transition in this regime, as proposed in some studies [35].
In order to continue with an accurate analysis of the critical SF-BG line in all regimes we study the superfluid correlations in the next section.

Correlation functions
We proceed further by studying the correlation functions, a method that turns out to yield higher accuracy in determining the transition line. Moreover, we double-check our results via another method for the determination of the MI-SF transition, namely the value of the Luttinger parameter K obtained from the hopping correlation functions [18] In the clean case the criterion for locating the BKT transition is given by K c = 1/2. Indeed, as critical exponents can be obtained with high numerical precision in TTN simulations [42], we expect this method to work well in the present scenario. Moreover, the PBC setting allows for an easy treatment of finite size effects in Eq. (4), simply by replacing the distance r with the effective distancer = crd(r), where crd(r) = N s /π sin(πr/N s ) is the chord function giving the effective distance of two sites j and j + r arranged on a regular polygon of N s edges. In Appendix A we demonstrate the validity of this procedure. Fig. 4 summarizes our numerical results, locating the transition at t c = 0.299 ± 0.002. Here, we set U = 1 for the sake of an easy comparison with Ref. [18]; see Tab. A1 for how this compares to the transition points obtained in the literature. We now once again turn back to the study of the disordered case. As already mentioned, this long-standing problem is unsolved in the low-interaction regime U < 2,  and even for the more easily tractable regime of higher interactions, results from different numerical methods are not in quantitative agreement [24,26,27,35].
We start by studying the hopping correlation functions of Eq. (4), whose critical exponent gives access to the Luttinger parameter K. The Giamarchi-Schulz criterion K c = 2/3 [23,47], where the overline indicates an average over a number of n r disorder realizations, results in the phase boundary shown in Fig. 5. We remark that while the K c = 2/3 criterion is known to be appropriate in the high U regime [35], there is no definite evidence that this is also the case on the left boundary of the SF lobe (low U regime). Therefore, the results of Fig. 5 require critical verification, which is presented in the following paragraphs.
As a first step into this direction we study the distributions of the Luttinger parameter K over disorder realizations. It has been shown [32] that in the SF phase the distributions display self-averaging, i.e. they become sharper and narrower for increasing system sizes. This is contrasted by a broadening and flattening in the BG phase [32,48]. Therefore, the change between the two different types of behaviors should signal the SF-BG transition. Fig. 6 shows our analysis for two different values of the disorder strength: the first one contains the crossing of the right phase boundary when moving along the line ∆ = 2, and the second one shows the crossing of the left phase boundary along the line ∆ = 3. We find that the K-distributions are well captured by log-normal distributions with probability density Here x = K −K min (where K min is the smallest measured K value among the n r different disorder realizations), since P (x < 0) = 0 for log-normal distributions. The parameters µ and σ determine the mean and the variance of the distribution. The results in Fig. 6 are consistent with the results in Fig. 5, meaning that we observe the crossover from self-averaging (∂[Var(K)]/∂N s < 0) to broadening (∂[Var(K)]/∂N s > 0) at locations compatible with the criterion K c = 2/3.
Finally, we propose a method for locating the SF-BG transition that requires minimal a priori knowledge on the nature of the transition, but merely relies on the assumption that the averaged hopping correlation functions C(r) (understood as an average of the correlation function C(r) in Eq. (4) over all disorder realizations for each distance r) decay algebraically in the SF phase while in the localized BG phase they decay exponentially. We can conveniently discriminate the different scaling behaviors by fitting C(r) in double logarithmic representation with a function of the form where ξ = log [crd(r)] and α is the curvature of the fit. Clearly, α = 0 indicates a straight line in log-log scale and therefore an algebraic decay, while α > 0 captures the downward bending of the exponential decay in the same plot. It turns out that the increase of α is quite sharp both in the low and high interaction regime, allowing us to apply this method throughout the numerically poorly explored U < 2 region. which is in very good agreement with the K c = 2/3 line.
Let us notice that these results seem to indicate that the SF lobe does not shrink due to rare events, as occurring in the mechanism leading to a "scratched XY universality class" [49]. This might well be due to the system sizes and amounts of disorder realizations we can reasonably access with our TTN method, but it would pose an interesting visibility question for realistic experimental setups that might suffer similar limitations. On the other hand, in the same Fig. 2  . The dashed gray lines are linear fits through the first five data points, serving as guides to the eye to discriminate the different scaling behaviors. The system size is N s = 128, the number of disorder realizations is n r = 100. Increasing n r does not result in changes bigger than the point size, the same applies for finite bond dimension effects. The extracted curvatures α as a function of U are shown in the right plots. The baselines are determined as the average of the data points close to zero; the threshold is used to locate the critical points.
phase boundary we obtained down to U = 0.4 . In conclusion, even though the increasing computational efforts needed to explore the regime of very small interactions U ≪ 1 prevent us from presenting converged and stable results in that region, we consider the phase diagram of Fig. 2 representative of the physics of the system.

Conclusions
In this manuscript we employed our recently introduced gauge-adaptive Tree Tensor Network algorithm [42] to study the disordered Bose-Hubbard model in one dimension.
Thanks to its natural ability to deal with periodic boundary conditions, we have computed the superfluid stiffness, the correlation functions and the scaling of entanglement entropy (see Appendix A) with a high degree of precision. We showed that criteria based on all these three quantities predict the location of the BKT transition for the clean (no disorder) case within an error of a few 10 −3 .
In order to identify the still partially undetermined, lobe-shaped, boundary between the superfluid phase and the Bose glass phase, we considered not only disorder-averaged quantities but also their statistical distributions. First, we showed that on the right side of the lobe (i.e., strongly interacting regime) the averaged stiffness exhibits a jump in the thermodynamic limit that can be fitted with the conventional BKT finite-size scaling (its height being a parameter to fit). Then, we determined the line corresponding to the averaged Luttinger parameter attaining a value K = 2/3, and observed that this is consistent with the line that discriminates between self-averaging and broadening regimes of its statistical distribution. The latter was indeed recently indicated [32] as a more reliable criterion for the left side of the lobe (i.e., weakly interacting regime), where a new universality class might take over the SF-BG transition. Finally, we put forward a criterion based on the decay of disorder-averaged correlation functions: the line determined by their change between algebraic and exponential decay also nicely agrees with the other established criteria.
The high compatibility between the results from the various analysis methods we employed is a strong proof of reliability of our numerical study. In conclusion, hardly any difference is detectable between the "self-averaging" criterion and the K = 2/3 one (or alternatively our averaged correlation decay), at least for what concerns system sizes and disorder repetitions we explored, which however are compatible with the conditions of current experimental verification. Finally, we stress that our data nicely match the analytically predicted critical line at tiny interactions [35], thus reconciling the two pictures.
As an interesting perspective, we mention here the possible extension of the TTN method to the time-dependent case which will allow to consider non-equilibrium situations that are out of reach with current numerical techniques. In this appendix we show that the TTN state representation combined with the PBC setting allows for an accurate determination of critical exponents. For the system sizes N s ≤ 256 used in this work, our computational resources allow to use bond dimensions of up to m = 150. At these bond dimensions the critical exponents fitted to the correlation functions have an accuracy of the order of 10 −3 , which is why no extrapolation in m is needed for our purposes. In Fig. A1 we demonstrate this for the clean system, similar results are obtained for the disordered case.
Finally, in order to complement the data on the BKT transition in the clean system, we use the entanglement entropy (a quantity that is easily accessible by means of tensor network methods [17]) as an alternative criterion to locate the phase transition. To this end, we numerically calculate the von Neumann entropy S(ℓ) for various bipartitions ℓ | N s − ℓ of the ring. Since the SF phase is known to be critical (with central charge c = 1), the entanglement entropy obeys [44] with a a non-universal constant. By fitting Eq. (A.1) (a being the only free parameter) for different values of U, we expect to get an increase in the fit's RMS error ∆ RMS starting from U c , where Eq. (A.1) is not valid any more because in the MI phase the von Neumann entropy S(ℓ) saturates for large ℓ (obeying an area law). As shown in Fig. A2, we find that the behavior of ∆ RMS as a function of U is well described by the form   [18] 0.297 ± 0.010 Kashurnikov and Svistunov [46] 0.304 ± 0.002 which is closely related to the opening of the energy gap E g ∼ exp −δ/ √ U − U c at the BKT transition [45], with some constants ∆ 0 , α, β, γ, δ. By fitting Eq. (A.2) to our data we obtain a numerical value of U c = 3.26 ± 0.13 for the transition point (equivalently, t c = 0.306 ± 0.011), see Fig. A2. Despite the fact that this result is compatible with the ones we obtained by means of the other methods, its precision is about one order of magnitude worse. This is due to the exponentially slow increase of ∆ RMS in Eq. (A.2), which prevents a more accurate determination of U c .
In Tab. A1 we summarize the results for the SF-MI transition points t c obtained from the various techniques that we employed (superfluid density in Section 3, Luttinger parameter in Section 4, entanglement entropy). The mutual agreement of the obtained values is of the order of some 10 −3 . SD: 0.14 ± 0.02 EV: 0.54 ± 0.01 n r = 100  Figure B1. Distributions of the Luttinger parameter K for two different system sizes N s and two different numbers of disorder realizations n r . In both plots the disorder strength is ∆ = 3.0 and the interaction is U = 2.2. The characteristics of the distributions (in particular their expectation value (EV) and standard deviation (SD)) remain within the statistical error when increasing n r from 100 to 500. Fitted lines are log-normal distributions.

Appendix B. Convergence of distributions with respect to the number of disorder realizations
In Fig. B1 we demonstrate that the distributions of the Luttinger parameter K can be meaningfully characterized with a number of disorder realizations of the order of n r = 100, even in the regime of relatively strong disorder ∆. This statement holds true for all system sizes N s used in this work. Moreover, it becomes apparent that the distributions can be faithfully described by the log-normal distributions defined in Eq. (5). We remark that the identification of rare events with very large K, as described in [32,35], would require system sizes and number of disorder realizations of the order of several thousand, which is not within the scope of the present investigation (see main text).
Appendix C. Comparison of phase diagram with data from the literature In this appendix we compare our data on the extension of the superfluid lobe (as summarized in Fig. 2) with results previously obtained in the literature, in particular with data from Prokof'ev et al. [26] and Rapsch et al. [27]. A direct comparison is shown in Fig. C1. As already mentioned in the main text, we observe good agreement with the Monte Carlo results reported in Ref. [26] (except for a slight deviation at the tip of the lobe). Ref. [26] Ref. [27] jump in ρ s K = 2/3 ∆ c ∼ U 3/4 SF BG U ∆ Figure C1. Comparison of the phase diagram from Fig. 2 with phase diagrams previously obtained in the literature.

Appendix D. Local boson occupation numbers
In order to get an estimate for the necessary local (on-site) dimension d, one has to consider the expectation values n j of the local boson number operators and their corresponding variances Var(n j ) = n 2 j − n j 2 . We plot these quantities in Fig. D1, both as a function of interaction strength U (with disorder strength ∆ fixed) and as a function of ∆ (U fixed). Adding two standard deviations to the highest occurring boson occupation expectation values along the ring, we find that up to d = 10 states are relevant for the lowest U values we simulated. By choosing d = 20 for this regime (as reported in the main text) we can be absolutely sure that errors due to the truncation of the local boson bases are negligible. We can easily afford such high local dimensions d because the impact on the computation time is relatively small in TTN simulations (e.g. for a N s = 128 simulation the runtime only increases by a factor of 1.3 when using d = 20 instead of d = 10).
Furthermore, Fig. D1 suggests that for constant U the maximally occurring occupation number increases linearly with the disorder strength ∆.