title: Ground state energy density, susceptibility, and Wilson ratio of a two-dimensional disordered quantum spin system
authors: J. -H. Peng D. -R. Tan F. -J. Jiang
text:  
  
  
  
  
  
  Department of Physics, National Taiwan Normal University,
  88, Sec.4, Ting-Chou Rd., Taipei 116, Taiwan
  
  Department of Physics, National Taiwan Normal University,
  88, Sec.4, Ting-Chou Rd., Taipei 116, Taiwan
  
  []fjjiang@ntnu.edu.tw
  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 of quenched disorder is investigated, using the first principles
    nonperturbative quantum Monte Carlo calculations (QMC). The employed
    disorder distribution has a tunable parameter p which can be considered
    as a measure of the corresponding randomness. In particular, when p=0,
    the disordered system becomes the clean one. Through a large scale QMC,
    the dynamic critical exponents z, the ground state energy densities E_0,
    as well as the Wilson ratios W of various p are determined with high
    precision. Interestingly, we find that the p dependence of z and W
    are likely to be complementary to each other. For instance, while the z
    of 0.4 ≤ p ≤ 0.9 match well among themselves and are statistically
    different from z=1 which corresponds to the clean system, the W for
    p < 0.7 are in reasonable good agreement with that of p=0.
    The technical subtlety of calculating these physical quantities for a
    disordered system is demonstrated as well. The results
    presented here are not only interesting from a theoretical perspective,
    but also can serve as benchmarks for future related studies.
    
           
  
  
  -0.25cm
  
  
  Ground state energy density, susceptibility, and Wilson ratio of a two-dimensional disordered quantum 
  spin system
      F.-J. Jiang
      December 30, 2023
  ==================================================================================================================
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  Spatial dimension two is extraordinary from a theoretical point of view.
  This is because
  according to the famous Mermin-Wagner theorem, for finite systems with
  short-range interactions, continuous symmetries cannot be
  broken spontaneously at any temperature T>0 <cit.>. As a result, for
  two-dimensional (2D) quantum spin antiferromagnets (AF), the associated
  studies have been focusing on certain exotic finite temperature
  properties of the systems.
  Particularly, several universal
  quantities are predicted and verified numerically. Such a temperature
  region where these unusual universal features exist
  is called the quantum critical regime (QCR) in the literature and
  has been explored in detail during the last few decades.
  <cit.>.
  
  For 2D quantum spin AF system, whenever QCR is mentioned, it typically
  refers to a finite temperature region. However, such an exotic regime
  extends to zero
  temperature at a quantum critical point (QCP). In addition, the (finite T)
  region above a QCP is where these profound
  characteristics can be uncovered the most clearly <cit.>.
  
  A physical observable, namely the spinwave velocity c plays an
  important role in those mentioned universal quantities of QCR for
  the 2D spin-1/2 antiferromagnets. For instance, the value of c
  without doubt has great impact on the determinations of two universal
  quantities of QCR, namely χ_u c^2/T ∼ 0.27185 and c/(Tξ) ∼ 1.04.
  Here χ_u and ξ are the uniform susceptibility and the correlation
  length, respectively <cit.>.
  In the phase with long-range antiferromagnetic order, c can be calculated
  efficiently using the spatial and the temporal winding numbers
  squared <cit.>.
  
  Considering a clean 2D spin-1/2 AF which comes with a given spatial
  arrangement of two types of antiferromagnetic couplings J' and J
  (J' > J),
  by tuning the ratio J'/J (i.e. the system is dimerized) a QCP may appear when J'/J exceeds a certain value
  (J'/J)_c. The dynamic critical exponent z associated with such a kind of QCP takes the value of 1.
  For a QCP g_c which is obtained by varying the associated parameter g,
  the physical quantity c
  scales as c ∝(g-g_c)^ν(z-1) close to g_c <cit.>, where ν is the correlation length exponent. As a
  result, c is a constant when the related z of a QCP is 1.
  For 2D quantum AF systems, the QCPs induced by dimerization introduce above
  belong exactly to this case. 
  When disorder is present, z > 1 and c is zero at g_c. Consequently, certain universal quantities of QCR cannot be
  calculated in a direct manner for disordered systems.
  
  While for a 2D disordered quantum spin antiferromagnet,
  certain quantities of QCR such as χ_u c^2/T cannot be directly accessed,
  yet some observables do not encounter the difficulty that c cannot be
  calculated with ease. One of them, namely
  the Wilson ratio W <cit.>, which will be defined later,
  is one of the main topics of
  our study presented here. In particular, we investigate the behavior of W
  with respect to the strength of randomness, which is controlled by a parameter p ≥ 0,
  of the employed disorder distribution. Here p=0 corresponds to the
  clean case. Apart from W, the dynamic critical
  exponents z as well as the ground state energy densities E_0 of several
  values of p considered in this study are determined as well.
  
  To carry out the proposed investigation, we have performed a large
  scale quantum Monte Carlo calculation (QMC). In addition, several
  p are considered and the simulations 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_0 grows
  monotonically with g_c(p) (hence p as well since g_c(p) ∝ p as shown in <cit.>), similar to that of the
  correlation length exponent
  ν <cit.>. Interestingly, the p dependence of z and W
  are likely to be complementary to each other. For instance, while the z
  of 0.4 ≤ p ≤ 0.9 match well among themselves and are statistically
  different from z=1 which corresponds to the clean system, the W for
  p < 0.7 are in reasonable good agreement with that of p=0 (W ∼ 0.1243).
  The subtlety of calculating these physical quantities
  for a disordered system is demonstrated here as well.
  Our investigation is important and interesting in itself
  from a theoretical perspective. In particular, the obtained outcomes can be
  used as benchmarks for future related studies.
  
  The rest of this paper is organized as follows. After the introduction,
  the studied model, the employed disorder distribution as well as the relevant
  observables are described. Following
  that we present our results. In particular, the numerical evidence for
  the mentioned complementary relation for z and W are
  demonstrated. Finally, a section concludes our study.
  
  
  
  
    
  
  
  § MICROSCOPIC MODELS AND OBSERVABLES
  
  
  
  The model investigated in our study has been described in detail in Refs. <cit.>.
  Here we briefly summarize certain technical perspectives of the considered
  system. The Hamiltonian of the investigated 2D
  disordered spin-1/2 herringbone Heisenberg model (on the square lattice)
  is given by
  
      H    =   ∑_⟨ ij ⟩J S⃗_i ·S⃗_j 
      + ∑_⟨ i'j' ⟩J' S⃗_i'·S⃗_j',
  
  where in Eq. (1) J (which are set to 1 here) and J' are the
  antiferromagnetic couplings (bonds) connecting nearest neighboring spins
  ⟨  ij ⟩ and ⟨  i'j' ⟩, respectively,
  and S⃗_i is the spin-1/2 operator at site i. In this study
  we use the convention J' > J.
  Fig. 1 is a cartoon representation of the considered model.
  The quenched disorder
  introduced into the system is based on the one employed in Refs. <cit.>.
  Specifically, for every bold bond in fig. <ref>, its
  antiferromagnetic strength J' takes the value of 1+(g-1)(1+p) 
  or 1+(g-1)(1-p) with equal probability. Here g>1 and 0 ≤ p ≤ 1.
  With the used conventions, the average and difference
  for these two types of bold bonds J' are given by g and 2p(g-1),
  respectively. In addition, p can be thought of as a measure for the disorder of the
  studied model as well.
  
  
  
  To perform the proposed calculations of determining the
  ground state energy density E_0, the dynamic critical exponent z, and the Wilson ratio W
  for the considered disordered system (with various p),
  the uniform susceptibility χ_u, the internal energy density E, 
  and the specific heat C_V (as functions of the temperature T or the inverse temperature β)
  are measured in our simulations.
  
  The uniform susceptibility χ_u is defined by
  
      χ_u = β/L^2⟨( ∑_iS_i^z)^2 ⟩,
   
  where β and L are the inverse temperature and linear box size
  used in the simulations, respectively.
  Furthermore, the internal energy density E and the specific heat C_V are given as 
  
      E = 1/L^2⟨ H ⟩, 
       
      C_V = ∂ E /∂ T.
  
  
  Using these observables, E_0, z, and W for various p of the studied
  disordered model can be determined with high precision.
  
  
  
  
  § 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 desired physical quantities, 
  we have carried out large-scale QMC
  using the stochastic series expansion (SSE) algorithm with very efficient 
  operator-loop update <cit.>. The simulations
  are done at the corresponding critical points g_c(p) for these chosen p.
  In addition, for every
  p, several hundred to
  few thousand randomness configurations with L=128 and (or) L=256 are generated.
  Each configuration is produced with its own random seed and then is used
  for all the calculations of the considered values of β.
  
  
  
   §.§ The strategy of calculating the Wilson ratio W
  
  
  In the framework of SSE,
  the quantities of internal energy density E and and specific heat C_V can
  be obtained by
  
      E    =    -1/L^2(⟨ n ⟩/β - 1/4∑_bJ_b ), 
       
        C_V    =   1/L^2(⟨ n^2⟩ - ⟨ n ⟩^2 - ⟨ n ⟩),
  
  respectively, where the summation is over all the bonds b and n is the number of nonidentity operators in the
  SSE operators sequence (operators string).
  
  Based on the large-N expansion of the relevant effective field theory,
  at the associated critical point it is predicted that for clean systems the (leading) low-T behavior of
  χ_u, E, and C_V are given by
  
      χ_u ∼1.0760/π c^2T,
            E ∼ E_0 + 2.8849/π c^2T^3, 
            C_V ∼8.6548/π c^2T^2,
    
  respectively, where the c appearing above is the spinwave velocity. With these
  leading T-dependence of χ_u, E and C_V, the Wilson ratio W can be
  expressed as 
  
      W = χ_uT/C_V∼ 0.1243.
  
  
  While C_V can be calculated directly from its definition 
  C_V = ∂ E/∂ T (or C_V = 1/L^2(⟨ n^2⟩ - ⟨ n ⟩^2 - ⟨ n ⟩)),
  as being shown in the literature, such a approach will lead to very noisy results at the
  region of low temperature <cit.>. 
  In addition, the fact that c is zero at the QCP of a disordered system
  prevents one from determining c directly.
  Motivated by the method outlined in Ref. <cit.>, 
  here we calculate W through the following procedures.
  
  Firstly, from the β-dependence of the internal energy density E, namely
  
  
      E(β) = E_0 + aβ^-1-2/z
  
  (here z is
  the dynamic critical exponent),
  one obtains a and z. Then the specific heat C_V, as a function of β,
  can be written as
  
  
      C_V = a(1+2/z)β^-2/z.
  
  
  Secondly, the β-dependence of χ_u is fitted
  to the expression 
  
      χ_u(β) = bβ^1-2/z.
   
  
  Finally, using Eqs. <ref> and <ref>,
  one arrives at the following formula for W
  
      W = b/a(1+2/z).
  
  
  In other words, instead of using C_V directly, here W is calculated through the
  coefficients z, a and b obtained
  from fitting the data of χ_u and E to their expected T-dependence ansatzes Eqs. <ref> and <ref>.
  
  
  
   §.§ 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. <ref>, <ref>, <ref>
  and <ref>.
  Both data of L=128 and L=256
  are put in these figures in order to demonstrate that the outcomes of L=256
  are (most likely) sufficient
  for size convergence.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ The analysis results associated with the clean model
  
  
  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. As a result,
  the following ansatz 
  
      χ_u = a_0 + aT + a_1T^2
  
  is considered to fit the χ_u data of p=0. Apart from that, the formula
  used to fit the data of E for p=0 is Eq. <ref> with z = 1 as well.
  
  By investigating the relevant data for L=128 and L=256, the finite size effect begins to
  appear when β > 16.0. Therefore the data of L=256 with β≤ 24.0 are used for
  the fits. We have additionally carried out fits using the L=256 data of β≤ 20.0
  and have found that these new results lead to a value of W which agrees quantitatively
  with that obtained using the L=256 data of β≤ 24.0.
  For each of data with a fixed range of β, in additional to taking care of finite
  size effect, the following procedures are adopted to calculate the corresponding W.
  
  Firstly, a bootstrap resampling (with respect to β) is conduct simultaneously
  for both χ_u and E. Secondly, fits
  for these obtained resampled data are performed. Finally W is determined using the
  outcomes of these fits.
  In these mentioned fits, Gaussian noises are considered as well. The above described steps are carried out for twenty thousand
  times and only those outcomes with both χ^2/DOF (DOF stands for degrees of freedom)
  of the two fits (for χ_u and E) being smaller than 3.0 are included as the candidate
  results of W.
  The resulting W, as well as its associated uncertainty, quoted for this set of data with that given fixed
  range of β are the mean and the standard deviation of these candidate results.
  
  We have conducted several calculations using various range of β, and each of these
  calculations comes with is own results
  (mean and uncertainty) for W. Moreover, to estimate the means and errors of the desired
  quantities appropriately,
  the weighted bootstrap resampling 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 resulting mean is given by
  
  
      ∑_i 1/σ_W_i^2W_i/∑_i1/σ_W_i^2.
  
  
  The reason for the use of above equation (called weighted mean in this study)
  is as follows. Notice that data with large standard deviations are less
  accurately determined than those with small standard deviations.
  As a result, those large standard deviation data should contribute
  less weight to the determination of the associated mean.
  
  After carrying out twenty thousand weighted bootstrap resampling steps, the resulting W
  is estimated to be W = 0.1238(3). The obtained W=0.1238(3) agrees very well with the
  theoretical prediction W=0.1243.
  This confirms the validity of the procedures introduced above for the calculation of W.
  
  The ground state energy density for the clean model is calculated by the same procedure
  and is given by E_0 = -1.022523(1).
  
  
  
   §.§ The results of the disordered model with various randomness strength p
  
  
  Since each generated configuration is used for all the simulations of the
  considered β, for a given set of p > 0 and L, the data themselves are correlated. Hence,
  to accurately estimate the associated errors for the coefficients in the fitting ansatzes,
  one should employ the correlated least χ^2 method for the analysis.
  However, the stability of the correlated least χ^2 method
  varies and depends on the quality of the data used for the fits. Moreover, biased outcomes may be reached
  if the associated covariance matrix for a given data set contains eigenvalues which have
  very small magnitude. Using the rule of thumb that the ansatz considered to fit correlated
  data should contain as few (to be determined) parameters as possible, we adopt
  the following approach to calculate E_0, z, and W for p > 0.
  
  First of all, for each β the bootstrap resampling method is performed for the raw data resulting from the generated disordered
  configurations.
  Second of all, these resampled data are used to calculate the disordered average of χ_u and E which are then
  considered for the relevant fits. Here the data employed for the fits of χ_u and E have different range of β.
  Indeed, as can be seen from figs. <ref>, <ref>, <ref>, for all the considered p > 0
  the associated E reaches its ground state value E_0 quickly, while this is not the case for χ_u. As a result, it is
  more appropriate to 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 fits of χ_u and
  E are done, the resulting results are then
  put back to calculate the associated correlated χ^2/DOF. Here a cut-off for the
  eigenvalues of the associated covariance matrix is imposed in order to avoid biased results.
  These introduced steps are
  performed for many times, and
  only those results which have correlated χ^2/DOF smaller than 3.0
  for both the fits of χ_u and E are considered for later calculations.
  Finally for each of the considered p, the above procedures
  have been applied to many sets consisting of various range of β. Each of these
  sets (The whole sets is denoted by S) has its own results (mean and standard deviation) of E_0, z,
  and W as well as the number of successful calculations. 
  
  For a considered p, the final quoted results of E_0, z, and W in this study are
  estimated by a bootstrap resampling procedure using the
  following formula to calculate the mean of every resampled data from S.
  
  
      ∑_i N_i O_i /σ_O_i^2/∑_i N_i/σ_O_i^2,
  
  where {O_i}, {σ_O_i} and {N_i} stand for the randomly picked outcomes in S, the associated
  standard deviations, as well as the related numbers of
  (successful calculated) results of these picked outcomes, respectively.
  Finally, such a resampling step is conducted for several thousand times, and the numerical values
  presented here for these considered
  physical quantities are the resulting means and standard deviations (estimated conservatively)
  of this procedure. The uncertainties of E_0 calculated by the described
  steps are much smaller than those of the original E_0 contained in S.
  Hence, for the data in S we have also calculated their associated weighted errors. 
  The dominant one of these two estimations, namely the standard deviations and the weighted errors, 
  are the final values quoted here.
  
  The E_0, z and W as functions of p calculated by the procedures introduced above
  are shown in figs. <ref>, <ref> ,<ref>.
  The related results for the clean model are shown in these figures as well for comparison.
  
  The E_0 as a function of p shows a monotonic behavior in magnitude, which is similar to
  that of the correlation length exponents ν obtained in Ref. <cit.>.
  
  Regarding the z presented in the figure, one observes that the magnitude
  of z increases with p until p reaches a specific p_c < 0.4.
  For p ≥ 0.4, all
  the calculated values of z lie between (around) 1.3 and (around) 1.4.
  If one takes into
  account the systematic errors due to the uncertainties of g_c(p), then the z for p ≥ 0.4
  are fairly close to each other. The solid and dashed lines in fig. <ref> represent the
  mean and standard deviation for all the values of z associated with p ≥ 0.4 (including both those of L=128 and L=256). These guided lines justify the claim
  made above.
  This phenomenon is consistent with the one shown in Ref. <cit.>,
  where the calculated z corresponding to various parameters take a universal value.
  We would like to point out that when conducting the determination of z from χ_u,
  the obtained results are somehow a little bit sensitive to the considered fitting range of χ_u.
  This motives the use of the resampling procedures described above. In conclusion, our analysis indicates that it
  is subtle to calculate the quantity z and a careful strategy is needed.
  
  Finally, in fig. <ref> we demonstrate the results of W as functions of p obtained from
  the analysis outlined previously. Intriguingly, similar to the scenario of z, for those W corresponding to p < 0.7,
  their values are more or less close to each other. The solid and dashed lines in the figure
  again stand for the mean and standard deviation for all the W with their associated p satisfy p < 0.7 (including both those of L=128 and L=256).
  Considering the impact resulting from the errors of g_c(p), the scenario that W take the same
  value (or at least values close to each other) for all the p such that
  p < 0.7 is probable. 
  Interestingly, the correlation length exponent ν is beginning to fulfill the Harris criterion
  when p > 0.8 and this is where the magnitude of W increases sharply. This observation
  implies that there may exist a relation between W and fulfillment of Harris criterion.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § DISCUSSIONS AND CONCLUSIONS
  
  
  
  
  
  
  
  
  
  
  
  In this study, we calculate the Wilson ratio W of a 2D spin-1/2 antiferromagnetic Heisenberg model
  with a specific quenched disorder, using the first principle nonperturbative quantum Monte Carlo
  simulations. The employed disorder distribution has a tunable parameter p which can be considered as a measure of
  randomness. The W of the clean case as well as that of p=0.2,0.4,0.5,0.6,0.7,0.8,0.9 are determined with
  high precision. The critical dynamic exponents z and the ground state energy densities E_0 are obtained as well.
  
  Remarkably, for the considered system with the employed 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 ≥ 0.4.
  This agrees with the outcomes determined in Ref. <cit.>. In addition, the calculated W for 0 < p < 0.7
  also have a trend of stay close to the result W ∼ 0.1243 of the clean model (p=0). In particular,
  the value of W begins to increase sharply when p is approaching p=0.9 where the Harris criterion is fulfilled.
  Considering the fact that with what conditions the Harris criterion is valid is still not known <cit.>,
  the results presented here may shed some light on setting up some
  useful guidelines to decide whether the Harris criterion is valid for a given disorder distribution.
  
  
  Apart from the subtlety of calculating z described previously, the determination of W is
  extremely non-trivial as well. Indeed, the W estimated here is based on Eq. <ref> which contains
  two constants a and b. Since a is a sub-leading coefficient in the
  associated ansatz, it is sensitive to the range of β used for the fits as well.
  Careful strategy and resampling procedure are conducted in this study in order to calculate
  W accurately.
  
  If the correlations among the data of various values of β are ignored, then
  the same resampling steps as well as the criterion of χ^2/DOF <3 (here the χ^2
  is the conventional uncorrelated χ^2, not the correlated χ^2 described in the main text) introduced in previous sections will lead to
  figs. <ref>, <ref>, and <ref>. Remarkably, while the outcomes of z shown in
  fig. <ref> are slightly different from those in fig. <ref>, the E_0 and W
  presented in figs. <ref> and <ref> agree very well with the ones demonstrated
  in figs. <ref> and <ref>. In particular, the trend claimed from the analysis associated with the correlated χ^2 regarding the p dependence
  of z and W, namely being complementary to each other, 
  is valid as well for the outcomes obtained using the conventional uncorrelated
  χ^2 (i.e. figs. <ref> and <ref>). This observation seems to reconfirm the conclusions resulting from investigating
  some lattice quantum chromodynamics data outlined in Refs. <cit.>.
  
  
  To summarize, the outcomes resulting from the investigations carried out here,
  specially the obtained numerical results of E_0, z, and W, are not only
  important accomplishments, but also can be considered as benchmarks for future
  related studies.
  
  
  
  0.5cm
  
  
  This study is partially supported by MOST of Taiwan.
  
   
   
   
  
   
   
   
   
   
   
   
   
   
   
   
   
   
  
  
  
  99
  
  Mer66
  N. D. Mermin and H. Wagner, Phys. Rev. Lett. 17, 1133 (1966).
    
  Hoh67
  P. C. Hohenberg, Phys. Rev. 158, 383 (1967).
  
  Col73
  Sidney Coleman, Communications in Mathematical Physics volume 31, pages 259–264 (1973).
  
  Gel01  
    Axel Gelfert and Wolfgang Nolting, J. Phys.: Condens. Matter 13 R505 (2001).
  
  Car96
    J. Cardy, Scaling and Renormalization in Statistical Physics, Cambridge
    University Press, Cambridge, UK, 1996.
    
  
  
  Sac11
  S. Sachdev, Quantum Phase Transitions  (Cambridge University Press, Cambridge, 2nd edition, 2011).
    
  Cha89
    S. Chakravarty, B. I. Halperin, and D. R. Nelson,
    Phys. Rev. B 39, 2344 (1989).
    
  
  Chu93
  A. V. Chubukov and S. Sachdev, Phys. Rev. Lett. 71, 169 (1993).
  
  Chu931
  A. V. Chubukov and S. Sachdev, Phys. Rev. Lett. 71, 2680 (1993).
  
  Chu94
  A. V. Chubukov, S. Sachdev, and J. Ye, Phys. Rev. B 49, 11919 (1994).
  
  San95
  A. W. Sandvik, A. V. Chubukov, and S. Sachdev, Phys. Rev. B 51, 16483 (1995)
  
  Tro96
  M. Troyer, H. Kantani, and K. Ueda, Phys. Rev. Lett. 76, 3822 (1996).
  
  Tro97 
  Matthias Troyer, Masatoshi Imada, and Kazuo Ueda, J. Phys. Soc. Jpn. 66, 2957 (1997).
  
  Tro98
  Jae-Kwon Kim and Matthias Troyer, Phys. Rev. Lett. 80, 2705 (1998).
  
  Kim99
  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).
  
  Kim00
  Y. J. Kim and R. J. Birgeneau, Phys. Rev. B 62, 6378 (2000).
  
  San11
  A. W. Sandvik, V. N. Kotov, and O. P. Sushkov, Phys.
  Rev. Lett. 106, 207203 (2011).
  
  Sen15
  A. Sen, H. Suwa, and A. W. Sandvik, Phys. Rev. B 92, 195145 (2015).
  
  
  Tan181
  	D.-R. Tan and F.-J. Jiang, Phys. Rev. B 98, 245111 (2018).
  
   Jia111
  F.-J. Jiang, Phys. Rev. B 83, 024419 (2011).        
  
  Jia112
  F.-J. Jiang and U.-J. Wiese,
  Phys. Rev. B 83, 155120 (2011).
  
  Pen20
  Jhao-Hong Peng, L.-W. Huang, D.-R. Tan, and F.-J. Jiang,
  Phys. Rev. B 101, 174404 (2020).
  
  Nvs14
  Nvsen Ma, Anders W. Sandvik, and Dao-Xin Yao, Phys. Rev. B 90, 104425 (2014).
  
  
  San99
  A. W. Sandvik, Phys. Rev. B 66, R14157 (1999).
  
  San10
  A. W. Sandvik, AIP Conf. Proc. 1297, 135 (AIP, New York, 2010).
  
  
  Har74
  A. B. Harris, J. Phys. C 7, 1671 (1974).
  
  Cha86
  J. T. Chayes, L. Chayes, D. S. Fisher, and T. Spencer, Phys. Rev.
  Lett. 57, 2999 (1986).
  
  Mot00
  O. Motrunich, S.C. Mau, D.A. Huse, and D.S. Fisher,
  Phys. Rev. B 61, 1160 (2000).
  
  San02d
  A. W. Sandvik, Phys. Rev. Lett. 89, 177201 (2002).
  
  Vaj02
  O. P. Vajk and M. Greven, Phys. Rev. Lett. 89, 177202 (2002).
  
  Skn04
  R. Sknepnek, T. Vojta, and M. Vojta, Phys. Rev. Lett. 93,
  097201 (2004).
  
  Yu05
  Rong Yu, Tommaso Roscilde, and Stephan Haas, Phys. Rev. Lett. 94,
  197204 (2005)
  
  San06d
  A. W. Sandvik, Phys. Rev. Lett. 96, 207201 (2006).
  
  Voj10
  T. Vojta, J. Low Temp. Phys. 161, 299 (2010).
  
  Yao10
  Dao-Xin Yao, Jonas Gustafsson, E. W. Carlson, and Anders W. Sandvik,
  Physical Review B, 82, 172409 (2010).
  
  
  Voj13
  Thomas Vojta, AIP Conference Proceedings 1550, 188 (2013).
  
  Voj14
    Thomas Vojta and José A. Hoyos, Phys. Rev. Lett. 112, 075702 (2014).
  
  Mic94
    C. Michael,
    Phys. Rev. D 49, 2616 (1994).
  
  Mic95
  C. Michael and A. McKerrell,
  Phys. Rev. D 51, 3745 (1995).
  
  
  
  
  
  
  
  
  
  
  