Green's function approach to the Bose-Hubbard model with disorder

We analyse the distinction between the three different ground states presented by a system of spinless bosons with short-range interactions submitted to a random potential using the disordered Bose-Hubbard model. The criteria for identifying the superfluid, the Mott-insulator, and the Bose-glass phases at finite temperatures are discussed for small values of the kinetic energy associated with the tunnelling of particles between potential wells. Field theoretical considerations are applied in order to construct a diagrammatic hopping expansion to the finite-temperature Green's function. By performing a summation of subsets of diagrams we are able to find the condition to the long-range correlations which leads to the phase boundary between superfluid and insulating phases. The perturbative expression to the local correlations allows us to calculate an approximation to the single-particle density of states of low-energy excitations in the presence of small hopping, which characterizes unambiguously the distinction between the Mott-insulator and the Bose-glass phases. We obtain the phase diagram for bounded on-site disorder. It is demonstrated that our analysis is capable of going beyond the mean-field theory results for the classification of these different ground states.


Introduction
Since the seminal paper of M. P. Fisher et al. [1], the study of interacting bosonic particles in random potentials has become an active field of research. A key aspect of this system is the interplay between localization and superfluidity produced by the combined effect of randomness and interaction. The experimental realization of Bose-Einstein condensation in ultracold atomic gases together with the precise control presented by optical lattice experiments provided a unique possibility to investigate fundamental questions about this problem [2,3,4]. The random potential can be experimentally achieved, for instance, with magnetic wire traps [5,6], where imperfections of the wire produce local disorder, or even by using speckle laser fields, where a diffuse laser front creates the random lattice [7,8]. Perhaps the most pronounced phenomenon exhibited by this system is the superfluid to insulator quantum phase transition.
The theoretical description of spinless bosons with short-range interactions moving in random external potentials is usually based on the disordered Bose-Hubbard Hamiltonian (BHH) whereâ † i andâ i are the bosonic creation and annihilation operators fulfilling the canonical commutation relations,n i =â † iâ i denotes the number operator, µ is the chemical potential, and the sum ij runs over nearest neighbours. In addition, the interaction between two particles at the same lattice site is parametrized by the energy U , the hopping parameter J corresponds to the kinetic energy associated with the tunnelling of a particle from a lattice site to one of its first neighbours, and the on-site energies i represent local imperfections which we assume to be uncorrelated at different sites and thus randomly spread over the lattice obeying some probability distribution p({ i }). A detailed analysis [9] demonstrates that the disorder generated,, for instance, by laser speckles would make all parameters random and depending on the speckle potential distribution. However, we will assume that (1) is valid as current experiments show that by choosing an appropriate holographic mask when imaging under a quantum gas microscope one can create arbitrary potential landscapes [10]. Additionally, it was shown in [11] that the speckle intensity can be customised to generate a wide variate of speckle patterns including a uniform distributions.
The competition between the system parameters gives rise to different phase transitions. If the hopping energy is much larger when compared to the interaction energy, the atoms can move without viscosity over the system's volume and the ground state is superfluid. In the opposite case, in which interactions dominate over the tunnelling energy, a finite number of atoms becomes localized around each potential minimum configuring a Mott-insulator state. The superfluid to Mott insulator transition was directly tested by observing the multiple matter wave interference pattern presented in absorption pictures of time-of-flight measurements taken for different lattice potential depths [12]. When disorder is introduced, a third phase intervenes between the latter two: the Bose-glass phase. This phase consists of rare superfluid regions inside an insulating background and for sufficiently strong disorder strength it can even destroy the Mott-insulator state. The superfluid-Bose glass transition has dynamically been probed using a quantum quench of disorder in an ultracold gas at non-zero temperature and measuring its excitations in an experiment by Meldgin et al. [13]. Although some attributes of the Bose-glass phase are well understood, detailed information concerning its phase boundary as well as on the nature of its low-lying excitations are still lacking both from experimental and theoretical points of view.
The analytic form of the BHH eigenstates and eigenenergies cannot be directly computed in general. Thus, with the aim of reaching predictions for the phase boundaries of the transitions, numerical methods such as Monte-Carlo simulations [14,15,16,17,18,19] and stochastic [20,21] as well as local [22] mean-field techniques have been applied, while analytic investigations have mostly been confined to the method of mean-field theory [23,24,25]. More recently, the application of field-theoretical approaches proved to be efficient in the analysis of the above described quantum phase transitions providing precise results when compared to Monte-Carlo simulations in the pure case [26,27,28,29]. The critical exponents for the superfluid to Bose-glass phase transition have been numerically calculated in [30]. Analytical information concerning the impact of temperature on the phase boundary between Mott insulator and Bose glass was obtained in [23]. In this study, the zero temperature characteristics of these two states were tested in a finite-temperature theory by analysing the single-particle density of states. Unlike the superfluid state, both insulating phases are distinguished by the absence of off-diagonal long-range order. However, the existence of a finite energy gap for particle-hole excitations in the Mott-insulator phase leads to a vanishing density of states at zero energy, in contrast to the Bose-glass phase which presents a gapless single-particle excitation spectrum and consequently a finite zero-energy density of states [1]. Even though such distinctions were tested on a finite-temperature theory, the influence of finite values of the tunnelling energy to their phase boundary remains to be studied. Therefore, it is important to investigate to which extent these definitions or similar ones still hold at finite temperatures and at least for small values of the hopping parameter. In this paper, we demonstrate that such an analysis can be performed by calculating corrections to the Green's function due to the hopping of particles. To this end, we develop a perturbative treatment to the BHH considering bounded onsite disorder. Using field theoretical considerations, similarly to [26,27,28,29], we construct a hopping parameter expansion to the two-point correlation function at finite temperatures. The phase boundary to the superfluid phase is identified by observing the divergence of the resummed expression of the correlation function, while the phase boundary distinguishing the Mott-insulator from the Bose-glass phase is computed by analysing the imaginary part of the same-site correlation function in real frequency space, which corresponds to the single-particle density of states.
In what follows, we first construct the perturbative expression to the Green's function in section 2 and then proceed the investigation to obtain the superfluid to insulator phase boundary in section 3. In section 4, we analyse the imaginary part of the local Green's function in order to calculate the first relevant correction to the singleparticle density of states. In section 5, we analyse the distribution of such a quantity to determine the Mott insulator to Bose glass phase boundary. A comparison of our results for the first Mott lobe with the numerical predictions of [20,21,22] for the 2D and 3D cases at both zero and finite temperatures is presented in section 6. In section 7, we then make a summary of our findings and conclude that this method is capable of going beyond mean-field theory for the analysis of the phase transition.

