Dilute Bose gas with correlated disorder: A Path Integral Monte Carlo study

We investigate the thermodynamic properties of a dilute Bose gas in a correlated random potential using exact path integral Monte Carlo methods. The study is carried out in continuous space and disorder is produced in the simulations by a 3D speckle pattern with tunable intensity and correlation length. We calculate the shift of the superfluid transition temperature due to disorder and we highlight the role of quantum localization by comparing the critical chemical potential with the classical percolation threshold. The equation of state of the gas is determined in the regime of strong disorder, where superfluidity is suppressed and the normal phase exists down to very low temperatures. We find a $T^2$ dependence of the energy in agreement with the expected behavior in the Bose glass phase. We also discuss the major role played by the disorder correlation length and we make contact with a Hartree-Fock mean-field approach that holds valid if the correlation length is very large. The density profiles are analyzed as a function of temperature and interaction strength. Effects of localization and the depletion of the order parameter are emphasized in the comparison between local condensate and total density. At very low temperature we find that the energy and the particle distribution of the gas are very well described by the T=0 Gross-Pitaevskii theory even in the regime of very strong disorder.


Introduction
The dirty boson problem has become a central and fascinating subject in condensed matter physics starting from the first theoretical investigations more than 20 years ago [1,2,3]. The interplay between quantum degeneracy, interactions and quenched disorder in a bosonic system gives rise to a rich scenario that exhibits new and peculiar features compared to the much older problem of the metal/insulator transition with electrons [4,5,6]. An important difference is that bosons can not rely on the Pauli pressure and repulsive interactions are crucial to avoid collapse in the lowest localized single-particle state. As a result perturbation schemes starting from the noninteracting model, that are most useful for fermions, are completely inappropriate in the case of bosons.
Theoretical investigations, including quantum Monte Carlo simulations, have mainly addressed the problem of bosons on a lattice with on site bound disorder, the so-called disordered Bose-Hubbard model. In this case commensurability, i.e. the integer ratio of the number of particles to the number of lattice sites, plays a major role allowing for a superfluid/insulator (of the Mott type) transition also in the absence of disorder that is purely driven by interaction effects. Furthermore, depending on the value of the interaction strength, disorder can act in favor of superfluidity, by randomizing the insulating state close to the Mott transition, or in opposition to it by localizing almost free particles into single-particle levels. The insulating phases occurring in the two regimes of strong and weak interactions are respectively often referred to as the Bose glass, when interactions suppress superfluidity, and the Anderson glass, when interactions compete with disorder enhancing superfluidity [7,8]. The disorder driven quantum phase transition occurring at T = 0 has been investigated in a series of numerical studies both at incommensurate and commensurate densities and in various dimensionalities [7,9,10,11,12,13,14,8]. The picture emerging from these studies, together with the crucial role of interactions to stabilize the system, is that superfluidity is lost for strong enough disorder, leading to a gapless normal phase different from the incompressible Mott insulator.
Random potentials in continuous systems have been considered using perturbative approaches based on the Bogoliubov theory [15,16,17,18,19]. These methods are reliable when both interactions and disorder are weak and allow for the determination of the effect of disorder on the thermodynamic properties, including the fraction of atoms in the condensate, the superfluid density and other thermodynamic functions. Exact numerical methods have also been applied both at zero [21] and at finite temperature [22,23]. In particular, the path-integral Monte Carlo simulations carried out at finite T addressed the problem of the elementary excitations [22] and of the transition temperature [23] of a Bose fluid in a random environment. In the case of the continuous-space liquid phase, disorder always acts against superfluidity, whereas interaction helps make the superfluid state more robust. For strong disorder and low temperatures one expects the system to enter an insulating phase (Bose glass) that is smoothly connected with the high-temperature normal phase existing when the disorder is weak.
On the experimental side a large body of work was devoted to 4 He adsorbed in porous media, such as Vycor glass and aerogels. These studies investigated the behavior of the heat capacity and of the superfluid response [24,25,26], as well as the dynamic structure factor [27,28] as a function of temperature and filling. A suppression of the λ transition is observed and the critical coverage for the onset of superfluidity is determined as a function of temperature, however, no clear evidence is found of a compressible Bose glass phase. More recently the dirty boson problem has been addressed using ultracold atoms, which offer unprecedented control and tunability of the disorder parameters and of the interaction strength. Transport and phase-coherence properties of an interacting gas in disordered optical potentials are investigated and an insulating state is reached by increasing the strength of disorder [29,30,31,32,33,34,35]. A large experimental effort has also been devoted to the suppression of diffusion for non-interacting particles (Anderson localization) [36,37].
In this work we report on a path-integral Monte Carlo (PIMC) study of an interacting Bose gas in the presence of correlated disorder produced by 3D optical speckles. This random potential is relevant for experiments and allows for an independent tuning of intensity and correlation length. By increasing the disorder strength, we find a sizable reduction of the superfluid transition temperature and the shift is larger for weaker interactions. We map out the normal to superfluid phase diagram, both in the chemical potential vs. disorder and in the density vs. disorder plane. For strong disorder and in the presence of small but finite interactions, the critical chemical potential varies linearly with the disorder intensity and is essentially independent of temperature and interaction strength, in agreement with the existence of a mobility edge separating localized from extended states. We also find that the critical chemical potential is much larger than the classical percolation threshold for the same random potential, implying that a major role is played by quantum localization effects. In the regime of strong disorder and for chemical potentials below the critical value, the equilibrium state is a highly degenerate normal gas. We investigate the thermodynamic properties of this phase, finding a T 2 dependence of the equation of state in agreement with the peculiar feature expected for the Bose glass phase. The effect of the disorder correlation length is discussed in detail and we show that a non-trivial behavior is obtained only when the correlation length is comparable with the mean interparticle distance. At T = 0 we also carry out calculations using the Gross-Pitaevskii (GP) equation and at finite T using a selfconsistent mean-field approach based on Hartree-Fock theory and on the local density approximation. The results of the GP equation for the ground-state energy and the spatial distribution of particles are accurate even in the regime of strong disorder with short-range correlations. This conclusion might be useful in view of investigating the structural properties of the Bose glass phase. We consider a system of N identical bosons of mass m, subject to the random field V dis and interacting with a repulsive pairwise potential. The Hamiltonian of the gas takes then the form: Interatomic interactions are modeled using a hard-sphere potential: where the diameter a coincides with the s-wave length of the two-body scattering problem. Furthermore, the system is in a cubic box of volume Ω = L 3 with periodic boundary conditions. The structure of the paper is as follows. In section 2 we introduce the random potential and its statistical properties. In section 3 we discuss classical percolation for the speckle potential and we estimate the percolation threshold in 3D. Some details of the PIMC numerical method are presented in section 4. In section 5 we report our results on the superfluid transition: shift of the critical temperature, critical chemical potential and critical density. Most of these results were already presented in a previous publication of some of us [38]. In section 6 we introduce a mean-field approach based on the GP equation at T = 0 and on a Hartree-Fock self-consistent theory at finite T and for long-correlated disorder. The low temperature thermodynamics is discussed in section 7, including the equation of state, the condensate and total density profiles and the behavior of superfluid density and condensate fraction as a function of temperature and interaction strength in the disordered phase. Finally, we draw our conclusions in section 8.

