Superfluid clusters, percolation and phase transitions in the disordered, two dimensional Bose-Hubbard model

The Bose glass (BG) phase is the Griffiths region of the disordered Bose Hubbard model (BHM), characterized by finite, quasi-superfluid clusters within a Mott insulating background. We propose to utilize this characterization to identify the complete zero-temperature phase diagram of the disordered BHM in $d\ge2$ dimensions by analyzing the geometric properties of what we call superfluid (SF) clusters, which are defined to be clusters of sites with non-integer expectation values for the local boson occupation number. The Mott insulator (MI) phase then is the region in the phase diagram where no SF clusters exist, and the SF phase the region, where SF clusters percolate - the BG phase is in between: SF clusters exist, but do not percolate. This definition is particularly useful in the context of local mean field (LMF, or Gutzwiller-Ansatz) calculations, where we show that an identification of the phases on the basis of global quantities like the averaged SF order parameter and the compressibility are misleading. We apply the SF cluster analysis to the LMF ground states of the two dimensional disordered BHM to produce its phase diagram and find a) an excellent agreement with the phase diagram predicted on the basis of quantum Monte Carlo simulations for the commensurate density $n=1$, and b) large differences to stochastic mean field and other mean field predictions for fixed disorder strength. The relation of the percolation transition of the SF clusters with the onset of non-vanishing SF stiffness indicating the BG to SF transition is discussed.


Introduction
The experimental proof of the Mott insulator (MI) to superfluid (SF) transition in ultracold atomic systems [1] opened a wide field of interesting research in this field. In particular the influence of disorder on a system of bosons in a regular (e.g. optical) lattice received much interest since the fundamental work of Fisher et al. [2]. Here the phase diagram and transitions for bosons in a disordered potential was analysed and the existence of a Bose glass (BG) phase was predicted. The BG represents a non-SF but, in contrast to the MI, compressible phase displaying an excitation spectrum with arbitrarily small excitation energies. The BG phase is the analogue of the Griffith regions occurring for instance in disordered magnets, whose physics is dominated by rare region effects due to arbitrarily large strongly coupled clusters [3,4,5].
Studies of the excitation spectrum of the disordered system in dependence of the disorder strength and time-of-flight measurements confirmed the predicted BG phase experimentally [6]. In addition to the well controllable optical lattice, disorder is introduced either by a non-commensurate periodic potential [6] or by speckle potentials [7,8]. A new view on the properties of ultra cold bosonic gases opened up as high resolution techniques allowed access to single-site detection recently [9,10]. This progress now yields a direct view on the population numbers within the different phases and the in-situ hopping dynamics of the bosons in their optical potential.
Theoretically the phase diagram of the disordered Bose Hubbard model (BHM) has been studied by various methods: The strong coupling expansion [11] is a perturbative method up to third order in the tunnelling rate yielding a prediction on the Mott lobes. The disordered BHM was widely studied by quantum Monte Carlo (QMC) methods in various dimension [12,13,14,15,4,16,17]. In addition, density-matrix renormalization group (DMRG) techniques were applied to 1D systems containing either quasi-periodic potentials [18,19,20] or a uniform distribution of disorder strength [21,22]. A frequently used alternative approach is the local mean field (LMF) approximation [23], which replaces the nearest neighbour hopping on the lattice by isolated bosonic degrees of freedom interacting via an effective mean field coupling with the neighbours. Based on the LMF approximation several numerical techniques, such as stochastic mean field (SMF) theory [24,25] and LMF theory [26,27,28], were proposed.
An intriguing question has been for a long time the potential existence of a direct MI-SF transition [12,13,14,15,4,16,17], which is now excluded by a rigorous theorem [13]. The occurrence of the BG phase intervening between the MI and SF is caused by Griffiths effects [29,2] due to arbitrarily large, but exponentially rare clusters of one phase within a background of another phase [30,31,3]. Since any exponentially rare event is hard to sample numerically, the existence of an intervening BG phase might have been eluded some studies, be it QMC [14,15,4,16,17], LMF [26,27,28] or SMF theory [25,24]. One might speculate that, rather than calculating spatially averaged quantities, a look at the aforementioned clusters in individual disorder realization itself would tell us more about the actual state the system is in. In this paper we propose a method to identify the different phases of the disordered BHM on the basis of geometric properties of what we call SF clusters, which are clusters of sites with non-integer boson occupation number. The MI phase then is the region in the phase diagram where no SF clusters exist, and the SF phase the region, where SF clusters percolate -the BG phase is in between: SF clusters exist, but do not percolate. We apply this criterion to results of LMF calculations and compare it with predictions of other methods: on one side SMF theory, where the individual phases are identified on the basis of spatially averaged LMF quantities like SF order parameter or compressibility, and on the other side QMC simulations, which are supposed to be exact up to statistical and extrapolation (L → ∞, T → 0) errors (which can, of course, be quite large). Our aim is to demonstrate that the use of averaged quantities in LMF theory leads to incorrect predictions and that the cluster analysis predicts the phase diagram in d ≥ 2 in very good agreement with the exact QMC results, even when applied to LMF data.
The paper is organized as follows: In section 2 we recapitulate the LMF and SMF theory for the disordered BHM. In section 3 we critically examine the use of the averaged SF order parameter ψ and the compressibility κ as indicators of the different phases of the disordered BHM in LMF and SMF theory and then introduce our new method to construct the phase diagram on the basis of an analysis of SF clusters. We apply this cluster analysis in section 4 to the 2d disordered BHM with commensurate filling and with fixed disorder strength and compare it with prediction from QMC simulations and from SMF theory. Section 5 concludes with a discussion of the equivalence of the predictions of the SF cluster analysis for the phase boundaries with the conventional definition of the MI-BG and BG-SF transition points.