Perturbation theory
We base our analysis on the single-particle Green's function as it characterizes the microscopic properties of the system. This function can be defined as the thermal average of the time-ordered product of the bosonic creation and annihilation operators in the Heisenberg representation whereT is the time-ordering operator and the Heisenberg representation for an arbitrary Schrödinger operatorÔ s is defined asÔ(τ ) = e τĤÔ s e −τĤ , with = 1. As one of our main interests is to describe the system at finite temperatures, we have used the Wick rotation t → −iτ to establish the imaginary-time formalism [31,32], where τ is the so-called imaginary time.
As an exact diagonalization of the BHH is not possible, the Green's function as well as other important quantities of the system may be calculated perturbatively. To this end, we first consider the BHH as belonging to a general class of Hamiltonians composed of a local term plus a hopping term where J ij is symmetric in i and j and J ii = 0. The BHH is recovered by settinĝ H 0 i = Un i (n i − 1)/2 − µ ini , where µ i = µ − i , and considering hopping only between first neighbouring sites. Following field-theoretic considerations [32,33,34], we then include a source term to the Hamiltonian with the intention of explicitly breaking any global symmetrieŝ Using the Dirac interaction picture, the initial value problem for the imaginary-time evolution operator takes the form whereĤ I (τ ) is the interaction picture representation of the hopping term plus the source termĤ Equation (5) has a solution which is given by the Dyson serieŝ Starting from the fully localized case, J ij = 0, the hopping-free partition function can be written as whereÛ with β = 1/k B T , where k B is the Boltzmann constant and T the temperature. Using the semi-group property of the imaginary-time evolution operator [28], we can express the full partition function as a power series in the hopping matrix elements given by With the same property it is possible to show that the Green's function can be calculated by taking functional derivatives of the full partition function with respect to the sources and then considering the limit where they vanish Thus, the problem is reduced to finding the expression of Z 0 and then taking its functional derivatives with respect to the sources in order to account for the hopping contributions. This calculation is simplified if we notice that in the local case the hopping-free partition function becomes a product of single-site contributions As a consequence, the free energy W 0 = −βF 0 = ln(Z 0 ) also becomes local. Hence, in the zero hopping case we have and we can expand each local term in the summation as a series in the sources where the functions W (2n) 0i (τ 1 , . . . , τ n ; τ 1 , . . . , τ n ) are the so-called local 2n-point correlation functions. We focus our analysis on the calculation of the 2-point correlations. For this purpose, it is convenient to use a diagrammatic notation similarly to [35]. To construct such a notation, we associate each 2n-point function W (2n) 0i with a vertex labelled with a site index i and containing n entering lines and n exiting lines which correspond to the imaginary-time variables τ n and τ n , respectively. Thus, the terms shown in the above equation are respectively represented by The notation can be shortened by suppressing both the vertex label when summing over all lattice sites, and each inward (outward) line label when multiplying by j i (j * i ) and integrating from 0 to β in the imaginary-time variable. Thereafter, we can rewrite (13) together with (14) as a sum of 1-vertex diagrams As our approach is applied only to the 2-point diagram, we choose not to show the subsequent terms in this expansion, which would consist of all the diagrams with an even number greater than 2 of external lines joined by a single vertex. The next diagram in the summation, the 4-point diagram, becomes important, for instance, in the effectiveaction approach of [26,27,28,29].
To compute the hopping corrections to the free energy, we must apply (10) to Z 0 = exp(W 0 ) and then take the logarithm of the result. In diagrammatic notation, the hopping matrix can be denoted by an internal line between two vertices, and the effect of the functional derivatives with respect to the sources j i and j * i on a local diagram is the introduction of an index i to its central vertex and the addition of an imaginary time variable τ to its inward or outward lines, respectively. The result of the full operator acting on products of diagrams is to generate different diagrams by joining an inward open line of one diagram with the outward open line of another. The linked-cluster theorem [36,37] assures that only the connected diagrams will contribute to W . Therefore, the hopping expansion to the free energy is given by where we have considered only tree-level corrections to the 2-point diagrams, which are the ones that will be essential to our analysis. Note that (18) differs from (16) as it contains all simple-chain diagrams with two external lines, where the order of each diagram in the hopping approximation is determined by the number of internal lines between the vertices. By using (18), to write Z[j, j * ] = exp(W [j, j * ]), and then applying (11) we obtain the Green's function which consists of a sum of all simple chain diagrams. With this expression we can now proceed to determine the phase boundaries for the quantum phase transitions.