Speckle potential
The random external field we consider is the one produced by 3D optical speckles. The local intensity is obtained from the following expression [39]: where V 0 is a positive constant and is the Fourier transform of the complex field ϕ(r), whose real and imaginary part are independent random variables sampled from a gaussian distribution with zero mean and unit variance. The function W (k) is a low-wavevector filter defined as  . (color online). Energy per particle of a classical non-interacting gas subject to the speckle potential. The homogeneous value E/N = 3k B T /2 is also shown as a reference.
The random potential in equation (3) is positive definite and the probability distribution of intensities is given by the normalized exponential law [40] If the system's size is large enough the above defined disordered potential is selfaveraging, i.e. spatial averages coincide with averages over different realizations. For a generic function f (V dis ) of the disorder intensity one can thus write the following identity According to this property, the spatial average and the corresponding root-meansquare displacement of the speckle potential are both determined by the energy scale The correlation length ℓ c of the random field is defined from the spatial autocorrelation function, as the length scale for which Γ(ℓ c /2) = Γ(0)/2.
The above equations characterizing the speckle intensity field in 3D can be straightforwardly generalized to 2D and 1D. In particular, in 1D the autocorrelation function Γ(x) takes the simple form [40] Γ and ℓ c = 0.88/Λ. In 3D we find that the autocorrelation function Γ(r) is well approximated by a gaussian (see figure 1) and we obtain numerically ℓ c = 1.1/Λ. It is important to notice that standard experimental realizations of optical speckles are 2D, i.e. the speckle pattern lies in the plane perpendicular to the propagation of the laser beam, featuring equal correlation lengths in the x and y planar directions and a much larger ℓ c in the z direction. We consider instead a 3D pattern, having the same correlation length in the three spatial directions. This random field can be realized, for example, by adding speckle patterns generated from different directions. The typical shape of the speckle potential V dis is shown in figure 2: typical wells have size ℓ c and depth V 0 . The energy 2 /mℓ 2 c , associated with the correlation length ℓ c , and V 0 provide the two relevant energy scales for the disorder potential. In particular, if V 0 ≫ 2 /mℓ 2 c the random potential is classical in nature, with typical wells that are deep enough to sustain many single-particle bound states. The opposite regime, V 0 ≪ 2 /mℓ 2 c , corresponds instead to quantum disorder, where typical wells of size ℓ c do not have bound states and these can be supported only by rare wells of size much larger than ℓ c or with depth much larger than V 0 .
The root-mean-square intensity V 0 and the correlation length ℓ c are the relevant parameters characterizing in general the various types of disorder. For example, the delta-correlated disorder, which has been considered in many theoretical investigations [15,18,19,20], is defined by the following autocorrelation function: where ∆V dis (r) = V dis (r) − V dis is the displacement from the average. By approximating the speckle Γ function (8) using a gaussian function, Γ(r) = V 2 0 e −r 2 /2σ 2 , with σ = ℓ c / √ 8 log 2 to recover the same half width at half maximum, one finds that the speckle field in the limit ℓ c → 0 reproduces a delta-correlated disorder with the strength κ given by In our simulations the length scale ℓ c is typically ∼ 100 times larger than the hardsphere diameter a, allowing for a wide range of disorder intensities where interaction effects are well described by the s-wave scattering length and the details of the interatomic potential are irrelevant. The typical box size used in the simulations ranges from L ∼ 20ℓ c to L ∼ 50ℓ c . An indication of self-averaging of disorder for these values of L is provided by figure 1, where we show the comparison between the autocorrelation function Γ averaged over many realizations of the random potential and the one corresponding to a single realization. The self-averaging property (7) allows one to calculate the thermodynamics of a classical non-interacting gas. For example, the average energy per particle obtained from the spatial average of the disordered potential over the Boltzmann factor yields the following simple result In figure 3 we compare the above analytical prediction with the results obtained from a direct spatial integration using a typical size of the simulation box. The good agreement found shows that in our simulations the self-averaging property is well satisfied for non-trivial functions of the disorder intensity.

