Glass to superfluid transition in dirty bosons on a lattice

We investigate the interplay between disorder and interactions in a Bose gas on a lattice in presence of randomly localized impurities. We compare the performance of two theoretical methods, the simple version of multi-orbital Hartree-Fock and the common Gross-Pitaevskii approach, showing how the former gives a very good approximation to the ground state in the limit of weak interactions, where the superfluid fraction is small. We further prove rigorously that for this class of disorder the fractal dimension of the ground state d* tends to the physical dimension in the thermodynamic limit. This allows us to introduce a quantity, the fractional occupation, which gives insightful information on the crossover from a Lifshits to a Bose glass. Finally, we compare temperature and interaction effects, highlighting similarities and intrinsic differences.


Introduction
The physics of disordered quantum systems is a very active field of research since the 1950s, when this topic received the fundamental contributions of Anderson and Mott [1]. Their works showed that, in presence of disorder, even ideal conductors undergo a phase transition towards an insulating state due to destructive quantum interference. An ideal playground for the investigation of quantum effects is offered by ultracold atoms, where a controllable amount of disorder may be implemented in many ways, e.g., by means of speckle potentials, bichromatic lattices with incommensurate frequencies, localized impurities, or site-resolved addressing in optical lattices. After a long search, Anderson localization of ultracold ideal gases was finally observed in one dimension (1D) [2], and very recently reported also in 3D [3].
At this point it is important to stress the fundamental difference between disordered fermionic and bosonic systems. For fermions the crucial role is played by the Fermi energy: the question whether a given system is a conductor or insulator depends on whether the single particle states close to Fermi level are localized or not. Bosons, in contrast, tend to condense in the lowest energy state (or states). The question whether they constitute a superfluid or an insulator reduces then to understanding if the lowenergy states are localized or not (for more extensive discussion see for instance [4,5]).  Figure 1. Schematic phase diagram for dirty bosons at zero temperature as a function of interactions (U ) and chemical potential (µ) containing glassy (G), superfluid (BEC), and Mott insulating (MI) phases. The continuous and dashed lines mark respectively the superfluid-insulator and glass-Mott transition, while the dotted line indicates a crossover from the Lifshits glass (heavily fragmented) to the Bose glass (extended) region. Outside of the BEC region, the one-body density matrix ρ(r, r ) decays as |r − r | → ∞. The Mott insulator lobes appear at sufficiently strong interactions, as the inclusion theorem guarantees that there may be no MI for U < V .
In the past decades there appeared an enormous literature on disordered systems, see, e.g., [4][5][6][7][8][9][10][11][12][13][14], although due to the complexity of the topic mathematically rigorous results are still scarce. A great deal of work has been done on random Schrödinger operators, which describe single particle behaviour, see [15] for an excellent survey. In the recent paper [16] two of us have provided rigorous tight bounds on the ground state energy, as well as approximations of the ground state and excited states wave-functions for the case of random impurity (Bernoulli) disorder in a 1D lattice. Since in this case the shape of the low-energy states is known, there are analytical ways to estimate interaction energy and obtain rigorous results also for interacting systems.
A challenging and still open problem is the understanding of interaction effects on disordered systems. A crucial step forward in our understanding of dirty bosonic systems was made by Giamarchi and Schulz [6] and Fisher et al. [7] who conjectured that, in presence of weak disorder, interactions give rise to an intermediate compressible but insulating phase, the Bose glass, in between the superfluid and the insulating stronglyinteracting phases (Tonks gas in continuous 1D systems, or Mott insulator in lattice systems). The disputed controversy, as to whether or not the insulating Mott phase is always surrounded by a glassy phase was finally settled in [13], who demonstrated, using the theorem of inclusions, that this is indeed the case for the dirty boson system with generic bounded disorder. Experimentally, this problem has been recently addressed in [17][18][19].
In the regime where disorder dominates over interactions, it was shown instead that the ground state of the quantum fluid may be described in terms of a qualitatively different glassy phase, called Lifshits (or Anderson) glass [10,20], characterized by exponentially localized and well separated "islands". The latter may be identified with the low-lying single-particle eigenstates residing in the regions of space where the potential is small. As repulsive interactions are strengthened, an increasing number of islands is populated, until their overlap becomes sufficient to establish phase coherence and transport between them. The gas as a whole then undergoes a phase transition towards an extended Bose-Einstein condensate (BEC). If the gas is confined in a lattice, for even larger interactions the repulsion between particles becomes so strong that on-site density fluctuations characteristic of the superfluid state become energetically unfavourable. In this situation the gas undergoes a further transition out of the BEC, first to the Bose glass phase and then into the incompressible Mott state. Between the Lifshits and the Bose glasses there is no phase transition, as both are gapless, insulating and compressible phases, but nonetheless the two states are qualitatively different, in the sense that in the former the gas is fragmented into independent and distant lowenergy islands (the Lifshits states), while the latter tends to extend over a large portion of the available volume. The phases discussed above are sketched in figure 1, and their most important properties listed in table 1.
In the present paper we consider disorder of Bernoulli type, i.e., randomlydistributed, identical and localized impurities, and we focus on the regime of weak to intermediate interactions in presence of such impurities. We first investigate the properties of the ideal gas, discussing the typical size and energy of the localized islands, and showing that the low-energy states indeed form a so-called Lifshits tail.
We further explore interaction effects on the ground state properties by comparing the performance of two theoretical methods. The first one is a simple version of the multi-orbital Hartree-Fock method (sMOHF), a variational method based on an expansion on single-particle states. This approximation is appropriate to describe the glassy regime where superfluidity is suppressed by disorder [10]. More elaborate versions of the Hartree-Fock method, that incorporate self-consistency and allow to describe fragmentation and quantum dynamics of interacting Bose systems can be found in [21] and references therein. The second method is the standard Gross-Pitaevskii (GP) approach, which generally yields an appropriate description for every interaction strength.
We calculate the ground state energy of the system, the associated superfluid fraction, the fractal dimension, and introduce a new parameter, the fractional occupation, allowing one to discern between the Lifshits and the Bose glass. Finally, we turn to the investigation of temperature effects in an ideal gas. Due to the increase of kinetic energy, temperature yields effects similar to those generated by repulsive interactions, in the sense that the gas occupies a larger number of islands, which increasingly overlap until the ground state covers a large amount of the available space. Nonetheless, we will point out that there are important differences between the two cases.
As an important result, in this paper we show that disordered systems with Bernoulli potential allow for analytical estimates which are very well supported by numerical simulations. These results provided us with intuitions on how to generalize the rigorous results of [16] to non-interacting 2D systems, and, at least partially, to interacting 1D and 2D systems. These rigorous results go beyond the scope of the present article, and will be published in a more appropriate mathematical physics journal. This paper is structured as follows. We introduce the disordered potential under study and the relevant Hamiltonian in Sec. 2. In Sec. 3 we identify the eigenstates of the ideal system, we numerically demonstrate the Lifshits tail behaviour, and we discuss the ground state energy scaling as a function of the size of the system. In Sec. 4 we introduce the sMOHF and GP methods, and in Sec. 5 we compute a number of relevant quantities for interacting systems. In Sec. 6 we compare the effects of interactions and non-zero temperature, and we present our conclusions in Sec. 7. Table 1. Summary of most common phases for dirty bosons in a lattice, and associated properties.