The model
Ultracold bosonic atoms in an two dimensional square optical lattice can be described by the BHMĤ where i = 1, . . . , M is the site index, M = L 2 the number of sites and L the lateral size of the square lattice. The chemical potential is described by µ, the inter particle repulsion by U and the tunnelling rate by J. The last sum runs over all nearest neighbour pairs (ij) of the underlying lattice. The operatorn i =â † iâi is the particle number operator of bosons on site i, which are annihilated and created by the operatorsâ i andâ † i . Moreover, on-site disorder is introduced by the parameter ǫ i , which is drawn randomly from a box distribution p(ǫ) = Θ (∆/2 − |ǫ|) /∆, where ∆ is the strength of the disorder.
In the SF regime tunnelling dominates the system. Thus, the ground state is a coherent state, which is an eigenstate of the tunnelling part of the Hamiltonian. The SF parameter ψ i = â i is the expectation value ofâ i evaluated in the ground state for T = 0, which is non-zero for a coherent state. An eigenstate of the diagonal part of the Hamiltonian on the other hand is a Fock-state, which is the ground state of the system for small tunnelling rates in the MI regime. The expectation value ofâ i in a Fock-state is zero in any case. Hence, the mean value of the SF parameters ψ = i ψ i /M is an order parameter for the SF phase.
The phase diagram of the pure system in dependence of JZ/U and µ/U consists of so called Mott lobes, in which the system is a MI with a fixed integer number of atoms per site. While in the MI regime the state is localized in real space, in the surrounding SF regime it is localized in k-space. When disorder is introduced, the Mott lobes shrink and a new phase, the BG phase, occurs. In order to distinguish all three phases the compressibility κ = n 2 − n 2 is necessary. Among the three phases only the MI is non-compressible and only in the SF phase the SF order parameter is nonzero. Consequently, the phase of the system can be identified by the SF order parameter ψ = i ψ i /M and the compressibility κ. While the SF order parameter is a measure for the coherence in the system, the compressibility describes the variance of the particle number per site.
In order to determine the phase diagram of the BHM, the ground state properties of the LMF Hamiltonian can be studied via SMF theory [24,25], which computes the PD self-consistently, or via LMF theory [23,28,27], which solves the coupled set of equations for the local SF parameter directly on the lattice. We will first describe the approximations made in LMF theory, followed by a discussion of the additional assumptions made in the SMF approach.

Local mean field theory
In LMF theory the tunnelling part of the Hamiltonian can be approximated viâ where terms of the form (â i − â i )(â † j − â † j ) are neglected [23]. The central quantities are the local SF order parameters which are defined as the expectation values of the annihilation operator at the individual site i in the ground state of the system. Because of the U(1)-symmetry they can be chosen to be positive and real, which leads tô Thus, the Hamiltonian can be decomposed in a sum of diagonal operators, whose tunnelling rate is replaced by an effective local rate Jη i , which depends on the local SF parameter of the neighbouring sites η i := j A ij ψ j , with A ij = 1 for i and j nearest neighbours on the square lattice with periodic boundary conditions and zero otherwise. This approximation reduces the full quantum problem to M quantum sites, which are coupled in a mean field way with a spatially varying coupling rate. In order to compute the phase diagram in LMF theory the coupled set of the selfconsistency equations is solved on a L×L lattice, where the expectation value is evaluated in the ground state of H i that itself depends an the local SF parameter ψ i . As a result of the decomposition (5) all states considered (in particular the ground state) are a direct product of individual single-site states. This means in particular that they are Gutzwiller states of the form with single-site states |φ i = ∞ n=0 c i n |n i given in particle number basis and |c i n | 2 describes the probability to find n bosons at site i and fulfils M i=1 |c i n | 2 = 1. In particular, the local SF order parameter is then given by ψ i = ∞ n=0 c i n−1 * c i n √ n i and the local boson number is represented by n i = ∞ n=0 |c i n | 2 n i . For the numerical implementation, starting from a random initial configuration for ψ i on the 2d lattice, the set (6) of equations is solved recursively. This involves solving the eigenvalue problem on each site and computing the expectation value of the annihilator in the numerically determined ground state. This is repeated until the averaged SF order parameter is determined with an accuracy of 10 −4 . In the disordered case the results are averaged over 200 different realizations of disorder, indicated by the brackets [. . .] av . Since we are working in a regime in which the maximum average particle number per site is three, it was numerically checked that it is sufficient to truncate the basis of the Hilbert space for each site at N = 10. With the solution found for the local SF parameters ψ i on the lattice, the ground state of Hamiltonian (5) is calculated numerically. Afterwards all desired expectation values and finally the compressibility withN = ini can directly be computed. The probability distribution (PD) is determined on the basis of the complete set of values of ψ i , additionally averaged over disorder realizations.

Stochastic mean field theory
The central idea of SMF theory is to circumvent the computation of all local order parameters ψ i by deriving a self-consistency equation for the probability distribution P (ψ) directly [24,25]. Additional approximations are of course necessary. The SF order parameter ψ = gs|â|gs , which is derived from the ground state of the full quantum Hamiltonian, can be determined by the ground state of the single-site Hamiltonian (5) in dependence of the stochastic variables ǫ and η. The parameter ǫ is then a stochastic variable drawn from the disorder distribution p (ǫ) and as a result ψ = gs|â|gs is a stochastic variable, drawn from the PD P (ψ), which must be determined selfconsistently. Since η is the sum of the SF parameters of the neighbouring sites it also is a stochastic variable drawn from Q (η). The problem of computing the ground state of the full quantum system for all lattice sites simultaneously is thereby replaced by analysing the ground state |gs (ǫ, η) of the site independent Hamiltonian as a function of ǫ and η. Thus, the probability distribution (PD) depends on the distribution Q (η) of the occurring values of η and the distributionP η (ψ) of the local SF parameters for given η. A direct analysis of gs (ǫ, η) |â|gs (ǫ, η) as a function of ǫ and η yields P η (ψ) = d dψ dǫp (ǫ) Θ (ψ − gs(ǫ, η)|â|gs(ǫ, η) ) .
Since η is the sum of the local SF parameters ψ of the neighbouring sites its distribution is given by where P Z (ψ 1 , . . . , ψ Z ) is the connected probability distribution function of the local order parameters ψ 1 , . . . , ψ Z of the Z neighbour sites of a single-site. Assuming that these Z local SF parameters are statistically independent equation (14) transforms into a convolution Since the assumption (15) implies the absence of correlations of the local SF parameters, one expects that it is not justified close to the phase boundaries, where the correlation length even diverges, when the transition is 2nd order. We examine the validity of this approximation in dependence of the system parameter JZ/U and µ/U in Appendix A. After determining the PD P (ψ) the SF order parameter ψ = dψP (ψ) ψ is given by the mean value of the distribution. The compressibility κ = [ N2 − N 2 ] av witĥ N = ini is computed. With these quantities at hand one can, on the basis of the underlying approximations, estimate the phase boundaries of the transitions between MI, BG and SF. They are shown in Figure 4 and 5 and will be discussed in the following.