Classical percolation
In this Section we investigate the problem of the conducting/insulating transition in a speckle potential from the point of view of classical percolation. The relevant question is to determine the mobility edge of a classical particle subject to the random field [41]. Given the disordered potential V dis (r), the fraction of space accessible to particles of energy ǫ is defined as the fractional volume where the reference energy exceeds the external field: The percolation threshold corresponds to the critical value Φ c of the fractional volume such that, if Φ(ǫ) > Φ c , there are infinitely extended volumes of allowed space and particles having energy ǫ can move across the whole system. In terms of energy, the value ǫ c determines the percolation threshold: Φ(ǫ c ) = Φ c . It corresponds to the classical mobility edge separating localized states with energy ǫ < ǫ c from delocalized ones with energy ǫ > ǫ c . In the case of speckles the function Φ(ǫ) can be simply expressed in terms of the disorder intensity V 0 using the property (7). One finds The determination of the percolation threshold for lattice or continuum models is in general a difficult numerical task and the study of percolation has become a mature branch of statistical physics (see e.g. [42]). Our discussion here is limited to the calculation of Φ c and of the corresponding thereshold energy ǫ c for the speckle potential. The role of dimensionality is crucial for this problem. In 1D the unbound nature of the disordered field implies that, approaching the thermodynamic limit, potential barriers occur higher than any finite energy ǫ and consequently ǫ c = +∞ and Φ c = 1. In 2D the percolation threshold of laser speckles was investigated experimentally [43] using photolithography on a conducting film obtaining the value Φ c = 0.407. This result was later confirmed by a numerical study [44]. To our knowledge there are no determination of Φ c for 3D speckle patterns.
We estimate the percolation threshold in the continuum by mapping the accessible and unaccessible regions of space on a finite grid. This is done by simply comparing the local value of the external field V dis (r i ) at the grid point r i with the reference energy ǫ. One then investigates the percolation of accessible grid points of the corresponding matrix. In figure 4 we show two typical configurations in the case of 2D speckles: panel a) corresponds to Φ = 0.30 well below the percolation threshold and panel b) to Φ = 0.50 where percolation, i.e. existence of an uninterrupted path of accessible points across the whole system in at least one of the spatial directions, has clearly occurred. Two issues have to be considered with care: i) the size L in units of the correlation length ℓ c must be increased in order to extrapolate to the thermodynamic limit and ii) for a given size L the grid points must be dense enough. In figure 5 we plot the fraction of configurations exhibiting percolation in 3D speckle patterns as a function of the accessible volume Φ and for different system sizes. An estimate of the threshold gives Φ c ≃ ǫ c /V 0 = 4(1) · 10 −4 . For comparison, we estimated the percolation threshold in 2D obtaining the value Φ c ≃ 0.4 in agreement with previous determinations [43,44]. The value found for Φ c in 3D is rather small, suggesting the presence of many long valleys where particles with a tiny fraction of the energy V 0 can freely move across the whole system. One should notice that a similar 3D continuum model, the percolation of voids between overlapping spheres (the so called "Swiss cheese" model), gives the much larger value Φ c ≃ 0.03 [45,46] for the percolation threshold.