The model
We study the properties of a bosonic gas on disordered 1D and 2D square lattices of linear dimension L. We consider disorder of the so-called Bernoulli type, i.e., the potential on each site is an independent random variable that assumes the value: This kind of potential may be ideally realized by using a two-component gas, and a species-selective lattice which deeply traps only one of the two components [19,22]. At sufficiently low energies, the component which is not trapped experiences s-wave collisions against localized scatterers. Alternatively, a Bernoulli disorder may be imprinted directly on the gas exploiting single-site addressable optical lattices [23]. The use of Bernoulli disorder is particularly appealing because its asymptotic properties, as discussed in the next sections, become apparent even for small lattices (e.g., with 50 D sites). As we show in Appendix A, the Bernoulli potential has actually optimal properties of convergence. Therefore, the Bernoulli disorder is suitable for numerical simulations and due to its simple form, it allows for various analytical estimates.
The Hamiltonian of the interacting system we consider is then where t denotes the tunnelling energy,ĉ † i is the creation operator of bosons on site i, and i, j denotes nearest neighbour pairs of sites. The interaction termÛ int will be discussed later in Section 4.

Properties of a non-interacting gas
Let us start by discussing the properties of a non-interacting system, and in particular look at its ground state energy, and the density of low energy excited states.

Ground state energy
In a large 1D system with Bernoulli disorder, the linear size L (1D) max of the largest island of contiguous zero-potential sites scales as L (1D) max ∝ log L [16]. The proof of this fact is based on the following simple observation. Let l 0 be the characteristic size of the largest island, then p l 0 estimates its occurrence. The number of islands of similar size is expected to be of order of L (strictly speaking L/ log L as we shall discuss below), so that the probability of occurring of any one of them is Lp l 0 → const as L → ∞, which indeed implies l 0 log L/| log(p)|. It is then clear that the ground state energy of the non-interacting bosonic gas in such a potential will be bounded above by the energy of the first harmonic, i.e., a half-sine function, of the largest island. A lower bound of the same order was proven in [16], showing that the ground state energy behaves asymptotically as π 2 /(log L) 2 when L → ∞.
Similarly, one can argue that in a large D-dimensional system the ground state will be occupying the largest spherical island of zero-potential sites since this shape minimizes the kinetic energy [24]. Its diameter can be shown to grow as L (D) max ∝ (log L) 1/D , and therefore we expect the energy of the ground state to scale as E 0 ∝ log L −2/D . Again, the proof is based on the similar argument as in the 1D case, except that one has to replace the characteristic length l 0 by the volume l D 0 in arbitrary dimension. The numerical confirmation of the scaling for the 2D case is shown in the inset of figure 2.
For comparison, we plot in the same inset the ground state energy obtained from a random-amplitude disorder (i.e., one where to each site corresponds a random potential V i , uniformly distributed in the interval [0, V ]). We see that the rate of convergence of the energy for a Bernoulli distribution is quicker than for a uniform distribution. Because of this, the potential with a Bernoulli distribution is ideal for ground state energy convergence as well as localization of low energy states. This may be quickly seen as follows, while further details are given in Appendix A. If the potential takes value zero with a positive probability p, one can define islands of zero potential as the natural locus of the ground eigenstate. In cases where the values of the potential are positive, even when they are arbitrarily close to zero, the low potential islands can still be defined, using a (volume-dependent) energy cut-off, but the contribution to energy from such islands will in this case be larger. Accordingly, the convergence of the ground state energy to zero will be faster in the former case. Among the distributions which assign a fixed probability p to the zero value of the potential, the optimal ones to work with are the Bernoulli distributions, in which some positive value V is assumed with probability 1 − p. While the rate of convergence of the ground state energy to zero is comparable for all such distributions, the advantage of the Bernoulli distribution is that it localizes the low energy states more clearly. The analysis of other potential distributions is more complicated, because it necessitates introducing an energy cutoff to define low potential islands. The Bernoulli potential is also easier to work with numerically, since it requires sampling fewer potential realizations.
Although we do not expect any of the properties discussed in the following to depend on the specific type of bounded disorder, we clearly see that the Bernoulli choice yields a much faster convergence of the ground state energy to the desired asymptotic behaviour.