Criterion for the phase transition
In this section we will first discuss the well known SF order parameter ψ and the compressibility κ, which are expected to indicate the phase transition: The ground state in the MI regime is a Fock-state, which is incompressible (κ = 0) and non-coherent ψ = 0, while conversely in the SF regime it is described by a coherent state ψ = 0, which is compressible (κ > 0). If disorder is introduced, those phases are separated by the BG phase, which is compressible (κ > 0), but not coherent ψ = 0. We will see that a precise prediction of the transition point on the basis of ψ or κ is not possible in LMF theory. Instead we will introduce an identification criterion of the different phases on the basis of the complete set of local occupation numbers.

SF order parameter and compressibility
In the ordered case the on-site energies ǫ i are zero and the lattice is homogeneous. The SF order parameter ψ clearly marks the location of the MI to SF phase transition as shown in Figure 1 on the left hand side for µ/U = 1.05 and 0.32. While the SF order parameter is zero for small tunnelling rates in the MI phase, it becomes non-zero and positive above a critical value of the tunnelling rate in the SF phase. The compressibility κ shows the same behaviour at the phase transition as the SF order parameter for both methods in the ordered case. Moreover, we analysed different lattice sizes L and the LMF results show no visible finite-size effects. In this way the phase transition in the ordered case can be determined very precisely, both within SMF and LMF theory. The resulting phase transitions agree perfectly with the perturbation predictions [32] where n denotes the mean number of particles per site and simultaneously counts the number of lobes. The calculation in [11] predicts that in the disordered case the upper (lower) part of the Mott lobes are shifted downwards (upwards) by ∆/2 but the shape remains unchanged.  Figure 1: Comparison of the LMF and SMF predictions for the average SF order parameter ψ and the compressibility κ for fixed chemical potential µ and varying tunnelling rate J. Left: Homogeneous case (∆ = 0), Right: disordered case with ∆/U = 0.6. Top row is for µ = 1.05, where the ordered system displays a MI SF transition and the disordered system a BG SF transition (κ > 0 for all values of J). Bottom row is for µ = 0.32, where the ordered system again displays a MI SF transition and the disordered system is expected to display MI, BG and SF phases (see section 3.2). For LMF theory the results for a 2d lattice with L = 100 (line), L = 50 (•), L = 10 (+) are depicted, which shows that finite-size effects can be neglected.
The situation for the disordered case is shown in Figure 1 on the right hand side. The SF order parameter is shown for µ/U = 1.05 and 0.32 as a function of the tunnelling rate for a disorder of ∆/U = 0.6. Whereas the SMF theory predicts a direct BG SF transition, at a critical value JZ/U ≈ 0.0241 and 1.0455, see 1 (b, d), above which ψ become non-zero the behaviour of ψ as predicted by LMF theory does not indicate a transition at all; it varies smoothly with the tunnelling rate J. This is not a finitesize effect as we have checked by examining different lattice sites, as shown in 1. The compressibility, which indicates the MI BG transition, displays the same behaviour.
It turns out that the reason for the failure of the average SF order parameter to predict the location of the BG SF boundary is the following: In the disordered case the value of the local SF parameter varies substantially from site to site due to the variation of the local potential of ǫ i . Close to the phase transition there are sites with zero local SF parameter and others, where the local SF parameter is still positive. This has been interpreted as an overestimation of the phase coherence in LMF description [24]. Our interpretation, however, is different: It is only the average SF order parameter ψ that overestimates the phase coherence. A closer look at the complete probability distribution P (ψ) of the local SF order parameters and their geometrical features provides an estimate of the SF regions in the phase diagram. Its prospects are discussed in the 4.3. In the next section we discuss, how a deeper understanding of the mechanisms driving the phase transitions and their location in the phase diagram can be obtained by studying the geometric characteristics of the spatial inhomogeneities of the local SF parameters ψ i and particle number per site n i .

