# Breakdown of self-averaging in the Bose glass

###### Abstract

We study the square-lattice Bose-Hubbard model with bounded random on-site energies at zero temperature. Starting from a dual representation obtained from a strong-coupling expansion around the atomic limit, we employ a real-space block decimation scheme. This approach is non-perturbative in the disorder and enables us to study the renormalization-group flow of the induced random-mass distribution. In both insulating phases, the Mott insulator and the Bose glass, the average mass diverges, signaling short range superfluid correlations. The relative variance of the mass distribution distinguishes the two phases, renormalizing to zero in the Mott insulator and diverging in the Bose glass. Negative mass values in the tail of the distribution indicate the presence of rare superfluid regions in the Bose glass. The breakdown of self-averaging is evidenced by the divergent relative variance and increasingly non-Gaussian distributions. We determine an explicit phase boundary between the Mott insulator and Bose glass.

###### pacs:

05.30.Jp, 64.70.P- 64.60.ae, 64.70.Tg,## I Introduction

Introducing quenched disorder into an otherwise pure system can lead to subtle and complex results and the disordered Bose-Hubbard (BH) model is no exception. While the clean BH model shows a relatively straightforward bosonic competition between repulsion and tunneling, the disordered model exhibits a new gapless insulating phase, the Bose glass (BG), the precise location of which has proven problematic from the outset.Fisher et al. (1989) With the advent of experimental methods that can engineer this model directly,Greiner et al. (2002) the problem has been inverted and this has sparked renewed interest in the role disorder plays in quantum systems.

In this Article, we develop a non-perturbative method to probe the nature of the transition between the localization-induced BG and the Mott insulator (MI) in the two-dimensional BH model with bounded potential disorder. The MI arises from on-site repulsions and hence dominates in the limit the hopping vanishes while the superfluid (SF) is the ground state in the opposite regime. It is in the difficult intermediate parameter space where the BG phase obtains. It has been argued by various authors Fisher et al. (1989); Kisker and Rieger (1997); Pollet et al. (2009); Gurarie et al. (2009); Iyer et al. (2012); Niederle and Rieger (2013) that the BG is a quantum Griffiths phase dominated by arbitrarily large SF regions that are, however, exponentially suppressed. Despite the abundance of numerical Scalettar et al. (1991); Krauth et al. (1991); Pai et al. (1996); Kisker and Rieger (1997); Sen et al. (2001); Lee et al. (2001); Prokof’ev and Svistunov (2004) and analytical Singh and Rokhsar (1992); Wu and Phillips (2008); Bissbort and Hofstetter (2009); Mukhopadhyay and Weichman (1996); Freericks and Monien (1996); Svistunov (1996); Herbut (1997, 1998); Weichman and Mukhopadhyay (2008); Krüger et al. (2009) work on the subject, it is only recently that several aspects of this model have been fully understood. This includes the confirmation that the BG always intervenes between MI and SF phases Pollet et al. (2009) (Fig. 1), the proof that the transition between the MI and the BG has to be of the Griffiths type,Gurarie et al. (2009) and the distinction between the MI and BG regarding whether fluctuations are self-averaging. Krüger et al. (2011)

As argued by Aharony and Harris,Aharony and Harris (1996) the breakdown of self-averaging can be identified from the renormalization-group (RG) flow of the relative variance of any extensive variable. If the relative variance does not renormalize to zero, the central limit theorem no longer applies and the system is not self-averaging. This concept has been used to characterize the phase transition between the MI and the BG Krüger et al. (2011) within a disorder averaged replica field theory. In both insulating phases, the mass of the theory diverges, signaling the presence of short-ranged SF correlations. In dimensions the variance of the mass distribution diverges as well and as a consequence, the breakdown of self-averaging can be readily understood as a competition between the spread of the distribution versus the shift of its average. In the MI, the shift dominates the spread leading to a vanishing relative variance. In the BG, the spread is faster and the relative variance diverges. This characterization of the MI/BG transition is consistent with the picture that the BG is dominated by rare SF regions since the negative mass values occur in the tail of the distribution. Whether the onset of this Griffiths instability Pollet et al. (2009); Gurarie et al. (2009) is correctly described by a perturbative RG calculation remains a central question.

