# Quantum Critical Scaling of Dirty Bosons in Two Dimensions

###### Abstract

We determine the dynamical critical exponent, , appearing at the Bose glass to superfluid transition in two dimensions by performing large scale numerical studies of two microscopically different quantum models within the universality class; The hard-core boson model and the quantum rotor (soft core) model, both subject to strong on-site disorder. By performing many simulations at different system size, , and inverse temperature, , close to the quantum critical point, the position of the critical point and the critical exponents, , and can be determined independently of any prior assumptions of the numerical value of . This is done by a careful scaling analysis close to the critical point with a particular focus on the temperature dependence of the scaling functions. For the hard-core boson model we find , and with a critical field of , while for the quantum rotor model we find , and with a critical hopping parameter of . In both cases do we find a correlation length exponent consistent with , saturating the bound as well as a value of significantly larger than previous studies, and for the quantum rotor model consistent with .

Most familiar quantum critical points (QCP’s) are characterized by Lorentz invariance implying a symmetry between correlations in space and time and consequently also between the respective correlation lengths . In turn, this implies that the dynamical critical exponent, defined through , is simply although and typically differ by an overall pre-factor. Systems for which are comparatively less common and possess a distinct anisotropy between space and time since and are not only different but now also scale differently. One model for which it is generally believed that is the Bose glass to super fluid (BG-SF) transition describing interacting bosons subject to disorder, the so called dirty-boson problem, modelled by the hamiltonian:

(1) |

Here , and are the boson creation and annihilation operators at site with the corresponding number operator. The parameters of the model are the hopping constant , Hubbard repulsion , and site-dependent chemical potential , inducing the disorder.

Experimental setups emulating dirty boson physics include optical lattices Greiner et al. (2002); *Pasienski10 adsorbed Helium in random media Crowell et al. (1997), Josephson-junction arrays van der Zant et al. (1992), thin-film superconductors Haviland et al. (1989); *Jaeger89; *Liu91; *Markovic98; *Lin2011 as well as quantum magnets such as doped-DTN Yu et al. (2010); *Rong12nature; *Rong12. For recent reviews see Ref. Zheludev and Roscilde (2013); Weichman (2008).

The dynamical critical exponent, , appearing at the BG-SF transition has proven exceedingly hard to determine. Initial theoretical work Fisher and Fisher (1988); *Fisher90; *Fisher89 argued that in any dimension. This has intriguing implications since it implies the absence of an upper critical dimension, although a scenario involving a discontinuous onset of mean field behavior has beeen proposed Dorogovtsev (1980); *Weichman89; *Mukhopadhyay96. Subsequent numerical studies Krauth et al. (1991); *Runge92; *Batrouni93; Sørensen et al. (1992); *Wallin94; Zhang et al. (1995); Kisker and Rieger (1997); Prokof’ev and Svistunov (2004); Pollet et al. (2009); *Gurarie09; *Carleo13; *Lin11 of the system in two dimensions were consistent with and the phase-diagram is well understood Bernardet et al. (2002); *Sengupta07; *Soyler11. More recently, the transition in three dimensions has been investigated both numerically Hitchcock and Sørensen (2006); *Gurarie09; *Yao14 and experimentally White et al. (2009), yielding evidence for . However, the majority of these numerical studies were not unbiased since implicit assumptions about are needed to fix the aspect ratio in the simulations.

The arguments leading to the equality starts with hyperscaling Stauffer et al. (1972); *Aharony74; *Hohenberg76 which states that the singular part of the free energy inside a correlation volume is a universal dimensionless number, With it follows that with a finite-size form: Fisher et al. (1989)

(2) |

Imposing a a phase gradient along one of the spatial directions will then give rise to a free energy difference where is the stiffness (superfluid density). Since must obey a form similar to Eq. (2) and since has dimension of inverse length implying , it follows that , with a finite-size scaling form of:

(3) |

