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

We consider a dilute gas of Bose particles described by the many-body Hamiltonian

$$\begin{equation} \hat{H}=\int \mathrm{d}\boldsymbol{r} \left[ \hat{\Psi}^\dagger(\boldsymbol{r}) \hat{H_0} \hat{\Psi}(\boldsymbol{r})+\frac{g}{2} \hat{\Psi}^\dagger(\boldsymbol{r})\hat{\Psi}^\dagger(\boldsymbol{r})\hat{\Psi}(\boldsymbol{r})\hat{\Psi}(\boldsymbol{r}) \right], \end{equation} \tag{ 1 }$$

where $\hat {\Psi }(\boldsymbol {r})$ is the bosonic field operator, g > 0 is the coupling constant parameterizing a repulsive contact interaction and $\hat {H}_0=-\frac {\hbar ^2}{2m}\nabla _{\boldsymbol {r}}^2+V(\boldsymbol {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 $\overline {V(\boldsymbol {r})V(\boldsymbol {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 three dimensions (3D) [67]. In this standard approach, the field operator $\hat {\Psi }(\boldsymbol {r})$ is split into a classical component Ψ₀(r) representing a condensate wave function with a well-defined phase and a field $\delta \hat {\Psi }(\boldsymbol {r})$ describing quantum fluctuations. An effective Hamiltonian is then derived from the expansion of $\hat {H}$ to second order in $\delta \hat {\Psi }$. In dimensions d ⩽ 2, however, the Mermin–Wagner–Hohenberg theorem [68, 69] rules out the presence of a condensate (i.e. long-range order [70, 71]) at any temperature T > 0, as well as T = 0 in 1D. At sufficiently 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

$$\begin{equation} G(\boldsymbol{r},\boldsymbol{r'})=\langle \hat{\Psi}^\dagger(\boldsymbol{r})\hat{\Psi}(\boldsymbol{r'}) \rangle \end{equation} \tag{ 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 divergences, and which captures the algebraic decay of correlation functions in quasicondensates [64, 74]. In this formulation the field operator reads $\hat {\Psi }(\boldsymbol {r})={\mathrm { e}}^{{\mathrm { i}}\hat {\theta }(\boldsymbol {r})}\sqrt {\hat {\rho }(\boldsymbol {r})}$, where $\hat {\theta }(\boldsymbol {r})$ is a phase operator that is safely defined in the high-density limit and $\hat {\rho }(\boldsymbol {r})$ is the density operator. The latter is split into a classical component ρ₀(r) corresponding to the mean-field condensate (or quasicondensate) density and a fluctuation $\delta \hat {\rho }(\boldsymbol {r})$. As in standard Bogoliubov theory, an effective Hamiltonian is derived from the expansion of $\hat {H}$ at leading order in the fluctuations. At low temperature, such an approach is valid wherever [64]

$$\begin{equation} \rho_0(\boldsymbol{r})\xi^d\gg1. \end{equation} \tag{ 3 }$$

In this expression, d is the spatial dimension and ξ is the healing length, defined here as

$$\begin{equation} \xi=\frac{\hbar}{\sqrt{mU}}, \end{equation} \tag{ 4 }$$

where $U=g\overline {\rho _0}$ is the average interaction energy. In a disordered system, the density ρ₀(r) may locally assume small values, but the regime of validity is recovered at identical U and local interaction energy U(r) = gρ₀(r) for sufficiently large $\overline {\rho _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 ($\overline {\rho _0}\to \infty$ and g → 0 with constant $U=g\overline {\rho _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 $\overline {\rho _0}\xi ^d\gg 1$ by at least one or two orders of magnitude (see e.g. [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 long-distance 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

$$\begin{equation} \overline{g_1}(\Vert\boldsymbol{r}-\boldsymbol{r'}\Vert)=\overline{g_1(\boldsymbol{r},\boldsymbol{r'})}, \end{equation} \tag{ 5 }$$

where g₁ is the reduced one-body density matrix

$$\begin{equation} g_1(\boldsymbol{r},\boldsymbol{r'})=\frac{G(\boldsymbol{r},\boldsymbol{r'})}{\sqrt{\rho(\boldsymbol{r})\rho(\boldsymbol{r'})}}, \end{equation} \tag{ 6 }$$

with ρ(r) the total gas density [64].

As anticipated above, G(r,r') and g₁(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 KPM.

In the ground state, the density ρ₀ obeys the Gross–Pitaevskii equation (GPE)

$$\begin{equation} [\hat{H}_0+g\rho_0(\boldsymbol{r})]\sqrt{\rho_0(\boldsymbol{r})}=\mu \sqrt{\rho_0(\boldsymbol{r})}, \end{equation} \tag{ 7 }$$

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 field $\hat {B}(\boldsymbol {r})=\delta \hat {\rho }(\boldsymbol {r})/(2\sqrt {\rho _0(\boldsymbol {r})})+\mathrm {i} \sqrt {\rho _0(\boldsymbol {r})}\hat {\theta }(\boldsymbol {r})$, which obeys the equation of motion [64, 76]

$$\begin{equation} \mathrm{i} \hbar \frac{\partial}{\partial t} \left(\begin{array}{@{}c@{}} \hat{B} \\ \hat{B}^\dagger \end{array}\right) = \mathcal{L}_{\mathrm{GP}} \left(\begin{array}{@{}c@{}} \hat{B} \\ \hat{B}^\dagger \end{array}\right). \end{equation} \tag{ 8 }$$

Here $\mathcal {L}_{\mathrm {GP}}$ is the standard Bogoliubov operator

$$\begin{equation} \mathcal{L}_{\mathrm{GP}} = \left(\begin{array}{@{}cc@{}} \hat{H}_{\mathrm{GP}}+g\rho_0(\boldsymbol{r})-\mu & g\rho_0(\boldsymbol{r}) \\ -g\rho_0(\boldsymbol{r}) & -\hat{H}_{\mathrm{GP}}-g\rho_0(\boldsymbol{r})+\mu \end{array}\right), \end{equation} \tag{ 9 }$$

with $\hat {H}_{\mathrm {GP}}=\hat {H_0}+g\rho _0(\boldsymbol {r})$. The field $\hat {B}$ admits the expansion [64, 77]

$$\begin{equation} \hat{B}(\boldsymbol{r})= \sum_j \left[u_j(\boldsymbol{r})\hat{b}_j+v_j^*(\boldsymbol{r})\hat{b}_j^\dagger\right] -\mathrm{i}\sqrt{N_0}\phi_0(\boldsymbol{r})\hat{Q}_{\mathrm{sb}}+\frac{\phi_a(\boldsymbol{r})}{\sqrt{N_0}}\hat{P}_{\mathrm{sb}}, \end{equation} \tag{ 10 }$$

where $\hat {b}_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)

$$\begin{equation} \mathcal{L}_{\mathrm{GP}} \left(\begin{array}{@{}c@{}} \vert u_j \rangle \\ \vert v_j \rangle \end{array}\right) =E_j \left(\begin{array}{@{}c@{}} \vert u_j \rangle \\ \vert v_j \rangle \end{array}\right) \quad(E_j>0). \end{equation} \tag{ 11 }$$

The operator $\mathcal {L}_{\mathrm {GP}}$ is non-Hermitian, but its eigenvalues are real in the ground state of Hamiltonian (1). As ${\mathcal {L}_{\mathrm {GP}}}^*=\mathcal {L}_{\mathrm {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 $\mathcal {L}_{\mathrm {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

$$\begin{equation} \langle u_j \vert u_{j'} \rangle-\langle v_j \vert v_{j'} \rangle =\int \mathrm{d}\boldsymbol{r} [u_j^*(\boldsymbol{r})u_{j'}(\boldsymbol{r})-v_j^*(\boldsymbol{r})v_{j'}(\boldsymbol{r})] =\delta_{j{j'}}. \end{equation} \tag{ 12 }$$

The operator $\mathcal {L}_{\mathrm {GP}}$ also has an eigenvector pertaining to the eigenvalue E = 0, namely (|ϕ₀〉,−|ϕ₀〉)^T, where $\phi _0(\boldsymbol {r})=\sqrt {\rho _0(\boldsymbol {r})/N_0}$ is the normalized ground state, with $N_0=\int \mathrm {d}\boldsymbol {r} \rho _0(\boldsymbol {r})$. As the set of eigenvectors of $\mathcal {L}_{\mathrm {GP}}$ is not complete, an eigenvector (|ϕ_a〉,|ϕ_a〉)^T of $\mathcal {L}_{\mathrm {GP}}^2$ with eigenvalue zero, such that $\phi _a(\boldsymbol {r})\in {\mathbb{R}}$ and 〈ϕ₀|ϕ_a〉 = 1/2, is added to the set to obtain the closure relation [77]

$$\begin{eqnarray} &&\fl {\mathbbm{1}}= \left(\begin{array}{@{}cc@{}}\vert \phi_0 \rangle\\ -\vert \phi_0 \rangle\end{array}\right) (\langle \phi_a \vert,-\langle \phi_a \vert) +\left(\begin{array}{@{}cc@{}}\vert \phi_a \rangle\\ \vert \phi_a \rangle\end{array}\right) (\langle \phi_0 \vert,\langle \phi_0 \vert) \nonumber\\ &&+\sum_{j} \left(\begin{array}{@{}cc@{}}\vert u_j \rangle\\ \vert v_j \rangle\end{array}\right) (\langle u_j \vert,-\langle v_j \vert) +\left(\begin{array}{@{}cc@{}}\vert v_j^* \rangle\\ \vert u_j^* \rangle\end{array}\right) (-\langle v_j^* \vert,\langle u_j^* \vert). \end{eqnarray} \tag{ 13 }$$

The $\hat {P}_{\mathrm {sb}}$ and $\hat {Q}_{\mathrm {sb}}$ terms in equation (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 ϕ₀(r). The field $\hat {\Lambda }(\boldsymbol {r})$ that describes the orthogonal fluctuations obeys an equation similar to equation (8):

$$\begin{equation} \mathrm{i} \hbar \frac{\partial}{\partial t} \left(\begin{array}{@{}c@{}} \hat{\Lambda} \\ \hat{\Lambda}^\dagger \end{array}\right)= \mathcal{L} \left(\begin{array}{@{}c@{}} \hat{\Lambda} \\ \hat{\Lambda}^\dagger \end{array}\right), \end{equation} \tag{ 14 }$$

where $\mathcal {L}_{\mathrm {GP}}$ has been replaced by [77]

$$\begin{equation} \mathcal{L} = \left(\begin{array}{@{}cc@{}} \hat{H}_{\mathrm{GP}}+g\hat{Q}\rho_0(\boldsymbol{r})\hat{Q}-\mu & g\hat{Q}\rho_0(\boldsymbol{r})\hat{Q} \\ -g\hat{Q}\rho_0(\boldsymbol{r})\hat{Q} & -\hat{H}_{\mathrm{GP}}-g\hat{Q}\rho_0(\boldsymbol{r})\hat{Q}+\mu \end{array}\right), \end{equation} \tag{ 15 }$$

and $\hat {Q}= {\mathbbm{1}}-\vert \phi _0 \rangle \langle \phi _0 \vert$ projects orthogonally to the ground state. Equation (10) is replaced by the modal expansion

$$\begin{equation} \hat{\Lambda}(\boldsymbol{r})={\sum}_j u^\perp_j(\boldsymbol{r})\hat{b}_j+{v^\perp_j}^*(\boldsymbol{r})\hat{b}_j^\dagger, \end{equation} \tag{ 16 }$$

where u^⊥_j(r) and v^⊥_j(r) are solutions of the modified BdGEs

$$\begin{equation} \mathcal{L} \left(\begin{array}{@{}c@{}} \vert u^\perp_j \rangle \\ \vert v^\perp_j \rangle \end{array}\right)= E_j \left(\begin{array}{@{}c@{}} \vert u^\perp_j \rangle \\ \vert v^\perp_j \rangle \end{array}\right) \quad(E_j>0). \end{equation} \tag{ 17 }$$

The operator $\mathcal {L}$ has the same spectrum as $\mathcal {L}_{\mathrm {GP}}$, and its positive and negative families of eigenvectors are simply obtained through the projections $\vert u^\perp _j \rangle =\hat {Q}\vert u_j \rangle$ and $\vert v^\perp _j \rangle =\hat {Q}\vert v_j \rangle$, which leave the biorthogonality relations (12) unaffected. Unlike $\mathcal {L}_{\mathrm {GP}}$, however, the operator $\mathcal {L}$ is diagonalizable. The zero eigenspace is spanned by the vectors (|ϕ₀〉,0)^T and (0,|ϕ₀〉)^T, so that the resolution of identity reads [77]

$$\begin{eqnarray} &&\fl {\mathbbm{1}}= \left(\begin{array}{@{}cc@{}}\vert \phi_0 \rangle\\ 0\end{array}\right) (\langle \phi_0 \vert,0) +\left(\begin{array}{@{}cc@{}}0\\ \vert \phi_0 \rangle\end{array}\right) (0,\langle \phi_0 \vert) \nonumber\\ &&+\sum_{j} \left(\begin{array}{@{}cc@{}}\vert u^\perp_j \rangle\\ \vert v^\perp_j \rangle\end{array}\right) (\langle u^\perp_j \vert,-\langle v^\perp_j \vert) +\left(\begin{array}{@{}cc@{}}\vert v^{\perp*}_j \rangle\\ \vert u^{\perp*}_j \rangle\end{array}\right) (-\langle v^{\perp*}_j \vert,\langle u^{\perp*}_j \vert). \end{eqnarray} \tag{ 18 }$$

As we shall see below, equations (7), (15) and (18) (or, equivalently, equations (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.

The expression of the one-body density matrix $G(\boldsymbol {r},\boldsymbol {r'})=\langle \hat {\Psi }^\dagger (\boldsymbol {r})\hat {\Psi }(\boldsymbol {r'})\rangle$ was derived in the density-phase formalism in [64]. At zero temperature, it can be cast into the form [54, 75]

$$\begin{equation} G(\boldsymbol{r},\boldsymbol{r'})=\sqrt{\rho(\boldsymbol{r})\rho(\boldsymbol{r'})}\exp\!\left[-\frac{1}{2}\sum_j \left|\frac{v^\perp_j(\boldsymbol{r})}{\sqrt{\rho_0(\boldsymbol{r})}}-\frac{v^\perp_j(\boldsymbol{r'})}{\sqrt{\rho_0(\boldsymbol{r'})}} \right|^2\right], \end{equation} \tag{ 19 }$$

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 $\overline {\rho _0}$ for any given interaction energy $U=g\overline {\rho _0}$. In this limit, one has ρ(r) ≃ ρ₀(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 equation (19) depend only on the product $U=g\overline {\rho _0}$ rather than on the coupling constant g and the average density $\overline {\rho _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 [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 equation (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 KPM [65], which circumvents the solution of the Bogoliubov eigenvalue problem and constitutes the main result of the present work.

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₀(x) = 1 and T₁(x) = x. Consider a Hermitian operator $\hat {X}$ with a discrete or continuous spectrum contained in I, and the correlation function

$$\begin{equation} f(\vert a \rangle,\vert b \rangle,x)=\langle a \vert\delta(\hat{X}-x)\vert b \rangle,\quad x\in I. \end{equation} \tag{ 20 }$$

The latter is a formal writing for

$$\begin{equation} f(\vert a \rangle,\vert b \rangle,x)= \sum_{j,\alpha}\delta(x_j-x) \langle a \vert x_{j,\alpha} \rangle\langle x_{j,\alpha} \vert b \rangle, \end{equation} \tag{ 21 }$$

where {|x_j,α〉} is an orthonormal eigenbasis of $\hat {X}$, which provides the spectral decomposition $\hat {X}=\sum _{j,\alpha } x_j \vert x_{j,\alpha } \rangle \langle x_{j,\alpha } \vert$, with x_j∈I, and the resolution of identity

$$\begin{equation} {\mathbbm{1}}=\sum_{j,\alpha} \vert x_{j,\alpha} \rangle\langle x_{j,\alpha} \vert. \end{equation} \tag{ 22 }$$

The above correlation function has the expansion

$$\begin{equation} f(\vert a \rangle,\vert b \rangle,x)=\frac{1}{\pi \sqrt{1-x^2}}\left[ \mu_0(\vert a \rangle,\vert b \rangle)+ 2 \sum_{n=1}^{\infty} \mu_n(\vert a \rangle,\vert b \rangle) T_n(x) \right], \end{equation} \tag{ 23 }$$

where the coefficient μ_n(|a〉,|b〉), called Chebyshev moment of order n, is defined as

$$\begin{equation} \mu_n(\vert a \rangle,\vert b \rangle)=\int_{-1}^1 f(\vert a \rangle,\vert b \rangle,x) T_n(x)\mathrm{ d}x. \end{equation} \tag{ 24 }$$

Owing to equations (21) and (22), the moments of f(|a〉,|b〉,x) are simply given by

$$\begin{equation} \mu_n(\vert a \rangle,\vert b \rangle) =\sum_{j,\alpha} T_n(x_j)\langle a \vert x_{j,\alpha} \rangle\langle x_{j,\alpha} \vert b \rangle =\langle a \vert T_n(\hat{X})\sum_{j,\alpha} \vert x_{j,\alpha} \rangle\langle x_{j,\alpha} \vert b \rangle, \end{equation} \tag{ 25 }$$

which boils down to the matrix element of a polynomial of $\hat {X}$:

$$\begin{equation} \mu_n(\vert a \rangle,\vert b \rangle)=\langle a \vert T_n(\hat{X}) \vert b \rangle. \end{equation} \tag{ 26 }$$

Instead of calculating $T_n(\hat {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, $\vert b_n \rangle =T_n(\hat {X})\vert b \rangle$ and $\vert b_{n-1} \rangle =T_{n-1}(\hat {X})\vert b \rangle$, with the initialization |b₀〉 = |b〉 and $\vert b_1 \rangle =\hat {X}\vert b \rangle$. Then, a single application of $\hat {X}$ (matrix–vector multiplication) yields the new vector $\vert b_{n+1} \rangle =2\hat {X}\vert b_n \rangle -\vert b_{n-1} \rangle$, 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

$$\begin{eqnarray} &&\fl f^{(N)}(\vert a \rangle,\vert b \rangle,x)=\frac{1}{\pi \sqrt{1-x^2}} \left[ g_0^{(N)}\mu_0(\vert a \rangle,\vert b \rangle) + 2 \sum_{n=1}^{N-1} g_n^{(N)}\mu_n(\vert a \rangle,\vert b \rangle) T_n(x) \right], \end{eqnarray} \tag{ 27 }$$

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').

Expression (19) for the one-body density matrix involves the eigenmodes of the Bogoliubov operator $\mathcal {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 $\mathcal {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 $\mathcal {L}\to \mathcal {L}/E_{\mathrm {max}}$.

To exhibit correlation functions akin to equations (20) and (21), we cast expression (19) into the following form:

$$\begin{equation} G(\boldsymbol{r},\boldsymbol{r'})=\sqrt{\rho(\boldsymbol{r})\rho(\boldsymbol{r'})}\exp\left[ -\frac{1}{2}\int_{0}^1 F(\boldsymbol{r},\boldsymbol{r'},\epsilon) \mathrm{d}\epsilon \right], \end{equation} \tag{ 28 }$$

with

$$\begin{eqnarray} F(\boldsymbol{r},\boldsymbol{r'},\epsilon)&= f(\boldsymbol{r},\boldsymbol{r},\epsilon)-f(\boldsymbol{r},\boldsymbol{r'},\epsilon)-f(\boldsymbol{r'},\boldsymbol{r},\epsilon)+f(\boldsymbol{r'},\boldsymbol{r'},\epsilon) \end{eqnarray} \tag{ 29 }$$

and

$$\begin{equation} f(\boldsymbol{r},\boldsymbol{r'},\epsilon)= -\sum_k \delta(\epsilon_k-\epsilon) \frac{\langle \boldsymbol{r} \vert w^\perp_{k,2} \rangle \langle \tilde{w}^{\perp}_{k,2} \vert \boldsymbol{r'} \rangle}{\sqrt{\rho_0(\boldsymbol{r})\rho_0(\boldsymbol{r'})}}, \end{equation} \tag{ 30 }$$

where |w^⊥_k,2〉 and $\vert \tilde {w}^{\perp }_{k,2} \rangle$ are the second components of the eigenvector (|w_k,1〉,|w_k,2〉)^T of $\mathcal {L}$ and its adjoint vector $(\vert \tilde {w}^{\perp }_{k,1} \rangle ,\vert \tilde {w}^{\perp }_{k,2} \rangle )^{\mathrm {T}}$, respectively. In equation (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 equation (28), the term f(r,r',) contributes to the exponent by

$$\begin{equation} \int_0^1 f(\boldsymbol{r},\boldsymbol{r'},\epsilon)\mathrm{ d}\epsilon= -\frac{1}{2 N_0} +\sum_j\frac{v^\perp_j(\boldsymbol{r}){v^\perp_j}^*(\boldsymbol{r'})}{\sqrt{\rho_0(\boldsymbol{r})\rho_0(\boldsymbol{r'})}}, \end{equation} \tag{ 31 }$$

where j enumerates the modes with E_j > 0. The −1/(2N₀) term, which stems from the zero eigenvector (0,|ϕ₀〉)^T, cancels out of the f sum in equations (28) and (29), so that equation (28) is indeed equivalent to equation (19).

The modes with E_k ⩽ 0 are included in sum (30) to use the resolution of identity (see equation (18))

$$\begin{equation} {\mathbbm{1}}=\sum_k (\vert w^\perp_{k,1} \rangle,\vert w^\perp_{k,2} \rangle)^{\mathrm{T}} (\langle \tilde{w}^{\perp}_{k,1} \vert,\langle \tilde{w}^{\perp}_{k,2} \vert). \end{equation} \tag{ 32 }$$

Indeed, rewriting equation (30) as

$$\begin{eqnarray} f(\boldsymbol{r},\boldsymbol{r'},\epsilon)&=-\sum_k \delta(\epsilon_k-\epsilon) \frac{ (0,\langle \boldsymbol{r} \vert) (\vert w^\perp_{k,1} \rangle,\vert w^\perp_{k,2} \rangle)^{\mathrm{T}} (\langle \tilde{w}^{\perp}_{k,1} \vert,\langle \tilde{w}^{\perp}_{k,2} \vert) (0,\vert \boldsymbol{r'} \rangle)^{\mathrm{T}}}{\sqrt{\rho_0(\boldsymbol{r})\rho_0(\boldsymbol{r'})}},\nonumber\\ \end{eqnarray} \tag{ 33 }$$

we find that the Chebyshev moments of f(r,r',) are

$$\begin{equation} \mu_n(\boldsymbol{r},\boldsymbol{r'})=- \frac{(0,\langle \boldsymbol{r} \vert)}{\sqrt{\rho_0(\boldsymbol{r})}} T_n\left(\frac{\mathcal{L}}{E_{\mathrm{max}}}\right) \frac{(0,\vert \boldsymbol{r'} \rangle)^{\mathrm{T}}}{\sqrt{\rho_0(\boldsymbol{r'})}}, \end{equation} \tag{ 34 }$$

and those of F(r,r',) follow as

$$\begin{equation} M_n(\boldsymbol{r},\boldsymbol{r'})=\mu_n(\boldsymbol{r},\boldsymbol{r})-\mu_n(\boldsymbol{r},\boldsymbol{r'})-\mu_n(\boldsymbol{r'},\boldsymbol{r})+\mu_n(\boldsymbol{r'},\boldsymbol{r'}). \end{equation} \tag{ 35 }$$

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

$$\begin{equation} \int_0^1 F(\boldsymbol{r},\boldsymbol{r'},\epsilon)\,\mathrm{ d}\epsilon= \frac{M_0(\boldsymbol{r},\boldsymbol{r'})}{2} +\frac{2}{\pi}\sum_{p=0}^{\infty} \frac{(-1)^p M_{2p+1}(\boldsymbol{r},\boldsymbol{r'})}{2p+1}. \end{equation} \tag{ 36 }$$

Note that the contributions of all even moments except M₀(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(\mathcal {L}/E_{\mathrm {max}})(0,\vert \boldsymbol {r'} \rangle )^{\mathrm {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

$$\begin{equation} \frac{G^{(N)}(\boldsymbol{r},\boldsymbol{r'})}{\sqrt{\rho(\boldsymbol{r})\rho(\boldsymbol{r'})}}= \exp\left[ -\frac{g_0^{(N)}}{4}M_0(\boldsymbol{r},\boldsymbol{r'}) -\sum_{p=0}^{\lfloor \frac{N}{2}-1\rfloor } \frac{(-1)^p g_{2p+1}^{(N)}}{(2p+1)\pi}M_{2p+1}(\boldsymbol{r},\boldsymbol{r'}) \right], \end{equation} \tag{ 37 }$$

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 equation (34) involve projections orthogonally to the ground state by means of $\hat {Q}= {\mathbbm{1}}-\vert \phi _0 \rangle \langle \phi _0 \vert$. Interestingly, the Bogoliubov operator $\mathcal {L}$ may be replaced by $\mathcal {L}_{\mathrm {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) − 〈ϕ₀|v_j〉ϕ₀(r), one easily sees that equation (19) also reads as

$$\begin{equation} G(\boldsymbol{r},\boldsymbol{r'})=\sqrt{\rho(\boldsymbol{r})\rho(\boldsymbol{r'})}\exp\left[-\frac{1}{2}\sum_j \left|\frac{v_j(\boldsymbol{r})}{\sqrt{\rho_0(\boldsymbol{r})}}-\frac{v_j(\boldsymbol{r'})}{\sqrt{\rho_0(\boldsymbol{r'})}} \right|^2\right]. \end{equation} \tag{ 38 }$$

The function F_GP(r,r',) is defined as the analogue of F(r,r',) in equations (28) and (29), with four terms of the form

$$\begin{equation} f_{\mathrm{GP}}(\boldsymbol{r},\boldsymbol{r'},\epsilon)=-\sum_{k'} \delta(\epsilon_{k'}-\epsilon) \frac{\langle \boldsymbol{r} \vert w_{k',2} \rangle \langle \tilde{w}_{k',2} \vert \boldsymbol{r'} \rangle}{\sqrt{\rho_0(\boldsymbol{r})\rho_0(\boldsymbol{r'})}}, \end{equation} \tag{ 39 }$$

where k' enumerates the eigenmodes of $\mathcal {L}_{\mathrm {GP}}$ with E_k' ≠ 0. Owing to the closure relation (13), the Chebyshev moments of f_GP(r,r',) read

$$\begin{equation} \mu_n^{\mathrm{GP}}(\boldsymbol{r},\boldsymbol{r'})= -\frac{(0,\langle \boldsymbol{r} \vert)}{\sqrt{\rho_0(\boldsymbol{r})}} T_n\left(\frac{\mathcal{L}_{\mathrm{GP}}}{E_{\mathrm{max}}}\right) \frac{(0,\vert \boldsymbol{r'} \rangle)^{\mathrm{T}}}{\sqrt{\rho_0(\boldsymbol{r'})}} +R_n^{(1)}(\boldsymbol{r},\boldsymbol{r'})+R_n^{(2)}(\boldsymbol{r},\boldsymbol{r'}), \end{equation} \tag{ 40 }$$

where

$$\begin{eqnarray} &&R_n^{(1)}(\boldsymbol{r},\boldsymbol{r'})= -\frac{\phi_a(\boldsymbol{r'})}{\sqrt{\rho_0(\boldsymbol{r})\rho_0(\boldsymbol{r'})}} (0,\langle \boldsymbol{r} \vert)\, T_n\left(\frac{\mathcal{L}_{\mathrm{GP}}}{E_{\mathrm{max}}}\right) \left(\begin{array}{@{}cc@{}}\vert \phi_0 \rangle\\ -\vert \phi_0 \rangle\end{array}\right), \end{eqnarray} \tag{ 41 }$$

$$\begin{eqnarray} &&R_n^{(2)}(\boldsymbol{r},\boldsymbol{r'})= +\frac{\phi_0(\boldsymbol{r'})}{\sqrt{\rho_0(\boldsymbol{r})\rho_0(\boldsymbol{r'})}} (0,\langle \boldsymbol{r} \vert)\, T_n\left(\frac{\mathcal{L}_{\mathrm{GP}}}{E_{\mathrm{max}}}\right) \left(\begin{array}{@{}cc@{}}\vert \phi_a \rangle\\ \vert \phi_a \rangle\end{array}\right). \end{eqnarray} \tag{ 42 }$$

Given that $\mathcal {L}_{\mathrm {GP}}(\vert \phi _0 \rangle ,-\vert \phi _0 \rangle )^{\mathrm {T}}=(0,0)^{\mathrm {T}}$ and $\mathcal {L}_{\mathrm {GP}}(\vert \phi _a \rangle ,\vert \phi _a \rangle )^{\mathrm {T}}=\alpha (\vert \phi _0 \rangle ,-\vert \phi _0 \rangle )^{\mathrm {T}}$ with constant α, as detailed in [77], we find

$$\begin{eqnarray} R_{2p}^{(1)}(\boldsymbol{r},\boldsymbol{r'})&=(-1)^{p}\frac{\phi_a(\boldsymbol{r'})}{\sqrt{N_0 \rho_0(\boldsymbol{r'})}}, \end{eqnarray}\noindent \tag{ 43 }$$

$$\begin{eqnarray} R_{2p+1}^{(1)}(\boldsymbol{r},\boldsymbol{r'})&=0, \end{eqnarray}\noindent \tag{ 44 }$$

$$\begin{eqnarray} R_{2p}^{(2)}(\boldsymbol{r},\boldsymbol{r'})&=(-1)^{p}\frac{\phi_a(\boldsymbol{r})}{\sqrt{N_0 \rho_0(\boldsymbol{r})}}, \end{eqnarray}\noindent \tag{ 45 }$$

$$\begin{eqnarray} R_{2p+1}^{(2)}(\boldsymbol{r},\boldsymbol{r'})&=(-1)^{p+1}\frac{(2p+1)\alpha}{N_0 E_{\mathrm{max}}}. \end{eqnarray} \tag{ 46 }$$

These R⁽ⁱ⁾_n(r,r') terms cancel out of the sum

$$\begin{equation} M_n^{\mathrm{GP}}(\boldsymbol{r},\boldsymbol{r'})= \mu_n^{\mathrm{GP}}(\boldsymbol{r},\boldsymbol{r})-\mu_n^{\mathrm{GP}}(\boldsymbol{r},\boldsymbol{r'})-\mu_n^{\mathrm{GP}}(\boldsymbol{r'},\boldsymbol{r})+\mu_n^{\mathrm{GP}}(\boldsymbol{r'},\boldsymbol{r'}). \end{equation} \tag{ 47 }$$

Hence, the comparison of equations (34) and (40) shows that $\mathcal {L}_{\mathrm {GP}}$ can be used instead of $\mathcal {L}$ in the Chebyshev iterations underlying equation (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 equation (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 operators $\mathcal {L}$ and $\mathcal {L}_{\mathrm {GP}}$ provide two other interesting examples in that context. Firstly, $\mathcal {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. Secondly, $\mathcal {L}_{\mathrm {GP}}$ is not diagonalizable, and yet its generalized eigenvectors (and their adjoints) can be used for the closure relation. In the derivation of equations (43)–(46), we took advantage of the fact that (|ϕ₀〉,−|ϕ₀〉)^T and (|ϕ_a〉,|ϕ_a〉)^T are generalized eigenvectors of low order, so that the action of $T_n(\mathcal {L}_{\mathrm {GP}})$ can be evaluated easily.

In [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 $\overline {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₁ is expected to decay at large distances with a power law given by [64, 86]

$$\begin{equation} g_1(r,r')\simeq\left(\frac{\mathrm{e}^{2-C}\xi}{4|r-r'|}\right)^{\frac{1}{2\pi\rho\xi}}, \quad |r-r'|\gg \xi, \end{equation} \tag{ 51 }$$

where $C=0.57721\ldots$ is Euler's constant, the density ρ can be approximated by ρ₀ within our Bogoliubov approach, and ξ is the healing length defined in equation (4). Figure 1 shows the result of a KPM calculation for a system of length L = 2¹⁶η, 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₁. Note the large system size achieved in this computation. In this homogeneous case, a single KPM iteration suffices for the computation of g₁(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 a few n_ℓ. This has to be compared to the storage space and computation time required for a full diagonalization of $\mathcal {L}_{\mathrm {GP}}$, which scale as n²_ℓ and n³_ℓ, respectively [65].

Zoom In Zoom Out Reset image size

For Δ > 0, the reduced one-body density matrix g₁(r,r') depends on the disorder configuration, and is no longer translation invariant. Figure 2 shows the behavior of g₁(0,r) for a single disorder configuration, and the results obtained with the KPM scheme for various moment numbers N (see equation (37)). When N is sufficiently high, the KPM results coincide with those obtained from complete diagonalization, and faithfully reproduce the spatial details of g₁(r,r'). As a general trend, we also observe that the KPM estimates of g₁(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₁, 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.

Zoom In Zoom Out Reset image size

After averaging over disorder, the one-body density matrix exhibits either a power-law or an exponential decay at 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

$$\begin{equation} g_1^L(r)=\frac{1}{L}\int_0^L g_1(r',r'+r)\,\mathrm{ d}r'. \end{equation} \tag{ 52 }$$

Figure 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.12E_c on the double logarithmic scale indicates a power-law decay of the averaged g₁, corresponding to the superfluid phase. For U = 0.3E_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 figure 3.

Zoom In Zoom Out Reset image size

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₁ behaves as a power law or an exponential. This procedure is put on a systematic footing by fitting the long-distance part of $\overline {g_1}(r)$ with both power laws and exponentials, and monitoring the fit quality. Figure 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 figure 4 are in good agreement with these scaling laws. Our findings thus reproduce the results of [52, 54] without the need for partial or complete diagonalization. This validates our approach.

Zoom In Zoom Out Reset image size

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₁(r,r') tends to a constant equal to the condensate fraction ρ₀/ρ for ∥r − r'∥ → ∞. To leading order in the strength of interactions, the condensate fraction is given by [59, 64, 66, 88]

$$\begin{equation} \frac{\rho_0}{\rho} \simeq1-\frac{g'}{8\pi} =1-\frac{1}{4\pi\rho_0 \xi^2}, \end{equation} \tag{ 53 }$$

where g' = 2mg/ℏ² 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 equation (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 [66] for the limit of weak interactions and weak disorder. To leading order in Δ/U, the depletion $\overline {\delta \rho }=\overline {\rho (r)}-\overline {\rho _0(r)}$ reads

$$\begin{equation} \overline{\delta\rho}= \delta\rho^{(0)}\left[ 1+\left(\frac{\Delta}{U}\right)^2\,h\left(\frac{\eta}{\xi}\right) \right], \end{equation} \tag{ 54 }$$

where δρ⁽⁰⁾ is the interaction-induced condensate depletion in a clean system, given by equation (53). The function h depends only on the ratio of the disorder correlation length η and the healing length ξ. Note that h differs from a similar function introduced in [66] by a trivial factor $\sqrt {2}$ in the argument due to a different definition of ξ.

The left panel of figure 5 shows the result of a KPM calculation of g₁(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₁ 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 figure 6, and stand in perfect agreement with equation (53).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The disordered case was examined with a similar procedure. The right panel of figure 5 displays the values of g₁(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 as 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 $\overline {g_1(\boldsymbol {r},\boldsymbol {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 GPE), 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 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₁(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₁(0,r) (i.e. the complement of the condensate fraction) as a function of the disorder strength. For weak disorder (Δ ≪ U), we indeed find 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 equation (54) and in $\overline {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 figure 7 thus suggest that higher orders in the weak-disorder expansion may reduce the condensate depletion.

Zoom In Zoom Out Reset image size