Paper • The following article is Open access

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

A E Niederle¹ and H Rieger¹

Published 30 July 2013 • © IOP Publishing and Deutsche Physikalische Gesellschaft
New Journal of Physics, Volume 15, July 2013 Citation A E Niederle and H Rieger 2013 New J. Phys. 15 075029 DOI 10.1088/1367-2630/15/7/075029

Download Article PDF

Article metrics

3558 Total downloads

Author e-mails

astrid@lusi.uni-sb.de

h.rieger@mx.uni-saarland.de

Author affiliations

¹ Theoretical Physics, Saarland University, D-66041 Saarbrücken, Germany

Dates

Received 18 April 2013
Published 30 July 2013

Buy this article in print

Journal RSS

Abstract

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 ⩾ 2 dimensions by analysing 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 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 inbetween: 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 such as the averaged SF order parameter and the compressibility is 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.

Export citation and abstract BibTeX RIS

Previous article in issue

Next article in issue

Content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

Supplementary data

1. 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 particular, the influence of disorder on a system of bosons in a regular (e.g. optical) lattice has received much interest since the fundamental work of Fisher et al [2]. Here the phase diagram and transitions for bosons in a disordered potential were analysed and the existence of a Bose glass (BG) phase was predicted. BG represents a non-SF but, in contrast to the MI, a 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–5].

Studies of the excitation spectrum of the disordered system as a function 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 ultracold bosonic gases has 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 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 dimensions [4, 12–17]. In addition, density-matrix renormalization group techniques were applied to one-dimensional systems containing either quasi-periodic potentials [18–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–28], were proposed.