Superfluid phase boundary
The evaluation of the diagrams introduced above involves an integration in the imaginary time variables. This process can be carried out by considering the transformation to the Matsubara frequency space with the bosonic Matsubara frequencies ω l = 2πl/β, where l ∈ Z. It turns out that working in the frequency space simplifies our calculations as the system that we are considering presents time-translation invariance. The diagrammatic expansion maintains the same form in Matsubara space. As shown in [33], the full correlation functions can be decomposed into connected correlation functions. In our theory, this decomposition for the 2-point function gives where and Θ(τ ) is the Heaviside step function. In the Matsubara representation, we get where δ ωω is the Kronecker delta and we have defined In order to deal with the random chemical potential, we consider the disorder to be frozen in time. Also, the magnitude of the local chemical potential at different lattice sites can be regarded as spatially uncorrelated, thus varying within a range where each value appears with a specific probability. This is equivalent to assume a lattice spacing which is much larger than the disorder correlation length. Accordingly, the local disorder i is assumed to be characterized by a probability distribution p( i ) of some kind. As a consequence, it becomes necessary to define a disorder ensemble average which has the following form The transition to the superfluid phase is characterized by diverging long-range correlations [32,33]. However, any finite order approximation in the hopping expansion (19) of the Green's function is a power series in J and therefore analytic. Hence, we must perform a summation of the infinite subset of chain diagrams, in the same manner as in [29]. This task becomes straightforward if we make the transformation of the disorder average of (19) to quasi-momentum space The above expression can be rewritten as with where J(k) = 2J D α=1 cos(k α a) is the D-dimensional lattice dispersion. This equation consists of a geometric series which can be directly evaluated to Note that the imaginary part of g i (ω l 1 ) in (24) vanishes for ω l 1 = 0. Furthermore, as phase transitions are governed by long-wavelength fluctuations [32,33], we must also set k = . In this situation we find that the Green's function diverges in J when where z = 2D is the lattice coordination number. This result for the phase boundary is exactly equal to the one obtained in [23] with mean-field theory. For the uniform disorder distribution the phase boundary becomes The result of the above equation for different disorder strengths can be observed in figure 1. There, we can see the phase boundary between superfluid and insulating phases at a finite temperature. Above the red line only the superfluid phase can exist. The region enclosed below this line corresponds to non-superfluid phases. We can see that, with an increase in the disorder strength, from figure 1(a) to 1(b), the phase boundary becomes smoother. However, no distinction between Mott insulator and Bose glass can be made in this diagram. As previously stated, such a distinction can be clarified by analysing the behaviour of the density of states for single excitations below the superfluid region, which is what we develop in the following.

