{'fraction_non_alphanumeric': 0.03767660910518053, 'fraction_numerical': 0.013036652788205583, 'mean_word_length': 4.497748592870544, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Bose-Hubbard model with random impurities: Multiband and nonlinear hopping effects

We investigate the phase diagrams of theoretical models describing bosonic atoms in a lattice in the presence of randomly localized impurities. By including multiband and nonlinear hopping effects we enrich the standard model containing only the chemical-potential disorder with the site-dependent hopping term. We compare the extension of the MI and the BG phase in both models using a combination of the local mean-field method and a Hartree-Fock-like procedure, as well as, the Gutzwiller-ansatz approach. We show analytical argument for the presence of triple points in the phase diagram of the model with chemical-potential disorder. These triple points however, cease to exists after the addition of the hopping disorder.


I. INTRODUCTION
Ultracold atoms in optical lattices serve now as a routine tool to study various lattice models derived from other areas of physics, such as condensed matter or high energy physics. They often enrich the original models with additional features accessible due to extreme controllability and versatility of possible experimental realizations (for recent reviews see [1][2][3]). This potential of cold atoms was recognized for the first time thanks to a seminal theoretical proposal of Jaksch et al. [4]. Soon followed the experimental demonstration of the superfluid-Mott insulator phase transition [5], and the intensive studies of the effects of disorder on cold-atom systems [6,7]. After a series of attempts [8][9][10][11][12][13][14] the Anderson localization for non-interacting (expanding) atoms was unambiguously demonstrated in [15], The experimental studies of the phase diagram of interacting bosonic atoms in a disordered potential revealed the long-debated gapless insulating Bose glass (BG) phase [16], stimulating at the same time the discussion on the ways to observe and detect the BG phase [17][18][19][20][21]. Various aspects of Anderson localization in cold atoms are reviewed in Ref. [22,23], while the importance of disorder studies in a broader context of quantum simulators is presented in [24].
Optical implementations of quenched disorder are unique in the sense that the disorder can be, in principle, controlled with high precision on demand. Various methods of creating disorder have been proposed. A speckle pattern collimated on the atomic sample [6] gives rise to a truly random intensity landscape, which follows the exponential distribution. Another proposed scheme [6,7] applies several (at least two) laser fields with different frequencies. For an appropriate choice of frequencies the resulting potential is quasi-periodic, and for a finite sample hardly distinguishable from a truly random case (see [25], and also [26][27][28] for recent results).
Yet another interesting way to create the disordered potential is to use the interactions between atoms. By pinning a secondary type of atoms in an optical potential we obtain the disorder with binary (or Bernoulli) distribution. Such proposal, originating from the work of Gavish and Castin [29], has been frequently discussed theoretically [30][31][32][33][34], and only recently implemented experimentally [35] (note that in the early papers [36] the impurities were mobile). This is the type of disorder we consider in this paper.
Specifically, we study the bosonic atoms in the optical lattice potential, which interact with immobile randomly distributed atoms of a secondary (fermionic) species. Such a situation is routinely described using the Bose-Hubbard model with the random potential i are the bosonic annihilation and creation operators, t is the tunneling, the interaction constant is denoted by U , µ is the chemical potential and i, j denotes the summation over the nearest neighbors. The parameter γ characterizes the strength of the disorder and ω i is a random variable with binary distribution (taking value 1 with probability p when the heavy background fermion is present at the ith site, and 0 with probability (1 − p), when no impurity is present). This Hamiltonian was already studied in several works [31][32][33], where the emergence of the Mott insulating (MI) phase with non-integer filling related to the impurity density was demonstrated, and the MI phase was shown to survive for arbitrarily strong disorder unlike in the continuous disorder case.
It has been shown recently that the simple Bose-Hubbard description for fermion-boson mixture may not be adequate for stronger interspecies interactions. The shift of the observed transition between the superfluid (SF) and the MI phases, observed experimentally in Refs. [36][37][38] (for the analogous effects with Bose-Bose mixtures and tightly trapped bosons, see [39,40]), could not be explained with this simple description. It has soon been realized that density-dependent tunneling terms as well as contributions coming from higher Bloch bands are necessary to describe the systems in question [41][42][43][44], whenever the interspecies interaction becomes strong enough. As a consequence, for bosons interacting with immobile fermionic impurities the disorder affects also the tunneling. While similar corrections could be also taken into account for boson-boson interactions, we do not include them to simplify the picture, assuming bosonboson interactions to be sufficiently weak (see a recent review [45] for discussion of different possible contributions). In this simplified picture we add to the Hamiltonian (1) the term T (ω i + ω j )]b † i b j yielding the tunneling dependent on the presence of the heavy fermion. The density induced tunneling coefficient T , proportional to boson-fermion interaction strength, seems at first independent of the standard tunneling t. However, this is not the case in the optical lattice potential, where both t and T depend on the potential depth (when proper Wannier functions are used to evaluate them), and are (see e.g. [45]) approximately proportional to each other for standard depths of optical lattices. Thus, we may assume T = αt obtaining the Hamiltonian: Here, α and γ depend on the interaction between the two species and, typically, |γ| > |α| [45]. They are of the same sign, and below we only discuss the case of positive α and γ. This corresponds to the repulsive boson-fermion interaction. Note that in this model a single local random variable related to the presence of the impurity, enters in the potential and the tunneling terms of the Hamiltonian of the system. For the sake of comparison we shall also study a different model in which the disorder in the tunneling is given by independent random variables Ω i = ω i (see also [46,47] for models with simultaneous potential and hopping disorder): The main goal of the present paper is to compare the phase diagrams of Hamiltonians H 1 , H 2 and H 3 . To this end we develop a method being a combination of a site-dependent decoupling mean-field method with a "Hartree-Fock-like" procedure [48]. Independently we use also the Gutzwiller-ansatz approach. Note that the type of disorder that we study does not fulfill the assumptions of the "theorem of inclusions", valid for continuously distributed bounded disorder [49]. We perform the analysis mostly in two dimensions (2D) as the mean-field approaches cannot be regarded as accurate in 1D. Our main results are: i) for the Hamiltonian H 1 there is a direct MI-SF transition at the tips of the Mott lobes in the thermodynamic limit; ii) for the Hamiltonian H 2 the Mott lobes are smaller, and the direct MI-SF transition at their tips disappears, although the region of BG is very narrow; iii) finally, for H 3 there is no direct MI-SF transition and the region of BG is wider in comparison to H 2 . One of the main aspects of this paper is also the use of a novel form of the mean-field theory combined with the simple Hartree-Fock approach, quite different from what has been proposed so far [50][51][52][53][54][55][56], and quite efficient in determining the phase boundary between the BG and the SF phase.
The paper is structured as follows. In section II we briefly recall the mean-field approaches applied in the studies of the disordered Bose-Hubbard model. In section III we introduce a new method being a combination of a local mean-field approach and a Hartree-Fock-like method. The standard Gutzwiller approach used later for comparison is presented in IV. In section V we apply the theory developed in section III to compare the phase diagrams of the bosonic atoms interacting with immobile impurities on a lattice described by the Hamiltonians H 1 , H 2 and H 3 and confront these results with the Gutzwiller approach. Finally, we conclude in section VI by summarizing the obtained results.