Lifshits tail
Since the disordered potential is chosen to be finite, the spectrum of the system is bounded from below. As a consequence, the spectrum is expected to exhibit a socalled Lifshits tail [10,25], in the sense that the cumulative density of states (cDOS) N ( ) = d D( ) behaves at low energies as The cDOS and the associated Lifshits tail for the considered 2D system were obtained numerically and are shown in figure 2. The Lifshits states lie on well-separated islands, and have almost degenerate energies since the islands' diameters are approximately the same.

Interacting case
In this Section we will present two common theoretical approaches used to describe a Bose gas with short-range interaction potential. The first one corresponds to a simplified version of Multi-Orbital Hartree-Fock (sMOHF) method, based on an expansion into single-particle eigenstates. As we will see, sMOHF is suitable to describe a weaklyinteracting system, whose ground state occupies few disconnected large islands of zero potential. For more sophisticated version of the Hartree-Fock approach applied to bosons, see [21]. The second method is the standard Gross-Pitaevskii (GP) equation, which describes the dynamics of a fully-coherent matter wave, and correctly describes also the regime of strong interactions, where large overlaps between the islands and self-interaction on each island play an important role. overlap between the islands. Before going into details, let us estimate when the sMOHF description based on non-interacting single particle states should be valid. As we have discussed, this regime corresponds to the Lifshits glass region. Let us assume that we haver filled islands of the characteristic size l 0 , and "volume" l D 0 , so that where the proportionality constant C depends, in general, weakly on p and the density ρ, since it results from the complex interplay between the kinetic and interaction energy. In the following we will ignore this dependence.
The Lifshits glass regime occurs then for an inter-particle interaction coupling g smaller than the characteristic value g ch , at which the kinetic energy is comparable to the interaction energy. An estimate of the characteristic coupling is easily obtained by equating the kinetic energy per particle 1/l 2 0 , with the interaction energy gN/rl D 0 . Using the above expression forr, we obtain This scaling provides a good estimate of the coupling at which the two energies depicted in figure 3 start to diverge.