We address this question within a real-space block decimation RG scheme on a square lattice. This approach is non-perturbative in the disorder and hence enables a study of the RG flow of the full random-mass distribution, a necessity for any definitive statement about Griffiths-type physics to be made. We reiterate that the onset of the Bose glass is well known to be mediatedPollet et al. (2009); Gurarie et al. (2009) by Griffiths rare-region physics and hence our conclusions are independent of whether the transition is studied from the Mott insulator or the superfluid. Our results confirm that the relative variance serves as the order parameter for the MI-to-BG transition. Determining the correlation length from the scale at which the relative variance becomes of order one, we extract a correlation-length exponent of about . This is close to the analytical value obtained within the perturbative 1-loop RG.Krüger et al. (2011) It has been argued Aharony and Harris (1996); Pázmándi et al. (1997) that a violation of the Harris-criterion bound for critical disordered systems Harris (1974); Chayes et al. (1986) is indicative of the lack of self-averaging. The absence of the central-limit theorem in the BG is further evidenced by an increasingly non-Gaussian shape of disorder distributions.

## Ii Model

Our starting point is the simplest form of the disordered BH model on a square lattice,

(1) |

which describes bosons tunneling with amplitude between nearest-neighbor sites and and interacting via an on-site repulsion . The bosonic raising and lowering operators are given by respectively where is the bosonic number operator. is the chemical potential shifted by the on-site disorder potential, . The random site energies are uncorrelated between different sites and uniformly distributed in the interval . From minimization of the energy in the atomic limit it is straightforward to see that for the phase diagram retains MI phases (see Fig. 1).

To facilitate a strong coupling expansion, we follow the standard procedure.Sachdev (2011); Sengupta and Dupuis (2005) After expressing the model by a coherent-state path integral in imaginary time, we decouple the hopping term by a Hubbard-Stratonovich transformation and trace over the original boson fields. We then perform a temporal gradient expansion to obtain the effective dual action on a lattice,

(2) | |||||

where denotes the lattice spacing and the sum over runs over the nearest neighbors of site . By construction, the complex fields correspond to the SF order parameter . In regions where , SF order is suppressed. The mass therefore corresponds to the local Mott gap. and are related to the single- and two-particle bosonic Green functions of the local on-site Hamiltonian, respectively, while the temporal gradient terms and are given by derivatives of the mass with respect to the chemical potential. We specialize to the first Mott lobe with bosons per site in which the coefficients are given by Sengupta and Dupuis (2005); Krüger et al. (2009)

(3a) | |||||

(3b) | |||||

(3c) | |||||

where is the coordination number of the square lattice. In the clean limit, , the mean-field phase boundary between the first MI lobe and the SF is obtained from . In the presence of disorder, , the coefficients , , , and depend on the disorder potential , which induces non-trivial disorder distributions of the coefficients. Note that initially the dual hopping amplitudes are uniform. We allow for a spatial dependence of the hopping since disorder will be induced under block decimation.

## Iii Block Decimation Real-Space RG

Equations (3) provide the initial conditions for our procedure where the random on-site energies are generated from a uniform distribution on the interval . In order to determine the phase for a given set of parameters , , , and of the disordered BH model (1), we derive a set of recursion equations using block decimation. We start with the discrete action (2) and eliminate short-range degrees of freedom by integrating out every other site, treating the quartic terms perturbatively to leading order. The remaining points form a new square lattice with lattice spacing and tilted from the original system. Recollecting the resulting terms and rescaling the action to look like the original, we find the RG recursion equations

(4a) | |||||

(4b) |

where is the static propagator and . The indices and correspond to nearest neighbors of site , whereas and correspond to the bonds adjacent to the bond connecting the sites and of the remaining lattice. The site is the common vertex of the bonds and (see Fig. 2). Since we are interested in the MI/BG transition at incommensurate filling, we can neglect any corrections to the coefficients , , and beyond dimensional scaling. Note that under the RG longer range couplings are generated which is a known problem of the block decimation method in . In the present case, however, locality is guaranteed since the mean of the mass distribution diverges in both insulating phase, leading to an exponential suppression of hopping amplitudes beyond nearest neighbors.