An intriguing question has been for a long time the potential existence of a direct MI–SF transition [4, 12–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 [2, 29] due to arbitrarily large, but exponentially rare, clusters of one phase within a background of another phase [3, 30, 31]. Since any exponentially rare event is hard to sample numerically, the existence of an intervening BG phase might have been eluded from some studies, be it QMC [4, 14–17], LMF [26–28] or SMF theory [24, 25]. 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 is the region where SF clusters percolate—the BG phase is in between: SF clusters exist, but do not percolate. We apply this criterion to the results of LMF calculations and compare it with the predictions of other methods: on one side we have SMF theory, where the individual phases are identified on the basis of spatially averaged LMF quantities such as SF order parameter or compressibility, and on the other side we have 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 theories for the disordered BHM. In section 3 we critically examine the use of the averaged SF order parameter $\skew3\bar{\psi }$ and the compressibility κ as indicators of the different phases of the disordered BHM in the LMF and SMF theories 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 two-dimensional (2D) disordered BHM with commensurate filling and with fixed disorder strength and compare 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.

2. The model

Ultracold bosonic atoms in a 2D square optical lattice can be described by the BHM

where i = 1,...,M is the site index, M = L² 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 〈i,j〉 of the underlying lattice. The operator $\hat {n} ^{\vphantom {\dagger }}_{i}=\hat {a} ^{\dagger }_{i}\hat {a} ^{\vphantom {\dagger }}_{i}$ is the particle number operator of bosons on site i, which are annihilated and created by the operators $\hat {a} ^{\vphantom {\dagger }}_{i}$ and $\hat {a} ^{\dagger }_{i}$ . Moreover, on-site disorder is introduced by the parameter _i, which is drawn randomly from a box distribution $p(\epsilon )=\Theta \left (\Delta /2-|\epsilon |\right )/\Delta$ , where Δ is the strength of the disorder.

In the SF regime tunnelling dominates the system and in the limit J → ∞ the ground state is a coherent state, which is an eigenstate of the tunnelling part of the Hamiltonian. The SF parameter $\psi _i=\langle \hat {a_i} \rangle$ is the expectation value of $\hat {a} ^{\vphantom {\dagger }}_{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 $\hat {a} ^{\vphantom {\dagger }}_{i}$ in a Fock state is zero in any case. Hence, the mean value of the SF parameters $\skew3\bar{\psi }=\sum _i \psi _i/M$ is an order parameter for the SF phase (in fact, $|\skew3\bar{\psi }|^2$ is proportional to the condensate fraction).

The phase diagram of the pure system as a function of JZ/U and μ/U consists of the so-called Mott lobes, in which the system is an 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 $\kappa =\langle \hat {n} {\vphantom {\dagger}}^2\rangle -\langle \hat{n} {\vphantom {\dagger}}\rangle ^2$ is necessary. Among the three phases, only the MI is non-compressible and only in the SF phase the SF order parameter is non-zero. Consequently, the phase of the system can be identified by the SF order parameter $\skew3\bar{\psi }=\sum _i \psi _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 probability distribution (PD) self-consistently, or via LMF theory [23, 27, 28], 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.

2.1. Local mean field theory

In LMF theory the tunnelling part of the Hamiltonian can be approximated via

where terms of the form $(\hat {a} ^{\vphantom {\dagger }}_{i}-\langle \hat {a} ^{\vphantom {\dagger }}_{i} \rangle ) (\hat {a} ^{\dagger }_{j}-\langle \hat {a} ^{\dagger }_{j} \rangle )$ are neglected [23]. The central quantities are the local SF 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 to

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 $\eta _i:=\sum _j A_{ij} \psi _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 self-consistency equations

is solved on an L × L lattice, where the expectation value is evaluated in the ground state of $\hat {H} ^{\vphantom {\dagger }}_{i}$ , that itself depends on 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 $|\phi _i\rangle = \sum _{n=0}^\infty c_n^i |n\rangle _i$ given in particle number basis, and |cⁱ_n|² describes the probability to find n bosons at site i and fulfils $\sum _{i=1}^M |c_n^i|^2=1$ . In particular, the local SF parameter is then given by $\psi _i=\sum _{n=0}^\infty {c_{n-1}^i}^* c_n^i \sqrt {n_i}$ and the local boson number is represented by $\langle \hat {n} ^{\vphantom {\dagger }}_{i} \rangle =\sum _{n=0}^\infty |c_n^i|^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⁻⁴. 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

with $\hat {N} {\vphantom {\dagger}} =\sum _i\hat {n} ^{\vphantom {\dagger }}_{i}$ can directly be computed. The PD

is determined on the basis of the complete set of values of ψ_i, additionally averaged over disorder realizations.

2.2. 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 PD P(ψ) directly [24, 25]. Additional approximations are of course necessary. The SF order parameter $\psi =\langle gs | \hat {a} {\vphantom {\dagger }}| gs\rangle$ , 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) as a function of the stochastic variables and η. The parameter is then a stochastic variable drawn from the disorder distribution $p\left (\epsilon \right )$ and as a result $\psi =\langle gs | \hat {a} {\vphantom {\dagger }} | gs\rangle$ is a stochastic variable, drawn from the PD $P \left (\psi \right )$ , which must be determined self-consistently. Since η is the sum of the SF parameters of the neighbouring sites it also is a stochastic variable drawn from $Q\left ( \eta \right )$ . 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\left (\epsilon , \eta \right )\rangle$ of the site-independent Hamiltonian

as a function of and η. Thus, the PD

depends on the distribution $Q\left ( \eta \right )$ of the occurring values of η and the distribution $\tilde P_\eta \left (\psi \right )$ of the local SF parameters for given η. A direct analysis of $\langle gs\left (\epsilon , \eta \right )|\hat {a} {\vphantom {\dagger }}|gs\left (\epsilon , \eta \right )\rangle$ as a function of and η yields

Since η is the sum of the local SF parameters ψ of the neighbouring sites, its distribution is given by

where $P_Z \left ( \psi _1,\ldots ,\psi _Z \right )$ is the connected PD function of the local order parameters ψ₁,...,ψ_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 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 second order. We examine the validity of this approximation as a function of the system parameter JZ/U and μ/U in appendix A.

After determining the PD P(ψ), the SF order parameter $\skew3\bar{\psi }=\int {\mathrm { d}}\psi P\left ( \psi \right ) \psi$ is given by the mean value of the distribution. The compressibility $\kappa =[\langle \hat {N^2} {\vphantom {\dagger }}\rangle -\langle \hat {N} {\vphantom {\dagger }}\rangle ^2]_{\mathrm { av}}$ with $\hat {N} {\vphantom {\dagger }} =\sum _i\hat {n} ^{\vphantom {\dagger }}_{i}$ 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 figures 4 and 5 and will be discussed in the following.

3. Criterion for phase transition

In this section we will first discuss the well-known SF order parameter $\skew3\bar{\psi }$ 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 $(\kappa =0)$ and non-coherent $\skew3\bar{\psi }=0$ , while conversely in the SF regime it is described by a coherent state $\skew3\bar{\psi }=0$ , which is compressible $(\kappa >0 )$ . If disorder is introduced, those phases are separated by the BG phase, which is compressible $(\kappa >0)$ , but not coherent $\skew3\bar{\psi }=0$ . We will see that a precise prediction of the transition point on the basis of $\skew3\bar{\psi }$ 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.

3.1. Superfluid order parameter and compressibility

In the ordered case, the on-site energies _i are zero and the lattice is homogeneous. The SF order parameter $\skew3\bar{\psi }$ clearly marks the location of the MI–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 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, within both SMF and LMF theories. The resulting phase transition agrees 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 is 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 $\skew3\bar{\psi}$ and the compressibility κ for fixed chemical potential μ/U and varying tunnelling rate JZ/U. Left: homogeneous case (Δ/U = 0); right: disordered case with Δ/U = 0.6. The top row is for μ/U = 1.05, where the ordered system displays a MI–SF transition and the disordered system a BG–SF transition (κ > 0 for all values of JZ/U). The bottom row is for μ/U = 0.32, where the ordered system again displays an 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.
Download figure:
Standard image High-resolution image

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 0.1092, see figures 1(b) and (d), above which $\skew3\bar{\psi }$ becomes non-zero, the behaviour of $\skew3\bar{\psi }$ as predicted by LMF theory does not indicate a transition at all, it varies smoothly with the tunnelling rate J. This is not a finite-size effect as we have checked by examining different lattice sites, as shown in figure 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 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 $\skew3\bar{\psi }$ that overestimates the phase coherence. A closer look at the complete PD P(ψ) of the local SF parameters and their geometrical features provides an estimate of the SF regions in the phase diagram. Its prospects are discussed in section 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 $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ .

3.2. 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 state 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ⁱ_n = |cⁱ_n|² (cf 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 $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ is an integer, whereas in the SF regime it is a 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 the MI regime one expects to encounter spatially inhomogeneous situations, in which some sites have sharp (integer) boson occupation numbers and others have fluctuating (non-integer) boson numbers. Introducing phase operator $\hat {\Phi }_i$ that is canonically conjugate to the boson number operators $\hat {n}_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–SF transition in the d-dimensional BHM (d ⩾ 2) coincides with the percolation transition of the connected clusters of lattice sites with a non-integer boson number expectation value $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ . 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 out in [26, 33, 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 sites with non-integer boson occupation number $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ as SF sites, and sites with integer $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ as MI sites. Analogously, we discriminate between 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 integers (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 an 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 = \langle gs|\hat {n} ^{\vphantom {\dagger }}_{i}|gs \rangle$ , 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⁻³ 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 particle number $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ has narrow peaks of width γ_n at integer values. The width decreases when we increase the number of iteration steps to solve the self-consistency 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 $\langle \hat {n} ^{\vphantom {\dagger }}_{i} \rangle$ and in the third the resulting discrete map S_i are 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 $\skew3\bar{\psi }$ 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 do not yet percolate. 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 coherent macroscopically, which means that all local SF parameters ψ_i are positive and compressible, as described above.

$\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle $ — **Figure 2.** Configurations of the local SF parameter ψ_i, the occupation number $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ and the discrete variable S_i (18) for a single realization of disorder for Δ/U = 0.6. The first row shows an example for MI (JZ/U = 0.0242,μ/U = 0.4394) followed by one for BG (JZ/U = 0.0182,μ/U = 1.0455) and SF (JZ/U = 0.141,μ/U = 1.0455). Note that blue marks the minimal value (zero in the left and right, one in the middle column) and green the maximal value (one in the left and right, two in the middle column). A video illustrating in more detail the evolution of the configurations with increasing tunnelling rate can be found in the supplementary data (video 1, available from stacks.iop.org/NJP/15/075029/mmedia).
Download figure:
Standard image High-resolution image

3.3. Percolation analysis

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 that we study the phase diagram as a function of the system parameters called x and y. Then, 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 ν is the critical exponent determining the divergence of the mean lateral cluster size at the transition. The scaling function $\tilde p(X)$ approaches zero for X ≪ 0 and one for X ≫ 0, which means that exactly at the transition x_c the curves for different system sizes should intersect (in the scaling limit). This intersection point, which 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 reformed 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 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 1 for Δ/U = 0.6, yielding J_cZ/U = 0.15 and 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.

**Figure 3.** The percolation probability p_Perc of the SF cluster (top) and finite-size scaling plot (bottom). The critical tunnelling rate according to the finite-size scaling (19) is given by J_cZ/U = 0.15 for μ/U = 0.439 (left) and J_cZ/U = 0.04 for μ/U = 1 (right); the critical exponent is ν = 1.33.
Download figure:
Standard image High-resolution image

4. Results

4.1. Commensurate filling—comparison with quantum Monte Carlo results

In this section, we determine the complete phase diagram for commensurate density as a function of Δ/2J and U/J, for which a prediction on the basis of QMC simulations is available [12]. We fix the particle density to $\bar{n}=\langle \sum _{i=1}^M \hat {n} ^{\vphantom {\dagger }}_{i} \rangle /M=1$ with an accuracy of 10⁻⁴ 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 $\bar{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 $\bar{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 the MI, BG and SF phases. 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 a 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 [12, 17, 36, 37].

**Figure 4.** Left: LMF cluster analysis phase diagram for commensurate boson density $\bar{n}=1$ determined with the discretized boson occupation number field S_i (18). The percolation transition of the SF sites (S_i = 1) occurs when crossing the red line, which indicates the BG–SF phase boundary. The blue line marks the boundary of the MI region, in which all sites are MI sites (S_i = 0). The black line indicates the MI–BG transition according to the perturbative result (17). Right: prediction for the phase diagram for commensurate boson density $\bar{n}=1$ based on the results of QMC simulations (data taken from [12]).
Download figure:
Standard image High-resolution image

**Figure 4.** Left: LMF cluster analysis phase diagram for commensurate boson density $\bar{n}=1$ determined with the discretized boson occupation number field S_i (18). The percolation transition of the SF sites (S_i = 1) occurs when crossing the red line, which indicates the BG–SF phase boundary. The blue line marks the boundary of the MI region, in which all sites are MI sites (S_i = 0). The black line indicates the MI–BG transition according to the perturbative result (17). Right: prediction for the phase diagram for commensurate boson density $\bar{n}=1$ based on the results of QMC simulations (data taken from [12]).
Download figure:
Standard image High-resolution image

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 BG–SF phase boundary with its characteristic nose and bumps clearly coincides. 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, cf figure 2 in [12]). Moreover, with our method we could also explore the weak interaction region, which is hardly accessible by QMC methods and we found a singular behaviour and a steep rise of Δ with U → 0, which is compatible with the analytically predicted behaviour $\Delta \propto \sqrt {U}$ [38]. Note that Δ rises steeply for U → 0; and in order to obtain data for Δ < 20, it would be necessary to compute the ground state of even larger system sizes. Finally, 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 the LMF theory.

Our result for the MI–BG transition line, which denotes the appearance of non-integer 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 an 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 with the theorems proven in [13], which exclude a direct MI–SF transition in any disordered system. Besides, we checked that the BG phase occurs for even higher values of U/J, which is in agreement with the results to be presented in the next section.

4.2. Fixed disorder—comparison with stochastic mean field 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 fixed disorder strength Δ/U = 0.6 and compare it with predictions of SMF theory. In this phase diagram the Mott lobes occur and the line given by $\bar{n}=1$ always passes the tip of the first one.

In section 2.2, we introduced SMF theory and already emphasized that SMF theory is based on the same approximation to the Hamiltonian as LMF theory, but it involves the additional approximation (15) on the distribution $P_Z \left ( \psi _1,\ldots ,\psi _Z \right )$ . 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 $\skew3\bar{\psi }$ as well as the compressibility κ computed within SMF theory are exactly zero in specific regions of the parameter space (cf figure 1), which one might want to identify with the MI and BG phases, as was done in [24, 25], cf 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 phase 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 that 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 $\skew3\bar{\psi }$ in SMF theory.

**Figure 5.** Comparison of LMF cluster analysis and SMF phase diagram for fixed disorder strength Δ/U = 0.6. Left: LMF cluster analysis phase diagram determined with the discretized boson occupation number field S_i (18). The percolation transition of the SF sites (S_i = 1) occurs when crossing the red line, which indicates the BG–SF phase boundary. The blue line marks the boundary of the MI region, in which all sites are MI sites (S_i = 0). The black line is the MI–BG transition according to the perturbative result given by (17). Right: SMF phase diagram determined by using the SF order parameter $\skew3\bar{\psi }$ and the compressibility κ [24, 25]. The red line indicates the critical tunnelling rate where the SF order parameter $\skew3\bar{\psi }$ becomes non-zero; the blue line indicates the critical tunnelling rate, where the compressibility κ becomes non-zero.
Download figure:
Standard image High-resolution image

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 $\bar{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 Δ/U = 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.

Table 1. Comparison of the parameters at the tip of the first Mott lobe of quantum models (1) with LMF cluster analysis results. The chemical potential μ/U for the QMC results of [12] is our LMF estimate for a density $\bar{n}=1$ and fixed values of JZ/U and Δ/2J. The BG–SF predictions of [12] were not obtained by QMC of the disordered BHM, but are based on gap data for the ordered BHM [39].

	MI–BG		BG–SF
	JZ/U	μ/U	JZ/U	μ/U	Δ/U
LMF (cluster analysis)	0.1115	0.4644	0.1509	0.4434	0.6
QMC results [12]	0.124	(0.4561)	0.2012	(0.4082)	0.6
Strong-coupling expansion [11]	0.1253	0.4345			0.6
LMF (cluster analysis)			0.0942	0.4868	1
QMC results [12]			0.1047	(0.4846)	1
LMF (cluster analysis)			0.0934	0.5043	2
QMC results [12]			0.1062	(0.4950)	2
Quantum rotors model [37]			0.112	0.375	2

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 [26, 27] to solve the self-consistency equation (3) and to calculate an 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 [26, 27], the LMF approximation, is identical to ours, the use 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 Gutzwiller projected variational techniques are 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 LMF 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 phases.

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.

4.3. Probability distribution of the local superfluid parameter

In this section, we take a closer look on the PD of the local SF parameter and describe its three different characteristic shapes. Moreover, we will discuss the occurring deviations of the PD in the LMF and SMF theories.

In LMF theory, for every realization of disorder, the self-consistent solution of (6) is determined. The resulting PDs of each are averaged over 200 different realizations of disorder and shown in figure 6 in blue for fixed chemical potential μ/U = 1.05 and disorder strength Δ/U = 0 and 0.6. In the ordered case, depicted in the first row, all local SF parameters are identical since all sites have the same on-site energy. The averaged order parameter $\skew3\bar{\psi }$ , 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 $\skew3\bar{\psi }$ . 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. 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. Although, in the MI regime the PD is a sharp delta function at ψ = 0 still, it becomes a broad distribution in the BG and SF phases. 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 a variety of positive values for the local SF parameter. Owing 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 call this distribution P_SF(ψ), since it represents sites, which we referred to as SF sites (S_i = 1) before:

**Figure 6.** PD P(ψ) of the local SF parameter at fixed chemical potential (μ/U = 1.05) for different tunnelling rates JZ/U, (top) in the ordered case (Δ/U = 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 $\skew3\bar{\psi }$ . 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.
Download figure:
Standard image High-resolution image

Starting with the same Hamiltonian (11) as LMF, SMF theory introduce the additional approximation (15), yielding 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 theories occur in the ordered $(\Delta /U=0)$ , they become visible in the disordered case. For Δ/U > 0 both distributions show roughly the same shape, but especially for small values of ψ, which are crucial for the determination of the phase transition, they differ significantly in the BG and SF phases.

In order to identify the reason for these deviations we have to take a closer look on 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 $\mathcal P_Z (\psi _i,\psi _j)$ and determined their deviation Δ_P in appendix A. As shown in figure A.1 the approximation (15) is best in the MI regime and the BG for very small tunnelling rates. But for increasing JZ/U it becomes worse and especially at the phase boundary it fails. In figure A.2 where P(ψ_i)P(ψ_j) and $\mathcal P_Z (\psi _i,\psi _j)$ are shown for different parameters; it is evident that they disagree especially for small values of the ψ. This disagreement is due to the presence of correlation of the local SF parameters ψ_i at different sites in the vicinity of the transition points, which are neglected in SMF theory and are therefore responsible for the deviation of the SMF PD from the LMF PD shown in figure 6.

5. 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 $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ are integers; 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 $\skew3\bar{\psi }$ and κ and all approaches to determine the LMF phase boundaries of the disordered system on the basis of the site and disorder averaged parameters, such as the SF order parameter $\skew3\bar{\psi }$ , the compressibility κ, the SF or condensate fraction, overestimate 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 by only a small amount.

The resulting cluster analysis phase diagram for a fixed commensurate density $\bar{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 remarkable, since other LMF approaches using averaged quantities, such as 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 because 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 on these data.

The questions that immediately arise in this context are: (i) is the MI–BG transition in the disordered BHM exactly where SF sites occur? And more interestingly: (ii) is the BG–SF transition actually identical to the percolation transition of SF clusters?

Concerning question (i), 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. $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ is non-integer. This is how we propose to identify the MI–BG transition.

Concerning question (ii), we argued in section 3.2 on the basis of the BHM to quantum rotor models, that the SF stiffness will always vanish as long as SF clusters within a MI background do not percolate. This BG situation then is a reminiscent of a d + 1-dimensional, 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 on 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 $\langle \hat {n} ^{\vphantom {\dagger }}_{i}\rangle$ (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 [3, 30, 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 disordered 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 superfluid parameter correlations

In this appendix, we check the statistical independence assumption underlying SMF theory. In addition to the LMF approximation (4) made in the tunnelling part in the Hamiltonian, SMF theory assumes that the local SF parameters ψ₁,...,ψ_Z of the Z neighbours of a chosen site i are uncorrelated and identically distributed, which is introduced by approximation (15). On the basis of LMF calculations, we want to test this approximation by comparing P(ψ_i)P(ψ_j) and $\mathcal P_Z\left (\psi _i,\psi _j\right )$ , of which some examples are shown in figure A.2. The function $P\left (\psi _i\right )P\left (\psi _j\right )$ is the product of the PD, describing the distribution of the local SF parameter as discussed in section 4.3. The function $\mathcal P_Z (\psi _i,\psi _j)$ is the PD of pairs $(\psi _i,\psi _j)$ , where i and j are neighbouring sites; it represents the probability of having a specific value for the pair $(\psi _i,\psi _j)$ . Exactly as $P\left (\psi _j\right )$ , $\mathcal P_Z (\psi _i,\psi _j)$ is computed for every realization of disorder and finally averaged. Both distributions should coincide if assumption (15) is valid. In figure A.1 the integral difference

of both distributions is shown in parameter space.

**Figure A.1.** The deviation Δ_P between P(ψ_i)P(ψ_j) and $\mathcal P_Z (\psi _i,\psi _j)$ is shown as a function 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 A.2 along a line μ/U = 1.0455 and at the tip of the Mott lobes, where the deviation is maximal.
Download figure:
Standard image High-resolution image

**Figure A.1.** The deviation Δ_P between P(ψ_i)P(ψ_j) and $\mathcal P_Z (\psi _i,\psi _j)$ is shown as a function 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 A.2 along a line μ/U = 1.0455 and at the tip of the Mott lobes, where the deviation is maximal.
Download figure:
Standard image High-resolution image

**Figure A.2.** The first column shows $\mathcal P_Z(\psi _i,\psi _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 to the black dots in figure A.1.
Download figure:
Standard image High-resolution image

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 A.2 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 A.1). In addition 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 $\mathcal P_Z (\psi _i \psi _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 figures A.2(d) and (h), which corresponds to the tip of the Mott lobe. These main differences in comparison with P(ψ_i)P(ψ_j) can be observed for all parameters shown in figure A.2 and mainly occur in the regime of small local SF parameters. Figure A.1 illustrates that assumption (15) made in the 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 be questioned.

The deviations occurring for small local 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 BG–SF phase transition as a percolation transition (cf sections 4.1 and 4.2).

Appendix B.: Characteristic shapes of the probability distribution

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 phases 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 sections 4.1 and 4.2.

Numerically, we identify two different benchmark of the histograms representing P(ψ) that we generate. The first one is the value of the histogram at the first bin, P₀, representing the potential delta function of P(ψ). The second characteristic point is the value of the histogram at the second bin, P₁, which is given by P₁ = P(ψ = δ_ψ) with δ_ψ being the bin size of the histogram representing P(ψ) (δ_ψ = 0.0025 in LMF and δ_ψ = 0.015 in SMF theory).

In figure B.1 the regions determined on the basis of these benchmarks are shown for fixed density $\bar{n}=1$ on the left and fixed disorder strength Δ/U = 0.6 on the right. The LMF results are depicted with circles and the SMF results with crosses. The blue curves enclose the regions in which P(ψ) is just a delta function at ψ = 0 (numerically P₀ > 0 and P₁ < 10⁻³), 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₀ < 10⁻⁶) 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 be either 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 PD of the local order parameter, P(ψ), but not the BG–SF boundary.

**Figure B.1.** Regions of the three characteristic shapes of P(ψ) for fixed density $\langle \hat {n} {\vphantom {\dagger }}\rangle =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 the 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 it is a superposition of the delta function at ψ = 0 and P_SF(ψ).
Download figure:
Standard image High-resolution image

The LMF and SMF theories 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 PD is purely given by P_SF(ψ). While these deviation are rather small for small disorder strengths and U/J > 23, they enlarge with increasing disorder. In appendix A, we tested the validity of 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 $\mathcal P_Z (\psi _i,\psi _j)$ . In figure A.1 its deviation Δ_P is shown as a function of μ/U and JZ/U for fixed disorder Δ/U = 0.6. In comparison with the diagram on the right, in figure B.1, 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 theories occur, since close to the BG–SF phase transition the approximation (15) is expected to be invalid.

Please wait… references are loading.

Video 1. (5.29 MB, MP4) Evolution of the configurations of the local superfluid parameter with increasing tunnelling rate.