Multi-orbital Hartree-Fock approach
In the sMOHF treatment, the wave function is expressed in the single particle eigenstates, and is taken to be in the product form where |0 is the vacuum of the system and a † k creates a particle in the kth non-interacting eigenstate (i.e., orbital). We consider repulsive interactions of strength g > 0, which yield an interaction energy given by where the sum runs only over terms conserving number of particles, i.e., such that k 1 +k 2 = k 3 +k 4 . The diagonal terms (â † kâ † kâ kâk ) contribute a factor g 2 k n k (n k −1)O k,k , where we have defined the overlaps O k,j = dx|ψ k (x)| 2 |ψ j (x)| 2 .
There are two other possibilities which conserve the number of particles: (k 4 = k 1 and k 3 = k 2 ), or (k 4 = k 2 and k 3 = k 1 ). These two terms, called direct and exchange terms, are equal for bosons, and their sum contributes a factor g 2 j =k 2n k n j O k,j . The average interaction energy reads then: The occupation probabilities n k yielding the ground state may now be found by minimizing the total energy subject to the constraints of normalization and positivity of all n k , k n k = N, ∀k, n k ≥ 0.
We present the analytic solution of this problem in Appendix B.

Gross-Pitaevskii approach
The properties of the system may be analysed also in the framework of the usual GP equation, which describes the behaviour of an interacting Bose gas in terms of a single-particle coherent wave-function. We find the ground state and the chemical potential µ by imaginary time evolution, applying an operator split method. The ground state energy E is related to µ by the formula

Comparison of sMOHF and GP approaches
In this Section, we present the numerical solutions of the sMOHF and GP equations, comparing the performance of the two methods. In particular, we discuss and compare the ground state energy, the superfluid fraction, and the fractal dimension. As we will see, in the regime of strong disorder (Lifshits regime) sMOHF provides ground state energies in accord with GP. Nonetheless, sMOHF has the advantage of being insensitive to convergence issues, and its solution is generally much faster than GP. The fact that in this regime both methods agree confirms the intuition that the particles populate low-lying eigenstates, the Lifshits islands. In the regime where interactions dominate over the disorder, we will show instead that GP provides energies which are perceptibly lower than sMOHF. Unless otherwise noted, the numerics presented in this section have been obtained for a system on a 2D lattice with 32 × 32 sites, N part = 10 4 particles, a Bernoulli potential with V = 5t or V = 50t, and we have set t = 2 /2m = k B = 1.
Where needed, we have performed appropriate averages to obtain results which are independent of the particular disorder configuration.