Identification of phases via local boson occupation number
The MI and SF phases can be discriminated by the boson number statistics at individual sites, as has also been demonstrated experimentally in [1]. The ground state in the extreme MI limit (J → 0) is a Fock states with a definite number of particles n at each site. In the extreme SF limit (U → 0), the ground state is a coherent state , in which the local boson number distribution is close to a Poissonian. Although, in the regime between these two extreme limits the ground state wave function can no longer be written as simple product states still the MI phase is characterized by a sharp, integer boson number per site and the SF phase by a fluctuating boson number per site, i.e. a non-vanishing variance of the boson number distribution p i n = |c i n | 2 (c.f. the expansion coefficient in the Gutzwiller wave function (7)). In other words in the MI regime the expectation value of the number of bosons per site n i is an integer, whereas in the SF regime it is non-integer.
Whereas in the ground state of the homogeneous BHM either all sites have an integer boson number (MI regime) or all sites have a non-integer boson number (SF regime) this is different in the disordered BHM. In particular, outside of the MI regime one expects to encounter spatially inhomogeneous situations, in which some sites have a sharp (integer) boson occupation numbers and others have fluctuating (non-integer) boson numbers. Introducing phase operators Φ that is canonically conjugate to the boson number operatorsn i the BHM can be mapped upon a Josephson junction array or more general to a quantum rotor model [2], in which superfluidity is indicated by long range order in these phases (d ≥ 2). Because of the Heisenberg uncertainty relation sites with sharp phases correspond to sites with fluctuating boson numbers, and connected clusters of sites with fluctuating boson numbers tend to have, roughly speaking, all the same phase. These clusters can therefore be identified with SF regions, although, true superfluidity only exists in the infinite system. Indeed, once these phase ordered clusters percolate, true superfluidity emerges, signified by a non-vanishing SF stiffness, which is the extra free energy cost to impose a uniform twist on the phases. Since such a uniform twist can be introduced by applying a certain twist at the boundary phases in one space direction, it is clear that the SF stiffness is zero as long as the clusters do not percolate: In the absence of long range order in the phases such a twist at the boundary does not cost energy.
On the basis of this qualitative picture we hypothesize that the BG to SF transition in the d-dimensional BHM (d ≥ 2) coincides with the percolation transition of connected clusters of lattice sites with a non-integer boson number expectation value n i . We expect this coincidence to hold as long as the SF phase displays true long range order, characterized within the phase description by a non-vanishing long distance limit of phase correlations -which means it should hold for d ≥ 2. The BG SF transition in one dimensional BHM might not be related to a percolation transition, since in d = 1 the SF phase has only quasi long range order (algebraically decaying correlations). We should note that the relation between the BG SF transition and percolation has already been pointed in [33,26,34], but has neither been used in a quantitative manner to determine phase boundaries nor checked against, for instance, Monte Carlo results.
In the following we denote the sites with non-integer boson occupation number n i as SF sites, and sites with integer n i as MI sites. Analogously, we discriminate SF clusters and MI clusters. Formally we map the boson occupation numbers to a discrete field S i that is set to S i = 1 for SF sites and S i = 0 for MI sites. Then we identify the different phases of the disordered BHM as follows: MI phase: S i = 0 for all sites i. All boson occupation numbers are integer (and identical), consequently the compressibility κ is zero.
SF phase: S i = 1 for a macroscopic fraction of sites, which form a percolating connected cluster. According to what we discussed above the percolating cluster has phase long range order and thus yields a non-vanishing SF stiffness (which is proportional to the SF density). BG phase: Characterized by a non-vanishing density of sites with S i = 1, none of the connected clusters formed by the SF sites percolates. The BG phase is thus characterized by isolated SF clusters within a MI sea. The phases of the isolated clusters are uncorrelated, hence phase long range order is lacking and the SF density vanishes (no SF order). Moreover, due to the number fluctuations on the SF sites the BG phase is compressible (κ > 0).
Within LMF theory the expectation values of the local boson occupation numbers are straightforward to calculate via n i = gs|n i |gs , where|gs is the ground state of the LMF Hamiltonian (5). For numerical reasons we introduce a threshold γ n into the definition of the discrete field where γ n = 5 · 10 −3 is chosen to serve as the cut-off in this algorithm. In the whole parameter range, where sites with integer particle number occur, the histogram of the mean particle number n i has narrow peaks of width γ n at integer values. The width decreases when we increase the number of iteration steps to solve the selfconsistency equations 6. The threshold parameter γ n introduced to identify MI sites (and complementarity SF sites) can be reduced by increasing the numerical effort without changing the final results.
In Figure 2 typical results for one realization of disorder for ∆/U = 0.6 and µ/U = 1.0455 are shown for three different values of the tunnelling rate JZ/U for the three phases . In the first row the local SF parameter ψ i , in the second the particle number per site n i and in the third the resulting discrete map S i is shown. In the MI regime all sites are occupied by the same integer number of particles (in this case one, since we are in the first Mott lobe). At the transition from the MI to the BG regime SF sites (S i = 1) with non-integer particle number occur. Because of these locally occurring SF sites ψ i > 0 the SF order parameter ψ is small but not zero in this regime. Since the SF islands are compressible, this phase has positive compressibility. In the BG phase the SF islands does not percolate, yet. They grow in number and size, until one of them finally percolates. The percolation represents the actual transition to the SF regime in parameter space. Just after the percolation the phase in the system is coherence macroscopically, which means that all local SF parameter ψ i are positive and compressible, as described above.  In this section we demonstrate how we determine numerically the percolation transition via a cluster analysis of the discrete map S i and finite size scaling [35]. Assume we study the phase diagram in dependence of the system parameter called x and y. Than, the percolation probability p Perc , i. e. the probability of having a percolating cluster, is given for fixed y as a function of x and will be determined for different system sizes L. The percolation probability is expected to obey the finite-size scaling form where x c is the percolation threshold, i.e. the value above which a percolating cluster exists with probability one, and ν the critical exponent determining the divergence of the mean lateral cluster size at the transition. The scaling functionp(X) approaches zero for X ≪ 1 and one for X ≫ 1, which means that exactly at the transition x c the curves for different system size should intersect (in the scaling limit). This intersection point, which we can easily be identified with the system sizes behaviour at hand, is thus a reliable indicator for the percolation transition. The cluster analysis of the discrete map S i is done for every disorder realization. Afterwards the results are averaged over 200 (L = 50, 100) and 2500 (L = 10) realizations of disorder. The percolation probability p Perc for this case is shown in Figure 3 (blue for MI and red for SF sites) for different system sizes as a function of the tunnelling rate JZ/U. Moreover, the finite size scaling analysis for the percolation transition at µ/U = 0.439 and µ/U = 1 for ∆/U = 0.6, yielding J c Z/U = 0.15, J c Z/U = 0.04 respectively and the critical exponent ν = 1.33 in both cases, is depicted. Thus, this transition is in the universality class of conventional 2d percolation [35]. We find the same universality class of the percolation transition for all parameter values that we studied.

