Mean-field phase diagram of disordered bosons in a lattice at non-zero temperature

Bosons in a periodic lattice with on-site disorder at low but non-zero temperature are considered within a mean-field theory. The criteria used for the definition of the superfluid, Mott insulator and Bose glass are analysed. Since the compressibility does never vanish at non-zero temperature, it can not be used as a general criterium. We show that the phases are unambiguously distinguished by the superfluid density and the density of states of the low-energy exitations. The phase diagram of the system is calculated. It is shown that even a tiny temperature leads to a significant shift of the boundary between the Bose glass and superfluid.


Introduction
The remarkable experimental control over ultracold atomic gases in optical lattices acquired in the last couple of years [1,2,3,4,5,6] has opened up completely new lines of interest in the field of Bose-Einstein condensation. One of these are ultracold atoms in optical lattices with disorder which can be created by several methods. One of the possibilities is to use a laser speckle field [7,8,9]. An alternative way of creating a disorder potential is the introduction of a tiny fraction of a second atomic species which are strongly localized on random sites [10,11,12]. Tuning to a Feshbach resonance, these random scatterers can even make the disorder very strong. The random potentials for atoms can be also created via the spatial fluctuations of the electric currents generating the magnetic wire traps [13,14,15] or with the aid of the incommensurate lattices [16,17]. Disordered lattices for ultracold rubidium atoms have been recently created by superimposing a regular periodic optical potential on the speckle field [18] and on the incommensurate lattice [19].
These very recent experimental developments prompt the theoretical question what the behavior of ultracold atomic gases in lattice potentials with disorder might be. Typical for the modern area of Bose-Einstein condensation (BEC) in particular and ultracold atomic gases in general is that such basic theoretical questions have been asked before in a different context. Indeed, there is an enormous body of literature on ultracold Fermi gases in disordered lattices referring, first of all, to electrons in amorphous solids, either in the normal or in the superconducting state [20]. The possibility to study fermionic atomic gases from this point of view is a very active and interesting field with which we shall not deal here, however. Suffice it to say that completely new questions appear like the BEC-BCS crossover [21,22]. The influence of disorder on this problem is still completely unknown.
Here, we shall focus exclusively on bosons in potentials with disorder. This is also known as the 'dirty boson problem'. It first came up in the pre-BEC area in the context of experimental investigations of superfluidity of 4 He in the random pores of Vycor. The surprising finding that, for sufficiently low coverage of the pores, the 4 He superfluidity would disappear even in an extrapolation to zero temperature [23,24,25] prompted many theoretical studies. These were based on Hartree-Fock theory [26,27], generalizations of the Bogoliubov and the Beliaev theory [28] for random potentials [29,30,31,32,33,34,35,36], field-theoretical considerations [37,38], meanfield theory [38,39], renormalization group theory [40,41,42,43,44] and quantum Monte Carlo simulations [45,46,47,48,49,50,51,52,53,54] as well as numerical diagonalizations [55]. The consensus which developed from these studies is that in a disordered lattice at temperature T = 0 two new phases of bosons may exist besides the superfluid phase which is theoretically defined by the presence of the off-diagonal long-range order: One is the Mott-insulator phase, which only exists at commensurate fillings of the lattice, and is distinguished by the absence of the off-diagonal long-range order, a non-zero energy gap and a vanishing compressibility. The other is the Boseglass phase, which is distinguished again by the absence of the off-diagonal long-range order, a non-vanishing density of states at zero energy, and a non-zero compressibility. A more recent suggestion towards identifying the Bose-glass phase has been made in Ref. [56].
While these operational definitions of a Mott-insulator phase and a Bose-glass phase at temperature T = 0 are precise and clear-cut, they run into the obvious difficulty that experiments and quantum Monte Carlo simulations are never performed at T = 0. It is therefore necessary to examine the extent to which these or similar definitions can be applied at least at small non-zero temperature. This is the goal of the present paper. In order to achieve our goal we have to investigate the low-lying states of a suitable model of strongly interacting bosons in a lattice with disorder. We shall choose for this purpose a Bose-Hubbard model with on-site disorder of bounded variation. For our purpose of defining and distinguishing the various phases at non-zero temperature it is sufficient to analyse the basic model within a mean-field theory which allows to detect all the phases and to see what happens if the system parameters are changed, although it does not provide precise conditions for various phase transitions encountered in the system. The paper is organized as follows. In Section 2 the model and the mean-field approach to its analysis are defined. Then, in Section 3, the phase boundary between the superfluid and the two non-superfluid phases is derived from the condition of vanishing off-diagonal long-range order. This is first done for the pure case without disorder, reproducing a well-known result, and then generalizing it to the case with disorder. In Section 4 follows the definition of the Bose-glass and Mott-insulator phases at non-zero temperature and the examination of the phase boundary between them. A common criterion for the distinction of these two phases at T = 0 is the compressibility. How well this quantity serves this purpose at finite temperature is therefore examined in Section 5. In Section 6 the paper ends with some final conclusions.