Energy
In figure 3 we compare the energies obtained through the sMOHF approach with those obtained from the GP equation as a function of the interaction strength g. For a coupling constant g smaller than a characteristic value g ch , the minimization of expression (8) and the GP equation yield the same values of energy proving that sMOHF method correctly describes the system in this limit. The characteristic value g ch is given by the interaction strength at which the energy obtained from the sMOHF starts to differs from the one obtained by the GP equation. The range of g for which sMOHF approach is an equivalent description shrinks with increasing p, as the scatterers become increasingly sparser, leaving large regions of zero potential. The range of applicability of sMOHF also shrinks with decreasing V . For values of the interaction g g ch , the disorder plays a negligible role and the energy per particle saturates to the analytic value E = E + gρ/2, where E is a constant that depends on V and p. This analytic result is recovered by the numerical solution of the GP equation, but not by the sMOHF approach.

Superfluid fraction
The superfluid fraction ρ sf may be defined in terms of the energy change of a periodic system in response to twisted boundary conditions along one direction. To calculate this quantity we solve the GP equation (10) with boundary conditions ψ i+L,j = e iΦ ψ i,j , and we obtain where E(Φ) is the energy per particle of a system with the total phase shift Φ. Our results are depicted in figure 4. In the regime where sMOHF provides a good description, the SF fraction is negligible because the localized single-particle eigenstates experience an exponentially small energy change due to the imposed phase shift. The GP equation instead yields a perceptibly lower energy than sMOHF when the superfluid fraction ρ sf becomes sufficiently large (ρ sf 0.1).