Following Ref. Fisher et al., 1989, an analoguous argument imposing a twist in the temporal direction leads to and assuming . On the other hand, if , the singular part of the compressibility must from Eq. (2) obey . Assuming that , so that diverges, it follows that which leads to . With this result would then contradict the (quantum) Harris inequality Chayes et al. (1986) invalidating the initial assumption. Hence, one must have . Finally, if one argues that for the disordered system cannot vanish at criticality the relation implies .

More recent theoretical work Weichman and Mukhopadhyay (2007); *Weichmanprb08; Weichman (2008) has questioned the arguments leading to the equality . In particular, in the presence of disorder breaking particle-hole symmetry, it was argued Weichman and Mukhopadhyay (2007); *Weichmanprb08 that is dominated by the analytical background and should not apply, invalidating the relation , leaving unconstrained. A recent numerical Priyadarshee et al. (2006) study of the hard-core version of Eq. (1) finds , by analyzing scaling behavior of quantities relatively far from criticality. However, in that study, relatively few disorder realizations () were used and the location of the QCP was not reliably determined. In more recent work, Meier et al Meier and Wallin (2012) performed a state of the art calculation of a soft-core version of Eq. (1) finding a significantly larger value of and . While the location of the QCP was determined with an impressive precision this latter study does not employ a fully quantum mechanical model but instead uses an effective classical model for which the temperature dependence, at the heart of the scaling with , is only approximately accounted for. It is for instance not possible to calculate the specific heat of the underlying quantum model using the representation of Meier and Wallin (2012).

At present, the value of at the dirty-boson QCP along with many of the other exponents most notably can therefore best be regarded as ill determined, at least for the fully quantum mechanical model. It is not known to what extent, if any, the relation is violated nor if the relation Chayes et al. (1986) is obeyed, although, as outline above it seems reasonable to expect if indeed the inequality is obeyed and, as shown in Fisher et al. (1989), in any dimension. Here we try to answer these questions by performing large-scale simulations on two fully quantum mechanical models; A hard-core boson model (HCB) modelled as a transverse field model and a soft-core quantum rotor model (QR), both of which are subject to strong on-site disorder. In all cases do we use disorder realizations over a large range of temperatures extending down to for linear system sizes .

Models: The first model we study, closely related to Eq. (1), is the quantum rotor (QR) model. It is defined in terms of conjugate phase and number operators satisfying on a lattice:

(4) |

where is the onsite repulsion, is the nearest neighbor tunneling amplitude and represents the uniformly distributed on-site disorder in the chemical potential. As before, . The disorder for a given disorder realization is not constrained and in all simulations we use , and tune through the BG-SF transtion varying at constant . In contrast to Eq. (1) can take negative as well as positive values and one can loosely associtate with deviations from the average filling in Eq. (1), . For convenience we study Eq. (4) using a link-current representation Witczak-Krempa et al. (2014) for which directed worm algorithms are available Alet and Sørensen (2003a); *Aletb. We use lattice ranging from to , with 50,000 disorder realizations for and 10,000 disorder realizations for in all cases with MCS per disorder realization. For the simulations of the QR model a temporal discretization of was used, sufficiently small that remaining discretization errors could be neglected.

The second model we consider is the hard-core limit of Eq. (1) where the boson occupation number is constrained to . It is therefore therefore equivalent to the following -model on a lattice in a random transverse field:

(5) |

where uniformly. In this case we tune through the transition by increasing the disorder strength, . Despite the representation as a spin model, we shall refer to this as the hard-core boson model (HCB). We use a directed loop version of the stochastic series expansion (SSE) Sandvik and Kurkijärvi (1991); *Olav02 to simulate this model. This technique does not have discretization errors and efficient directed algorithms Syljuåsen and Sandvik (2002); Alet et al. (2005) are available. For the SSE calculations, we use a beta-doubling scheme Sandvik (2002) that allows us to very quickly equilibrate at large values. For each temperature in the beta-doubling scheme, we average over Monte Carlo sweeps (MCS), with each sweep consisting of one diagonal update and directed loop updates. is set during the equilibration phase so that on average vertices are visited, where is the number of non-trivial operators in the SSE string. In contrast to the QR model, we use here a microcanonical ensemble for the disorder by constraining each disorder realization to have exactly . This facilitates the analysis and is not believed to affect the results Dhar and Young (2003); *Monthus04. We use at least disorder realizations per data point, a large improvement over Priyadarshee et al. (2006).