II. BRIEF SURVEY OF MEAN-FIELD APPROACHES FOR THE DISORDERED BOSE-HUBBARD MODEL
The phase diagram of the disordered Bose-Hubbard model has been studied by a number of methods including the quantum Monte Carlo [57][58][59][60][61], renormalization group [62,63], density-matrix renormalization group techniques [64][65][66], tensor networks-based algorithms, or various mean-field approaches [6,50,67,68]. In this work we propose an extension of the local mean-field method, thus let us first briefly review the mean-field approaches used earlier.
The local mean-field method was introduced in Ref. [50], and further developed in [52], where the boundary of the Mott lobe was linked to the stability of the zero solution of the self-consistency equations, which is then studied through linearization of those equations. Moreover the authors suggest that the presence of the BG surrounding the MI phase could be inferred from the spectral properties of the random matrix, which appear in the linearized problem. In [53] the authors generalize the inhomogeneous site-dependent mean-field to clusters, which allows them to include the (short-range) spatial correlations. Sheshadri and co-workers [51] proposed to study the Bose-Hubbard model with disordered potential using the inhomogeneous generalization of the site-decoupling mean-field method, i.e., the local mean-field theory . There the hopping term is decoupled as b j yielding single site Hamiltonians coupled to neighbouring sites only through SF amplitudes ψ j = b j . The authors then diagonalize the local Hamiltonians in the occupation-number basis and determine self-consistently the values of ψ j minimizing the energy. The BG-SF transition is characterized through the percolation of sites with non-zero SF parameter.
A description equivalent to the local mean-field theory may be achieved by minimizing the average energy over a variational manifold composed of products of single site wave vectors. Indeed, under such assumptions a i a † j = a i a † j . This numerical ansatz is called the Gutzwiller ansatz.
Yet another approach, the stochastic mean-field method, was proposed by Bissbort and Hofstetter in [55] (see also [56]). There, the starting point is also the decoupling of the tunneling term, however instead of choosing a different mean-field parameter for each site, the authors consider the probability distribution P (ψ), reducing the description to effectively single-site problem. The probability distribution of ψ is then found self-consistently to ensure the compatibility of P (ψ) with the distribution on the neighbouring sites.
Finally, in a recent work Niederle and Rieger [54] compare the results obtained with the local mean-field theory and the stochastic mean-field method with the quantum Monte-Carlo results. They conclude that the identification of different phases based on averaged quantities obtained through the mean-field approaches is misleading. Instead, the authors propose to distinguish the phases through the presence and percolation of the SF clusters finding an excellent agreement with the quantum Monte-Carlo studies.
In what follows we will take the latter approach, i.e., study the percolation of the SF clusters; we propose, however, a novel method to determine the distribution of the "superfluid particles" in the lattice. To that end we combine the mean field approach with Hartree-Focklike method mixing different mean field modes. In this way, at a little numerical effort, we can go beyond the standard mean field approach.