Fractal dimension
As we have discussed in the previous sections, the ground state of a weakly-interacting fluid is localized on one or a few Lifshits islands. When interactions start to play an important role, the extension of the ground state increases generally until the gas occupies all available space. In order to provide a quantitative measure of the extension of the ground state, we calculate its fractal dimension d * [1], and introduce here the concept of fractional occupation c.
The fractal dimension is defined as a minimum d * such that Here P = 1/ dx|ψ 0 (x)| 4 = 1/O 0,0 is the so-called "participation number" of the ground state [1]. For an extended state such as a plane wave one finds that P is proportional to the volume of the system, and therefore the fractal dimension equals the Euclidean dimension, d * = D. For a state which is instead localized in a volume Q, the participation number behaves as Q, and the corresponding fractal dimension vanishes if Q grows slower than any power of L. For instance, the ground state of the non-interacting system is localized on an island of volume Q ∼ log L which for all α grows slower then L α . So, in general, one has 0 ≤ d * ≤ D.
Analysing the results obtained from both sMOHF and GP approaches, we find that d * = D for any g > 0. This may be seen as follows. Assuming that d * is bounded strictly below D, we see that the interaction energy (gN/2) dx |ψ 0 (x)| 4 will diverge. Physically, there is not enough space for the particles to keep the interaction energy bounded unless they spread out through a non-zero proportion of the whole space. show this, in figure 5 we plot log L P as a function of 1/ log L for a 1D system and a 2D system (inset).
The quantity c, which takes values between 0 and 1, can be interpreted as the fraction of space occupied by the ground state, and we will therefore refer to it as the fractional occupation of the ground state. We show our results for the fractional occupation in figure 6.
For strong disorder (V t) the fractional occupation c displays at intermediate interactions a plateau at ∼ p. This coincides with a similar plateau in the plot of the superfluid fraction (cf. figures 4b and 6b). In 1D the explanation for such behaviour of c and ρ sf is the following. For small g the wave function fills only a part of zero-potential space. Indeed, for g = 0, we have a linear operator and therefore the ground state is approximately the wave function on the longest island of zero potential, yielding a nearly zero c. For g small such that gN/2 is approximately 1/(log L) 2 , the ground state spreads to many of the long islands: the kinetic energy increases a little (on the order 1/(log L) 3 ) while interaction energy decreases by a factor 1/(# of islands). As g grows to be on the order of 1, it supports itself on all sites of zero potential save for the shortest islands, so c ≈ p. This is where the graphs level off in (Fig 6), where the ground state does not have large norm on short islands of zero potential and sites of V potential.
For large interactions one may use the Thomas-Fermi Approximation, cf. [5]. Suppose we choose an ansatz such that the wave function equals m √ pL on sites of zero potential and √ 1−m 2 √ qL , then energy optimizing m would be m 2 = p + pqV gρ and 1 − m 2 = q − pqV gρ . If gρ ≤ pV , then the above ansatz is invalid, and the ground state stays exclusively on sites of zero potential. If we compare the energy of a wave function equal to 1 √ pL on sites of zero potential and zero elsewhere to a wave function that is equal to 1 √ L everywhere, and ignore the kinetic energy terms (which we can control by another parameter), we see that the first wave function provides a better variational ansatz if gρ 2 ≤ pV . This means that for gρ 2pV , we are in an area where c ≈ p, that is, the kinetic energy and interaction energy are balancing, but any significant support of the ground state on the sites occupied by scatterers is too costly. Since the wave function does not spread significantly in this regime, also the superfluid fraction remains constant.

Comparison of temperature and interaction effects
In this section we turn to the analysis of the effects of non-zero temperature in a noninteracting system. We will see that those remind very much effects of interactions at zero temperature, and yet are quite different.

Single-particle density distribution
For free boson systems, the effects of temperature may be studied by expanding arbitrary states in the basis of non-interacting eigenstates, in analogy with what we did to treat the interacting case. In this approximation, the average density of the gas reads An approximation of a typical wave-function is then given by where e iφ k are random phases, and w k is a Gaussian random variable with distribution P(w k ) = n −1 k e −w k /n k , such that w k P = n k . The occupation factors depend on the statistics of the particles. Identical bosons follow the Bose-Einstein distribution Distinguishable (or classical) particles would instead populate the non-interacting eigenstates following a Boltzmann distribution, n k = e −(E k −µ)/T . At zero temperature, only the lowest energy state will have a non-zero population, i.e., n k = N δ 0,k . Non-zero temperatures will modify the occupation probabilities, redistributing population over the higher energy states, but we expect that for T t only the lowest energy states (the Lifshits states) will be populated.
We note here that in an interacting 2D Bose gas of finite extension one generally expects a crossover between two qualitatively different regimes. At very low temperatures (T < T C ) phase correlations in the gas decay algebraically with distance. For temperatures larger than the critical value T C instead, phase correlations decay exponentially with distance. This switch in the decay of correlations can be traced back to the dissociation of bound vortex-antivortex pairs at T > T C . The crossover becomes a real phase transition in an infinite system, and the underlying mechanism goes under the name of Berezinskii-Kosterlitz-Thouless (BKT) transition [26].
Our approach to finite temperatures is rather crude, and is not well suited to describe BKT physics, which requires considering both temperature and interactions at the same time. However, for the range of parameters considered in this manuscript, our a) 15 10  Figure 7. Comparison of temperature and interaction effects. First two rows: thermally averaged single particle density with respectively Boltzmann (top) and Bose (centre) statistics; from left to right, T /t =0, 0.1, 0.2, 0.5. Bottom row: single particle density of an interacting gas; from left to right, gρ/t = 0, 0.1, 2, 10. The red dots indicate the sites occupied by the disordered scatterers. The average density is ρ = 0.9particles/site. model provides a qualitative, although simplified, explanation of the complementary roles payed by temperature and interactions. Moreover, a similar scenario to the one considered here is present in three dimensions, where a true BEC exists at low temperatures.