PIMC method
The partition function Z of a bosonic system with inverse temperature β = (k B T ) −1 is defined as the trace over all states of the density matrixρ = e −βĤ properly symmetrized. The partition function satisfies the convolution equation where τ = β/M , R collectively denotes the position vectors R = (r 1 , r 2 , ..., r N ), P R denotes the position vectors with permuted labels P R = (r P (1) , r P (2) , ..., r P (N ) ) and the sum extends over the N ! permutations of N particles. The calculation of the partition function in equation (16) can be mapped to a classical-like simulation of polymeric chains with a number of beads M equal to the number of terms of the convolution integral. In a PIMC calculation, one makes use of suitable approximations for the density matrix ρ(R, R ′ , τ ) at the higher temperature 1/τ in equation (16) and performs the multidimensional integration over R, R 2 ,...,R M as well as the sum over permutations P by Monte Carlo sampling [51]. The statistical expectation value of a given operator O(R), is calculated by generating stochastically a set of configurations {R i }, sampled from a probability density proportional to the symmetrized density matrix, and then by averaging over the set of values {O(R i )}. An approximation for the high temperature density matrix, which is particularly well suited for studies of dilute gases, is based on the pair-product ansatz [51] In the above equation ρ 1 is the single-particle ideal-gas density matrix and ρ rel is the two-body density matrix of the interacting system, which depends on the relative coordinates r ij = r i − r j and r ′ ij = r ′ i − r ′ j , divided by the corresponding ideal-gas term The two-body density matrix at the inverse temperature τ , ρ rel (r, r ′ , τ ), can be calculated for a given spherical potential V (r) using the partial-wave decomposition where P ℓ (x) is the Legendre polynomial of order ℓ and θ is the angle between r and r ′ . The functions R k,ℓ (r) are solutions of the radial Schrödinger equation with the asymptotic behavior holding for distances r ≫ R 0 , where R 0 is the range of the potential. The phase shift δ ℓ of the partial wave of order ℓ is determined from the solution of equation (22) for the given interatomic potential V (r). For the hard-sphere potential (2) a simple analytical approximation of the hightemperature two-body density matrix due to Cao and Berne [52] has been proven to be highly accurate [53]. The Cao-Berne density matrix ρ CB rel (r, r ′ , τ ) is obtained using the large momentum expansion of the hard-sphere phase shift δ ℓ ≃ −ka + ℓπ/2 and the large momentum expansion of the solutions of the Schödinger equation (22) R k,ℓ (r) ≃ 2/π sin[k(r − a)]/r. This yields the result Simulations are based on the worm algorithm [54], which allows for an efficient sampling of permutation cycles. In this scheme one samples both diagonal configurations, contributing to averages of the type (17) where all paths are closed, and off-diagonal configurations where one path is open. These latter configurations contribute to the one-body density matrix (OBDM) defined as where r P (1) = r ′ 1 . The long-range behavior of the OBDM determines the condensate density In our simulations the largest displacement of the OBDM we consider is |r−r ′ | = L/2. If the size L is large enough the number N 0 of condensate particles can be written as where r ′ is fixed by the constraint |r − r ′ | = L/2 and we perform an average over the solid angle. The quantity under the integral corresponds to the local condensate density at position r, which could be highly non uniform in the presence of a random potential.
Beside the condensate density n 0 , in the present study we consider also the superfluid density ρ s . The superfluid component is the part of the fluid that remains at rest when an infinitely slow movement is applied to the walls that contain the system. In the path-integral formalism, the superfluid fraction of a fluid contained in a box with periodic boundary conditions can be related [51] to the fluctuations of the winding number via the equation The winding number W is defined as: It is a topological property of the configuration. It counts the net number of paths that cross any plane perpendicular to one axis. The worm algorithm is particularly suitable to perform ergodic random walks that span all possible winding numbers since it extends the configurations space by including configurations with an open path. Only the Monte Carlo moves that modify the open path can efficiently change the winding number. We perform calculations both in the canonical (at fixed density n) and in the grand-canonical ensemble (at fixed chemical potential µ) [54]. We supplement the worm algorithm with two additional Monte Carlo updates that change the particle number N . The first update adds one particle to the system by placing a closed path at a randomly selected position. The second update erases a randomly selected closed path. The acceptance probability of the first (second) update is fixed by the fugacity e βµ (by its inverse), by the change in the interaction energy due to the path to be inserted (erased) and by the factor ΩC N +1 ( N ΩC ) that takes into account the density change and the normalization of the free particle propagator C ≡ 2π 2 β/m − 3 2 .