III. LOCAL MEAN-FIELD APPROACH COMBINED WITH A HARTREE-FOCK-LIKE METHOD
The proposed method is a compromise between the local mean-field description which, based on the productstate description of the system, may not capture correctly the long-range correlations, and the resource-demanding numerical approaches. We also make an attempt to include the spatial correlations on top of the simple local mean-field description.
Since the Hamiltonians H 1 and H 2 can be considered as special cases of the Hamiltonian H 3 , in what follows we will concentrate on the latter.

A. Standard mean-field approach -local Hamiltonian
Like in the routine mean-field approach [50] we begin by decoupling the Hamiltonian (3) using the standard and introducing a local mean-field parameter ψ i = b i . As a result we obtain: where we also express all the parameters in the units of U , i.e.,t = t/U,μ = µ/U,γ = γ/U .

B. Standard mean-field approach -energy minimization and self-consistency
The ground state, or more generally Gibbs energy of the local mean-field Hamiltonian is a highly nonlinear function of the local mean fields ψ j 's. The next step in the standard approach is to find a minimum of the energy under the constraint that b j = ψ j . This is in general a complicated task, but as long as we are interested in finding the boundaries of the MI phase, the analysis can be restricted to small values of ψ j 's, where the energy is a quadratic form of ψ j 's. The solutions of the self-consistency equations can be then obtained via a perturbative expansion up to first order int: with the random matrix The MI phase corresponds to sucht,μ that the system admits only a trivial solution (the energy has the minimum at ψ j = 0 for all j). Clearly, this occurs if and only if det(tR − 1) = 0, in other words, whenever where λ(·) denotes the spectrum of a matrix. Oncet exceeds the critical value (the condition (8) is violated) we enter a phase with at least one unstable mode that attains a non-zero value, determined by the full nonlinear dependence of the energy on ψ i 's; obviously, the linear theory predicts only the instability and, formally, a value of the amplitude of the corresponding unstable mode tending to infinity.
C. Non-standard mean-field approach -populating unstable modes The idea here is a simple one -we consider all modes that are unstable (i.e. these that violate the inequality (8)). Note that for finite systems the spectrum of the matrix is discrete, hence only for specific values oft different eigenvectors of R become the solutions of Eq. (6), i.e., a fixed point of the linear maptR. Hence, in general it is more physical to consider the vectors ψ belonging to the unstable manifold of the maptR, rather than the individual solutions of the self-consistency equation (6). We denote by Q(t) the set of indices for the eigenvectors oftR corresponding to the eigenvalues larger then one. These are the modes that we expect to be populated.
Consequently, we define new modes, a k , a † k , corresponding to right and left eigenvectors of R, ψ (k) , ψ (k) , and related to the original modes as and express the initial (not decoupled) Hamiltonian (3) in terms of these operators. Our aim is to minimize the energy with respect to the population {n k } of the new modes. Taking the ground state in the form |g.s. = where Q(t) is the set of indices defined before, and The number of particles is adjusted for each value oft to match the chemical potentialμ in the definition of R.
Knowing the population of the modes, we determine the distribution of the number of particles n i in the lattice as: As long as the regions of non-zero (i.e. of absolute value greater then a threshold value) number of particles are disconnected we identify the phase as the BG.
The boundary of this phase and the transition to the SF phase is given byt for which the sites with non-zero n i begin to percolate. This condition seems analogous to the approach of Niederle and Rieger [54], however, there is an important difference. While in [54] the percolation of mean-field occupations of sites is directly taken as a superfluid border, in our approach we rebuild these occupations from the Hartree-Fock-like procedure discussed above.