Non-zero temperature
As temperature is raised, the cloud gets to occupy more and more islands; when T V it will populate all unoccupied sites, and for T V it will occupy even sites inhabited by scatterers.
The averaged densities in presence of positive non-zero temperature are shown in figure 7a) and b). While a classical cloud quickly spreads over various low energy states, a bosonic gas at low temperatures tends to condense in a single or few Lifshits states. Interaction effects in a bosonic cloud at T = 0, instead, yield at first sight a state which very much resembles the Boltzmann case. Even for moderate repulsive interactions, the cloud spreads over various Lifshits states. Effects of interactions at T = 0 are shown in figure 7c. The occupation probabilities n k for the three cases analysed above (Boltzmann, Bose, and interacting) are shown in figure 8. Since the lowest-energy single-particle states correspond to the localized and well separated Lifshits states, a population distribution which is narrow in energy corresponds to a ground state localized on one or very few well-separated islands, the Lifshits glass. As the temperature increases the Bose liquid tends to have a rather narrow population distribution, while the Boltzmann and interacting distributions quickly spread. Nonetheless, a few important differences may be noted comparing the two latter cases, as maybe seen in Figs. 7a and 7c. While the population distribution for the Boltzmann case smoothly decreases, with a long tail, the distribution for the interacting case has a local minimum at very low k, while at larger k goes through a maximum and then quickly drops to zero. These features can be understood in the following way. The eigenstates with lowest energies tend to have higher inverse participations (overlaps), i.e., the Lifshits states are rather concentrated on a single or few islands while states at intermediate energies are delocalized over the whole system. For an interacting system is therefore preferable to occupy states indexed by states with intermediate k values, as their larger spatial extension helps to reduce the interaction energy. At very high energies (k 1) instead the states have again very high O kk , since they become completely localized on single sites, and their populations abruptly drop to zero, as it is too costly energetically to put particles with repulsive interactions in small regions of space.

Conclusions
In this paper we have considered a bosonic lattice gas in the presence of Bernoulli disorder, given by randomly-localized identical scatterers. We have shown that in this case it is possible to provide precise analytical estimates even for the interacting case.
We have compared two theoretical schemes, the simplified multi-orbital Hartree-Fock and the Gross-Pitaevskii approaches, showing how the first is very accurate in the glassy regime of strong disorder, but it fails when interactions bring the system into an extended state and the superfluid fraction reaches values ρ sf 0.1. Further, we have shown analytically that the fractal dimension for this kind of potential tends to the physical dimension of the system. As a result, we have introduced a quantity termed fractional occupation, which characterizes the typical extension of the system. When the latter becomes of order p, i.e., the fraction of physical space where scatterers are absent, the system crosses over from the Lifshits to the Bose glass. Finally, we have discussed the similarities and differences between effects due to interaction and temperature. This question is of fundamental importance for ongoing experiments, since in a laboratory the two effects are often difficult to distinguish. We hope that our results will provide new insights in the complex route towards understanding the interplay between disorder, interactions, and temperature in quantum mechanical systems.
For a specific and related probability P ( ), the typical largest island composed by sites where the potential is less than has size L max ≈ log (P ( )(1 − P ( ))L) log P ( ) (A.2) Then the energy of the ground state is expected to be approximately (| log P ( )|) 2/D (log(P ( )(1 − P ( ))L)) 2/D + mγ γ+1 2(γ + 1)(a + m γ) (A.3) For this to be minimized in the large system limit, the order of each term must be approximately the same, so must scale with L as: For the Bernoulli distribution, a = p and m = 0, so the rate of log L −2/D is recovered. For the uniform distribution, a = 0 and m = 1 V and has rate of convergence ( log log L log L ) 2/D which is slower than the Bernoulli distribution. In fact, any distribution with a = 0, which means the potential can equal exactly zero with positive probability, must have its associated ground state energy converge to zero faster than for distributions with a = 0. Considering distributions with a > 0 and comparing the cases m = 0 and m > 0, the ground state energy in the former case will be larger than in the latter case, but asymptotically they will converge to zero at the same rate. More importantly, in the case where m = 0, the longest island is clearly defined whereas its definition involves an arbitrary (or artificial) cutoff. in the case where m > 0.

Appendix B. Solution of the sMOHF equations
The occupation probabilities n k may be found by minimizing the total energy Now, γ and µ should be chosen such that the following conditions, named KKT [27], are fulfilled: The KKT conditions allow for solving the optimization problem in the case in which some of the constraints are in the form of inequalities. These conditions lead to a nonlinear system for γ and µ.