Hamiltonian
We consider a system of spinless bosons in a homogeneous infinitely extended lattice of dimension d = 1, 2, 3 described by the Bose-Hubbard Hamiltonian (in units ofh = 1) where J is the tunneling matrix element, U is the on-site interaction constant, and µ is the chemical potential. In this work we assume that the random on-site energies ǫ i at different sites are uncorrelated and equally distributed with the probability density p(ǫ). We introduce the superfluid order parameter ψ = a i , where . . . = Tr [. . . exp(−βH)] / Tr [exp(−βH)] and . . . = Π i +∞ −∞ . . . p(ǫ i )dǫ i denote quantummechanical and disorder averaging, respectively, and β = 1/(kT ). Making use of the decoupling mean-field approximation in the hopping term [57,58,59] which is valid for sufficiently high-dimensional systems, we obtain the following on-site Hamiltonian where we have omitted the site index. The phase diagram of the system can be obtained in the following manner. First of all one has to calculate the disorder-averaged free energy of the system F (ψ) corresponding to the Hamiltonian (2). Then minimizing it with respect to ψ to determine ψ = ψ m one can distinguish the superfluid (ψ m = 0) and non-superfluid (ψ m = 0) regions of the parameter space. Applying a small phase gradient to the atomic matter-field operator and calculating the corresponding correction to the free energy in a manner similar to Ref. [60], one can show that the superfluid density following from the Hamiltonian (2) equals |ψ m | 2 . In the non-superfluid region, one has to work out the disorder average of the static superfluid susceptibility χ or the density of states ρ(ω) for the single-particle excitations [38]. In the region where ρ(ω) = 0 in the interval 0 ≤ ω < ω g we have, by the definition we apply, the Mott-insulator phase with the energy gap ω g . On the other hand, again by definition, the Bose-glass phase occurs when lim ω→0 ρ(ω) = 0 which corresponds to a divergent superfluid susceptibility [38].
The form of the mean-field Hamiltonian (2) implies that, in our approximation, the properties of the Mott-insulator phase as well as the Bose-glass phase, where ψ vanishes, do not depend on the tunneling matrix element J. This is consistent with the fact that the boundary between these phases occurs only for small values of J. The transition to superfluidity, where ψ m = 0 starts to appear, does depend on J also in our approximation.

Boundary between the superfluid and non-superfluid phases
In order to calculate the free energy of the system, one has to solve the eigenvalue problem for the Hamiltonian (2). This can be done exactly by means of numerical calculations. However, the boundary between the superfluid and non-superfluid phases can also be determined with high accuracy treating the first term in the Hamiltonian (2) as a perturbation. The free energy can only depend on |ψ| 2 since a change of the phase in ψ can be undone by the unitary transformation a → exp(−iϕ)a, a † → exp(iϕ)a † . Indeed, the calculations show that the result has the following structure: The explicit form of a 4 as well as a 0 and a 2 for T = 0 was obtained in Ref. [59]. The generalization to T = 0 and the average over disorder needed here is straightforward.
Since a 4 turns out to be always positive and a 2 can be either positive or negative, the superfluid/non-superfluid transition is of second order. The equation a 2 = 0 determines the phase boundary which is given by being the partition function without hopping.