IV. GUTZWILLER-ANSATZ
The results obtained following the approach discussed above will be in the next Section confronted with the standard Gutzwiller approach. In the latter the minimization of the energy functional E[ψ] = ψ|H|ψ over the product states of the form and subject to the normalization condition ψ|ψ = 1 is performed. One needs to minimize over the expansion coefficients f i (n), with |n i denoting the Fock states at site i. The numerical minimization of such a nonlinear problem is simple for a homogeneous system, as minimization variables f i (n) without loss of generality may be considered site-independent. For low densities one limits possible occupations to, say, n max = 5 making it a five-parameter minimization (taking into account the normalization). In such case the standard conjugategradient minimization algorithms always converge then to the global minimum. In contrast, in the presence of disorder the number of minimized parameters increases to n max L where L is the number of sites. More importantly, the energy landscape of the energy functional E[ψ] may contain plenty of local minima in the presence of disorder. To reach the (hopefully) global ground state additional precautions have to be made (starting from different initial conditions, perturbing the found minima to check whether they are the local or global ones etc.).
We have used a 2D lattice which contained M × M = 40 × 40 lattice sites (L = 1600) with n max = 4. Periodic boundary conditions have been assumed.
In the MI phase the mean field solution yields a Fock state at each site with a vanishing variance of the occupation number Thus a convenient criterion for the disappearance of MI is that σ 2 ≡ max σ 2 i exceeds a given threshold value s m . Of course the obtained results depend to some extent on the value of s m . Practically the dependence is quite small, we found s m = 0.001 leads to an almost perfect agreement between the Gutzwiller ansatz prediction and the eigenvalue condition (8).
Having the distribution of σ 2 i one could define the BG-SF transition as a border at which non-zero σ 2 i percolate. However, we employ another method calculating the classical property, the superfluid fraction (SF), ρ s which should vanish in BG phase. It is obtained using the "boost" method [6,7,70] transferring the system to the moving frame by making the tunneling amplitude complex. Explicitly for tunneling along the x-axis we change t → t exp(±iϕ) for tunnelings (with sign corresponding to the direction of tunneling and ϕ being a small angle). This corresponds to the presence of a constant flux proportional to ϕ and the SF is obtained as [6,7,70] ρ s = (E(ϕ) − E(0))/N tϕ 2 where N is the total number of atoms and E(ϕ) is the ground state energy at a given value of ϕ. In practice the parabolic dependence of E(ϕ) on ϕ is tested to extract a reliable SF. Let us, however, mention that the SF calculated in this manner is not entirely correct in the presence of the density dependent tunneling (for a detailed discussion of superfluid fraction definition in various situations see [71]). We believe, however, that since in our case the density dependent tunneling terms are small, the onset of the superfluidity may be well estimated by the traditional approach.
Before discussing the results let us stress that the SF fraction introduced in this way has little in common with the "averaged SF order parameter" criticized in [54]. The results obtained for the superfluid border using a proper superfluid fraction reproduce, in fact, the percolation border of [54] with quite a good accuracy. They are shown below together with the HF percolation threshold.