Superfluid transition
The effect of disorder is to suppress both the superfluid and the condensate density. In figure 6 we illustrate the behavior of these two quantities when the disorder strength is increased. The value of ρ s /ρ and n 0 /n in the absence of disorder is determined by the temperature T and by the value of the gas parameter na 3 . The figure shows a linear decrease of the superfluid and condensate components with increasing disorder in the classical regime V 0 > 2 /mℓ 2 c . An important question is to understand whether disorder also affects the superfluid transition temperature. The results are shown in figure 7. The transition temperature T c is expressed in units of the critical temperature of the non-interacting gas with ζ(3/2) ≃ 2.612. In these calculations the value of the scattering length and of the disorder correlation length are kept fixed and for the latter we choose the value ℓ c = 0.6n −1/3 , such that there is typically one particle in each sphere of radius ℓ c : n4πℓ 3 c /3 ≃ 1. We show results corresponding to two different densities. The reported values of T c in the absence of disorder are taken from Ref. [55]. At the density corresponding to the gas parameter na 3 = 10 −4 we find no appreciable change of the transition temperature compared to clean systems by increasing the disorder strength up to V 0 ∼ 2 /mℓ 2 c . For larger intensities we find a sizable shift that is well described again by a linear dependence in V 0 . For a given strength V 0 the reduction of the transition temperature is enhanced for smaller values of the gas parameter, consistently with the instability of the ideal Bose gas in the presence of disorder [16]. The value of T c is extracted from the results  of the superfluid fraction ρ s /ρ (ρ = mn is the total mass density), corresponding to systems with different particle number N , using the scaling ansatz Here, t = (T − T c )/T c is the reduced temperature, ν is the critical exponent of the correlation length ξ(t) ∼ t −ν , and f (x) is a universal analytic function, which allows for a linear expansion around x = 0. The validity of the scaling behavior (31) is shown in figure 8, where the effect of different realizations of the random potential is also shown. The quantity N (1+η)/3 n 0 /n, involving the condensate fraction n 0 /n and the correlation function critical exponent η = 0.038 of the XY-model universality class, is also expected to obey a scaling relation of the form (31). For all reported disorder µ/µ C = 1.05 µ/µ C = 0.975 strengths V 0 , the extracted value of the critical exponent ν is compatible with the result ν = 0.67 corresponding to clean systems [55]. It is worth noting that the values of T c , obtained from the scaling law of the superfluid ρ s /ρ and of the condensate fraction n 0 /n, coincide within our statistical uncertainty (see figure 7). The presence of disorder reduces the quantum delocalization of particles occupying the deepest wells of the potential and, consequently, their contribution to the superfluid behavior. Superfluidity takes place when the degeneracy condition is met for the effectively smaller density of "delocalized" particles, resulting in a suppressed value of T c . In Ref. [18] the shift δT c = T c − T 0 c of the superfluid transition temperature is calculated using a perturbative approach for the δ-correlated disorder ∆V dis (r)∆V dis (r ′ ) = κδ(r − r ′ ), where ∆V dis (r) = V dis (r) − V dis . The T c shift is found to be quadratic in κ, implying for our speckle potential that δT c /T 0 , where we used a gaussian fit to the radial dependence of the autocorrelation function Γ and considered the limit ℓ c → 0. We report this prediction in figure 7 (we also add the interaction contribution not accounted for by Ref. [18], so that in the clean case an exact result is reproduced). Our data in the regime of very weak disorder do not have enough precision to allow for a quantitative comparison and diverge from the theory before δT c /T 0 c becomes appreciable. The effect of disorder on the critical temperature of a hard-sphere gas was also investigated using PIMC methods in Ref. [23] where, however, no significant reduction of T c was reported. For stronger intensities of disorder, the calculation of T c becomes increasingly difficult, since the dependence on the realization gets more important (see figure 8) and larger systems are needed in order to have a satisfactory self-averaging of the random potential.
In figure 9 we report results for the critical chemical potential µ c obtained from calculations carried out in the grand-canonical ensemble. A small change of µ around µ c translates into a drastic change in the long-range behavior of the OBDM (see figure 10): for µ < µ c the OBDM decays to zero and corresponds to a normal phase, for µ > µ c the OBDM reaches a constant value characteristic of the superfluid state. If interactions are small but finite, we also find that the value of µ c is essentially insensitive to a change of temperature and of interaction strength. For weak disorder, this result is accompanied by a very small critical density (see figure 11) and corresponds to a renormalization of µ c due to disorder in an extremely dilute gas. For strong disorder, it is instead consistent with the picture of a mobility edge, separating localized single-particle states from extended ones, that depends only on the parameters of the random potential. In this latter regime we find a linear dependence of µ c as a function of V 0 , in agreement with the qualitative T = 0 prediction of Refs. [57,56] in the case of classical disorder. The figure also shows the classical percolation threshold µ = ǫ c , whose value for the speckle potential has been determined in section 3. One should notice that in the whole range of disorder intensities the critical chemical potential is significantly larger than ǫ c as a consequence of quantum localization effects. In fact, in terms of a mobility edge, classical particles with energy larger than ǫ c can freely move across the entire system, while in the quantum world extended states appear only for significantly larger energies bound by the inequality ǫ > µ c .
To conclude the study of the critical behavior, we analyze the dependence of the critical density n c on the intensity of the random potential. The calculations are carried out in the canonical ensemble at fixed temperature and scattering length. The method used to determine n c is shown in the inset of figure 11. For a given value of V 0 one increases the density and calculates the superfluid ρ s /ρ and the condensate fraction n 0 /n. The results are then fitted by a power-law dependence ρ s /ρ ∼ (n− n c ) ν and n 0 /n ∼ (n − n c ) ν(1+η) for n > n c , where the proportionality coefficients are expected to be non-universal parameters. In the inset of figure 11 we show the results corresponding to a configuration without disorder (V 0 = 0) and with strong disorder (V 0 = 6.4 2 /mℓ 2 c ). The reported values are averaged over a few realizations of the random potential and their scatter gives an idea of the relevance of this effect. For the small value of the scattering length used here, the critical density at V 0 = 0 coincides with the non-interacting result n 0 c = ζ(3/2)(mk B T /2π 2 ) 3/2 , while for the large V 0 one finds that n c is about a factor of eight greater than n 0 c . More comprehensive results are shown in figure 11 where n c is estimated from the superfluid fraction, which is less sensitive to finite-size effects. The results clearly show an increase of the critical density as a function of V 0 , from the non-interacting degenerate density n 0 c up to values ∼ 20 times larger. It is also worth noticing that for strong disorder, an increase of the scattering length a is accompanied by a decrease of n c resulting in a constant value of the critical chemical potential (see figure 9).