Commensurate filling -Comparison with QMC results
In this section we determine the complete phase diagram for commensurate density in dependence of ∆/2J and U/J, for which a prediction on the basis of QMC simulations is available [12]. We fix the particle density to n = M i=1ni /M = 1 with an accuracy of 10 −4 by adjusting the chemical potential for each point (∆/2J, U/J) in the phase diagram that we study. Outside of the Mott lobes this result is unique, whereas in the MI regime the chemical potential is fixed to the middle of the MI gap. In the µ/∆ versus JZ/U representation, where the Mott lobes are visible and which we will discuss in section 4.2, this n = 1 line always passes the tip of the first Mott lobe. In the ∆/2J versus U/J parameter space the corresponding line for fixed ∆ is a straight line through the origin with slope ∆/2U.
With the chemical potential that fixes the density n to one we compute the ground state of the LMF Hamiltonian (5) and determine the discretized boson number field S i (18), which we use to identify MI, BG and SF phase. The resulting phase diagram is shown in Figure 4 on the left. As expected [13] the SF region is completely surrounded by the BG phase (except at ∆ = 0). Its boundary has some characteristic features: It extends in a slight bump up to quite large disorder strength (up to ∆/2J ∼ 75) and in a pronounced nose up to the interaction strength U/J ∼ 52. This nose gives rise to a re-entrant behaviour: Moving vertically from a point within the MI phase, which has long range positional order, one enters first the BG phase, which is disordered and then, upon further increasing the disorder strength, enters the SF phase, which has off-diagonal long range order. Weak disorder thus supports superfluidity in the BHM, as has been observed before [17,12,36,37]. Remarkably, our prediction on the basis of a cluster analysis of LMF ground states agrees very well with the results of QMC simulations [12] shown for comparison in Figure  4 on the right. The shape of the SF-BG phase boundary with its characteristic nose and bumps clearly coincide. The quantitative agreement is very good, too, regarding the substantial error bars of the QMC data in the large disorder and large interaction regime (the QMC estimate for the extreme value of ∆ in the bump is ∆/2J ∼ 72 ± 4 and of inter particle interaction in the nose U/J = 49 ± 3, c.f Fig. 2 in [12]). Moreover, with our method we could also explore the weak interaction region, which is hardly accessible by QMC methods, and found a singular behaviour of ∆ with U → 0, which is compatible with the analytically predicted behaviour ∆ ∝ √ U [38]. We conclude that the percolation criterion that we introduced in section 3.2 to locate the SF-BG transition produces remarkably accurate predictions even in LMF theory.
Our result for the MI-BG transition line, which denotes the appearance of noninteger boson occupation numbers and thus SF sites, agrees well with the perturbative result 17, shown in Figure 4 on the left. Moreover, they agree with the line ∆ = E g/2 obtained using the gap data from [39], shown in Figure 4 on the right.
In passing, we note that for weak disorder the MI clusters percolate close to the BG-SF transition line, whereas for stronger disorder they percolate deeper inside the BG phase. Whereas for weak disorder the individual sites of a MI cluster in the BG phase all have the same integer occupation number, this is in general not the case any more for strong disorder: the integer occupation number of MI clusters can vary from site to site. Finally, we note that SMF theory as described in section 3.1 predicts a direct MI-SF transition along the lower border of the SF region in the parameter range shown in Figure 4. The characteristic BG region for small disorder strength is absent in this parameter range, which is in contradiction to the theorems proven in [13], which exclude a direct MI-SF transition in any disordered systems. Besides, we checked that the BG phase occurs for even higher values of U/J, which is in agreement with results to be presented in the next section.

Fixed Disorder -Comparison with SMF theory
After we have seen in the last section that our method to determine the phase diagram of the 2d disordered BHM leads to results that agree very well with QMC predictions, we determine in this section the (µ/U-JZ/U) phase diagram for a fixed disorder strength ∆/J = 0.6 and compare it with predictions of SMF theory. In this phase diagram the Mott lobes occur and the line given by n = 1 always passes the tip of the first one.  (17). Right: SMF phase diagram determined by using the SF order parameter ψ and the compressibility κ [24,25]. The red line indicates the critical tunnelling rate J where the SF order parameter ψ becomes non-zero, the blue line the critical tunnelling rate, where the compressibility κ becomes non-zero.
the additional approximation (15) on the distribution P Z (ψ 1 , . . . , ψ Z ). The validity of this restriction fails close to the phase transitions as we show in Appendix A. Despite or perhaps because of this approximation the SF order parameter ψ as well as the compressibility κ computed within SMF theory are exactly zero in specific regions of the parameter space (c.f. Figure 1), which one might want to identify with MI and BG phase, as was done in [24,25], c.f. section 3.1. The LMF cluster analysis and SMF [24,25] phase diagrams for fixed disorder strength ∆/U = 0.6 are shown in Figure  5. One immediately observes substantial differences: Firstly, in LMF theory the BG phase always separates the MI from the SF phase. The intervening BG phase is actually predicted to be quite large even at the tip of the Mott lobes, not just a "thin sliver" [2]. SMF theory, however, predicts a direct MI-SF transition, in contradiction to [13]. Secondly, large differences in the critical tunnelling rate for the BG-SF transition occur especially in the region around µ = 1. Assume we fix the chemical potential there. In this case the SMF theory predicts the phase transition at JZ/U = 0.0241. The percolation of the SF cluster, however, takes place at JZ/U = 0.0430. Thus, significant changes of the system in this case occur for values of the tunnelling rate twice as large as predicted by ψ in SMF theory. A direct comparison of our results with the QMC data of [12] is not possible here, since the latter are obtained for the canonical ensemble, where the chemical potential is absent. However, we can take our LMF estimate for the value of µ that fixes the particle density at n = 1 for fixed U/J and ∆/2J to obtain an approximate comparison -see Table 1, where we also show the prediction of the strong coupling expansion [11] for the MI-BG transition for ∆ = 0.6. One observes deviations of the QMC and strong coupling predictions from our LMF cluster analysis results at the tip of the Mott lobe, but a good agreement for stronger disorder, ∆/U = 2. At this disorder strength also a QMC prediction for the quantum rotor model exists [37], which differs by 25% from the predictions for the BHM, ours and the one from QMC. We also note that the tapered shape of the Mott lobe predicted by the strong coupling expansion [11] agrees well with our result of the LMF cluster analysis shown in Figure 5.
In addition to our LMF cluster analysis and the SMF theory discussed above a number of other approximative methods have been applied to calculate the phase diagram of the BHM at fixed disorder. In the zero temperature mean field phase diagram of [40] the BG phase is completely absent, which might be true in infinite dimensions, but certainly not in d = 2 or 3.
A LMF theory has been used in [27,26] to solve the self consistency equations (3) and to calculate a LMF expression for the stiffness or SF fraction and the compressibility. Using these two observables the (µ/U-JZ/U)-phase diagram is then determined, which displayed a round shape of the Mott lobes, a direct MI-SF transition for small disorder and an intervening BG phase at larger disorder. It should be noted that although the starting point of the calculation in [27,26], the LMF approximation, is identical to ours, the usage of a different criterion to identify the phases leads to a phase diagram that differs significantly from the one predicted by us.
A multi-site LMF theory is introduced in [41], where every plaque of two by two sites is treated quantum, which keeps the spatial correlation therein. Instead of single sites these plaques are coupled in a LMF way (analogous to section 2.1). The Mott lobe is determined for both, the single-site and multi-site LMF theory, on the basis of the condensate fraction. In agreement with [27] and SMF theory it shows a round shape at the tip. The multi-site LMF theory predicts a larger MI region than the single-site LMF theory. Note that the condensate fraction smoothly approaches zero, analogous to our observations on the SF order parameter and the compressibility made in section 3.1; a linear fit is used to determine the transition point.
In [42] the so-called the Gutzwiller projected variational techniques is introduced in order to determine a canonical transformation of the quantum Hamiltonian, which requires the truncation of the hopping term. Thus, it is possible to minimize the expectation value of the transformed Hamiltonian in Gutzwiller type local mean field states with respect to its variational parameters. Finally, the SF stiffness and the compressibility yield the phase diagram, which shows a remarkably narrow BG region between the MI and the SF phase.
In all mean field calculations mentioned the tip of the Mott lobe is predicted for far higher values of the tunnelling rate as our results based on the LMF cluster analysis, and results of QMC or strong coupling expansion methods, as listed in Table 1.