Single-particle density of states
According to [1], for vanishing hopping, the density of states of the low-lying excitations in the well-localized regime of the Bose-glass phase is constant at zero excitation energy by virtue of the continuous distribution of the random potential. It was further stated that this situation should be sustained when the hopping parameter is made slightly positive resulting in a single-particle density of states also constant at zero energy. In light of these arguments, we now consider the hopping corrections of our perturbation theory to the local correlation function from which we can obtain the density of states of single excitations.
In the same site case, i = j, the first-order term in (19) vanishes due to the fact that there should be no possibility of a particle hopping from a site to itself, J ii = 0, and therefore the lowest relevant contribution to the local Green's function is of second order, resulting in the following diagrammatic expression Note that by using local correlations we restrict our calculation to be valid only inside the non-superfluid part of the phase diagram. In Matsubara space, (33) can be written as which explicitly reads with a (1) n,m (µ i , µ j ) = where we show only those terms which will be important to our analysis. We have chosen to represent the double poles implicit in (34) as derivatives of first-order poles. This representation simplifies the process of analytic continuation which must be taken in order to go from the Matsubara frequencies to the real frequency domain. Such an analytic continuation can be expressed by the transformation iω l → ω ±iη, with η → 0 + . Following [31], the density of states can be obtained by using ρ i (ω) = − 1 π Im(G i (ω)). Hence, in the second-order hopping expansion it gives By taking the disorder ensemble average of the density of states we get where which, in the case where ω = 0, is essentially the same integration as (30). The derivatives that appear in the second-order correction of (39) are a product of our perturbation theory. In the exact solution, the Green's function should only present simple poles in Matsubara space, as is demonstrate, for instance, in [38,Chapter 9] or [39,Chapter 3]. Therefore, these anomalies in the expression of the density of states are related to the double poles implicit in (34). In order to deal with such a problem, one must renormalize the location of these poles. For that reason, we propose the following transformations where the initial situation is recovered by setting λ = 1. The idea is to use these transformation to expand the density of states in λ and then determine the values of the renormalized frequencies that eliminate the derivatives that appear in (39).
In the Poincaré-Lindstedt method, the renormalization of the frequency aims to eliminate secular terms in a response that should be periodic [40,41,42,43]. If we take a Fourier transform in time from the Poincaré-Lindstedt method, purely periodic terms will appear as simple poles, while secular terms will appear as higher order poles. Analogously, we can understand the delta functions as the imaginary part of the simple poles and its derivatives as the imaginary part of higher order poles. In this fashion, we can interpret the renormalization of the poles in the Green's function approach as the Poincaré-Lindstedt method in the Fourier space. Therefore, the renormalization of the delta functions is equivalent to the previous method considering only the imaginary parts.
Using the transformations of (41), the first-order in the λ expansion for the density of states gives Setting λ = 1, we now choose the values of Ω − n and Ω + n such that they cancel the derivatives of the disorder distribution in (42) thus finding the following expressions The result of this process is that the relevant correction to the density of states in the second hopping order expansion consists of a shift in the argument of the disorder where we have considered hopping processes only between first neighbouring sites. Note that (45) constitutes a resummation of (39) and consequently both equations are equivalent up to second order in J. With this result, we have therefore computed the corrections due to the hopping of particles to the expression of the density of states calculated in [23]. Now the task is to solve the integral of (40). By choosing a uniform disorder distribution of the form of (31) such an integration is reduced to As explained earlier, the Mott insulator and the Bose glass can be distinguished from one another by the fact that, in the case of ω = 0, the single-particle density of states is zero in the Mott phase and finite in the Bose-glass phase. The result of (45) allows us to construct a phase diagram by distinguishing on a density plot the regions below the superfluid line (32) where ρ(0, µ) vanishes from the ones where it is finite, which can be observed in figure 2 for homogeneous disorder. As argued in [1], in the case where ∆ < U , there should always be a region, namely µ ∈ [U (n − 1) + ∆/2, U n − ∆/2] for J = 0, where the average number of particles that minimizes the local energy f n (µ i ) given in (22) becomes fixed at an integer value, which represents the Mott insulating state. When J is made slightly positive, this situation is sustained and both Mott insulator and Bose glass are guaranteed to appear at different regions of the insulating part of the phase diagram. As can be seen in figure 2(a), the Boseglass phase, represented by the grey regions, emerges between the Mott lobes which correspond to the black regions. In the case where ∆ ≥ U , the average number of particles per site never sticks to an integer value and the Bose-glass phase suppresses the Mott states dominating the insulating part of the phase diagram, which is shown in figure 2(b). In the two cases, above the red line, only the superfluid phase exists. The result shown in figure 2 demonstrates that the single-particle density of states distinguishes unambiguously the Mott insulator from the Bose-glass state, however it gives no detail about the exact phase boundary between the two. In the situation where the three phases coexist, our perturbation theory should be valid to determine the such a phase boundary. This calculation will now be presented.