Mean-field approach
A simple description of the thermodynamic properties of disordered systems can be provided in terms of a mean-field approach. At T = 0 this approach is based on the solution of the Gross-Pitaevskii (GP) equation for the order parameter in the random external field and it yields quantitatively reliable results for both the chemical potential and the density profiles. At finite temperature the mean-field theory can be efficiently applied in the case of random potentials with exceedingly long-range correlations where the local density approximation holds valid.

Zero temperature
The relevant equation is the stationary GP equation for the order parameter ψ(r) in the presence of the random potential where g = 4π 2 a/m is the coupling constant. Within the GP approach one does not distinguish between order parameter and particle density and the following identity holds valid: n(r) = |ψ(r)| 2 . The GP equation can be obtained from the following energy functional [47] using the variational ansatz that corresponds to finding configurations that minimize the energy functional (34) with the normalization constraint d 3 r|ψ(r)| 2 = N . The solutions of the GP equation are obtained numerically by discretizing the wave function ψ on a 3D box of size L (the number of grid points ranges from 64 3 to 128 3 ). Then, the energy functional E[ψ] is minimized by using a conjugate gradient algorithm as described in [48,49], for a given realization of the potential and a given particle density n = N/L 3 . The numerical solution yields the value of the chemical potential µ as well as the spatial particle distribution. Averages over disorder are obtained by repeating the calculation for various realizations of the random field.
It is useful to rewrite the GP equation (32) in terms of the dimensionless variables v dis (r) = V dis (r)/V 0 ,r = r/ℓ c ,Ṽ 0 = V 0 /( 2 /mℓ 2 c ) andμ = µ/( 2 /mℓ 2 c ). One finds where we have used the wavefunction φ(r) = ψ(r)/ √ n rescaled in terms of the average density n and the healing length ξ = 1/ √ 8πna. From the above equation two regimes can be investigated analytically. On the one hand, if the correlation length ℓ c is much smaller than the length 2 /mV 0 fixed by the disorder strength or equivalently for weak enough disorder (Ṽ 0 ≪ 1), the order parameter is almost uniform |φ(r)| 2 ≃ 1 and the effect of the random potential on the chemical potential is just a shift of the mean-field energy On the other hand, if ℓ c ≫ ξ, one can neglect the kinetic energy term in the GP equation (35) and use the Thomas-Fermi approximation µ = gn(r) yielding the following result for the local particle density Here Θ(x) is the Heaviside function: Θ(x) = 1 if x > 0 and zero otherwise. The normalization condition n = 1/L 3 drn(r) determines the chemical potential µ in terms of the average density n. By using the self-averaging property (7) one obtains the equation relating the chemical potential to the disorder strength V 0 . i) If V 0 ≪ gn then µ = gn + V 0 and the disorder acts as a small shift of the pure interaction term, similarly to the short-ℓ c regime of equation (36).
ii) If V 0 ≫ gn one finds to the lowest order corresponding to an energy per particle E/N = 2µ/3.
The variation of µ with the disorder correlation length is shown in figure 12, where we report the results of the GP equation for the fixed value na 3 = 10 −6 of the gas parameter and disorder strength V 0 = 10k B T 0 c . In this case ξn 1/3 ∼ 1 and the figure clearly shows the two limiting regimes of short and long correlation length. For the same disorder strength, in figure 13, we show instead the density dependence of the chemical potential for a fixed value of ℓ c . The equations of state corresponding to large and small correlation lengths are also shown as a reference.
In the regime ℓ c ≫ ξ one can make use of the Thomas-Fermi approximation for the order parameter which is expected to become more and more accurate as ℓ c increases. Within this approximation one can make contact with the classical percolation problem of section 3 by noticing that, if µ is larger than the threshold energy ǫ c , the condensate density represented by equation (38) is different from zero on a percolating path ensuring thus the superfluid behavior of the system. For large V 0 the chemical potential increases as √ V 0 according to equation (40), while ǫ c is proportional to V 0 . The small value of the ratio ǫ c /V 0 implies though that the quantum phase transition to the insulating state takes place only at extremely large disorder intensities.