Pure case
In the pure case we have p(ǫ) = δ(ǫ), so we get from (4) for the phase boundary [61] 2dJ In the zero-temperature limit this equation reduces to the well-known result for the boundary between the superfluid and Mott-insulator [58,59] Here n denotes the positive integer at which E m (µ) is minimal with respect to m. This fixes n as the smallest integer larger than or equal to µ/U.

Disorder with homogeneous distribution
In the following we choose for simplicity a homogeneous disorder distribution in the interval ǫ ∈ [−∆/2, ∆/2] so we have for the phase boundary (4) 2dJ ∆ We discuss first the special case T = 0 and consider ∆ < U. It is assumed that µ ∈ [U(n − 1), Un], n = 1, 2, . . .. Eq. (9) then gives if µ ∈ [U(n − 1) + ∆/2, Un − ∆/2], otherwise 2dJ = 0. In the limit ∆ → 0, Eq. (10) reduces to (7). The phase boundary at T = 0 following from Eq. (10) is shown in Fig. 1a for typical parameter values. It has a lobe structure and the size of the lobes decreases with increasing ∆. If T = 0 but ∆ > U, we obtain 2dJ = 0 as the transition line for any value of µ, i.e., the lobes disappear and the superfluid phase appears as soon as J is turned on [38].
In the case of non-zero temperature, Eq. (9) gives always non-vanishing values of 2dJ. The boundary between the superfluid and non-superfluid phases for small T and different values of ∆ is shown in Figs. 1b,c,d. The plots obtained from Eq. (9) and that obtained by numerically diagonalizing the Hamiltonian (2) and minimizing the free energy with respect to ψ are indistinguishable. If the temperature increases, the boundary between the superfluid and non-superfluid phases goes upwards and the size of the non-superfluid region grows. The presence of even a small temperature (kT /U = 0.01) changes rather strongly the phase boundary. This is in contrast to the pure case where values of kT /U ∼ 1 are required in order to get a noticeable shift of the boundary between the Mott-insulator and the superfluid phases [62]. The superfluid density |ψ m | 2 obtained by the numerical diagonalization for the same parameters as in Fig. 1b and the values of J indicated by dotted lines in Fig. 1b is plotted in Fig. 2.

Boundary between the Mott-insulator and Bose-glass phases
In the Mott-insulator as well as in the Bose-glass phase the superfluid order parameter ψ = ψ m vanishes, which implies that J disappears from the mean-field Hamiltonian (2). Therefore, the properties of these two phases do not depend on J in the mean-field approximation. This is related to the fact that the two phases are localized and, therefore, the dependence on J should be weak. However, more accurate calculations beyond the mean-field theory should give some dependence on J, in particular close to the boundary with the superfluid phase. The Mott-insulator phase is characterized by   the gap in the excitation spectrum and it has a finite superfluid susceptibility. The Boseglass phase has no gap and the superfluid susceptibility diverges. All this is directly related to the properties of the Green's functions and the density of states.

