Ground state energy density, susceptibility, and Wilson ratio of a two-dimensional disordered quantum spin system
GGround state energy density, susceptibility, and Wilson ratio of a two-dimensionaldisordered quantum spin system
J.-H. Peng, D.-R. Tan, and F.-J. Jiang ∗ Department of Physics, National Taiwan Normal University, 88, Sec.4, Ting-Chou Rd., Taipei 116, Taiwan
A two-dimensional (2D) spin-1/2 antiferromagnetic Heisenberg model with a specific kind ofquenched disorder is investigated, using the first principles nonperturbative quantum Monte Carlocalculations (QMC). The employed disorder distribution has a tunable parameter p which can beconsidered as a measure of the corresponding randomness. In particular, when p = 0, the disorderedsystem becomes the clean one. Through a large scale QMC, the dynamic critical exponents z , theground state energy densities E , as well as the Wilson ratios W of various p are determined with highprecision. Interestingly, we find that the p dependence of z and W are likely to be complementaryto each other. For instance, while the z of 0 . ≤ p ≤ . z = 1 which corresponds to the clean system, the W for p < . p = 0. The technical subtlety of calculating these physicalquantities for a disordered system is demonstrated as well. The results presented here are notonly interesting from a theoretical perspective, but also can serve as benchmarks for future relatedstudies. PACS numbers:
I. INTRODUCTION
Spatial dimension two is extraordinary from a theo-retical point of view. This is because according to thefamous Mermin-Wagner theorem, for finite systems withshort-range interactions, continuous symmetries cannotbe broken spontaneously at any temperature
T > T ) region above a QCP is wherethese profound characteristics can be uncovered the mostclearly [14, 19].A physical observable, namely the spinwave velocity c plays an important role in those mentioned universalquantities of QCR for the 2D spin-1/2 antiferromagnets.For instance, the value of c without doubt has great im-pact on the determinations of two universal quantitiesof QCR, namely χ u c /T ∼ . c/ ( T ξ ) ∼ . χ u and ξ are the uniform susceptibility and the cor-relation length, respectively [10, 14, 19]. In the phasewith long-range antiferromagnetic order, c can be calcu-lated efficiently using the spatial and the temporal wind- ∗ [email protected] ing numbers squared [18, 20, 21].Considering a clean 2D spin-1/2 AF which comes witha given spatial arrangement of two types of antiferromag-netic couplings J (cid:48) and J ( J (cid:48) > J ), by tuning the ratio J (cid:48) /J (i.e. the system is dimerized) a QCP may appearwhen J (cid:48) /J exceeds a certain value ( J (cid:48) /J ) c . The dynamiccritical exponent z associated with such a kind of QCPtakes the value of 1. For a QCP g c which is obtained byvarying the associated parameter g , the physical quan-tity c scales as c ∝ ( g − g c ) ν ( z − close to g c [10, 13],where ν is the correlation length exponent. As a result, c is a constant when the related z of a QCP is 1. For 2Dquantum AF systems, the QCPs induced by dimerizationintroduce above belong exactly to this case. When dis-order is present, z > c is zero at g c . Consequently,certain universal quantities of QCR cannot be calculatedin a direct manner for disordered systems.While for a 2D disordered quantum spin antiferromag-net, certain quantities of QCR such as χ u c /T cannot bedirectly accessed, yet some observables do not encounterthe difficulty that c cannot be calculated with ease. Oneof them, namely the Wilson ratio W [10, 17, 18], whichwill be defined later, is one of the main topics of our studypresented here. In particular, we investigate the behaviorof W with respect to the strength of randomness, whichis controlled by a parameter p ≥
0, of the employed dis-order distribution. Here p = 0 corresponds to the cleancase. Apart from W , the dynamic critical exponents z as well as the ground state energy densities E of severalvalues of p considered in this study are determined aswell.To carry out the proposed investigation, we have per-formed a large scale quantum Monte Carlo calculation(QMC). In addition, several p are considered and thesimulations are done at the corresponding critical point g c ( p ) of each studied p . Based on our numerical results,we find the magnitude of E grows monotonically with a r X i v : . [ c ond - m a t . d i s - nn ] A ug g c ( p ) (hence p as well since g c ( p ) ∝ p as shown in [22]),similar to that of the correlation length exponent ν [22].Interestingly, the p dependence of z and W are likely tobe complementary to each other. For instance, while the z of 0 . ≤ p ≤ . z = 1 which correspondsto the clean system, the W for p < . p = 0 ( W ∼ . z and W are demonstrated. Finally,a section concludes our study. II. MICROSCOPIC MODELS ANDOBSERVABLES
The model investigated in our study has been de-scribed in detail in Refs. [22, 23]. Here we briefly summa-rize certain technical perspectives of the considered sys-tem. The Hamiltonian of the investigated 2D disorderedspin-1/2 herringbone Heisenberg model (on the squarelattice) is given by H = (cid:88) (cid:104) ij (cid:105) J (cid:126)S i · (cid:126)S j + (cid:88) (cid:104) i (cid:48) j (cid:48) (cid:105) J (cid:48) (cid:126)S i (cid:48) · (cid:126)S j (cid:48) , (1)where in Eq. (1) J (which are set to 1 here) and J (cid:48) arethe antiferromagnetic couplings (bonds) connecting near-est neighboring spins (cid:104) ij (cid:105) and (cid:104) i (cid:48) j (cid:48) (cid:105) , respectively, and (cid:126)S i is the spin-1/2 operator at site i . In this study we usethe convention J (cid:48) > J . Fig. 1 is a cartoon represen-tation of the considered model. The quenched disorderintroduced into the system is based on the one employedin Refs. [22, 23]. Specifically, for every bold bond infig. 1, its antiferromagnetic strength J (cid:48) takes the value of1+( g − p ) or 1+( g − − p ) with equal probability.Here g > ≤ p ≤
1. With the used conventions,the average and difference for these two types of boldbonds J (cid:48) are given by g and 2 p ( g − p can be thought of as a measure for the disor-der of the studied model as well.To perform the proposed calculations of determiningthe ground state energy density E , the dynamic criticalexponent z , and the Wilson ratio W for the considereddisordered system (with various p ), the uniform suscepti-bility χ u , the internal energy density E , and the specificheat C V (as functions of the temperature T or the inversetemperature β ) are measured in our simulations. JJ’
FIG. 1: The 2D dimerized spin-1/2 herringbone Heisenbergmodel on the square lattice investigated in this study. Theantiferromagnetic coupling strengths for the thick and thinbonds are J (cid:48) and J , respectively. The uniform susceptibility χ u is defined by χ u = βL (cid:42) (cid:32)(cid:88) i S zi (cid:33) (cid:43) , (2)where β and L are the inverse temperature and linear boxsize used in the simulations, respectively. Furthermore,the internal energy density E and the specific heat C V are given as E = 1 L (cid:104) H (cid:105) , (3) C V = ∂E∂T . (4)Using these observables, E , z , and W for various p of the studied disordered model can be determined withhigh precision. III. THE NUMERICAL RESULTS
For each of the considered values of p = 0.0, 0.2, 0.4,0.5, 0.6, 0.7, 0.8, and 0.9, to calculate the associated de-sired physical quantities, we have carried out large-scaleQMC using the stochastic series expansion (SSE) algo-rithm with very efficient operator-loop update [24, 25].The simulations are done at the corresponding criticalpoints g c ( p ) for these chosen p . In addition, for every p , several hundred to few thousand randomness config-urations with L = 128 and (or) L = 256 are generated.Each configuration is produced with its own random seedand then is used for all the calculations of the consideredvalues of β . A. The strategy of calculating the Wilson ratio W In the framework of SSE, the quantities of internal en-ergy density E and and specific heat C V can be obtainedby E = − L (cid:32) (cid:104) n (cid:105) /β − (cid:88) b J b (cid:33) , (5) C V = 1 L (cid:0) (cid:104) n (cid:105) − (cid:104) n (cid:105) − (cid:104) n (cid:105) (cid:1) , (6)respectively, where the summation is over all the bonds b and n is the number of nonidentity operators in the SSEoperators sequence (operators string).Based on the large- N expansion of the relevant effec-tive field theory, at the associated critical point it is pre-dicted that for clean systems the (leading) low- T behav-ior of χ u , E , and C V are given by χ u ∼ . πc T, (7) E ∼ E + 2 . πc T , (8) C V ∼ . πc T , (9)respectively, where the c appearing above is the spinwavevelocity. With these leading T -dependence of χ u , E and C V , the Wilson ratio W can be expressed as W = χ u TC V ∼ . . (10)While C V can be calculated directly from its definition C V = ∂E/∂T (or C V = L (cid:0) (cid:104) n (cid:105) − (cid:104) n (cid:105) − (cid:104) n (cid:105) (cid:1) ), as be-ing shown in the literature, such a approach will lead tovery noisy results at the region of low temperature [18].In addition, the fact that c is zero at the QCP of a dis-ordered system prevents one from determining c directly.Motivated by the method outlined in Ref. [18], here wecalculate W through the following procedures.Firstly, from the β -dependence of the internal energydensity E , namely E ( β ) = E + aβ − − /z (11)(here z is the dynamic critical exponent), one obtains a and z . Then the specific heat C V , as a function of β , canbe written as C V = a (1 + 2 /z ) β − /z . (12)Secondly, the β -dependence of χ u is fitted to the ex-pression χ u ( β ) = bβ − /z . (13)Finally, using Eqs. 12 and 13, one arrives at the fol-lowing formula for WW = ba (1 + 2 /z ) . (14)In other words, instead of using C V directly, here W is calculated through the coefficients z , a and b obtainedfrom fitting the data of χ u and E to their expected T -dependence ansatzes Eqs. 11 and 13. u L = 128 L = 256 - E L = 128 L = 256 FIG. 2: χ u (top) and minus (internal) energy density E (bot-tom) as functions of β for p = 0 .
0. The shown errors are thecorresponding mean errors.
B. The obtained χ u and E from simulations The obtained data of χ u and − E for p = 0 .
0, 0.4, 0.6,and 0.9 are depicted in figs. 2, 3, 4 and 5. Both data of L = 128 and L = 256 are put in these figures in orderto demonstrate that the outcomes of L = 256 are (mostlikely) sufficient for size convergence. C. The analysis results associated with the cleanmodel
For the clean model p = 0, it is well known that z = 1.Hence, z will be fixed to 1 in our analysis for p = 0. Asa result, the following ansatz χ u = a + aT + a T (15)is considered to fit the χ u data of p = 0. Apart fromthat, the formula used to fit the data of E for p = 0 isEq. 11 with z = 1 as well.By investigating the relevant data for L = 128 and L = 256, the finite size effect begins to appear when β > .
0. Therefore the data of L = 256 with β ≤ . L = 256 data of β ≤ . W which agrees u - E L = 128 L = 256 FIG. 3: χ u (top) and minus (internal) energy density E (bot-tom) as functions of β for p = 0 .
4. The shown errors are thecorresponding mean errors. quantitatively with that obtained using the L = 256 dataof β ≤ .
0. For each of data with a fixed range of β , inadditional to taking care of finite size effect, the followingprocedures are adopted to calculate the corresponding W .Firstly, a bootstrap resampling (with respect to β ) isconduct simultaneously for both χ u and E . Secondly,fits for these obtained resampled data are performed. Fi-nally W is determined using the outcomes of these fits.In these mentioned fits, Gaussian noises are consideredas well. The above described steps are carried out fortwenty thousand times and only those outcomes withboth χ /DOF (DOF stands for degrees of freedom) ofthe two fits (for χ u and E ) being smaller than 3.0 areincluded as the candidate results of W . The resulting W , as well as its associated uncertainty, quoted for thisset of data with that given fixed range of β are the meanand the standard deviation of these candidate results.We have conducted several calculations using variousrange of β , and each of these calculations comes withis own results (mean and uncertainty) for W . More-over, to estimate the means and errors of the desiredquantities appropriately, the weighted bootstrap resam-pling method is applied to all the mentioned results of W . Specifically, for every randomly generated data set( W j , σ W j ) obtained using the bootstrap procedure ( σ W j is the standard deviation associated with W j ), the result- u L = 128 L = 256 - E L = 128 L = 256 FIG. 4: χ u (top) and minus (internal) energy density E (bot-tom) as functions of β for p = 0 .
6. The shown errors are thecorresponding mean errors. ing mean is given by (cid:80) i σ Wi W i (cid:80) i σ Wi . (16)The reason for the use of above equation (calledweighted mean in this study) is as follows. Notice thatdata with large standard deviations are less accuratelydetermined than those with small standard deviations.As a result, those large standard deviation data shouldcontribute less weight to the determination of the asso-ciated mean.After carrying out twenty thousand weighted boot-strap resampling steps, the resulting W is estimated tobe W = 0 . W = 0 . W = 0 . W .The ground state energy density for the clean model iscalculated by the same procedure and is given by E = − . u L = 128 L = 256 - E L = 128 L = 256 FIG. 5: χ u (top) and minus (internal) energy density E (bot-tom) as functions of β for p = 0 .
9. The shown errors are thecorresponding mean errors.
D. The results of the disordered model withvarious randomness strength p Since each generated configuration is used for all thesimulations of the considered β , for a given set of p > L , the data themselves are correlated. Hence, to ac-curately estimate the associated errors for the coefficientsin the fitting ansatzes, one should employ the correlatedleast χ method for the analysis. However, the stabilityof the correlated least χ method varies and depends onthe quality of the data used for the fits. Moreover, biasedoutcomes may be reached if the associated covariance ma-trix for a given data set contains eigenvalues which havevery small magnitude. Using the rule of thumb that theansatz considered to fit correlated data should contain asfew (to be determined) parameters as possible, we adoptthe following approach to calculate E , z , and W for p > β the bootstrap resamplingmethod is performed for the raw data resulting from thegenerated disordered configurations. Second of all, theseresampled data are used to calculate the disordered av-erage of χ u and E which are then considered for the rel-evant fits. Here the data employed for the fits of χ u and E have different range of β . Indeed, as can be seen fromfigs. 3, 4, 5, for all the considered p > E reaches its ground state value E quickly, while this is not the case for χ u . As a result, it is more appropriateto use different range of β for the fits of χ u and E .After carrying out the fit of χ u , the obtained result of z is employed as an input for the fit of E . When both fitsof χ u and E are done, the resulting results are then putback to calculate the associated correlated χ /DOF. Herea cut-off for the eigenvalues of the associated covariancematrix is imposed in order to avoid biased results. Theseintroduced steps are performed for many times, and onlythose results which have correlated χ /DOF smaller than3.0 for both the fits of χ u and E are considered for latercalculations. Finally for each of the considered p , theabove procedures have been applied to many sets con-sisting of various range of β . Each of these sets (Thewhole sets is denoted by S ) has its own results (meanand standard deviation) of E , z , and W as well as thenumber of successful calculations.For a considered p , the final quoted results of E , z , and W in this study are estimated by a bootstrap resamplingprocedure using the following formula to calculate themean of every resampled data from S . (cid:80) i N i O i /σ O i (cid:80) i N i /σ O i , (17)where { O i } , { σ O i } and { N i } stand for the randomlypicked outcomes in S , the associated standard devia-tions, as well as the related numbers of (successful cal-culated) results of these picked outcomes, respectively.Finally, such a resampling step is conducted for severalthousand times, and the numerical values presented herefor these considered physical quantities are the result-ing means and standard deviations (estimated conserva-tively) of this procedure. The uncertainties of E calcu-lated by the described steps are much smaller than thoseof the original E contained in S . Hence, for the data in S we have also calculated their associated weighted er-rors. The dominant one of these two estimations, namelythe standard deviations and the weighted errors, are thefinal values quoted here.The E , z and W as functions of p calculated by theprocedures introduced above are shown in figs. 6, 7 ,8.The related results for the clean model are shown in thesefigures as well for comparison.The E as a function of p shows a monotonic behaviorin magnitude, which is similar to that of the correlationlength exponents ν obtained in Ref. [22].Regarding the z presented in the figure, one observesthat the magnitude of z increases with p until p reachesa specific p c < .
4. For p ≥ .
4, all the calculated val-ues of z lie between (around) 1.3 and (around) 1.4. Ifone takes into account the systematic errors due to theuncertainties of g c ( p ), then the z for p ≥ . z associated with p ≥ . L = 128 and L = 256). These guided lines justifythe claim made above. This phenomenon is consistent p E L=256L=128
FIG. 6: - E as functions of p . The results are obtained fromthe analysis using the correlated χ described in the maintext. The solid squares and solid circles are for L = 256and L = 128, respectively. For some values of p , the L =128 results contains those corresponding to g c , the lower andupper bounds of g c . with the one shown in Ref. [31], where the calculated z corresponding to various parameters take a universalvalue. We would like to point out that when conduct-ing the determination of z from χ u , the obtained resultsare somehow a little bit sensitive to the considered fit-ting range of χ u . This motives the use of the resamplingprocedures described above. In conclusion, our analysisindicates that it is subtle to calculate the quantity z anda careful strategy is needed.Finally, in fig. 8 we demonstrate the results of W asfunctions of p obtained from the analysis outlined previ-ously. Intriguingly, similar to the scenario of z , for those W corresponding to p < .
7, their values are more orless close to each other. The solid and dashed lines inthe figure again stand for the mean and standard devia-tion for all the W with their associated p satisfy p < . L = 128 and L = 256). Consid-ering the impact resulting from the errors of g c ( p ), thescenario that W take the same value (or at least valuesclose to each other) for all the p such that p < . ν is beginning to fulfill the Harris criterion when p > . W increases sharply.This observation implies that there may exist a relationbetween W and fulfillment of Harris criterion. IV. DISCUSSIONS AND CONCLUSIONS
In this study, we calculate the Wilson ratio W of a2D spin-1/2 antiferromagnetic Heisenberg model with aspecific quenched disorder, using the first principle non-perturbative quantum Monte Carlo simulations. Theemployed disorder distribution has a tunable parame-ter p which can be considered as a measure of ran-domness. The W of the clean case as well as that of p = 0 . , . , . , . , . , . , . z and the p z L = 256 L = 128 FIG. 7: z as functions of p . The results are obtained fromthe analysis using the correlated χ described in the maintext. The solid squares and solid circles are for L = 256and L = 128, respectively. For some values of p , the L =128 results contains those corresponding to g c , the lower andupper bounds of g c . The solid and dashed lines in the figurerepresent the mean and standard deviation for all the z suchthat their associated p satisfy p ≥ . L = 128 and L = 256). p W L = 256 L = 128 FIG. 8: W as functions of p . The results are obtained fromthe analysis using the correlated χ described in the maintext. The solid squares and solid circles are for L = 256and L = 128, respectively. For some values of p , the L =128 results contains those corresponding to g c , the lower andupper bounds of g c . The solid and dashed lines in the figurerepresent the mean and standard deviation for all the W suchthat their corresponding p satisfy p < . L = 128 and L = 256). ground state energy densities E are obtained as well.Remarkably, for the considered system with the em-ployed quenched disorder, the p dependence of W and z seems to be complementary to each other. The obtained z are likely to take a universal value for p ≥ .
4. Thisagrees with the outcomes determined in Ref. [31]. In ad-dition, the calculated W for 0 < p < . W ∼ . p = 0). In particular, the value of W begins toincrease sharply when p is approaching p = 0 . p E L=256L=128
FIG. 9: - E as functions of p . The results are obtained fromthe analysis using the conventional uncorrelated χ . The solidsquares and solid circles are for L = 256 and L = 128, respec-tively. For some values of p , the L = 128 results containsthose corresponding to g c , the lower and upper bounds of g c . p z L = 256 L = 128 FIG. 10: z as functions of p . The results are obtained fromthe analysis using the conventional uncorrelated χ . The solidsquares and solid circles are for L = 256 and L = 128, respec-tively. For some values of p , the L = 128 results containsthose corresponding to g c , the lower and upper bounds of g c .The solid and dashed lines in the figure represent the meanand standard deviation for all the z such that their associ-ated p satisfy p ≥ . L = 128 and L = 256). some light on setting up some useful guidelines to decidewhether the Harris criterion is valid for a given disorderdistribution.Apart from the subtlety of calculating z described pre-viously, the determination of W is extremely non-trivialas well. Indeed, the W estimated here is based on Eq. 14which contains two constants a and b . Since a is a sub-leading coefficient in the associated ansatz, it is sensitiveto the range of β used for the fits as well. Careful strategyand resampling procedure are conducted in this study inorder to calculate W accurately. If the correlations among the data of various values of β are ignored, then the same resampling steps as well asthe criterion of χ / DOF < χ is the conven-tional uncorrelated χ , not the correlated χ described inthe main text) introduced in previous sections will leadto figs. 9, 10, and 11. Remarkably, while the outcomes p W L = 256 L = 128 FIG. 11: W as functions of p . The results are obtained fromthe analysis using the conventional uncorrelated χ . The solidsquares and solid circles are for L = 256 and L = 128, respec-tively. For some values of p , the L = 128 results containsthose corresponding to g c , the lower and upper bounds of g c .The solid and dashed lines in the figure represent the meanand standard deviation for all the W such that their corre-sponding p satisfy p < . L = 128and L = 256). of z shown in fig. 10 are slightly different from those infig. 7, the E and W presented in figs. 9 and 11 agreevery well with the ones demonstrated in figs. 6 and 8.In particular, the trend claimed from the analysis associ-ated with the correlated χ regarding the p dependenceof z and W , namely being complementary to each other,is valid as well for the outcomes obtained using the con-ventional uncorrelated χ (i.e. figs. 10 and 11). Thisobservation seems to reconfirm the conclusions resultingfrom investigating some lattice quantum chromodynam-ics data outlined in Refs. [38, 39].To summarize, the outcomes resulting from the investi-gations carried out here, specially the obtained numericalresults of E , z , and W , are not only important accom-plishments, but also can be considered as benchmarks forfuture related studies.This study is partially supported by MOST of Taiwan. [1] N. D. Mermin and H. Wagner, Phys. Rev. Lett. , 1133(1966). [2] P. C. Hohenberg, Phys. Rev. , 383 (1967).[3] Sidney Coleman, Communications in Mathematical Physics volume 31, pages 259264 (1973).[4] Axel Gelfert and Wolfgang Nolting, J. Phys.: Condens.Matter 13 R505 (2001).[5] J. Cardy,
Scaling and Renormalization in StatisticalPhysics , Cambridge University Press, Cambridge, UK,1996.[6] S. Sachdev,
Quantum Phase Transitions (CambridgeUniversity Press, Cambridge, 2nd edition, 2011).[7] S. Chakravarty, B. I. Halperin, and D. R. Nelson, Phys.Rev. B , 2344 (1989).[8] A. V. Chubukov and S. Sachdev, Phys. Rev. Lett. ,169 (1993).[9] A. V. Chubukov and S. Sachdev, Phys. Rev. Lett. ,2680 (1993).[10] A. V. Chubukov, S. Sachdev, and J. Ye, Phys. Rev. B , 11919 (1994).[11] A. W. Sandvik, A. V. Chubukov, and S. Sachdev, Phys.Rev. B , 16483 (1995)[12] M. Troyer, H. Kantani, and K. Ueda, Phys. Rev. Lett. , 3822 (1996).[13] Matthias Troyer, Masatoshi Imada, and Kazuo Ueda, J.Phys. Soc. Jpn. 66, 2957 (1997).[14] Jae-Kwon Kim and Matthias Troyer, Phys. Rev. Lett. , 2705 (1998).[15] Y. J. Kim, R. J. Birgeneau, M. A. Kastner, Y. S. Lee,Y. Endoh, G. Shirane, and K. Yamada, Phys. Rev. B 60,3294 (1999).[16] Y. J. Kim and R. J. Birgeneau, Phys. Rev. B , 6378(2000).[17] A. W. Sandvik, V. N. Kotov, and O. P. Sushkov, Phys.Rev. Lett. , 207203 (2011).[18] A. Sen, H. Suwa, and A. W. Sandvik, Phys. Rev. B ,195145 (2015).[19] D.-R. Tan and F.-J. Jiang, Phys. Rev. B , 245111(2018). [20] F.-J. Jiang, Phys. Rev. B , 024419 (2011).[21] F.-J. Jiang and U.-J. Wiese, Phys. Rev. B , 155120(2011).[22] Jhao-Hong Peng, L.-W. Huang, D.-R. Tan, and F.-J.Jiang, Phys. Rev. B , 174404 (2020).[23] Nvsen Ma, Anders W. Sandvik, and Dao-Xin Yao, Phys.Rev. B , 104425 (2014).[24] A. W. Sandvik, Phys. Rev. B , R14157 (1999).[25] A. W. Sandvik, AIP Conf. Proc. 1297, 135 (AIP, NewYork, 2010).[26] A. B. Harris, J. Phys. C , 1671 (1974).[27] J. T. Chayes, L. Chayes, D. S. Fisher, and T. Spencer,Phys. Rev. Lett. , 2999 (1986).[28] O. Motrunich, S.C. Mau, D.A. Huse, and D.S. Fisher,Phys. Rev. B , 1160 (2000).[29] A. W. Sandvik, Phys. Rev. Lett. , 177201 (2002).[30] O. P. Vajk and M. Greven, Phys. Rev. Lett. , 177202(2002).[31] R. Sknepnek, T. Vojta, and M. Vojta, Phys. Rev. Lett. , 097201 (2004).[32] Rong Yu, Tommaso Roscilde, and Stephan Haas, Phys.Rev. Lett. , 197204 (2005)[33] A. W. Sandvik, Phys. Rev. Lett. , 207201 (2006).[34] T. Vojta, J. Low Temp. Phys. 161, 299 (2010).[35] Dao-Xin Yao, Jonas Gustafsson, E. W. Carlson, and An-ders W. Sandvik, Physical Review B, , 172409 (2010).[36] Thomas Vojta, AIP Conference Proceedings 1550, 188(2013).[37] Thomas Vojta and Jos´e A. Hoyos, Phys. Rev. Lett. ,075702 (2014).[38] C. Michael, Phys. Rev. D , 2616 (1994).[39] C. Michael and A. McKerrell, Phys. Rev. D51