Finite temperature
At T = 0 one should combine the GP equation for the condensate with a proper description of the thermally excited states in the random potential. The theory becomes simple if the disorder is correlated over large distances, as one can apply the local density approximation (LDA) within standard mean-field techniques suitable for dilute gases. At low temperatures the validity condition requires ℓ c to be much larger than the healing length ξ, while at higher temperatures the correlation length must exceed the thermal wavelength λ T .
We use the self-consistent Hartree-Fock (HF) scheme within LDA. This meanfield approach is based on the following expression for the elementary excitations of the system in terms of their momentum p and position r [50] ǫ(p, r) = p 2 2m + V dis (r) − µ + 2gn(r) .
The thermal density of atoms n T (r) is obtained from the momentum integral of the distribution of elementary excitations The condensate density n 0 (r) is determined by the finite-T extension of the GP equation (38) in the Thomas-Fermi approximation The sum of thermal and condensate density gives the total density n(r), which must satisfy the overall normalization providing the closure relation of the mean-field equations. The total energy E of the system can be calculated from the following integral E = dp dr (2π ) 3 p 2 /2m e ǫ(p,r)/kBT − 1 + drV dis (r)n(r) The HF scheme defined by the above equations neglects the quantum depletion of the condensate and at T = 0 yields n 0 (r) = n(r), in agreement with section 6.1. Furthermore, it neglects the contribution of collective modes (phonons) to thermodynamics since all excitations are single particle. This approximation is known to be very accurate in dilute non-uniform systems both at high and low temperatures [50]. In particular, at low temperatures, one can neglect the thermal density contribution to the expression (41) of the elementary excitations and one finds the simple spectrum ǫ(p, r) = p 2 /2m + |V dis (r) − µ(T = 0)|. The term in modulus vanishes at the condensate boundaries and thermal excitations accumulate at these minima of the effective potential. The single-particle excitations close to these minima are the dominant ones at low temperature, being more important than the phonons of the bulk condensate. In fact, one can show that at low energy the single-particle excitations have the following density of states g(ǫ) = 4 3 proportional to ǫ 3/2 in contrast to g(ǫ) ∝ ǫ 2 of phononic excitations. A similar situation occurs in harmonically trapped condensates [50]. The above semiclassical approach provides an estimate of the temperature at which Bose-Einstein condensation sets in locally in some deep well. This temperature is defined as the point where the local density n(r), corresponding to some deep well in the random field, reaches the critical value n(r)λ 3 T = ζ(3/2) ≃ 2.612. By neglecting interactions one finds the following implicit equation for the temperature T * n 2π 2 The temperature T * is always larger than the temperature T 0 c of the occurrence of Bose-Einstein in non-interacting clean systems. In particular, for large disorder strength one finds . This effect comes from the reduced available volume and the corresponding higher local particle density. In the presence of weak interactions, local Bose-Einstein condensation sets in at a temperature slightly lower that T * , because density is reduced in the wells of the random field due to repulsion and a lower temperature is needed to reach the critical value. We would like to stress that the temperature scale T * corresponds to the appearance of local condensates at the minima of V dis (r) and should not be confused with the critical temperature T c at which extended superfluidity sets in. Within the above semiclassical approach this latter temperature corresponds to the chemical potential reaching the percolation threshold of the effective potential V dis (r)+ 2gn T (r) where, according to equation (43), the condensate density n 0 (r) is different from zero on a percolating path. For weakly-interacting systems though, because of the small value of the percolating volume fraction Φ c , the temperatures T c and T * are very close, unless for extremely large disorder intensities. As an example we consider the configuration shown in figure 14 corresponding to na 3 = 10 −6 and V 0 = 10k B T 0 c in the regime of extremely long-range correlation length ℓ c . The value of the temperature T * obtained from equation (47) is given by T * = 3.63T 0 c . The self-consistent solution of the HF equations yields a temperature T BEC for the local onset of Bose-Einstein condensation in the range 3.5T 0 c < T BEC < T * . The transition temperature T c where the condensate density percolates is found to be in the range 3.4T 0 c < T c < T BEC . The HF equation of state and the corresponding transition temperature are shown in figure 14 for a fixed density of the gas and for large disorder strength. In the figure is also reported the equation of state corresponding to the regime of a very short correlation length ℓ c , where the energy is merely shifted by the average disorder intensity V 0 from the value of the clean system. This result is consistent with the T = 0 prediction (32) and the corresponding transition temperature coincides with T c in the absence of disorder.