Green's function
The bosonic single-particle Green's function G(t) is defined as [63] where a(t) is the annihilation operator in the Heisenberg representation. Straightforward calculations lead to the result where µ ′ = µ + ǫ denotes the random local chemical potential. One can easily show that the imaginary-time Green's functions satisfy the periodicity condition G > (τ + β) = G < (τ ), where τ = it is the imaginary time.
At T = 0 the imaginary-time Green's function takes the form where n ± is the smallest integer greater than or equal to (µ ± ∆/2)/U. If n + = n − which corresponds to µ ∈ [U(n − 1) + ∆/2, Un − ∆/2] for ∆ < U, the first term in Eq. (13) vanishes and the superfluid susceptibility defined by the integral over the realtime Green's function χ = ∞ 0 G > (t)dt is a finite quantity. This means that we have not the Bose-glass phase, i.e., we are in the Mott-insulator phase (see Fig. 1a). If n + > n − which corresponds to µ ∈ [U(n − 1), U(n − 1) + ∆/2] ∪ [Un − ∆/2, Un] for ∆ < U or arbitrary µ for ∆ > U, the first term in Eq. (13) survives and renders χ divergent which is the distinguishing property of the Bose-glass phase. In the case ∆ > U, the lobes in Fig. 1a disappear completely which means that the Mott-insulator phase is destroyed by the disorder [38].
At non-zero temperature, it is more difficult to analyze the structure of the disorder averaged Green's function. Expanding Eq. (12) for large but finite values of β shows that the Green's function has a similar structure as Eq. (13) but the explicit expressions become very long and we do not display them. Typical τ -dependences of G > in the Mott-insulator as well as in the Bose-glass phase are shown for different temperatures in Fig. 3. Due to the different scales of τ , Figs. 3a and 3b indicate that the Mottinsulator and the Bose-glass phase are characterized by an exponential and algebraic decay of G > (τ ), respectively. These analytical results agree qualitatively with Monte Carlo simulations [47,48,49,52,54].
Since even at finite temperature G > (τ ) decays only like 1/τ , the superfluid susceptibility still diverges logarithmically. Vice versa, the divergent superfluid suceptibility χ has the consequence that the density of states for the single-particle excitations at zero energy does not vanish in contrast to the case of finite χ [38]. In the next section we will see that the density of states is easier to analyze at finite temperature than the full Green's function.

Density of states
The density of states for the single-particle excitations can be determined in terms of the Fourier transformed single-particle Green's functionG(ω) = +∞ −∞ d t exp (iωt) G(t) as ρ(ω, µ) = − 1 π ImG(ω) [28]. The Fourier transformation of Eqs. (11), (12) gives the density of states for the pure case The two δ-functions correspond to the hole and particle excitations, respectively. After the disorder averaging we obtain This disorder averaged density of states is plotted in Fig. 4a for the Mott-insulator phase with the energy gap ω g = Un − ∆/2 − µ and in Fig. 4b for the Bose-glass phase. Since E m+1 (Um) = E m (Um), we get finally For the homogeneous disorder distribution (8) the summation in Eq. (16) is restricted by m = n − , . . . , n ′ + , where n − is the smallest integer greater than or equal to (µ − ∆/2)/U and n ′ + is the greatest integer less than or equal to (µ + ∆/2)/U. If ∆ < U, Eq. (16) takes the form Simple analytical expressions (17), (18) show that the lines µ = U(n − 1) + ∆/2 and µ = Un − ∆/2 determine the boundaries between the Mott-insulator and Bose-glass phases at arbitrary temperature (see Fig. 1). In the case ∆ > U of strong disorder, ρ(0, µ) does not vanish for finite temperature and the Mott-insulator phase does not exist.

Compressibility
The compressibility of the system is defined as κ(µ) = −∂ 2 F (µ)/∂µ 2 , where F (µ) = −ln Z(µ)/β. In a non-superfluid phase Z(µ) = Z 0 (µ) is given by Eq. (5). Partial integration over ǫ gives For the homogeneous distribution (8) we get where is the mean particle number per lattice site in the pure case. For the Bose-glass phase, the compressibility (20) does not vanish. One can easily show that If βU ≫ 1, the compressibility (20) expanded for small temperatures has the form where δ(µ) = E n − min(E n−1 , E n+1 ) is the energy difference between the first excited state and the ground state in the pure case (cf. (5)), and α is some finite constant. This equation shows that the Mott-insulator phase, which occurs for n + = n − , has an exponentially small compressibility at non-zero temperature, in contrast to the Boseglass phase.
The dependence of the compressibility on µ for small temperature is shown in Figs. 6,7. Since the compressibility does not vanish at non-zero temperature and is a continuous function of the system parameters, it can not be used as a criterion to distinguish between different phases. Thus, we deduce that the transitions are better defined in terms of the superfluid order parameter ψ and the density of states.

