Superfluid-insulator transition in weakly interacting disordered Bose gases: a kernel polynomial approach

An iterative scheme based on the kernel polynomial method is devised for the efficient computation of the one-body density matrix of weakly interacting Bose gases within Bogoliubov theory. This scheme is used to analyze the coherence properties of disordered bosons in one and two dimensions. In the one-dimensional geometry, we examine the quantum phase transition between superfluid and Bose glass at weak interactions, and we recover the scaling of the phase boundary that was characterized using a direct spectral approach by Fontanesi et al. [Phys. Rev. A 81, 053603 (2010)]. The kernel polynomial scheme is also used to study the disorder-induced condensate depletion in the two-dimensional geometry. Our approach paves the way for an analysis of coherence properties of Bose gases across the superfluid-insulator transition in two and three dimensions.


Introduction
More than 20 years ago, the discovery that superfluidity may be suppressed in 4 He adsorbed on porous media [1][2][3][4] triggered investigations into the conducting and insulating phases of interacting bosons in quenched disorder. In this effort to understand what is now known as the dirty-boson problem [5], most studies focused on the zerotemperature quantum phases, with a variety of approaches including Luttinger-liquid theory [6,7], general scaling arguments [8,9], Bogoliubov theory [10][11][12][13], strongcoupling expansions [14], as well as numerical calculations with Monte-Carlo [15][16][17] and DMRG algorithms [18]. The picture that emerged revealed a rich interplay of bosonic statistics, disorder, repulsive interactions, and commensurability effects in the presence of a lattice. The hallmark of this interplay is the restriction of the superfluid phase to a regime of moderate interactions and weak disorder, surrounded both at weak and strong interactions (or strong disorder) by a compressible gapless insulator called Bose glass [6,9,15,16]. This picture holds for bosons both in the continuum and on a lattice, with the difference that in the case of lattice Bose gases, commensurate fillings give rise to an additional incompressible Mott-insulator phase at strong coupling, provided the disorder is bounded [9]. The Bose glass then intervenes between the superfluid phase and the Mott insulator [9,[17][18][19][20], a feature of the commensurate lattice case that may nevertheless prove difficult to observe in experiments [20,21]. At incommensurate fillings or weak interactions, on the other hand, the lattice case qualitatively resembles the continuous case [6,18,20]. Additionally, in the presence of special symmetries, the lattice Bose gas may exhibit insulating phases of another kind, such as an incompressible and gapless Mott glass [22][23][24]. In spite of an important body of available results, however, the characterization of the generic Bose-glass phase and the superfluid-insulator transition remains to a large extent an open problem. This appears to be the case even for the one-dimensional geometry in view of recent works [25][26][27][28], which highlight the challenges related to the connection of the weakly and strongly interacting regimes and the extension of the ground-state phase diagram to finite temperatures.
The superfluid-insulator transition of disordered bosons attracted a renewed interest in the context of ultracold atoms [29], due to the high degree of control over disorder, interactions and confining potentials achieved in these systems [30][31][32]. While the pioneering experiments with disordered bosons [33][34][35][36][37] aimed at observing Anderson localization in the noninteracting limit [38][39][40][41], more recent ones have provided first results towards a quantitative characterization of the phase diagram for nonvanishing interactions [42][43][44][45][46][47][48]. With this experimental activity, theoretical investigations also turned to the weakly interacting regime, which hitherto had been only poorly characterized. The scaling of the superfluid-insulator phase boundary as a function of the strength of disorder and interactions was established at the mean-field level [49][50][51][52][53], and shown to depend in an essential way on the microscopic disorder correlations. For the 1D geometry, the fragmentation mechanisms driving the transition were analyzed by means of Bogoliubov theory [54], while universal features of the transition and manybody corrections at intermediate disorder strengths were worked out with real-space renormalization group (RG) techniques [25,55]. In the latter approach, making contact with experiments is a challenging task [55], and the method itself is not generalized to higher dimensions in a straightforward way [24].
To date, the details of the superfluid to Bose-glass transition in dimension d > 1 are not well known, and the mechanisms driving the transition are expected to be more complex than in 1D. The notion of weak links and fragmentation, for instance, involves connectivity in higher dimensions i.e. percolation [24,51,56]. Moreover, while providing a natural step towards higher dimensions, the 2D case is particularly interesting in several respect. First, it stands for the marginal dimension of Anderson localization at the single-particle level, for the orthogonal and unitary Wigner-Dyson universality classes [57,58]. Naively, one would therefore expect this geometry to be very sensitive to the introduction of interactions. Second, the clean (disorder-free) weakly interacting system forms a true condensate at T = 0, and an algebraic superfluid for 0 < T < T BKT , where T BKT is the temperature of the Berezinskii-Kosterlitz-Thouless transition [59]. An outstanding question in this respect is how the zero-temperature superfluid-Bose glass transition connects to the clean BKT transition as disorder and temperature are varied. Experiments are ongoing in this regime to characterize the properties of the disordered Bose gas in 2D [60,61].
The features of the 2D superfluid to Bose glass transition beyond the mean field have been addressed recently. In Ref. [62], the Lifshitz-tail physics associated with the deep insulating regime was analyzed by means of a multi-orbital Hartree-Fock method based on a set of low-lying single-particle states [62]. A real-space RG approach was devised for the 2D dirty-boson problem in Ref. [24], and applied to the particle-holesymmetric case, where the insulating phase is an incompressible Mott glass instead of a Bose glass. The analysis also emphasized the possible limitations of the strong-disorder RG in the 2D study. In Ref. [63] the T = 0 phase diagram was studied by means of a weak-disorder expansion of Bogoliubov theory, valid far away from the transition, and quantum Monte-Carlo calculations. Upon finite-size scaling of superfluid fractions obtained numerically, the authors concluded in favor of a smooth crossover between the superfluid and insulating phases. However, the ad-hoc scaling law used in the analysis may deserve a careful examination, as the system sizes used in the numerics cover a relatively modest range. Hence, the development of a method reaching beyond the mean field and affording large system sizes appears highly desirable in order to fill the existing gaps in the understanding of the two-dimensional case.
In this contribution, we present a numerical scheme for the efficient computation of the one-body density matrix of a weakly interacting Bose gas in the framework of Bogoliubov theory, appropriately extended to account correctly for diverging phase fluctuations in low dimensions [64]. The asymptotic behavior of the one-body density matrix determines the superfluid or insulating behavior of the disordered Bose gas, and thereby allows a discussion of the phase diagram. Our scheme accomodates arbitrary disorder strengths, it has no intrinsic limitation in dimensionality, and it admits a straightforward extension to nonvanishing temperatures within the range of validity of Bogoliubov theory. The underlying density-phase representation allows for an accurate description of condensate, quasicondensate and insulating phases in the limit of large densities for any fixed interaction energy. The key feature of our approach is that it is based on an iterative scheme called the kernel polynomial method (KPM) [65], which allows the computation of correlation functions in large systems.
The article is organized as follows. Section 2 provides a reminder on Bogoliubov theory in the density-phase formulation and on the form of the one-body density matrix within that framework. The KPM scheme for the computation of the one-body density matrix is detailed in section 3. In section 4, we validate our scheme by applying it to the case of disordered bosons at T = 0 in 1D and 2D, and by comparing our results with existing literature. In the 1D geometry, we analyze the destruction of quasi-long-range order by disorder, and recover the phase diagram obtained through a direct approach in Refs. [52,54]. In the 2D case, we compute the condensate depletion induced by both interactions and disorder, and compare our findings with those of Ref. [66]. In section 5, we conclude and discuss extensions of the present work.
2. One-body density matrix of weakly interacting Bose gases 2.1. Long-range and quasi-long-range order We consider a dilute gas of Bose particles described by the many-body Hamiltonian whereΨ(r) is the bosonic field operator, g > 0 is the coupling constant parameterizing a repulsive contact interaction, andĤ 0 = − 2 2m ∇ 2 r + V (r) is the single-particle Hamiltonian. In the following, the external potential V is a homogeneous random potential with zero mean, root-mean-square amplitude ∆, and a correlation length η defined as the spatial width of the two-point correlation function V (r)V (r ). Here and in the following, the bar denotes a statistical average over the disorder configurations, while the brackets . indicate quantum-mechanical expectation values.
In the weakly interacting regime and close to the ground state, the properties of the dilute Bose gas are well described by standard Bogoliubov theory in 3D [67]. In this standard approach, the field operatorΨ(r) is split into a classical component Ψ 0 (r) representing a condensate wave function with well-defined phase and a field δΨ(r) describing quantum fluctuations. An effective Hamiltonian is then derived from the expansion ofĤ to second order in δΨ. In dimensions d ≤ 2, however, the Mermin-Wagner-Hohenberg theorem [68,69] rules out the presence of a condensate (i.e. longrange order [70,71]) at any temperature T > 0, as well as T = 0 in 1D. At sufficienty low temperatures in 2D and T = 0 in 1D, the Bose gas nevertheless forms a quasicondensate with a power-law decay of the one-body density matrix G(r, r ) = Ψ † (r)Ψ(r ) (2) at large distances ( r − r → +∞), and exhibits superfluidity [26,59]. While the absence of a true condensate in those cases precludes a straightforward application of standard Bogoliubov theory, the reformulation of the latter in a density-phase representation [72,73] leads to a theory that is free of divergencies, and which captures the algebraic decay of correlation functions in quasicondensates [64,74]. In this formulation the field operator writesΨ(r) = e iθ(r) ρ(r), whereθ(r) is a phase operator that is safely defined in the high-density limit, andρ(r) is the density operator. The latter is split into a classical component ρ 0 (r) corresponding to the mean-field condensate (or quasicondensate) density and a fluctuation δρ(r). As in standard Bogoliubov theory, an effective Hamiltonian is derived from the expansion ofĤ at leading order in the fluctuations. At low temperature, such an approach is valid wherever [64] ρ 0 (r)ξ d 1.
In this expression, d is the spatial dimension and ξ is the healing length, defined here as where U = gρ 0 is the average interaction energy. In a disordered system, the density ρ 0 (r) may locally assume small values, but the regime of validity is recovered at identical U and local interaction energy U (r) = gρ 0 (r) for sufficiently large ρ 0 (i.e. small g). It is also worth noting that, although the regions of low density determine the physics of the superfluid-insulator transition, the latter is also observed within Bogoliubov theory in the limit of asymptotically large average densities (ρ 0 → ∞ and g → 0 with constant U = gρ 0 ), where the theory is expected to be exact [54,75]. Many-body effects beyond Bogoliubov theory arise as subleading terms in an expansion in terms of the inverse density [55,64]. Current experiments with weakly interacting disordered Bose gases fulfill inequality ρ 0 ξ d 1 by a least one or two orders of magnitude (see e.g. Refs. [44,60]), so that the spatial regions where Bogoliubov theory breaks down may be neglected in a good approximation.
While the relation between superfluidity and condensation or quasicondensation is rather subtle, the presence of a compressible superfluid appears both necessary and sufficient for the existence of a condensate or quasicondensate [31]. Thus, the longdistance behavior of the one-body density matrix (2) allows a distinction between superfluid and insulating phases. In the clean (interacting) 1D system, G(r, r ) decays algebraically at T = 0. In the clean 2D system at T = 0, the one-body density matrix exhibits a plateau at long distances, which characterizes a true condensate. For 0 < T < T BKT , the 2D system is also an algebraic superfluid with a power-law decay of correlations. In 3D, finally, the system is a true Bose-Einstein condensate below the critical temperature T BEC . All these phases can be distinguished from the normal phases found at higher temperatures, which exhibit an exponential decay of G(r, r ). Similarly, the complete suppression of superfluidity by disorder at the phase transition to the Bose glass phase coincides with the destruction of long-range order or quasi-long-range order, i.e. with the emergence of an exponential decay of the one-body density matrix [52]. Since the ground state of the interacting Bose gas is globally extended and G(r, r ) depends on the disorder configuration, the disordered phases can be characterized by the long-distance behavior of the statistical average where g 1 is the reduced one-body density matrix with ρ(r) the total gas density [64]. As anticipated above, G(r, r ) and g 1 (r, r ) can be calculated in a density-phase formulation of Bogoliubov theory [64]. Sections 2.2 and 2.3 provide a reminder on both the number-conserving and nonconserving variations of such a theory. These results are used in the following sections for the computation of the one-body matrix with the kernel polynomial method.

Bogoliubov theory in number-conserving and nonconserving approaches
In the ground state, the density ρ 0 obeys the Gross-Pitaevskii equation (GPE) where µ corresponds to the chemical potential in a grand-canonical description of the system. The quantum fluctuations and elementary (quasiparticle) excitations of the Bose gas are described by the fieldB(r) = δρ(r)/(2 ρ 0 (r)) + i ρ 0 (r)θ(r), which obeys the equation of motion [64,76] i ∂ ∂t Here L GP is the standard Bogoliubov operator withĤ GP =Ĥ 0 + gρ 0 (r). The fieldB admits the expansion [64,77] whereb j is a bosonic quasiparticle operator, and the two last terms are explained below. The mode functions u j (r) = r|u j and v j (r) = r|v j are given by the solutions of the usual Bogoliubov-de Gennes equations (BdGEs) The operator L GP is non-Hermitian, but its eigenvalues are real in the ground state of Hamiltonian (1). As L GP * = L GP , the components u j (r) = r|u j and v j (r) = r|v j can be chosen as real-valued functions. We nevertheless consider the more general case of complex u j and v j . For each eigenvector (|u j , |v j ) T with eigenvalue E j > 0, the operator L GP also has an eigenvector (|v * j , |u * j ) T with eigenvalue −E j < 0. The adjoint vectors of these positive and negative eigenvectors are (|u j , −|v j ) T and (−|v * j , |u * j ) T , respectively [78]. The biorthogonality of both positive and negative solutions of the BdGEs is thus expressed by the well-known relation The operator L GP also has an eigenvector pertaining to the eigenvalue E = 0, As the set of eigenvectors of L GP is not complete, an eigenvector (|φ a , |φ a ) T of L 2 GP with eigenvalue zero, such that φ a (r) ∈ R and φ 0 |φ a = 1/2, is added to the set to obtain the closure relation [77] TheP sb andQ sb terms in Eq. (10) account for fluctuations in the particle number and are responsible for the phase diffusion of condensates in symmetry-breaking approaches [79]. These terms do not arise in number-conserving approaches [77,80,81], which retain only fluctuations that are orthogonal to the ground state φ 0 (r). The fieldΛ(r) that describes the orthogonal fluctuations obeys an equation similar to Eq. (8): where L GP has been replaced by [77] andQ = 1 − |φ 0 φ 0 | projects orthogonally to the ground state. Equation (10) is replaced by the modal expansion where u ⊥ j (r) and v ⊥ j (r) are solution of the modified BdGEs The operator L has the same spectrum as L GP , and its positive and negative families of eigenvectors are simply obtained through the projections |u ⊥ j =Q|u j and |v ⊥ j =Q|v j , which leave the biorthogonality relations (12) unaffected. Unlike L GP , however, the operator L is diagonalizable. The zero eigenspace is spanned by the vectors (|φ 0 , 0) T and (0, |φ 0 ) T , so that the resolution of identity writes [77] As we shall see below, Eqs. (7), (15) and (18) [or, equivalently, Eqs. (9) and (13)] are all that is needed for the efficient calculation of the one-body density matrix of disordered Bose gases in the weakly interacting regime.

One-body density matrix in the density-phase representation
The expression of the one-body density matrix G(r, r ) = Ψ † (r)Ψ(r ) was derived in the density-phase formalism in Ref. [64]. At zero temperature, it can be cast into the form [54,75] G(r, r ) = ρ(r)ρ(r ) exp where j enumerates the Bogoliubov modes with E j > 0. This expression is valid in the limit of small density and phase fluctuations, which is realized at large average density ρ 0 for any given interaction energy U = gρ 0 . In this limit, one has ρ(r) ρ 0 (r), which in the presence of a true condensate amounts to a small condensate depletion (see section 4.2). Note also in this respect that, owing to the form of the GPE (7) and BdGEs (17), the v ⊥ j (r) numerators in the exponent of Eq. (19) depend only on the product U = gρ 0 rather than on the coupling constant g and the average density ρ 0 independently. Because of the denominators, the above exponent thus admits a simple scaling as a function of density for fixed interaction energy U .
Remarkably, expression (19) accurately describes weakly interacting Bose gases in any dimension. In particular, it is not plagued by divergences in low dimensions, and captures the power-law decay of the one-body density matrix of quasicondensates in 1D Bose gases at T = 0 [64,74,82]. This expression was also used in Ref. [54] in conjunction with a numerical diagonalization of the Bogoliubov operator to analyze the destruction of quasi-long-range order by disorder in the 1D geometry. While Eq. (19) involves all Bogoliubov modes, the sum is typically dominated by the modes of the low-energy phonon regime. However, even with a restriction to low-energy modes, the calculation of the one-body density matrix G(r, r ) through complete or partial diagonalization of the Bogoliubov operator becomes prohibitive for large system sizes. In the following section, we present an alternative scheme, based on the kernel polynomial method [65], which circumvents the solution of the Bogoliubov eigenvalue problem and constitutes the main result of the present work.

Kernel polynomial scheme for the one-body density matrix
The kernel polynomial method (KPM) [65,83,84] is an iterative numerical scheme for the computation of correlation functions. The KPM bypasses the spectral decomposition of the operators involved in those correlation functions, which may be numerically intractable. The KPM technique and some applications have been recently reviewed in Ref. [65]. In section 3.1 we introduce the elementary aspects of KPM that are relevant to the present study. In section 3.2, we show how a KPM scheme can be devised to compute the one-body density matrix G(r, r ) on the basis of Eq. (19).

Basics of the kernel polynomial method
The core ingredient of KPM is the expansion of correlation functions on Chebyshev polynomials of the first kind. The latter are orthogonal polynomials on the interval I = [−1, 1], defined by the recurrence relation T n+1 (x) = 2xT n (x) − T n−1 (x) with T 0 (x) = 1 and T 1 (x) = x. Consider a Hermitian operatorX with a discrete or continuous spectrum contained in I, and the correlation function The latter is a formal writing for where {|x j,α } is an orthonormal eigenbasis ofX, which provides the spectral decompositionX = j,α x j |x j,α x j,α |, with x j ∈ I, and the resolution of identity The above correlation function has the expansion where the coefficient µ n (|a , |b ), called Chebyshev moment of order n, is defined as Owing to Eqs. (21) and (22), the moments of f (|a , |b , x) are simply given by which boils down to the matrix element of a polynomial ofX: Instead of calculating T n (X) for each new n index, which is computationally costly, one takes advantage of the recurrence relation between Chebyshev polynomials and keeps track of two vectors, |b n = T n (X)|b and |b n−1 = T n−1 (X)|b , with the initialization |b 0 = |b and |b 1 =X|b . Then, a single application ofX (matrix-vector multiplication) yields the new vector |b n+1 = 2X|b n − |b n−1 , along with the next Chebyshev moment µ n+1 (|a , |b ) = a|b n+1 . In practice, the expansion (23) is truncated at some finite order N , and f (|a , |b , x) is approximated by where the g (N ) n factors are introduced to damp the Gibbs oscillations caused by the truncation [65]. These factors amount to the action of a convolution kernel on f (|a , |b , x), whence the name of KPM. Finally, the functional dependence of f (N ) (|a , |b , x) on x is usually efficiently computed for a set of points x k ∈ I by using a discrete cosine transform [65].
In summary, the KPM offers a simple iterative scheme for the calculation of correlation functions akin to expression (20), and avoids the numerical complexity associated with the spectral representation (21). Let us now examine how this method can be applied to the Bogoliubov operator for the computation of G(r, r ).

Calculation of the one-body density matrix
Expression (19) for the one-body density matrix involves the eigenmodes of the Bogoliubov operator L, the spectrum of which lies on the real line. As all subsequent numerical calculations are carried out with a finite-difference scheme and finite-size systems, the spectrum of L has a compact support [−E max , E max ], where E max depends on the hopping term t associated with the Laplacian in the finite-difference scheme, the strength of interactions and the disorder configuration V (r). In the calculations presented in section 4, the upper bound E max is obtained by solving the sparse Bogoliubov eigenvalue problem (11) for the largest eigenvalue with a Lanczos method. Taking a slightly larger E max to ensure good KPM convergence at the spectrum boundaries [65], the spectrum is mapped to I by the rescaling L → L/E max .
To exhibit correlation functions akin to Eqs. (20) and (21), we cast expression (19) into the following form: with and where |w ⊥ k,2 and |w ⊥ k,2 are the second components of the eigenvector (|w k,1 , |w k,2 ) T of L and its adjoint vector (|w ⊥ k,1 , |w ⊥ k,2 ) T , respectively. In Eq. (30), the index k runs over the positive (E k > 0), null (E k = 0) and negative (E k < 0) families of eigenvectors, and k = E k /E max . Because of the integration boundaries in Eq. (28), the term f (r, r , ) contributes to the exponent by where j enumerates the modes with E j > 0. The −1/(2N 0 ) term, which stems from the zero eigenvector (0, |φ 0 ) T , cancels out of the f sum in Eqs. (28) and (29) Indeed, rewriting Eq. (30) as we find that the Chebyshev moments of f (r, r , ) are and those of F (r, r , ) follow as M n (r, r ) = µ n (r, r) − µ n (r, r ) − µ n (r , r) + µ n (r , r ).
Finally, there is no need for a Chebyshev inversion by discrete cosine transform, as the expansion (23) can be integrated analytically on [0, 1], and we obtain Note that the contributions of all even moments except M 0 (r, r ) are integrated out on [0, 1]. The moments µ n (r, r ) are calculated iteratively following the recurrence scheme outlined in section 3.1, for the four (r, r ) pairs in M n (r, r ). This requires only two Chebyshev sequences as, for instance, T n (L/E max )(0, |r ) T may be used to compute µ n (r, r ) and µ n (r , r ). When the Chebyshev iterations are truncated at order N , the reduced one-body density matrix is approximated by where the g (N ) n are convolution kernel factors (see section 3.1). We found the standard Jackson kernel [65] to be suitable in this scheme.
The Chebyshev iterations based on Eq. (34) involve projections orthogonally to the ground state by means ofQ = 1 − |φ 0 φ 0 |. Interestingly, the Bogoliubov operator L may be replaced by L GP in the iterations, so that projections are not necessary. This provides a further simplification of the KPM calculation of the one-body density matrix. Indeed, upon expansion with v ⊥ j (r) = v j (r) − φ 0 |v j φ 0 (r), one easily sees that Eq. (19) also writes as The function F GP (r, r , ) is defined as the analogue of F (r, r , ) in Eqs. (28) and (29), with four terms of the form where k enumerates the eigenmodes of L GP with E k = 0. Owing to the closure relation (13), the Chebyshev moments of f GP (r, r , ) read where Given that L GP (|φ 0 , −|φ 0 ) T = (0, 0) T and L GP (|φ a , |φ a ) T = α(|φ 0 , −|φ 0 ) T with constant α, as detailed in Ref. [77], we find 2p (r, r ) = (−1) p φ a (r) N 0 ρ 0 (r) These R (i) n (r, r ) terms cancel out of the sum M GP n (r, r ) = µ GP n (r, r) − µ GP n (r, r ) − µ GP n (r , r) + µ GP n (r , r ).
Hence, the comparison of Eqs. (34) and (40) shows that L GP can be used instead of L in the Chebyshev iterations underlying Eq. (37). The expressions (34), (35) and (37) are the main results of this section. They provide an iterative scheme for the calculation of the reduced one-body density matrix of a weakly interacting Bose gas, once the solution of the GPE (7) is given. In the following section, this scheme is applied in various geometries, with and without disorder.
Remarkably, the above approach can be extended in a straightforward way to nonzero temperatures within the framework of Bogoliubov theory. In addition to the v j (r) terms, the thermal G(r, r ) also contains contributions from the Bogoliubov components u j (r), which can be calculated through additional KPM iterations. As the u j and v j terms are weighted by Bose-Einstein occupation factors, the integration in Eq. (36) no longer has a simple analytical solution. However, this analytical step may be replaced by a discrete cosine transform at low computational expense.
Our results also show that the operator involved in the correlation function and the KPM iterations need not be Hermitian. Some illustrations of this fact can be found in the literature, with special cases such as the computation of retarded Green's functions [65] or the solution of generalized Hermitian eigenvalue problems [85]. The Bogoliubov operator L and L GP provide two other interesting examples in that context. First, L is diagonalizable, albeit non-Hermitian, and its eigenvalues are real. As a consequence, the eigenvectors and their adjoints form a complete biorthogonal set that can be used for the closure relation, and there is no need for a twofold KPM iteration to retrieve the spectral information on a compact set of the complex plane. Second, L GP is not diagonalizable, and yet its generalized eigenvectors (and their adjoints) can be used for the closure relation. In the derivation of Eqs. (43) to (46), we took advantage of the fact that (|φ 0 , −|φ 0 ) T and (|φ a , |φ a ) T are generalized eigenvectors of low order, so that the action of T n (L GP ) can be evaluated easily.

Application to disordered bosons
We now employ the KPM scheme introduced above to analyze the asymptotic behavior of the one-body density matrix of disordered bosons in 1D and 2D. While the approach of the previous sections is general, we consider here a Gaussian random potential V with Gaussian correlation function The numerical procedure for the computation of g 1 (r, r ) is the following. The continuum problem (1) is represented on a finite volume L d in a finite-difference scheme with lattice spacing = L/n . To emulate the continuum limit, the hopping term t = 2 /(2m 2 ) is chosen to be much larger than all the other energy scales of the problem. In all the calculations presented here, the correlation length of the disorder is taken to be η = 4 , which is sufficient for our purposes. The correlation length η and the associated energy are used as reference scales throughout this section, even in the absence of disorder (∆ = 0). For each configuration V (r), the ground-state solution ρ 0 (r) of the GPE (7) and the corresponding chemical potential µ are computed through imaginarytime propagation with a standard Crank-Nicolson scheme, at fixed particle number N 0 = drρ 0 (r). We denote by the average mean-field interaction energy. Periodic boundary conditions are imposed on the GPE and on the BdGEs (11). The Bogoliubov operator L GP is set up on the basis of V (r) and ρ 0 (r). Then, g 1 (r, r ) may be calculated with the KPM iteration detailed in the last section for any (r, r ) pair.

Superfluid to Bose glass transition in one dimension
In Refs. [52,54], the destruction of quasi-long-range order was used as a signature of the superfluid to Bose glass transition at T = 0 in 1D, and this criterion was used to draw the quantum phase diagram on the basis of the asymptotic behavior of the (averaged) reduced one-body density matrix g 1 (|r − r |). Here, we use the 1D setting as a testbed for the KPM approach described above. In the absence of disorder, g 1 is expected to decay at large distances with a power law given by [64,86]: where C = 0.57721 . . . is Euler's constant, the density ρ can be approximated by ρ 0 within our Bogoliubov approach, and ξ is the healing length defined in Eq. (4). Figure 1 shows the result of a KPM calculation for a system of length L = 2 16 η, and the excellent agreement obtained with the power law (51) for |r − r | ξ, i.e. in its regime of validity. The regrowth observed at large |r − r | is due to the periodic boundary conditions, and does not affect significantly the data for |r − r | L/4. In all the subsequent analyses, we retain only this range for determining the asymptotic behavior of g 1 . Note the large system size achieved in this computation. In this homogeneous case, a single KPM iteration suffices for the computation of g 1 (r, r ). The number of moments N required to resolve all individual Bogoliubov modes in the phonon regime and achieve a good convergence of the KPM result grows linearly with the system size n = L/ . The required storage space is also of a few n . This has to be compared to the storage space and computation time required for a full diagonalization of L GP , which scale as n 2 and n 3 respectively [65]. For ∆ > 0, the reduced one-body density matrix g 1 (r, r ) depends on the disorder configuration, and is no longer translation invariant. Fig. 2 shows the behavior of g 1 (0, r) for a single disorder configuration, and the results obtained with the KPM scheme for various moment numbers N [see Eq. (37)]. When N is sufficiently high, the KPM results coincide with those obtained from complete diagonalization, and faithfully reproduce the spatial details of g 1 (r, r ). As a general trend, we also observe that the KPM estimates of g 1 (r, r ) converge from above with increasing N . This can be attributed to the fact that low N values imply poor spectral resolution, and hence are not able to resolve the small energy scales associated with the long-distance decay of g 1 , such as for example the vanishing energy separation that can arise when Bogoliubov modes are strongly localized in different spatial regions. While the number of moments required for convergence scales linearly with the system size in the clean case, we also observed that this scaling is slightly faster than linear in the disordered case.
After averaging over disorder, the one-body density matrix exhibits either a powerlaw or an exponential decay a large distances, depending on the strength of interactions and disorder. For long-enough systems a similar behavior may already be observed qualitatively with a single disorder configuration in the spatial average Fig. 3 shows spatial averages computed with the KPM scheme for two sets of parameters in a system of length L = 512 η. We display here the range r L/4, and the crossover to the long-distance behavior is visible on both panels. The linear behavior found for U = 1.12 E c on the double logarithmic scale indicates a power-law decay of the averaged g 1 , corresponding to the superfluid phase. For U = 0.3 E c , on the other hand, we find an exponential decay indicating a Bose glass. For the same system size, the convergence of the KPM result in the insulating phase requires a higher number of moments than in the superfluid phase. This may be attributed to the increase of the Bogoliubov density of states near zero energy and to a reduced level repulsion. Indeed, as the disorder increases, the system turns progressively into a collection of superfluid puddles separated by high potential barriers. As a consequence, the efficiency of the KPM scheme is reduced deep in the insulating regime, where the number of moments required for convergence typically becomes very large. In the superfluid regime and in the parameter range of interest around the superfluid-insulator transition, however, our KPM approach converges quickly and outperforms complete diagonalization in terms of memory usage and computation time already for system sizes as moderate as those used in Fig. 3.
The quantum phase diagram of the weakly interacting regime can be drawn by varying ∆ and U , and determining for each parameter set whether the asymptotic part of the disorder-averaged g 1 behaves as a power law or an exponential. This procedure is put on a systematic footing by fitting the long-distance part of g 1 (r) with both power laws and exponentials, and monitoring the fit quality. Fig. 4 shows the phase diagram obtained with the KPM technique. The blue circles and green squares have been identified as belonging to the superfluid and Bose glass phases, respectively. The black triangles cannot be attributed to one phase or the other with the accuracy of the data, and are assumed to lie on the phase boundary. The red solid lines are fits to the black triangles in the regimes U E c and U E c . In the white-noise regime U E c (i.e. ξ η), the phase boundary is expected to scale as ∆/E c ∼ (U/E c ) 3/4 , while in the Thomas-Fermi regime U E c (i.e. ξ η) the critical ∆ is expected to grow linearly with U [51,52,54,87]. The fits of Fig. 4 are in good agreement with these scaling laws. Our findings thus reproduce the results of Refs. [52,54] without the need for partial or complete diagonalization. This validates our approach.

Condensate depletion in two dimensions
In the absence of disorder, weakly interacting bosons form a true condensate at T = 0 in 2D. In this case, the reduced one-body density matrix g 1 (r, r ) tends to a constant equal to the condensate fraction ρ 0 /ρ for r−r → ∞. To leading order in the strength of interactions, the condensate fraction is given by [59,64,66,88] where g = 2mg/ 2 and g is the 2D coupling constant. While scattering theory shows that this constant actually depends logarithmically on the 2D density ρ and the 3D scattering length a [59], we take it as a given parameter of Hamiltonian (1) in 2D, even for the strongly inhomogeneous case. According to Eq. (53), interactions reduce the condensate fraction. In the presence of disorder, the condensate fraction is expected to be further reduced (but nonzero) in the superfluid regime, and completely suppressed in the Bose glass phase, due to the (exponential) decay of the one-body density matrix.
The KPM algorithm introduced above is expected to reduce significantly the computational cost of a study of the superfluid-insulator transition on the basis of the one-body density matrix. A fully-fledged analysis of this transition nevertheless lies beyond the scope of the present paper, and is left for future work. Here, we restrict ourselves to the superfluid regime and consider the calculation of the disorder-induced condensate depletion as a benchmark for the KPM technique. This depletion has been calculated analytically in Ref. [66] for the limit of weak interactions and weak disorder. To leading order in ∆/U , the depletion δρ = ρ(r) − ρ 0 (r) reads where δρ (0) is the interaction-induced condensate depletion in a clean system, given by Eq. (53). The function h depends only on the ratio of the disorder correlation length η and the healing length ξ. Note that h differs from the similar function introduced in Ref. [66] by a trivial factor √ 2 in the argument due to a different definition of ξ. The left panel of Fig. 5 shows the result of a KPM calculation of g 1 (0, r) in a clean 2D system. Starting from the origin, the one-body density matrix drops and reaches a plateau beyond a few healing lengths. The regrowth of g 1 on the system edges is due to the periodic boundary conditions, as in the 1D case. The condensate fraction is extracted from the value assumed at the center of the system. With this procedure, we studied the dependence of the condensate fraction on the interaction strength. The numerical results are plotted in Fig. 6, and stand in perfect agreement with Eq. (53).
The disordered case was examined with a similar procedure. The right panel of Fig. 5 displays the values of g 1 (0, r) obtained for a single disorder configuration with ∆/U = 0.125. For this value the effect of disorder is weak, and the one-body density matrix still exhibits a plateau, at roughly the same level than the clean case shown in the left panel. As a definition of the averaged condensate fraction in the disordered case, we use the asymptotic value of the disorder-averaged reduced one-body density matrix g 1 (r, r ) at large distances. This definition follows the Penrose-Onsager criterion for Bose-Einstein condensation, based on off-diagonal long-range order [70], and agrees with the definition of the condensate fraction in the clean case. In the presence of disorder, this definition yields a condensate fraction equal to one at the mean-field level (where the Bose gas is described solely by the Gross-Pitaevskii equation), and properly takes into account the role of condensate deformation [66]. Because of this deformation of the condensate in the presence of an inhomogeneous potential, the condensate depletion cannot simply be associated to the fraction of atoms with momenta k = 0. More generally, in the presence of disorder, the superfluid component should be characterized by spatial inhomogeneity. Hence, the superfluid fraction is not expected to be related in a straightforward way to the fraction of atoms with k = 0. Even if the two quantities might vanish simultaneously at the superfluid-insulator phase boundary, the precise relation between the condensate fraction, the superfluid fraction and the fraction of  Figure 6. Condensate fraction ρ 0 /ρ versus interaction strength for a 2D system. The interaction strength is parametrized by the reduced coupling constant g = 2mg/ 2 , where g is the 2D coupling constant. The data points are the values extracted from numerical calculations of the one-body density matrix for particle densities N 0 /L 2 = 160 η −2 (blue circles; system size 64η × 64η) and N 0 /L 2 = 10 η −2 (green crosses; system size 256η × 256η), and variable interactions energies U = gN 0 /L 2 . The solid red line is a fit to the data with N 0 /L 2 = 160 η −2 . The fitted slope agrees with the factor 1/(8π) predicted by Eq. (53) within 2.5%. The other set of data is not fitted for the sake of clarity. The data for N 0 /L 2 = 160 η −2 and g = 0.2 corresponds the plateau in the first panel of Fig. 5. atoms with vanishing momentum in the presence of disorder remains an open issue. With the above definition, the statistical average and fluctuations of the condensate fraction are evaluated by calculating g 1 (0, r), with r at the system center, for several disorder configurations. Figure 7 shows the average and the standard deviation of the fractional depletion 1 − lim r→∞ g 1 (0, r) (i.e. the complement of the condensate fraction) as a function of the disorder strength. For weak disorder (∆ U ), we indeed find a good agreement with the theoretical prediction (54), as shown in the inset. For ∆ 1.3 E c (i.e. ∆ 0.2 U ), however, the numerical averages clearly lie below the theoretical curve. While differences due to the averages used in Eq. (54) and in g 1 are not excluded, this discrepancy is likely to be due to the breakdown of leading-order perturbation theory in the disorder amplitude. The results of Fig. 7 thus suggest that higher orders in the weak-disorder expansion may reduce the condensate depletion.

Summary and outlook
We have developed an iterative scheme, based on the kernel polynomial method, for the efficient computation of the one-body density matrix of weakly-interacting Bose gases in the framework of Bogoliubov theory. Such a scheme is relevant for regimes of strong disorder, which cannot be tackled analytically. The scheme was applied to the case of disordered bosons at T = 0 in one and two dimensions. In the one-dimensional case, we characterized the superfluid-insulator phase transition on the basis of the long-range behavior of the one-body density matrix, and successfully reproduced the results of Refs. [52,54] with a low computational overhead. In the two-dimensional geometry, we analyzed the quantum depletion induced by interaction and disorder in the superfluid regime, and found a good agreement with results available for the weakly disordered regime. These case studies validate our approach and suggest that it may be used to study the coherence properties of weakly interacting Bose systems for system sizes that remain hardly tractable with other numerical techniques. This feature is particularly interesting for investigations into the superfluid-insulator transition in higher dimensions. As outlined here, our approach is also easily extended to regimes of low but nonzero temperatures, which are relevant to ongoing experiments with ultracold atomic gases.