In the following denotes the disorder average while denotes the thermal average. In simulations of both models two independent replicas of each disorder realization are simulated in parallel so that combined averages may be correctly estimated as .

Observables: Our main focus is the scaling behaviour of the superfluid stiffness, , for which the finite-size scaling form Eq. (3) was derived. For both models we measure as:

(6) |

where and are the winding numbers in the spatial directions. (For the HCB model Eq. (6) is multiplied by to yield .) From Eq (6) it follows that has a particularly attractive scaling form when , which we may write:

(7) |

where we define (QR model) and (HCB model). We also make extensive use of the correlation functions, defined as for the QR model and as for the HCB model.

Results, QR: A large number of independent simulations of Eq. (4) were carried out at different and close to the QCP. Since we expect to approach zero in an exponential manner as is increased at fixed and since is likely exponentially suppressed in the insulating phase it seems reasonable to approximate the function in Eq. (7) as with . If the temperature dependence is carefully mapped out Sup () one indeed sees that has a clear exponential dependence. As a first step, we then assume . We can then fit all 142 data points to this form determining the coefficients along with and . The results are shown in Fig. 1 with a scaling plot using the scaling variable . A more refined analysis Sup () shows that the temperature dependence likely involves a correction term . The correction term is here proportional to and disappears as tends to zero. It is straight forward to fit all our data to this form which yields identical estimates for , along with . Estimating the AIC (Akaike Information Criterion) for the two forms heavily favors the latter.

With a reliable estimate of we can fix the scaling argument by appropriately selecting for each . If we then study the Binder cumulant we see that at fixed it should follow a simplified form of Eq. (7), . As shown in Fig. 2, lines for different will then cross at . This is indeed the case, confirming our previous estimates.

Our results for the correlation functions for the QR models are shown in Fig. 3 for a lattice at for a range of temperatures. Asymptotically, one expects Fisher et al. (1989) and . Clearly, drops off much faster than confirming that . However, pronounced finite temperature effects are visible in arising because the limit has not yet been reached and we have not been able to reliable determine the power-law describing . However, from we determine = and hence using our previous estimate .

For the QR model we have also verified that the compressibility, , remains finite and independent of throughout the transition, consistent with . Furthermore, a direct evaluation of directly at for fixed , expected from Eq. (7) to scale as , yields consistent with our previous results.

Results, HCB: Due to the hard-core constraint number fluctuations are dramatically suppressed in the HCB model. Combined with the very effective beta-doubling scheme we can effectively reach much lower temperatures with the HCB model relative to the QR model. Hence, we use a simplified form of Eq. (3):

(8) |

suppressing the temperature dependence. We have extensively verified that this is permissible for the system sizes used Sup () and that our data appear independent of temperature at to within numerical precision. We then fit our data for at to an expansion of in Eq. (8) to second order: obtaining the estimates: , , . The result of this fit is shown in Fig. 4, where our simulation data (unfilled markers) is superimposed with the quadratic fit (dotted lines). In panel of Fig. 4, we show the crossing of the scaling function , at . By including all results for it is also possible to perform an identical analysis to the one performed for the QR model in Fig. 1. Such an analysis similar results for , and .

As was the case for the QR model the correlation functions show a pronounced temperature dependence as shown in Fig. 5(a) for . However, as we lower the temperature, reaches a stable power-law form at for all studied showing that the regime is reached for all and confirming that the temperature dependence can be neglected in Eq. (8). To determine the anomalous dimension we then fit the results in Fig. 5(a) for , to a power-law form with = as shown in Fig. 5(b). Using our earlier estimate of , we obtain in reasonable agreement with the result for the QR model.

For the HCB model we have also calculated the compressibility, . It remains roughly constant and independent of through the transition.

Conclusion: Our results for for both models studied indicate clearly that is satisfied as an equality. For the dynamical critical exponent , describing the BG-SF transition, we find a value that is significantly larger than previous estimates. While there is a slight disagreement in the estimate of for the two models we studied it seems possible that indeed . In light of this it now seems particularly interesting to focus attention on the transition in .