The probability distribution of the local SF parameter
In this section we discuss the probability distribution (PD) of the local SF order parameter P (ψ). In the ordered case ∆ = 0, depicted in the first row, all local SF parameter are identical since all sites have the same on-site energy. The averaged order parameter ψ, depicted as a blue cross, is identical to each local SF parameter ψ i and the variance of this value is zero. Therefore, the PD P (ψ) is a delta function at the values of ψ. Within the MI region the local SF parameter is zero everywhere and thus the PD is a delta function at ψ = 0. In the SF regime still the PD is a delta function but at positive ψ, which increases with the tunnelling strength.  Figure 6: Probability distribution P (ψ) of the local SF parameter at fixed chemical potential (µ/U = 1.05) for different tunnelling rates J, (top) in the ordered case (∆ = 0) and (bottom) in the disordered case (∆/U = 0.6) as predicted by the LMF (blue) and SMF (red) theory. The crosses (×) at the ψ-axis represent the mean of the PD, which is the average SF order parameter ψ. In the ordered case the PD is a delta function. With disorder all local SF parameters are zero for the MI; the PD in the BG phase is a superposition of a delta function at ψ = 0 and a SF tail; the SF phase is characterized by a broad distribution of positive non-zero local SF parameters.
This situation changes if disorder is introduced, since then the on-site energy is different on every site resulting in a variety of different values of the local SF parameter ψ i . In the MI regime the PD is a sharp delta function at ψ = 0 still and becomes a broad distribution in the BG and SF phase. In the BG phase sites with zero local SF parameter (corresponding to MI sites with S i = 0) coexist with sites, which have non-zero local SF parameter (corresponding to SF sites with S i = 1) and where called SF islands before. In the SF regime the PD is a broad distribution representing the variety of positive values for the local SF parameter. Due to this characteristic behaviour the PD can be written as a superposition of a delta function at ψ = 0 and a broad distribution representing the values ψ > 0: We denote this distribution P SF (ψ), since it represents sites, which we refereed to as SF sites (S i = 1) before: Starting with the same Hamiltonian (11) as LMF, SMF theory introduces the additional approximation (15) yielding in a self-consistent equation for the PD itself, which is given by equation (12). The resulting distributions are shown in Figure 6 in red. Whereas no deviations between LMF and SMF theory occur in the ordered case (∆ = 0), they become visible in the disordered case. For ∆ > 0 oth distribution have similar shapea, but especially for small values of ψ, which are crucial for the determination of the phase transition, they differ significantly in the BG and SF phase.
In order to identify the reason for these deviations we scrutinized the validity of the additional SMF restriction (15). We checked its validity by comparing the LMF results for the product of the PD of two different sites P (ψ i ) P (ψ j ) with the pair PD P Z (ψ i , ψ j ) and determined their deviation ∆ P in Appendix A. As shown in Figure A1 the approximation (15) is best in the MI regime and the BG for very small tunnelling rates. But for increasing JZ/U is becomes worse and especially at the phase boundary it fails. In Figure A2 where P (ψ i ) P (ψ j ) are P Z (ψ i , ψ j ) are shown for different parameters it is visible that they disagree especially for small values of the ψ. This disagreement is due to the presence of correlations of the local SF parameters ψ i at different sites in the vicinity of the transition points, which are neglected in SMF theory.