Since the gradient terms are renormalized under block decimation the effective local mass should be defined relative to the kinetic hopping terms,

(5) |

As a consistency check, we evaluate our recursion relations in the clean limit. For , , and , the RG equations (4) reduce to and , leading to the recursion relation for the effective mass. This indeed correctly describes the mean-field transition between the MI and the SF at .

## Iv Results

In the following, we integrate the RG equations (4) numerically for the inhomogeneous system obtained for one particular disorder realization and keep track of the values of the coupling constants on each lattice site. This allows us to extract the mass distribution at each iteration step of the block decimation. We vary while keeping the disorder fixed at for chemical potential values and several lattice constants . Since Eq. (2) has only nearest neighbor or on-site terms, we can ignore the boundary points after each iteration without affecting the overall distribution. The major limitation to this method is the necessity of finite size lattices. Estimates of any diverging quantities near the critical point must take into account finite system size effects. In the data given below we use an initial square grid of points with side length sites. This side length is not a power of two because each decimation step concludes by throwing away points affected by the boundary.

As expected for insulating phases, in both the MI and the BG, the mean of the mass distribution increases exponentially under the RG signaling short-range SF correlations. To distinguish the behavior in the MI and the BG, we normalize to a mean of unity and analyze the evolution of . The variance of this rescaled distribution corresponds to the relative variance of the mass distribution, which should serve as order parameter.Krüger et al. (2011)

In Fig. 3a, the RG flow of in the MI phase is shown. Note that the initial mass distribution is asymmetric and non-Gaussian due to the functional dependence (3a) on the uniformly distributed on-site energies . As a consequence of the relatively small hopping value , the bulk of the sites begin well above and continue together towards larger values under repeated iteration. Increasing the number of iterations results in a distribution well described by a Gaussian. Further, the width of the distribution narrows, indicating a vanishing relative variance. This demonstrates that in the MI disorder is irrelevant and the system is self-averaging.

The situation changes dramatically as we increase the value of and enter the BG phase (see Fig. 3). While the overall shape of the initial distribution looks quite similar to the one in the MI, now a large fraction of the initial sites lies close to or below . With enough sites close to the transition point, some regions of the system take much longer to flow to large positive values than others. This results in a drastic spread of values upon repeated iteration of the RG, leading to a divergence of the relative variance. In addition, the distribution develops a non-Gaussian form under successive iterations, indicating a violation of the central limit theorem. This demonstrates a breakdown of self-averaging in the BG.

Our results show that the relative variance of the random-mass distribution serves as an order parameter for the transition between the MI and the BG. We can therefore determine a correlation length in the BG from the average number of iterations it takes before the relative variance becomes of order unity. For the parameters used in Fig. 3b this happens between 6 and 7 iterations. A more precise value is obtained by averaging over several disorder realizations. Note that once the relative variance becomes of order unity, the left tail of the distribution pushes through zero. Therefore, the correlation length corresponds to the typical distance between SF droplets in the system.

We generate a phase diagram for the transition between the MI and BG by estimating values for such that the data for several initial conditions collapse onto a single curve. This method, and consequently the resulting phase diagram, can only predict the relative value of between different values of and not the absolute location of . Therefore, we set to be a small constant value for initial conditions in which the BG is suppressed and plot the results for intermediate values in Fig. 4. Other omputational methods have produced a phase boundary on the same order of magnitude with the same characteristic shape Niederle and Rieger (2013).

At the transition to the MI, the correlation length diverges as a power law, . To extract the correlation-length exponent , we vary the hopping slightly above the transition point for fixed values of the lattice spacing and several values of the chemical potential and extract the correlation length as described above by averaging over several disorder realizations. In the following, we use , and average over 10-20 disorder realizations. Note that the value of the transition point is non-universal. With the values of in Fig. 4, we plot vs for various values of . The data collapse onto a single curve with an asymptotic value of as shown in Fig. 5. Values of the chemical potential near commensurate density require larger system sizes to see universal behavior as was found in previous momentum-shell work.Krüger et al. (2011) This is indicated in Fig. 5 where the mid-range values of are further from the universal asymptotic fit than those at the extrema. As a result, the system sizes used here of are too small to obtain enough data to estimate for filling values slightly less than commensurate.