During the final stages of writing this manuscript we became aware of Ref. Álvarez Zúñiga et al., 2014 which for the HCB model reach conclusions similar to ours.

## I Acknowledgements

###### Acknowledgements.

We would like to thank Fabien Alet and Rong Yu for useful discussions. This work was supported by NSERC and made possible by the facilities of the Shared Hierarchical Academic Research Computing Network (SHARCNET:www.sharcnet.ca) and Compute/Calcul Canada.## References

- Greiner et al. (2002) M. Greiner, O. Mandel, T. Esslinger, T. W. Hansch, and I. Bloch, Nature 415, 39 (2002).
- Pasienski et al. (2010) M. Pasienski, D. McKay, M. White, and B. DeMarco, Nature Physics 6, 677 (2010).
- Crowell et al. (1997) P. A. Crowell, F. W. Van Keuls, and J. D. Reppy, Phys. Rev. B 55, 12620 (1997).
- van der Zant et al. (1992) H. S. J. van der Zant, F. C. Fritschy, W. J. Elion, L. J. Geerligs, and J. E. Mooij, Phys. Rev. Lett. 69, 2971 (1992).
- Haviland et al. (1989) D. B. Haviland, Y. Liu, and A. M. Goldman, Phys. Rev. Lett. 62, 2180 (1989).
- Jaeger et al. (1989) H. M. Jaeger, D. B. Haviland, B. G. Orr, and A. M. Goldman, Phys. Rev. B 40, 182 (1989).
- Liu et al. (1991) Y. Liu, K. A. McGreer, B. Nease, D. B. Haviland, G. Martinez, J. W. Halley, and A. M. Goldman, Phys. Rev. Lett. 67, 2068 (1991).
- Marković et al. (1998) N. Marković, C. Christiansen, and A. M. Goldman, Physical Review Letters 81, 5217 (1998).
- Lin and Goldman (2011) Y.-H. Lin and A. M. Goldman, Physical Review Letters 106, 127003 (2011).
- Yu et al. (2010) R. Yu, S. Haas, and T. Roscilde, EPL (Europhysics Letters) 89, 10009 (2010).
- Yu et al. (2012a) R. Yu, L. Yin, N. S. Sullivan, J. S. Xia, C. Huan, A. Paduan-Filho, O. J. N. F., S. Haas, A. Steppke, C. F. Miclea, F. Weickert, R. Movshovich, E.-D. Mun, B. L. Scott, V. S. Zapf, and T. Roscilde, Nature 489, 379 (2012a).
- Yu et al. (2012b) R. Yu, C. Miclea, F. Weickert, R. Movshovich, A. Paduan-Filho, V. Zapf, and T. Roscilde, Phys. Rev. B 86, 134421 (2012b).
- Zheludev and Roscilde (2013) A. Zheludev and T. Roscilde, Comptes Rendus Physique 14, 740 (2013).
- Weichman (2008) P. Weichman, Modern Physics Letters B 22, 2623 (2008).
- Fisher and Fisher (1988) D. S. Fisher and M. P. A. Fisher, Phys. Rev. Lett. 61, 1847 (1988).
- Fisher et al. (1990) M. P. A. Fisher, G. Grinstein, and S. M. Girvin, Phys. Rev. Lett. 64, 587 (1990).
- Fisher et al. (1989) M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
- Dorogovtsev (1980) S. Dorogovtsev, Physics Letters A 76, 169 (1980).
- Weichman and Kim (1989) P. B. Weichman and K. Kim, Phys. Rev. B 40, 813 (1989).
- Mukhopadhyay and Weichman (1996) R. Mukhopadhyay and P. B. Weichman, Physical Review Letters 76, 2977 (1996).
- Krauth et al. (1991) W. Krauth, N. Trivedi, and D. Ceperley, Phys. Rev. Lett. 67, 2307 (1991).
- Runge (1992) K. J. Runge, Phys. Rev. B 45, 13136 (1992).
- Batrouni et al. (1993) G. G. Batrouni, B. Larson, R. T. Scalettar, J. Tobochnik, and J. Wang, Phys. Rev. B 48, 9628 (1993).
- Sørensen et al. (1992) E. Sørensen, M. Wallin, S. M. Girvin, and A. P. Young, Physical Review Letters 69, 828 (1992).
- Wallin et al. (1994) M. Wallin, E. S. Sørensen, S. M. Girvin, and A. P. Young, Phys. Rev. B 49, 12115 (1994).
- Zhang et al. (1995) S. Zhang, N. Kawashima, J. Carlson, and J. Gubernatis, Phys. Rev. Lett. 74, 1500 (1995).
- Kisker and Rieger (1997) J. Kisker and H. Rieger, Phys. Rev. B 55, R11981 (1997).
- Prokof’ev and Svistunov (2004) N. V. Prokof’ev and B. Svistunov, Physical Review Letters 92, 015703 (2004).
- Pollet et al. (2009) L. Pollet, N. V. Prokof’ev, B. Svistunov, and M. Troyer, Physical Review Letters 103, 140402 (2009).
- Gurarie et al. (2009) V. Gurarie, L. Pollet, N. V. Prokof’ev, B. V. Svistunov, and M. Troyer, Physical Review B 80, 214519 (2009).
- Carleo et al. (2013) G. Carleo, G. Boéris, M. Holzmann, and L. Sanchez-Palencia, Physical Review Letters 111, 050406 (2013).
- Lin et al. (2011) F. Lin, E. Sørensen, and D. Ceperley, Phys. Rev. B 84, 094507 (2011).
- Bernardet et al. (2002) K. Bernardet, G. Batrouni, and M. Troyer, Phys. Rev. B 66, 054520 (2002).
- Sengupta and Haas (2007) P. Sengupta and S. Haas, Phys. Rev. Lett. 99, 050403 (2007).
- Söyler et al. (2011) Ş. G. Söyler, M. Kiselev, N. V. Prokof’ev, and B. V. Svistunov, Physical Review Letters 107, 185301 (2011).
- Hitchcock and Sørensen (2006) P. Hitchcock and E. Sørensen, Physical Review B 73, 174523 (2006).
- Yao et al. (2014) Z. Yao, K. P. C. da Costa, M. Kiselev, and N. V. Prokof’ev, Physical Review Letters 112, 225301 (2014).
- White et al. (2009) M. White, M. Pasienski, D. McKay, S. Q. Zhou, D. Ceperley, and B. DeMarco, Phys. Rev. Lett. 102, 055301 (2009).
- Stauffer et al. (1972) D. Stauffer, M. Ferer, and M. Wortis, Phys. Rev. Lett. 29, 345 (1972).
- Aharony (1974) A. Aharony, Phys. Rev. B 9, 2107 (1974).
- Hohenberg et al. (1976) P. C. Hohenberg, A. Aharony, B. I. Halperin, and E. D. Siggia, Phys. Rev. B 13, 2986 (1976).
- Chayes et al. (1986) J. T. Chayes, L. Chayes, D. Fisher, and T. Spencer, Physical Review Letters 57, 2999 (1986).
- Weichman and Mukhopadhyay (2007) P. B. Weichman and R. Mukhopadhyay, Phys. Rev. Lett. 98, 245701 (2007).
- Weichman and Mukhopadhyay (2008) P. Weichman and R. Mukhopadhyay, Physical Review B (2008).
- Priyadarshee et al. (2006) A. Priyadarshee, S. Chandrasekharan, J.-W. Lee, and H. U. Baranger, Phys. Rev. Lett. 97, 115703 (2006).
- Meier and Wallin (2012) H. Meier and M. Wallin, Phys. Rev. Lett. 108, 055701 (2012).
- Witczak-Krempa et al. (2014) W. Witczak-Krempa, E. S. Sørensen, and S. Sachdev, Nature Physics 10, 361 (2014).
- Alet and Sørensen (2003a) F. Alet and E. S. Sørensen, Phys. Rev. E 67, 015701 (2003a).
- Alet and Sørensen (2003b) F. Alet and E. S. Sørensen, Phys. Rev. E 68, 026702 (2003b).
- Sandvik and Kurkijärvi (1991) A. W. Sandvik and J. Kurkijärvi, Phys. Rev. B 43, 5950 (1991).
- Syljuåsen and Sandvik (2002) O. F. Syljuåsen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002).
- Alet et al. (2005) F. Alet, S. Wessel, and M. Troyer, Phys. Rev. E 71, 036706 (2005).
- Sandvik (2002) A. W. Sandvik, Phys. Rev. B 66, 024418 (2002).
- Dhar and Young (2003) A. Dhar and A. Young, Physical review. B, Condensed matter 68, 134441 (2003).
- Monthus (2004) C. Monthus, Physical Review B 69, 054431 (2004).
- (56) See supplementary material.
- Álvarez Zúñiga et al. (2014) J. P. Álvarez Zúñiga, D. J. Luitz, G. Lemarié, and N. Laflorencie, preprint (2014), arXiv:1412.5595 .