Mott insulator to Bose glass phase boundary
We can calculate the analytical phase boundary between Mott insulator and Bose glass by analysing the distributions in (45). For the uniform disorder distribution of (31), the Bose-glass region, where ρ(0, µ) = 0, is limited by the following inequalities where we have defined As explained in the last section, for sufficiently weak disorder, ∆ < U , the Mott states should always appear in the insulating part of the phase diagram. Therefore, based on (47), the Mott lobes would correspond to the regions satisfying the following condition to the hopping parameter What stands out in this result is the fact that the location of any Mott lobe can be obtained by plugging in the above equation the integer n which minimizes the local energy. The phase boundaries of (32) and (49) are shown in figure 3 for a uniform disorder distribution. According to our analysis, the region inside the dashed black lines corresponds to the Mott lobes, where ρ(0, µ) = 0, and outside of this region, where ρ(0, µ) = 0, we would have the Bose-glass phase. With our approximation, which corresponds to calculating hopping-dependent corrections to the density of states analysed in [23], we have found an equation for the phase boundary between the two insulating phases valid for small values of J. Our finding shows that our Green's function approach is capable of going beyond mean-field theory for the prediction to the phase boundary between the insulating phases, demonstrating that the single-particle density of states still serves to unambiguously distinguish the Bose-glass and Mott-insulator phases at finite temperatures and when the tunnelling energy is made slightly positive. However, above the continuous red line, which corresponds to the result of (32), only the superfluid phase should exist. As we can see, the Mott lobes of (49) trespass this limit, allowing a direct transition between Mott insulator and superfluid phase. It has been analytically proved that such a transition should not occur in the presence of any bounded disorder [44]. Our result may be explained by the fact that we have considered only the first relevant hopping correction to the single-particle density of states. We expect that by considering higher order approximations to this quantity, the triple critical points, represented by blue dots located at both sides of each Mott lobe, would join below the red curve, such that the Bose glass would always intervene between the Mott-insulator and superfluid phases, never allowing a direct transition between the latter two. Further analysis of the energy behaviour inside that region could also lead to a determination of an accurate phase boundary. An effective-action approach, similar to the one of [26,27,28], should also be of use in this respect, where an Edwards-Anderson like order parameter, such as the one suggested in [45] and used in [22,46,47,48], could be applied to identify the Bose-glass phase, allowing also to obtain information on the nature its collective excitations. The investigation of these questions characterizes a natural progression of this work.
We now turn our attention to the comparison between our results and numerical data from the literature.