Conclusion
In this paper we have introduced a new criterion to identify the different phases of the disordered BHM in d ≥ 2 on the basis of the complete set of local boson occupation numbers {n i } of each sample and applied it to the ground states calculated using the LMF approximation. In the MI phase all n i are integer, in the BG phase some of them are non-integer and form SF clusters in a MI background and in the SF phase at least one of these clusters percolates. The emergence of SF clusters, with an average lateral size that is expected to be of the order of the SF correlation length, have a finite density, which gives rise to a non-vanishing mean of the average SF parameters although the system is not SF. The latter happens only when these SF clusters percolate, which is the hallmark of the BG-SF transition. Moreover, the SF clusters have a fluctuating boson occupation number resulting in a small but non-vanishing compressibility. Thus their appearance is the indicator of the MI-BG transition, i.e. from the incompressible (κ = 0) to the compressible (κ > 0) phase. Consequently, the BG phase displays arbitrarily small but non-vanishing values for ψ and κ and all approaches to determine the LMF phase boundaries of the disordered system on the basis of the site and disorder averaged parameters, like the SF order parameter ψ, the compressibility κ, the SF or condensate fraction, overestimates the SF and MI phases substantially and are doomed to fail: the putative phase boundaries move systematically and substantially when increasing or decreasing the threshold only by a small amount.
The resulting cluster analysis phase diagram for a fixed commensurate density n = 1 is in excellent agreement with the prediction of QMC simulations, not only qualitatively in reproducing the characteristic shape of the SF region in the (∆-U)-diagram, but also quantitatively within the numerical error bars. This is remarkably, since other LMF approaches using averaged quantities, like the mean SF order parameter or the compressibility as indicators predict much larger MI regions in the phase diagram or even fail to identify the BG transition, since the used indicator varies smoothly at the expected phase transition. Small deviations between QMC calculations and our LMF cluster analysis might be due to the fact that the local occupation numbers calculated by using the LMF approximation deviate in some regions of the phase diagram from the exact expectation values. Obviously it would be desirable to calculate the latter by QMC simulations and to perform the cluster analysis we propose to these data.
The questions that immediately arises in this context is: 1) Is the MI-BG transition in the disordered BHM exactly where SF sites occur? And more interestingly 2) Is the BG-SF transition actually identical with the percolation transition of SF clusters?
Concerning question 1): Although not proven rigorously, the MI-BG transition is supposed to occur, where the gap E g/2 , i.e. the energy for particle-hole excitations, of the pure, ordered BHM is equal to the disorder strength ∆ [2]. It seems plausible that when this happens individual sites or small clusters will occur, where the addition or removal of a particle does not cost energy and thus the local boson occupation number fluctuates, i.e. n i is non-integer. This is how we propose to identify the MI-BG transition.
Concerning question 2) we argued in section 3.2 in the basis of the BHM to quantum rotor models that the SF stiffness will always vanish as long as SF clusters with a MI background do not percolate. This BG situation then is reminiscent of a d + 1dimensional, classical XY model with columnar disorder, in a state with (quasi) phase ordered finite clusters in a phase disordered background. Application of a phase twist at the system boundaries will only cost a macroscopic amount of energy when the ordered regions actually percolate -in the SF region. Note that this argument is based upon the existence of true long-range order in the SF phase of the pure, ordered BHM, thus we expect it to be valid for d ≥ 2.
A complementary picture is based on the path-integral computation of the SF density [43], which is used in QMC simulations to identify SF order. The SF density or stiffness is proportional to the mean-square of the winding number of boson world lines in the path integral representation. When on average a finite fraction of boson world lines wrap around the whole system, the mean-square winding number is positive and the system is SF. To wrap around the whole system (with periodic boundary conditions), a boson world line, on its way through imaginary time, has to move along a path that traverses the whole system, thus attributing particle number fluctuations to the individual sites of this path. These sites will consequently attain non-integer expectation values for the boson occupation numbers n i (since for some time the boson was there and for some time not), thus in the end there must be at least one percolating SF cluster in the system.
It should be noted that other quantum phase transitions of disordered systems are naturally percolation transitions, too: The critical point of the random transverse Ising model is governed by an infinite randomness fixed point (in d ≥ 1 dimensions [30,3,44,45]), which signals the percolation of strongly coupled clusters that away from criticality constitute the Griffiths phase. The percolation transition that we observe in our calculations falls into the universality class of a conventional, 2d site percolation -which means it does not carry the signature of the critical properties of the proper BG-SF transition. This is most probably a consequence of the LMF approximation that we use, since it does not properly account for spatial correlations -if applied to the exact ground state one would expect the critical exponents of the percolation transition to be related to the critical exponents of the BG-SF transition.
In addition to providing an intuitive picture and a deeper understanding of the underlying physics of the phase transitions in the BHM the cluster analysis may serve as a reliable tool to locate the transitions in situations, in which the application of other criteria to discriminate the different phase might lead to erroneous predictions -as for instance in LMF theories. Finally, since experiments recently reached the regime of single site detection [9,10] and are now able to observe the particle numbers at each site, an experimental application of the cluster analysis that we propose appears in reach.