Low temperature thermodynamics
The energy per particle obtained from PIMC simulations is reported in figure 14 as a function of temperature and for random potentials with different correlation lengths in the range of the mean interparticle distance. In particular, for the value ℓ c = 0.6n −1/3 we also indicate the value of the superfluid transition temperature. The equation of state and the value of T c corresponding to the limiting regimes of exceedingly long and short correlation lengths are also shown as a reference. Three important remarks about this figure are worth stressing. a) At fixed particle density and for a given    disorder strength the largest suppression of T c is achieved with ℓ c on the order of the interparticle distance. b) At low temperature the energetics of the system is well described by the GP equation even when the disorder is strong and n 1/3 ℓ c ∼ 1. c) For random potentials of large strength and with spatial correlations in the range n 1/3 ℓ c ∼ 1, a window of temperatures opens up where the system is normal (i.e. not superfluid), though highly degenerate (nλ 3 T > 1). The properties of the "exotic" normal phase displayed at low temperatures for strong disorder are worth investigating. We find that the equation of state is well fitted by a quadratic temperature dependence as shown in figure 14. This T 2 -law for the energy, and consequently linear dependence of the specific heat, is a remarkable feature of the phase. In a clean superfluid, or in a dirty one with short-range disorder correlations (see figure 14), the energy at low temperatures exhibits the T 4 -law typical of the thermal excitation of phonons. In the opposite regime of long-range disorder correlations one can apply the HF-LDA approach described in section 6.2, which is expected to become more and more accurate as ℓ c increases. Within this approach the system consists of large condensate lakes that may or may not be connected through a percolating path, the relevant excitations at low temperature are of single-particle nature and are localized at the boundary of the condensate lakes. These excitations contribute to the energy with a T 7/2 -law ‡. A Bose glass is predicted to exhibit a non-vanishing density of states at zero energy [3], which results in a T 2 -law for the low temperature equation of state. We interpret the quadratic dependence found in the low temperature normal phase as an evidence of the Bose glass phase.
Another important output of PIMC simulations with the random potential is the spatial particle distribution and the distribution of particles contributing to the condensate density. In figure 15 we show the results of the particle density n(r) for a given realization of the random potential. Starting from a high temperature distribution spread over the entire system, as temperature is decreased, the density becomes more and more peaked in correspondence of the minima of the external potential. Around T ≃ 0.4T 0 c the system turns superfluid (see figure 19) and the density changes only slightly down to the lowest temperatures. Remarkably, the comparison with the findings of the GP equation for the same realization of the random field is rather good both for the position and the relative intensity of the peaks.
The comparison between particle density and condensate density profiles is reported in figure 16 for two temperatures, above and below T c . The result at the lowest temperature shows that the condensate density follows the particle distribution, but does not exhibit its pronounced peaks. This behavior is clearly represented in figure 17, where we show the profiles along a cut through the plane of figure 16. Notice that the particles contributing to the condensate are only a small fraction (∼ 20%) of the total number of particles (see figure 19). Finally, in figure 18 we report the density profiles as a function of the interaction strength. The figure clearly shows the effect of particle localization in the deepest wells of the random field as the value of the gas parameter is reduced.
The behavior of the superfluid and condensate fraction in the proximity of the phase transition with and without disorder are shown in figure 19. It is interesting to notice the large depletion of ρ s and n 0 even at the lowest temperatures and the fact that, in the regime of strong disorder, the superfluid component becomes significantly ‡ In the insulating Bose glass phase the asymptotic low-temperature law is expected to be ∼ T 2 independent of the value of the correlation length. However, for increasing ℓc, the range of temperatures where the asymptotic law applies is suppressed.
smaller than the condensate one. Similar results are reported in figure 20 for a fixed temperature and for varying values of the interaction strength. The figure clearly shows that by making interactions weaker the system will eventually turn normal and that the difference between superfluid and condensate fraction disappears with strong enough interactions.
To complete the picture we would like to mention the limit of weak disorder (V 0 mℓ 2 c / 2 ) ≪ 1 and weak interactions na 3 → 0 considered in Ref. [56]. It explains how interactions stabilize the superfluid phase starting from the localized phase of the non-interacting gas. At zero temperature the boundary between the insulating (Bose glass) and superfluid phases is set (up to logarithmic corrections) by percolation at the Larkin energy scale V 4 0 m 3 ℓ 6 c / 6 which has to be compared with the chemical potential of the weakly interacting gas gn. This leads to the relation for the critical density (na) c ∼ V 4 0 , i.e. for small V 0 only a tiny density of particles is required to off-set the localization effects. In practice, one observes extremely robust superfluidity of the weakly-interacting Bose gas even in the presence of relatively large disorder, see Figs. 7,19,20; for na ≫ (na) c the critical temperature is essentially the same as in the ideal Bose gas.

Conclusions
We find that in a quantum degenerate bosonic gas a random potential is most efficient in suppressing superfluidity if it is correlated over length scales comparable with the mean interparticle distance. However, for the typical diluteness conditions of ultracold gases, disorder intensities significantly larger than the energy scale k B T 0 c set by the BEC transition of the ideal gas are required for a significant reduction of the superfluid critical temperature and stronger interactions make the superfluid state more robust. In the regime of weak interactions and strong disorder, the superfluid transition turns out to be well characterized by the existence of a mobility edge, separating localized from extended states, that is largely independent of temperature and interaction strength. This picture is similar to the percolation threshold of classical particles, but we find that quantum localization effects drive the system normal in a large region of energies where classical states would be extended. Furthermore, most of the particles are localized in the deepest wells of the random potential and only a small proportion contributes to the extended condensate state. The effective density of these delocalized particles is much smoother due to the screening of the external field from the other particles and sets the critical density of superfluidity which however starts at significantly lower temperatures compared to clean systems. For larger disorder intensities an "exotic" normal phase appears in the degenerate regime, even though we can not reach T = 0. This phase is characterized by a peculiar T 2 dependence of the equation of state, that is markedly different from the T 4 law of homogeneous superfluids and from the T 7/2 law found for large condensate lakes within LDA and is in agreement with the predictions for the Bose glass phase. Remarkably, some aspects of this phase, such as the T = 0 equation of state and the spatial distribution of particles can be correctly described using the mean-field GP theory.
FerMix", was supported by funds from the CNR and the EC Sixth Framework Programme. NP acknowledges support from NSF grant PHY-0653183. Calculations have been performed on the HPC facility Wiglaf at the Physics Department of the University of Trento and on the BEN cluster at ECT * in Trento.