## Appendix A Supplementary material

### a.1 HCB simulation details

A central assumption made in our analysis of the HCB results was that a simplified scaling function, Eq. (8), for . In the main text, we showed the convergence of the equal-time spatial correlation function: in Fig. 5 at the critical point using data from the beta-doubling procedure. For completeness, we show the convergence of the stiffness in Fig. 6 for a lattice. Unlike the system, we need to go up to before finite temperature effects are eliminated. The data indicates that for the simulation of even larger system sizes, one may need to go up to to validate the use of the simplified scaling form.

Another important point that is often overlooked, is the convergence of observable data with the number of disorder realizations, used. For the SSE parameters chosen, we find that it was necessary to use at least disorder samples before disorder fluctuations are reasonably small. This is demonstrated in Fig. 7 for the largest lattice size: . For reliable data with controlled errorbars, all the HCB data points were averaged over least independent disorder realizations. This is contrast to earlier studies Priyadarshee et al. (2006) most relevant to our work, where only disorder realizations were used. We note that it is quite unlikely that self-averaging applies in this model and increasing the lattice size does therefore not decrease the number of disorder realizations needed.

## Appendix B Temperature dependence of

Insight into the temperature dependence of can be gained by first studying for the QR model without disorder. A model for which it is know that . Results for for a lattice are shown in Fig. 8 for . For we expect to go zero in an exponential manner while for it should approach a constant. The simplest ansatz is therefore:

(9) |

where we have tentatively included an dependence. However, as is clearly evident in Fig. 8, has a maximum close to . The presence of this maximum signals that there are likely two contributions to describing the and regimes. (Although we note that the existence of two terms does not imply a maximum.) We therefore assume the presence of a correction term proportional to with . Such a term will therefore disappear in the zero temperature limit. For our final ansatz we therefore take:

(10) |

This is the form used in Fig. 8 and it gives an essentially perfect fit over the entire range of the figure. We can immediately generalize to a scaling form for :

(11) |

A fit to this form yields in relative good agreement with similar correction terms used in Ref. Witczak-Krempa et al. (2014).

We now turn to the QR model in the presence of disorder. In this case we assume that the scaling variable generalizes to . We therefore expect to find:

(12) |

In Fig. 9 results are shown for for two different lattice sizes Assuming in this case we plot the results against demonstrating that the results fall on a single curve with only slight deviations from a straight line. It is perhaps surprising that it is the variable that appears as opposed to but this can very clearly be verified from the simulations. Performing a fit to the ansatz we find exceedingly good agreement with the expected form with a correction exponent , close to the value for the model without disorder. An inspection of our results for for this system size shows that in this case there is no maximum in versus .

We have performed a similar analysis of as a function of for the QR model again clearly confirming the overall exponential dependence and the presence of the two terms.

## Appendix C Crossing of for the QR model.

With the and determined from the fit in Fig. 1 plotted for different at a fixed aspect ratio should also cross at the critical point . This is indeed the case and is shown Fig. 10. Since this data is essentially already shown in Fig 1(a) we have in the main text opted to show the crossing using the related quantity shown in Fig. 2.