Comparison with numerical results
By considering the results of (49) only below the curve of (32), we can locate the region corresponding to any Mott lobe in the phase diagram. This makes a direct comparison between our results and numerical data from the literature possible. To this end, we first consider our prediction of the phase boundary in the zero-temperature limit, which corresponds to β → ∞. This results in for the phase boundary with the superfluid phase, and for the Mott insulator to Bose glass phase boundary, where we have defined and n 0 minimizes the energy inside each Mott lobe. Using the above expressions, we are now able to compare our result with the numerical data of [22], where an Edwards-Anderson order parameter was used in order to characterize the Bose-glass phase for a 2D square lattice at zero temperature. For the 3D case, we use the data from [20,21] for zero and finite temperatures, respectively. Such a comparison can be observed in figure 4.
As is demonstrated, our results compare quite well with the numerical data from the literature. In figure 4(a), the numerical points lie mostly inside the red continuous line indicating that our prediction overestimates the first Mott lobe in the 2D case. However, in 4(b) and 4(c), we see a remarkable agreement of our prediction with numerical data for the 3D phase boundary for both zero and finite temperatures. We notice in figure 4(c)  (50) and (51), the results of our Green's function approach, while the blue dots indicate the numerical predictions of [22] for the 2D case and [20,21] for the 3D case. Figures (a) and (b) show the 2D and 3D zero-temperature phase boundary with ∆/U = 0.6 and z = 4 and z = 6, respectively. Figure (c) presents the finite temperature 3D phase boundary with ∆/U = 0.5, z = 6, and k B T /U = 0.03.
that the introduction of temperature leads to a smoother curve for the phase boundary to the superfluid phase, characterized by (32). As a consequence, in the points where this solution meets the Mott insulator to Bose glass phase boundary of (49), kinks emerge, locating exactly the tri-critical points where the three phases coexist. Despite the discrepancy in the comparison for the 2D case, the relative deviation between our prediction and the numerical data for the tip of the Mott lobe in the three cases represents an error of less than 2%. Therefore, we can conclude that our results are in significant accordance with the numerical calculations.

Summary and Conclusions
We have presented an analytic approach for perturbatively calculating the finitetemperature Green's function of the disordered Bose-Hubbard model as well as the analysis of this quantity in order to provide a distinction for the three possible ground states of the system. By summing up a subset of the contributions in the hopping expansion of the 2-point Green's function we were able to reproduce the results of mean-field theory for the phase boundary between the superfluid and insulating phases obtained in [23]. A renormalization method was employed allowing us to compute the first relevant correction to single-particle density of states due to the hopping of particles, which made it possible to construct a phase diagram with unambiguous distinction between Mott insulator and Bose-glass, thus confirming that this quantity still serves to differentiate both phases even for slightly positive values of the kinetic energy. Our results compared well with the numerical results of [22] for the 2D zero-temperature phase boundary and show noticeable agreement with the 3D numerical results of [20,21] both for zero and finite temperatures, providing an error of less than 2% for the first Mott lobe tip. However, our result for the phase diagram predicts the possibility of a direct transition between Mott-insulator and superfluid phases, which should not occur in the presence of any bounded disorder, as was proven in [44]. Despite this limitation, it was made clear that this approach is able of going beyond mean-field theory for the distinction of the two insulating phases and we expect that the inclusion of higher order corrections in the hopping expansion of the density of states should lead to a more accurate prediction to the superfluid to Bose-glass phase transition. Further studies regarding these questions should be the subject of future investigation.