V. RESULTS
Consider first the MI phase and determination of its borders (MI lobes). We study the 1D systems of different lengths L = 100, 1000, 10000 and the 2D system of size L = 40, 50, 80 with open boundary conditions. This allows for semi-analytical expressions for the MI borders in the thermodynamic limit. In this limit the properties of the system do not depend on the density of the impurities [42] thus we conveniently choose p = 0.5. We also choose γ = 0.3, α = 0.1 in the model. Those parameters are, as mentioned above, determined by the boson-fermion interaction strength. For too strong interactions one would have to consider higher bands [42,43] not included in our model.
First, we find the boundary of the MI phase analyzing the spectrum of the matrix R. In the 1D case R is a tridiagonal matrix with the following random offdiagonal upper (X + i ) and lower (X − i ) elements and zeros elsewhere. As discussed by Mering and Fleishhauer [31] (see also [33]) in the thermodynamic limit the border of the fully incompressible Mott phase may be determined from non-random the situation, i.e., assuming that all sites ω i are identical. It is due to the fact that it is the tunneling which kills MI and it is optimal at resonance, i.e., between identical sites. This argument holds also for our model with density induced tunneling terms. Then the random variables X ± i take the same value for each i which allows us to estimate the spectrum of R using the formula for the tridiagonal Toeplitz matrix [72].
Now for large L the condition (8) takes the form Analogously for 2D we obtain where X ± i,j , Y ± i,j are given by expressions similar to (14) (see Appendix A for a more detailed derivation). In Fig.  1 we compare the theoretical boundary of the Mott lobes for the Hamiltonians H 1 , H 2 , H 3 with the one obtained for finite systems of sizes L = 100, 1000, 10000 (1D) and L = 40, 80, 100 (2D). Note that for the binary disorder Mott lobes with non-integer average filling equal to m + (1 − p), m = 0, . . . (number of free sites) appear between the standard Mott lobes [33] as for m <μ < m +γ the local number of bosons is site-dependent, i.e., m <μ <n i < 1 +μ < 1 + m +γ, for ω i = 0 (18) m −γ <μ −γ <n i < 1 +μ −γ < 1 + m, for ω i = 1, which gives n i = 1 + m in the first case and n i = m in the second. In order to assess whether in any of the three analyzed cases a direct MI-SF transition is possible or if the system always passes through an intermediate BG phase we analyse the properties of the spectrum of R. As indicated in Eqs. (16) and (17) the most relevant eigenvalues stemming from homogeneous region are X + i X − i for any i in that region. For large enough region the cosine term in (15) may be approximated as unity. It turns out that for the Hamiltonian H 1 for each choice of the parameter γ there exists µ for which the variance of the random variables X + i X − i , Var(Λ) vanishes. This value of µ corresponds to the tip of the MI lobe. In other words, at this point, the matrix R is exactly Toeplitz and not random. We compare the variances of the largest eigenvalues of R, corresponding to large homogeneous regions, for the three Hamiltonians in 1D and 2D in Fig. 2. Clearly, for H 2 and H 3 the variance of the eigenvalues never becomes zero. Moreover, the minimum variance of the eigenvalues for H 3 is much larger than for H 2 . This leads us to the conclusion that H 1 may have triple points and a direct MI-SF phase transition whenever Var(Λ) vanishes. The addition of the disorder in the tunneling term in H 2 removes the triple points. In this case at least a thin layer of the BG phase surrounds the Mott lobe everywhere. In the last situation, in which the disorder in the tunneling is independent from the random chemical potential, the area of the BG should be much wider despite the same strength of the disorder.
We now determine the boundary of the BG by finding the distribution of n i in the lattice as in Eq. (12). We compute the energy and the number of particles in the SF clusters using the redefined modes and the information about their population. The total number of particles which is adjusted to match the chemical potential provides an estimate of the condensate fraction. In Fig. 3 we show a typical distribution of the SF clusters in the regime of the BG and in the point of the transition to the SF phase in 2D. , the regions of the lattice with non-zero local occupation by the "superfluid" particles, in the systems described by Hamiltonians H1, H2, H3. In the first column the clusters do not percolate, hence the phase is identified as the Bose-glass. In the second column, the clusters begin to percolate and the SF phase emerges.
In Fig. 4 we compare the phase diagrams for the Hamiltonians H 1 , H 2 , H 3 obtained using the percolation approach with the BG-SF border coming from the Gutzwiller ansatz as explained in the previous sections. Observe that the percolation approach gives a consistently larger BG region, notably for H 3 , although in all the cases we find that the BG regions are quite small in 2D for our choice of parameters. This is in contrast to the 1D situation with prominent BG regions [33] for H 1 and also for H 2 and H 3 (not shown). We do not present the results for 1D since the mean field does not give reasonable quantitative predictions in this case.