## V Discussion and Conclusion

We used a real-space block decimation method to characterize the MI-to-BG transition. By analyzing the RG flow of the induced random-mass distribution in the dual field theory, we have demonstrated that the transition is characterized by a breakdown of self-averaging. The associated correlation length corresponds to the typical separation of rare SF regions. Note that while the transition is not related to spontaneous symmetry breaking, the divergence of the correlation length defined here has been shown to be related to a vanishing compressibility.Krüger et al. (2011) Our work provides an explicit confirmation that the instability of the MI towards the formation of the BG is of the Griffith type, as previously argued by other authors.Pollet et al. (2009); Gurarie et al. (2009)

The method employed here enables us to study the RG flow of entire disorder distributions, whereas the perturbative 1-loop momentum-shell RG based on the disorder averaged replica theory is restricted to the mean and the variance of the random mass distribution.Krüger et al. (2011) Both approaches, however, show that the relative variance diverges in the BG. This is an important result as it demonstrates that contrary to the general belief, the onset of Griffiths instabilities is captured in perturbative RG. The comparison of the correlation-length exponents obtained by the two methods suggests that corrections beyond 1-loop order are small.

Griffiths phases can been classified based on the comparison of the dimension of defects (or rare regions) with the lower critical dimension.Vojta and Schmalian (2004); Vojta (2006); Hoyos and Vojta (2008); Mohan et al. (2010) Since the SF droplets are one-dimensional – rod-like in imaginary time and zero dimensional in space – the rare regions are below the lower critical dimension of the problem, leading to weak Griffiths singularities characterized by an essential singularity in the free energy.cla ()

Our results show that the MI/BG transition is characterized by a fixed point with finite relative variance. This is consistent with strong disorder RG calculations that show that the transition between the SF and the disordered insulator is governed by a finite disorder fixed point.Altman et al. (2004, 2010); Iyer et al. (2012) While the strong disorder RG approach becomes asymptotically exact in the limit of infinite disorder, it might produce unphysical results in the regime of weak disorder.Iyer et al. (2012) Another difference to the block decimation scheme is that the lattice coordination is not preserved in dimensions . It would be interesting to systematically compare the evolution of disorder distributions for the two complementary methods.

The real-space RG approach presented here has a wide range of future applications. It can be used to study the effects of spatial correlations in the disorder and entails the search for self-similar disorder characterized by scale invariant distribution functions. Finally, the method is not restricted to disorder distributions, but can be used to study other inhomogeneities such as the so-called wedding cake structures of alternating MI and SF regions found in optical-lattice systems.Gemelke et al. (2009)

Acknowledgment: The authors benefited from stimulating discussions with B. Brinkman, A. V. Chubukov, I. F. Herbut, J. A. Hoyos, J. Keeling, R. Moessner, and H. Rieger. This work was supported by EPSRC under grant code EP/I 004831/1 and by NSF-DMR-1104909.

## References