Appendix A. Local SF order parameter correlations
In this appendix we check the statistical independence assumption underlying SMF theory. Additional to the LMF approximation (4) made in the tunnelling part in the Hamiltonian, SMF theory assumes that the local SF parameters ψ 1 , . . . , ψ Z of the Z neighbours of a chosen site i are uncorrelated and identically distributed which is introduced by the approximation (15). On the basis of LMF calculations we want to test this approximation by comparing P (ψ i ) P (ψ j ) and P Z (ψ i , ψ j ), of which some examples are shown in Figure A2. The function P (ψ i ) P (ψ j ) is the product of the PD, describing the distribution of the local SF parameter as discussed in section 4.3. The function P Z (ψ i , ψ j ) is the PD of pairs (ψ i , ψ j ), where i and j are neighbouring sites. It represents the probability of having a specific value for the pair (ψ i , ψ j ). Exactly as P (ψ j ) the P Z (ψ i , ψ j ) is computed for every realization of disorder and finally averaged. Both distributions should coincide if the assumption (15) is valid. In Figure A1 the integral difference of both distributions is shown in parameter space.
In the MI region, where the P (ψ) is a delta function at ψ = 0 and for very small tunnelling rate JZ/U the deviations are small, whereas they are significant in the region of the phase transition and in the SF regime. For illustration both PDs are shown in Figure A2 along a line of µ/U = 1.0455 and at the tip of the Mott lobe, where the deviation ∆ P reaches its maximal value (corresponding to the black dots in Figure  A1). Additional to the fact that all distributions are symmetric naturally, P (ψ i ) P (ψ j ) shows a rectangular symmetry, which intrinsically follows from the fact that it is a product of the same PD P (ψ). The PD P Z (ψ i ψ j ) containing further information of the occurring pairs shows systematic deviations. Whereas the values on the diagonal are reproduced quite well, the off-diagonal contributions are squeezed to the diagonal. This is especially pronounced in Figure A2 (d) and (h), which corresponds to the tip of the Mott lobe. These mean differences in comparison with P (ψ i ) P (ψ j ) can be observed for all parameters shown in Figure A2 and mainly occur in the regime of small local SF parameters. Figure A1 illustrates that the assumption (15) made in SMF theory, is well fulfilled in the MI regime but becomes worse in the region of the phase transition. Whether this theory predicts the phase transition reliably in this regime, is therefore to question.
The deviations occurring for small SF parameters in the limit of small ψ as shown in Figure 6 have also been discussed in [24]. In this work the authors concluded that LMF theory overestimates the phase coherence in the BG regime. But this is also true for SMF theory, since it is based on the same approximation of the tunnelling term of the Hamiltonian. In this paper we resolved this apparent problem by interpreting the SF-BG phase transition as a percolation transition (c.f. section 4.1 and 4.2).  Figure A1: The deviation ∆ P between P (ψ i ) P (ψ j ) and P Z (ψ i , ψ j ) is shown in dependence of the system parameter. In the MI regime and for very small tunnelling rates the deviations are small, whereas they grow at the phase transitions and in the SF regime. The black dots make the parameters used in Figure A2 along a line µ/U = 1.0455 and at the tip of the Mott lobes, where the deviation is maximal.  Figure A2: The first column shows P Z (ψ i , ψ j ) and the second P (ψ i ) P (ψ j ) in the disorder case (∆/U = 0.6). In the first row the parameters are given by JZ/U = 0.0283, µ/U = 1.0455 followed by JZ/U = 0.0586, µ/U = 1.0455 and JZ/U = 0.1414, µ/U = 1.0455 and JZ/U = 0.1434, µ/U = 0.4394 in the last row corresponding the black dots in Figure A1.
Appendix B. The characteristic shapes of the PD In section 4.3 we discussed three different shapes of the PD in the disordered case, which are depicted at the bottom of Figure 6. In the MI phase the PD is given by a delta function at ψ = 0, whereas in the BG and SF phase a broad distribution occurs, which means that P (ψ) can be represented by a superposition of a delta function at ψ = 0 and a continuous part P SF (ψ) caused by SF sites with ψ > 0 as defined in equation (20), and in the following denoted as SF distribution. Here, we will identify regions of the three different shapes of P (ψ) in parameter space and discuss their connection to the phase transitions determined in section 4.1 and 4.2.  Figure B1: Regions of the three characteristic shape of P (ψ) for fixed density n = 1 (left) and fixed disorder ∆/U = 0.6 (right). The black line represents the MI-BG transition according to the perturbative result given by (17). Inside of the blue line P (ψ) consists only of a delta function at ψ = 0. Within the region bounded by the red line on left and to the right of the red line on the right P (ψ) only has a continuous part P SF (ψ). In other parts of the parameter space is is a superposition of the delta function at ψ = 0 and the P SF (ψ).
Numerically we identify two different benchmarks of the histograms representing P (ψ) that we generate: The first one is the value of the histogram at the first bin, P 0 , representing the potential delta peak of P (ψ). The second characteristic point is the value of the histogram at the second bin, P 1 , which is given by P 1 = P (ψ = δ ψ ) with δ ψ being the bin size of the histogram representing P (ψ) (δ ψ = 0.0025 in LMF and δ ψ = 0.015 in SMF theory).
In Figure B1 the regions determined on the basis of these benchmarks are shown for fixed density n = 1 on the left and fixed disorder strength ∆/U = 0.6 on the right. The LMF results are depicted with circles, SMF results with crosses. The blue curves enclose the regions in which P (ψ) is just a delta function at ψ = 0 (numerically: P 0 > 0 and P 1 < 10 −3 ), i.e. the system contains only MI sites and is therefore in the MI phase. The red curves delimit the regions, in which the delta function part of P (ψ) vanishes (numerically P 0 < 10 −6 ), which means that all sites are SF sites. In the remaining part of the phase diagram P (ψ) consists of a superposition of a delta function at ψ = 0 and a continuous part P SF (ψ), i.e. the system has MI and SF sites. In this region the system can either be in the BG phase, or, if the SF sites percolate, in the SF phase. Therefore only the MI-BG phase boundary can be extracted from the characteristics of the probability distribution of the local order parameter, P (ψ), but not the BG-SF boundary.
LMF and SMF theory coincide very well at the blue line, which describes the occurrence of the SF distribution. Deviations are visible at the red line, where the delta function at ψ = 0 disappears and the PD is purely given by P SF (ψ). While these discrepancies are rather small for small disorder strengths and U/J > 23, they enlarge with increasing disorder. In Appendix A we tested the validity the SMF approximation (15) by comparing the LMF results for the product of the PD of two different sites P (ψ i ) P (ψ j ) with the pair PD P Z (ψ i , ψ j ). In Figure A1 its deviation ∆ P is shown in dependence of µ/U and JZ/U for fixed disorder ∆/U = 0.6. In comparison with the diagram on the right in Figure B1 it is obvious that in the region of the blue line the SMF assumption (15) for P (ψ) is valid. However, at the red curve deviations between LMF and SMF theory occur, since close to the BG-SF phase transition the approximation (15) is expected to be invalid.