VI. SUMMARY AND DISCUSSION
We studied the Bose-Hubbard model with binary disorder obtained by pinning a secondary (fermionic) type of atoms in the optical potential. We revealed the following differences between the models describing such system: i) the model with disorder entering only in the potential term admits a direct MI-SF transition, which is not in contradiction with the theorem of inclusions; ii) in the model with disorder given by the same random variable affecting the tunneling and the chemical potential the transition to the SF phase goes always through an intermediate Bose-glass phase; iii) the disorder in the tunneling and in the potential given by independent random variables makes the intermediate Bose-glass region a bit thicker then the correlated disorder of point ii). Still the BG region in 2D for the given type and strength of disorder yields only a rather narrow slip around the Mott lobes.
Moreover we introduced a new type of mean-field approach for disordered systems, which combines the local mean-field approach with a simple "Hartree-Fock-like" procedure. Its advantage stems from the fact that maintaining the simplicity of the local mean-field approach it allows one to bring some spatial correlations into the description.  ) give the MI-BG border, red line (triangles) corresponds to Gutzwiller ansatz BG-SF border obtained using superfluid fraction threshold, green line (squares) estimates the same border using HF percolation approach. Note that the latter is consistently higher than the SF prediction. framework of the Polish Innovation Economy Operational Program (POIG.02.01.00-12-023/08).
Appendix A: Estimation of the spectrum of the random matrix R We study the 1D case as the first step in the iterative procedure which allows us to estimate the spectrum of the matrix R for any dimension.
The structure of the matrix R for 1D is the following: where and µ − γω i <n i < 1 + µ − γω i . The elements in the off-diagonal bands form a set of identical random variables, therefore to estimate the spectrum we treat them as identical elements in the tridiagonal Toeplitz matrix of dimension L × L and evaluate the spectrum as Note that the random variable Λ k does not depend on the choice of i, as the random variables X ± i have identical distributions. The index is written explicitly only to specify how many independent ω's in (A2) should be taken into account. The maximum of Λ k is achieved for k = L and is equal to the maximum value the random variable 2 X + i X − i+1 can take. Having the spectrum (A3) one can also estimate the variance of the maximal eigenvalues of R, as In the 2D case the random matrix R has the following block structure: where R 1D i is the matrix (A1) for the ith row of the lattice, with X ± i,j = [1 + α(Ω i,j + Ω i,j±1 )] (A7) The off-diagonal block-band contains diagonal matrices of the form We again estimate the eigenvalues of R 2D using the formula for the spectrum of the L × L Toeplitz matrix analogously to the 1D case, only now the bands consist of of identical blocks of size L × L. The random variables estimating the spectrum then read: Note that the random variable Λ k,l again does not depend on the choice of i, as the random variables X ± i and Y ± i have identical distributions and once more we write it explicitly only to specify independent ωs that should be taken for the calculation. The maximum of Λ k,l is achieved for k = L, l = L and is equal to the maximum value of the sum 2 X + i,j X − i,j+1 + 2 Y + i,j Y − i+1,j . The variance of the eigenvalues of R in the 2D case can be estimated similarly to the 1D case.
The procedure can be further iterated for higher dimensions.