- Fisher et al. (1989) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
- Greiner et al. (2002) M. Greiner, O. Mandel, T. Hänsch, and I. Bloch, Nature 419, 51 (2002).
- Kisker and Rieger (1997) J. Kisker and H. Rieger, Phys. Rev. B 55, R11981 (1997).
- Pollet et al. (2009) L. Pollet, N. V. Prokof’ev, B. V. Svistunov, and M. Troyer, Phys. Rev. Lett. 103, 140402 (2009).
- Gurarie et al. (2009) V. Gurarie, L. Pollet, N. V. Prokofev, B. V. Svistunov, and M. Troyer, Phys. Rev. B 80, 214519 (2009).
- Iyer et al. (2012) S. Iyer, D. Pekker, and G. Rafael, Phys. Rev. B 85, 094202 (2012).
- Niederle and Rieger (2013) A. E. Niederle and H. Rieger, New Journal of Physics 15, 075029 (2013).
- Scalettar et al. (1991) R. T. Scalettar, G. G. Batrouni, and G. T. Zimanyi, Phys. Rev. Lett. 66, 3144 (1991).
- Krauth et al. (1991) W. Krauth, N. Trivedi, and D. Ceperley, Phys. Rev. Lett. 67, 2307 (1991).
- Pai et al. (1996) R. V. Pai, R. Pandit, H. R. Krishnamurthy, and S. Ramasesha, Phys. Rev. Lett. 76, 2937 (1996).
- Sen et al. (2001) P. Sen, N. Trivedi, and D. M. Ceperley, Phys. Rev. Lett. 86, 4092 (2001).
- Lee et al. (2001) J.-W. Lee, M.-C. Cha, and D. Kim, Phys. Rev. Lett. 87, 247006 (2001).
- Prokof’ev and Svistunov (2004) N. Prokof’ev and B. Svistunov, Phys. Rev. Lett. 92, 015703 (2004).
- Singh and Rokhsar (1992) K. G. Singh and D. S. Rokhsar, Phys. Rev. B 46, 3002 (1992).
- Wu and Phillips (2008) J. Wu and P. Phillips, Phys. Rev. B 78, 014515 (2008).
- Bissbort and Hofstetter (2009) U. Bissbort and W. Hofstetter, Europhys. Lett. 86, 50007 (2009).
- Mukhopadhyay and Weichman (1996) R. Mukhopadhyay and P. B. Weichman, Phys. Rev. Lett. 76, 2977 (1996).
- Freericks and Monien (1996) J. K. Freericks and H. Monien, Phys. Rev. B 53, 2691 (1996).
- Svistunov (1996) B. V. Svistunov, Phys. Rev. B 54, 16131 (1996).
- Herbut (1997) I. F. Herbut, Phys. Rev. Lett. 79, 3502 (1997).
- Herbut (1998) I. F. Herbut, Phys. Rev. B 57, 13729 (1998).
- Weichman and Mukhopadhyay (2008) P. B. Weichman and R. Mukhopadhyay, Phys. Rev. B 77, 214516 (2008).
- Krüger et al. (2009) F. Krüger, J. Wu, and P. Phillips, Phys. Rev. B 80, 094526 (2009).
- Krüger et al. (2011) F. Krüger, S. Hong, and P. Phillips, Phys. Rev. B 84, 115118 (2011).
- Aharony and Harris (1996) A. Aharony and A. B. Harris, Phys. Rev. Lett. 77, 3700 (1996).
- Pázmándi et al. (1997) F. Pázmándi, R. T. Scalettar, and G. T. Zimányi, Phys. Rev. Lett. 79, 5130 (1997).
- Harris (1974) A. B. Harris, J. Phys. C. 7, 1671 (1974).
- Chayes et al. (1986) J. T. Chayes, L. Chayes, D. S. Fisher, and T. Spencer, Phys. Rev. Lett. 57, 2999 (1986).
- Sachdev (2011) S. Sachdev, Quantum phase transitions (Cambridge University Press, Cambridge New York, 2011), ISBN 0521514681.
- Sengupta and Dupuis (2005) K. Sengupta and N. Dupuis, Phys. Rev. A 71, 033629 (2005).
- Vojta and Schmalian (2004) T. Vojta and J. Schmalian, Phys. Rev. B 72, 045438 (2004).
- Vojta (2006) T. Vojta, J. Phys. A 39, R143 (2006).
- Hoyos and Vojta (2008) J. A. Hoyos and T. Vojta, Phys. Rev. Lett. 100, 240601 (2008).
- Mohan et al. (2010) P. Mohan, R. Narayanan, and T. Vojta, Phys. Rev. B 81, 144407 (2010).
- (35) Note that such weak Griffiths singularities are sometimes called ‘classical’, see e.g. Ref. [Vojta and Schmalian, 2004]. Although this case is indeed generic to classical systems with point disorder, the terminology ‘classical’ is misleading since the BG is a phase.
- Altman et al. (2004) E. Altman, Y. Kafri, A. Polkovnikov, and G. Rafael, Phys. Rev. Lett. 93, 150402 (2004).
- Altman et al. (2010) E. Altman, Y. Kafri, A. Polkovnikov, and G. Rafael, Phys. Rev. B 81, 174528 (2010).
- Gemelke et al. (2009) N. Gemelke, X. Zhang, C.-L. Hung, and C. Chin, Nature 460, 995 (2009).