Conclusions
The Mott-insulator phase and the Bose-glass phase at vanishing temperature can be defined either by their thermodynamic properties or by the spectral properties of their low-lying excitations. Both characterizations are, of course, closely related. Both phases  Fig. 1b). The results are obtained by numerically diagonalizing the Hamiltonian (2) and calculating the second derivative of the free energy with respect to µ at ψ = ψ m . are non-superfluid, i.e., the corresponding Goldstone modes, the phonons, are absent. In the case of the Mott-insulator phase the spectral characterization by an energy gap implies a vanishing compressibility and vanishing particle-number fluctuations at T = 0. In the Bose-glass phase the non-vanishing density of states at zero energy implies a nonvanishing compressibility. These features allow a sharp distinction between the two phases at zero temperature.
However, at non-vanishing temperatures, the characterization of the Mott-insulator and Bose-glass phases by their thermodynamic properties is no longer sharp -the Mottinsulator phase has an exponentially small but finite compressibility which corresponds to non-vanishing fluctuations of the particle number density. Still, as we have pointed out in this paper, the characteristic spectral features remain present also at T > 0 and can therefore be used for a sharp definition and distinction between these low temperature phases. We employed this possibility to calculate finite-temperature phase diagrams within a Bose-Hubbard model with on-site disorder within the mean-field approximation. For experiments with optical lattices it is usually easiest to change system parameters like the tunnelling amplitude J at fixed temperature, i.e., phase diagrams of the format of Fig. 1, where system parameters are used as variables, are most natural from this point of view. However, from a thermodynamic point of view, it may be more natural to give the phase diagram in the (µ, T )-plane. This is done in Fig. 8 for the two values of J marked in Fig. 1b as dotted lines. On the high temperature side, the phase diagram of Fig. 8 is incomplete, because there the transition to the normal gas phase must occur, which we have not considered in the present work. For both the Mott-insulator and the Bose-glass phase this transition would be sharp, if the energy gap or the finite density of states would start to appear suddenly at a critical temperature. Alternatively, the transition could also take the form of a smooth crossover. For the Mott-insulator phase the crossover could occur at the temperature kT ≈ ω g , where ω g is the energy-gap for thermal excitations [57,62]. For the Bose-glass phase this would happen at the temperature where the density of states starts to be dominated by the normal Bose gas.
In our calculations, we find direct transitions from the superfluid phase either to the Mott-insulator or to the Bose-glas phase under variation of J, depending on the values of µ and the disorder strength ∆ (Fig. 1). The character of the transition from the superfluid to the Bose-glas phase is in agreement with Monte Carlo calculations. The direct transition from the superfluid to the Mott-insulator phase in the case of weak disorder ∆ < U is also in agreement with the results obtained by the pathintegral Monte Carlo techniques [45] and by the Monte Carlo simulations based on the J-current model, which neglects amplitude fluctuations of the bosonic quantum fields, for µ/U = 0.5 [48,49,51,52] and for some finite interval of µ near the tip of the Mottinsulator lobe [50]. However, more recent Monte Carlo investigations of the J-current model [53], still neglecting the amplitude fluctuations, show that the superfluid and the Mott-insulator regions on the phase diagram are separated by a narrow region occupied by the Bose-glas phase if the number of lattice sites is large enough. Thus, different Monte Carlo techniques give different results on the character of the superfluid -Mottinsulator transition in the presence of disorder. Other methods lead also to mutually conflicting results concerning this point [64]. Here we have presented the results of the mean-field approach.
Experimentally, the Bose-glass phase may not be easy to identify with ultracold atoms in a suitably disordered lattice. In fact, it might best be identifiable indirectly by the absence of properties which are present in the competing phases for T → 0, like the absence of a macroscopic wave function and the absence of an energy gap or of incompressibility [19]. The finite density of states at ω → 0 would show up in a specific heat proportional to T for T → 0 and in a logarithmically diverging susceptibility χ. It would certainly be of great interest if a way could be found to measure ρ(0) directly.