Stochastic approximation of dynamical exponent at quantum critical point

# Stochastic approximation of dynamical exponent at quantum critical point

Shinya Yasuda Department of Applied Physics, University of Tokyo, Tokyo 113-8656, Japan    Hidemaro Suwa Department of Physics, University of Tokyo, Tokyo 113-0033, Japan    Synge Todo Department of Physics, University of Tokyo, Tokyo 113-0033, Japan Institute for Solid State Physics, University of Tokyo, Kashiwa, 277-8581, Japan
July 4, 2019
###### Abstract

We have developed a unified finite-size scaling method for quantum phase transitions that requires no prior knowledge of the dynamical exponent . During a quantum Monte Carlo simulation, the temperature is automatically tuned by the Robbins-Monro stochastic approximation method, being proportional to the lowest gap of the finite-size system. The dynamical exponent is estimated in a straightforward way from the system-size dependence of the temperature. As a demonstration of our novel method, the two-dimensional quantum model in uniform and staggered magnetic fields is investigated in the combination of the world-line quantum Monte Carlo worm algorithm. In the absence of the uniform magnetic field, we obtain the fully consistent result with the Lorentz invariance at the quantum critical point, , i.e., the three-dimensional classical universality class. Under a finite uniform magnetic field, on the other hand, the dynamical exponent becomes two, and the mean-field universality with effective dimension governs the quantum phase transition.

quantum spin system, quantum phase transition, quantum Monte Carlo, Robbins-Monro algorithm, anisotropy
###### pacs:
05.10.Ln, 05.30.Rt, 64.60.F-, 75.10.Jm

## I Introduction

Recent enhancement of the computational power has enabled us to simulate larger-scale systems with higher precision than ever before. In particular, with the help of the recent development of simulation algorithms for strongly correlated quantum systems, a number of simulations have been performed to elucidate the novel nature of quantum phase transitions, in which many-body physics plays an essential role Sachdev (1999); Avella and Mancini (2013). Quantum phase transitions occur at absolute zero temperature, triggered by quantum fluctuations. Through the quantum-classical mapping, a quantum phase transition in dimensions, if it is of second order, can be generally described by the critical theory as the temperature-driven phase transition in a -dimensional classical system with the same symmetries, where is the so-called dynamical exponent Suzuki (1976); Vojta (2003).

A world-line quantum Monte Carlo (WLQMC) method is one of the most powerful tools for investigating quantum critical phenomena without any bias or approximation Suzuki (1994); Landau and Binder (2005). A quantum system in dimensions is mapped to a classical system in dimensions in the WLQMC method. The system length along the additional direction, the imaginary-time direction, is given by the inverse temperature, .

When one performs a WLQMC simulation to investigate quantum criticality, the choice of for each system size is essential. The reason is that the quantum critical system can be extremely anisotropic even if the interactions are isotropic in real space. While the correlation length in the real-space direction diverges as , that in the imaginary-time direction does as , where is the coupling constant that controls quantum fluctuations, the quantum critical point, and the critical exponent. If the dynamical exponent is one, the space-time isotropy is kept aside from a scale factor, or the velocity of low-energy excitation. In the meanwhile, there are phase transitions with a dynamical exponent larger than one. The Bose-Hubbard model with randomness that exhibits quantum criticality with has been extensively investigated analytically Fisher et al. (1989); Singh and Rokhsar (1992); Herbut (1997); Weichman and Mukhopadhyay (2007), numerically Sørensen et al. (1992); Zhang et al. (1995); Prokofâev and Svistunov (2004); Priyadarshee et al. (2006); Meier and Wallin (2012); Álvarez Zúñiga et al. (2015) as well as experimentally Hüvonen et al. (2012); Yu et al. (2012).

Let us review the renormalization group and the scaling theory near a quantum critical point. Consider the scale transformation with a certain length scale, . A physical quantity, denoted as , is generally transformed as

 F(g−gc,L−1,β−1)=byFF(b1/ν(g−gc),bL−1,bzβ−1)=LyF~F(L1/ν(g−gc),Lzβ−1), (1)

where is the scaling dimension of the quantity under consideration. In the second line of Eq. (1), we chose and introduced . This equation has several unknown constants, , , , , and the scaling function itself.

In order to determine the constants in the finite-size scaling ansatz (1), one had to repeat simulations densely in the three-dimensional parameter space , and perform a multi-parameter finite-size scaling analysis as in Refs. 19 and 20. It typically requires considerable computational resources to scan the multi-dimensional parameter space. Instead of the exhaustive scanning, simulations with some assumed were performed in most previous studies. The consistency was checked after the calculation as in Ref. 21. This approach, however, would be awfully inefficient in the case without knowledge of the value of in advance. Another approach for estimation was to focus on the temperature dependence of the correlation length, at and  Priyadarshee et al. (2006). After the correlation length in the thermodynamic limit was extrapolated at each temperature, the low-temperature asymptotic behavior was analyzed. This two-step procedure requires additional computational cost, and possibly introduces some uncontrollable systematic error from the extrapolations, even if the location of quantum critical point, , is known. In the meantime, the winding numbers of the world-lines in space and time directions were exploited in Refs. 22 and 23. A parameter, or , was interpolated so that the winding number squared in each direction averagely took the same value. However, such an interpolation is again non-trivial in a multi-parameter space and multi conditions.

One of the most effective strategies to overcome the difficulty of a multi-parameter scaling is to introduce an auto-tuning technique. A number of auto-tuning techniques have already been used in numerical simulations in the field of statistical physics. For example, the invaded cluster algorithm Machta et al. (1995) or the probability-changing cluster algorithm Tomita and Okabe (2001) can automatically locate the critical point. The Wang-Landau algorithm Wang and Landau (2001) enables us to directly estimate the density of states of a system.

In the present paper, we employ the stochastic approximation method. Recently in Ref. 27, a method to automatically optimize the aspect ratio of a quantum system was proposed for analyzing quantum criticality under strong spatial anisotropy. The relative correlation length, , where and are the correlation length and the system size in direction, respectively (, , or for two-dimensional systems and ), was adjusted for making the system virtually isotropic as . Figure 1 schematically illustrates virtually anisotropic and isotropic one-dimensional quantum systems. In Ref. 27, the stochastic approximation scheme was applied to the staggered dimer antiferromagnetic Heisenberg model, and the universality class of the quantum critical point was successfully identified in spite of the existence of strong finite-size corrections that easily lead a naive finite-size scaling analysis to an incorrect conclusion Wenzel and Janke (2009); Wenzel et al. (2008). Note that the tuning method using the correlation length Yasuda and Todo (2013) is applicable to general systems, while the method based on the winding number Jiang (2011, 2012) works only for systems with U(1) symmetry.

The aim of the present paper is to propose a unified finite-size scaling method based on the stochastic approximation technique for quantum criticality with general . The relevant critical exponents including the dynamical exponent will be obtained simultaneously without any prior knowledge or assumption of the values. We will demonstrate our approach for the two-dimensional quantum model in uniform and staggered magnetic fields along direction (in spin space). We will clarify that the dynamical exponent becomes two under a finite uniform magnetic field, while it does one in the absence of a uniform field.

This paper is organized as follows: Sec. II introduces the scaling ansatz of the correlation lengths for space-time anisotropic systems, the Robbins-Monro stochastic approximation method, and its convergence property. It is also discussed how the stochastic approximation method is applied to the present finite-size scaling analysis. In Sec. III, the model considered in the present paper and the WLQMC method are introduced. The numerical results are shown in Sec. IV. Finally our study is concluded in Sec. V. The technical details are reported in Appendices A and B.

## Ii Stochastic approximation method

### ii.1 Conditions for realizing space-time isotropy

As noted in Sec. I, for the system with , one should pay attention to the space-time aspect ratio when considering the finite-size scaling analysis. In this section, we explain conditions to realize a virtually isotropic system during a simulation. For simplicity, we assume that the model considered hereafter has no anisotropy in real space. Generalization to systems with spatial anisotropy is straightforward.

Let us start by choosing in Eq. (1). In this case, is one, i.e.,

 ξ(g−gc,L−1,β−1)=L~ξ((g−gc)L1/ν,Lz/β). (2)

Another choice, , yields

 ξτ(g−gc,L−1,β−1)=Lyτ~ξτ((g−gc)L1/ν,Lz/β), (3)

where we set . At the quantum critical point, , the correlation length in the imaginary-time direction exhibits the power law, , in the limit of . By the definition of the dynamical exponent, one finds . Then, dividing both sides of Eq. (3) by yields

 ξτ(g−gc,L−1,β−1)/β=~ξ′τ((g−gc)L1/ν,Lz/β), (4)

where .

Here, let us introduce two conditions,

 ξ/L=R, (5)

where is an arbitrarily chosen constant, and

 ξτ/β=Rτ, (6)

where is another constant. Assume that conditions (5) and (6) are both satisfied by tuning and in simulating systems with different system sizes. If this is the case, and are kept constant even though they are different functions. Meanwhile, the set of arguments of in Eq. (2) and that of in Eq. (4) are the same as each other. That is, the different functions and sharing the same arguments are kept constant at different system sizes. This means that each of the arguments should be constant if the functions have some reasonable monotonicity. The monotonicity of the scaling functions is expected to hold near a generic critical point and supported by our numerical calculation shown below. Then Eqs. (5) and (6) provide solutions, and , for each system size:

 gc(L)−gc∝L−1/ν (7)

and

 β(L)∝Lz. (8)

Thus, the coupling constant, , automatically converges to the critical point as increases. Moreover, the dynamical exponent can be simultaneously estimated from the asymptotic dependence of the inverse temperature, .

### ii.2 Robbins-Monro stochastic approximation method

In this section, we introduce an iteration procedure to fulfill the conditions proposed in the previous section. Our task is to solve the system of nonlinear equations, and , with respect to and for given , , and . The solution cannot be obtained by standard iterative methods for nonlinear equations, such as the Newton-Raphson method. It is because and have statistical errors coming from the Monte Carlo sampling that make the conventional methods unstable. We thus employ the stochastic approximation method explained below.

Let us see a concrete example of the stochastic approximation. For simplicity, assume that is already set to . We estimate the optimal that satisfies the relation . The solution of this equation is denoted as . First, one runs a short Monte Carlo simulation with a trial parameter and measures the correlation length, then calculates . Next, one updates the parameter, , by using the Robbins-Monro type update procedure Robbins and Monro (1951); Bishop (2006)

 β(n+1)=β(n)−pnA(β(n)) (9)

with and repeats the above until converges to a certain value with increasing . Here, is a (constant) parameter that determines the gain of the feedback. Regardless of the choice of the gain, it is proved that converges to in with probability one Robbins and Monro (1951); Albert and Gardner (1970).

As explained in Appendix A, the mean of the probability distribution of at the -th step (denoted as ) converges as , where is the derivative of at , and the sign of is chosen as the same with . For , the variance of at the -th step (denoted as ) is evaluated as . For , on the other hand, , where is the asymptotic variance and is the statistical error resulting from a Monte Carlo estimation of . Here we should set to minimize the variance (see the detailed discussion in Appendix A). By this choice of , it is also guaranteed that the systematic error of decreases faster than the statistical (standard) error. In actual simulations, one needs to perform some (10 at least) independent stochastic approximation processes to estimate error bars of and physical quantities. The number of steps of each approximation process has to be large enough ( typically) for the systematic error to become negligible in comparison to the statistical error.

The present stochastic approximation method can be extended to multi-dimensional problems in a straightforward way. Below, we will apply the method to the quantum phase transition of two-dimensional model in uniform and staggered magnetic fields in order to demonstrate the efficiency of the present approach and clarify the quantum phase transitions.

## Iii Model and quantum Monte Carlo method

### iii.1 S=1/2 quantum Xy model in uniform and staggered magnetic fields

The Hamiltonian of the two-dimensional quantum model in uniform and staggered magnetic fields is defined as follows:

 H=−12∑⟨j,k⟩(S+jS−k+S+kS−j)−∑j[hu+hs(−1)σ(j)]Szj, (10)

where () is the -component raising (lowering) operator at site , denotes a pair of nearest-neighboring spins, and () is the amplitude of the uniform (staggered) magnetic field. Here we consider the square lattice of linear extent with the periodic boundary conditions, and the lattice is bipartite with even . If site belongs to one of the sublattices, takes zero, otherwise .

This model can be mapped to the hard-core boson model with the uniform and the staggered chemical potentials Lieb et al. (1961). The phase diagram of the model consists of several phases Hen et al. (2010); Hen and Rigol (2009): (i) the disordered phase that corresponds to the insulating or pinning phase in the boson model, (ii) the -plane ferromagnetic phase with non-zero transverse magnetization, or the compressible superfluid phase, and (iii) the fully-polarized phase along , or the empty (fully-occupied) phase. We will fix to some value and change across the phase boundary. When is smaller than the saturation field, , a phase transition from the ferromagnetic phase to the disordered phase occurs as increases. If is larger than , an additional phase transition from the fully-polarized phase to the ferromagnetic phase occurs. When , the particle-hole symmetry holds and the phase transition is known to belong to the three-dimensional (3D-) universality class, i.e., . Phase transitions different from the 3D- universality with , on the other hand, are expected for  Hen et al. (2010); Fisher et al. (1989).

### iii.2 World-line quantum Monte Carlo worm algorithm

In order to simulate the system described by Hamiltonian (10), we used the worm (directed-loop) algorithm Prokof’ev et al. (1998); Syljuasen and Sandvik (2002); Kawashima and Harada (2004) with the continuous-time path-integral representation. In the continuous-time representation, we introduce imaginary time as

 Z=Tre−βH=Tr[exp(−∫β0dτH)], (11)

where is the partition function. The continuous-time formulation was adopted because of the convenience for calculating the Fourier component of the imaginary-time correlation function, which we will exploit to calculate . Expanding the exponential in the r.h.s. of Eq. (11), we insert the identity, , between the operators, where is a complete basis set of the Hilbert space. We then obtain

 Z= 1+∞∑n=1∑(ϕ1,…,ϕn)∫β0dτ1⋯∫βτn−1dτn ×n∏ℓ=1⟨ϕℓ|(−H)|ϕℓ+1⟩, (12)

where . In our WLQMC simulation, a state in the basis set is the direct product of the eigenstate of the local operator (up or down). The Hamiltonian (10) conserves total of the system and thus a space-time configuration forms continuous lines of up spins (or down spins), i.e., the world-lines.

One can consider the integrand in Eq. (12) as a weight (probability measure) of each world-line configuration. In order to make a simulation efficient, the second (site) term in Hamiltonian (10) is included in the bond term as

 14∑⟨j,k⟩[hu(Szj+Szk)+hs(−1)σ(j)(Szj−Szk)], (13)

where the factor comes from the coordination number of the square lattice. The matrix elements of the combined bond term are expressed as

 ⎡⎢ ⎢ ⎢ ⎢⎣hu/40000hs(−1)σ(j)/41/2001/2−hs(−1)σ(j)/40000−hu/4⎤⎥ ⎥ ⎥ ⎥⎦|↑↑⟩|↑↓⟩|↓↑⟩|↓↓⟩. (14)

A constant larger than or equal to needs to be added to the diagonal elements for ensuring the non-negativity of the weights.

In the worm algorithm, extended world-line configurations are introduced. The configurations to sample in the Monte Carlo method consist of the original world-lines and the world-lines with a pair of kinks, points of discontinuity. Such a pair is called a worm, and each of discontinuity head or tail. In the present spin model, the worm is represented by the pair of spin ladder operators, (, ) or (,), each of which is defined at a space-time point.

The whole update process of the world-line configuration consists of the diagonal update and the worm update Kawashima and Harada (2004). In the former, the diagonal bond operators in the Hamiltonian (10) are inserted into or removed from a world-line configuration according to the diagonal elements. In the latter, first a worm, i.e., a pair of the raising and lowering operators, is inserted at a randomly chosen space-time point, and either operator is chosen as the head. The order of the ladder operators is uniquely determined in the case with since the local degree of freedom is binary (up or down). The worm head then moves along the imaginary-time direction until it arrives at a bond operator. At the operator, the worm head is scattered and its moving direction and/or sitting site may be changed stochastically according to the matrix elements (14). In Fig. 2, an example of the worm-scattering process is illustrated. We choose transition probabilities so as to minimize the bounce probability (see Fig. 2), by breaking both the detailed balance of each worm-scattering process and even that of the whole Monte Carlo dynamics Suwa and Todo (2010). This scattering process is repeated until the worm head reaches back its own tail. Then the head and tail destroy each other. The worm is inserted at several times in each Monte Carlo step.

We will investigate the phase transition between the -plane ferromagnetic phase and the disordered phase. The order parameter is the transverse (off-diagonal) magnetization in or direction. Although it is non-trivial to measure off-diagonal correlation in a WLQMC simulation, one can efficiently calculate the structure factor

 S0=1L2⟨∑j,k(SxjSxk+SyjSyk)⟩, (15)

the transverse susceptibility

 χ=1L2β⟨∑j,k∫β0dτ1∫β0dτ2S+j(τ1)S−k(τ2)⟩, (16)

and the Fourier component of the (spatial and temporal) correlation functions, exploiting the virtue of the worm-update process Prokof’ev et al. (1998). Here the symbols are defined as follows: and . The correlation lengths in , , and directions are then calculated from the Fourier components by the second-moment method Cooper et al. (1982); Todo and Kato (2001). The detail of the measurement of these quantities is explained in Appendix B.

At the critical point, the structure factor and the susceptibility exhibit the following power-law behavior:

 S0(L)∝Lθ, (17)

and

 χ(L)∝Lγ/ν, (18)

respectively. Here we introduce , where is not the inverse temperature but the critical exponent of the order parameter. The exponent of the susceptibility is denoted by . Note that in Eqs. (8), (17), and (18), one can use the quantities evaluated at , the solution of Eqs. (5) and (6) for each system size , instead of the true critical point , as both give the same exponent.

## Iv Numerical results

### iv.1 For hu=0

First we discuss the case without the uniform magnetic field. We performed WLQMC simulations for system sizes , and obtained the optimal inverse temperature and staggered magnetic field for each by solving Eqs. (5) and (6) using the stochastic approximation. The optimal inverse temperature ensures that the system is virtually isotropic, and the optimal staggered magnetic field gives an estimate of the critical point. The target relative correlation lengths, and , were chosen as , , or 0.7, where is the approximate value of at the critical point of the three-dimensional classical model Campostrini et al. (2001). Note that a particular choice of and does not introduce any bias to the final estimates; it affects only the speed of convergence to the thermodynamic limit as seen below.

The critical strength of the staggered magnetic field, , can be estimated based on the asymptotic form (7), i.e.,

 hsc(L)=hsc+c(R)L−1/ν, (19)

where is a certain constant which depends only on (, here). The system-size dependence of is shown in Fig. 3. We extrapolated the critical point from the finite-size data by using the ansatz Eq. (19) and assuming the same for three different values of . Seven parameters, and (, ) for each , were determined by the least squares fitting. (Note that also is a fitting parameter in our analysis.) In order to estimate the statistical error, we performed the following bootstrap procedure. From several hundreds Robbins-Monro runs, the physical quantities and the parameters [ and ] were obtained with some statistical error bars. Then the data were resampled from the Gaussian distribution with the estimated mean and variance. The generated samples were fitted for large-enough system-size data that minimize , where the asymptotic form (19) would approximate the plots well (we used the data for and obtained ). This procedure was repeated 4000 times, which yielded as many fitted functions. By taking the average of the function values, we obtained a whole shape of that should be asymptotically accurate, which is shown in Fig. 3 as the solid line for each . Finally, the bootstrap estimation gave the critical point for , where the number in the parenthesis denotes the standard error of the estimation in the last digit(s). It is consistent with the previous report,  Priyadarshee et al. (2006), but more precise by an order of magnitude.

It should be noted that in Fig. 3 the fastest convergence is achieved by the choice of , i.e., the leading correction seems at . Meanwhile, the leading correction is likely in the order of at and , where  Campostrini et al. (2001) is the critical exponent of the correlation length for the 3D- universality class. These findings indicate that this phase transition belongs to the 3D- universality class.

In Fig. 4, we present the system-size dependence of the inverse temperature, the static structure factor, and the transverse susceptibility. Assuming the asymptotic forms [Eqs. (8), (17), and (18)], we conclude , , and . These results are fully consistent with the scenario of the 3D- universality class, and  Campostrini et al. (2001), with . The final estimates for , , and , which were evaluated from the asymptotic behavior of the different quantities, indeed satisfy the scaling relation

 γ/ν=θ+z (20)

within the error bar.

### iv.2 For hu=0.5

Next, we discuss the critical point and the critical exponents for the case with finite uniform magnetic field . In this case we used , , and for . Following the same procedure with the case for , we obtained , the quantum critical point. In the fitting procedure, the data with were used ().

The system-size dependence of the physical quantities is shown in Fig. 6. In comparison to the case with shown in Fig. 4, larger corrections to scaling are seen, especially for . To cope with the strong finite-size corrections, we took the following procedure: Assume we have data points at system sizes , , , , where . First, we construct triads consisting of the data with three consecutive system sizes as , , , and , defining as the average system size of each triad. Next, we fit each triad with a simple power function, i.e., with and the fitting parameters depending on . The error of is estimated by the bootstrap method as explained above. Then, we fit with a quadratic function of , and extrapolate the critical exponent in the thermodynamic limit. In our fitting procedure of the exponent, we assume that the fitting function should be monotonic.

The extrapolation results are shown in Fig. 7. For the dynamical exponent, we estimated ; the effective dimension of the imaginary-time axis changes from one to two by the introduction of uniform magnetic field . As for the other critical exponents, and were obtained. These values coincide with the mean-field exponents, i.e., and . This is consistent with the result for the dynamical exponent, , by which the effective dimension of the critical theory becomes four, the upper critical dimension. Thus we have demonstrated that our finite-size scaling method enables us to extract the dynamical exponent successfully without any prior knowledge of the value of . The universality of the quantum phase transition without particle-hole symmetry belongs to the mean-field universality class. This is consistent with the discussion on the Bose-Hubbard model in Ref. 7.

## V Summary and discussions

In the present paper, we have presented the unified finite-size scaling method that works well regardless of the value of the dynamical exponent. During the WLQMC simulation, the system size in the imaginary-time direction in the path-integral representation is adjusted automatically so as to satisfy the conditions, and , based on the Robbins-Monro stochastic approximation. This auto-tuning procedure guarantees that the coupling constant converges to the critical point and the inverse temperature is proportional to for large enough .

We then applied the method to the two-dimensional quantum model in uniform and staggered magnetic fields. The correlation lengths were measured by the worm algorithm based on the continuous imaginary-time representation. In the absence of the uniform magnetic field, , our numerical results are consistent with the 3D- universality class. This system can be mapped into the half-filled hard-core boson system. The Lorentz invariance, , reflects the particle-hole symmetry at half-filling. In the case with , we have concluded that the dynamical exponent changes to two and the other exponents take the mean-field values. This result of the mean-field universality is consistently explained by the conclusion that the dimension of the effective field theory is four; , the upper critical dimension. Our conclusion agrees on the discussion in Ref. 7, in which the authors claimed that the fourth-order term in the effective action becomes irrelevant.

The method proposed in the present paper is applicable also to models with randomness. For example, the method will be effective for systems whose dynamical exponent depends on the parameters, such as the random Ising model in random transverse field Pich et al. (1998); Rieger and Kawashima (1999), where the dynamical exponent may take an irrational value or even becomes infinite Fisher (1995). It is thus extremely difficult to analyze the properties of quantum criticality by the conventional strategies. By using our method, one does not need any assumption about the dynamical exponent, which should be quite effective for the systematic investigation of the randomness-driven quantum critical point causing extreme space-time anisotropy.

## Acknowledgments

The simulation code used in the present study has been developed based on an open-source implementation of the worm algorithm WOR () using the ALPS Library Bauer and et al. (2011); ALP () and the BCL (Balance Condition Library) Suwa and Todo (2010); BCL (). The authors acknowledge the support by KAKENHI (No. 23540438, 26400384) from JSPS, the Grand Challenge to Next-Generation Integrated Nanoscience, Development and Application of Advanced High-Performance Supercomputer Project from MEXT, Japan, the HPCI Strategic Programs for Innovative Research (SPIRE) from MEXT, Japan, and the Computational Materials Science Initiative (CMSI). S.Y. acknowledges the financial support from Advanced Leading Graduate Course for Photon Science (ALPS).

## Appendix A Convergence by the Robbins-Monro algorithm

In this appendix, the time evolution of the probability distribution function driven by the Robbins-Monro iteration [Eq. (9)] is discussed. Let us consider the situation in which the physical quantity is obtained by Monte Carlo simulation and the distribution function of is given by a Gaussian (normal) distribution written as

 P(A) =N(f(β),σ2) =1√2πσ2exp[−(A−f(β))22σ2]. (21)

Also, we assume that can be expanded as

 f(β)≈a(β−βc), (22)

near , the zero of . Then, the asymptotic recursion relation between the distribution functions of and is written as

 Pn+1(β(n+1))∝∫dβ(n)Pn(β(n)) ×exp{−12σ2[np(β(n)−β(n+1))−a(β(n)−βc)]2}. (23)

In our procedure, we start from an initial condition . Here we assume that the distribution function of will be approximated by a Gaussian for large enough . Then the recursion relations for the mean, , and the variance, , of are obtained as

 μn+1=(1−apn)μn+apnβc (24) σ2n+1=σ2(pn)2+σ2n(1−apn)2, (25)

respectively. Eq. (24) can be rewritten as

 μn+1−βc=(1−apn)(μn−βc). (26)

The absolute value of in Eq. (26) is less than 1 for sufficiently large . Then for . Similarly, it can be proved that for .

Next, let us assume the leading term of as , where is an unknown constant and is the asymptotic variance of . From the approximation for , the lowest-order terms in Eq. (25) are evaluated as

 −s2αsnαs+1≈σ2p2n2−2s2apnαs+1. (27)

When , dominates over . Then and

 σ2n∼1n2apfor ap<12. (28)

When , on the other hand,

 σ2n≈s2n=1nσ2p22ap−1for ap>12. (29)

There is no solution under the assumption as for the case where , but it is expected that only some correction from will appear.

Similar discussion holds for the mean of distribution, . Assuming , with some constants and , we obtain

 βc+knαm(1−αmn)=(1−apn)(βc+knαm)+apnβc (30)

from Eq. (24). This results in , and thus we have

 μn−βc∼1napfor ap>0. (31)

According to the asymptotic forms (28), (29), and (31), let us discuss the dependence of the final error on the number of Monte Carlo steps. It is obvious from Eqs. (28) and (29) that it is better to set the gain large enough so that is satisfied. Otherwise the convergence of the variance becomes slower. We will call one iteration of the Robbins-Monro feedback process [Eq. (9)] a Robbins-Monro step (RMS), and suppose that a whole calculation consists of  RMSs and each RMS has Monte Carlo updates. The total computational cost is proportional to . From Eq. (29), the variance of the estimate after RMSs is given by for . Using , where is the asymptotic variance of the Monte Carlo estimation, we obtain

 σNR∼s2MCNtot. (32)

This means that the asymptotic variance depends only on the total number of Monte Carlo updates .

Next, let us discuss the optimal choice of the gain . For , the convergence of is faster than that of . Thus, we should minimize in Eq. (29); then we derive

 p=1a. (33)

Fig. 8 shows the dependence of and calculated from an exemplary case:

 f(β)=exp(β)−1, (34)

where and . While the convergence of the mean becomes faster monotonically as increases, the variance takes a minimum value at and increases again for larger . In the simulation presented in the main text, we performed short preparatory calculations for small systems in order to roughly estimate , and set for succeeding long runs. We performed several hundreds of Robbins-Monro iterations, each of which has worm updates.

## Appendix B Measurement of off-diagonal correlation in the worm algorithm

We explain the way to measure off-diagonal correlation in the worm algorithm. Let us begin with measuring the static structure factor defined in Eq. (15). This quantity is easily rewritten by the spin raising and lowering operators. We then need to evaluate the thermal average,

 (35)

The numerator and denominator in Eq. (35) correspond to the partition function of extended world-line configurations with a worm and those of the original world-line configurations, respectively. In other words, Eq. (35) can be read as the frequency of events that the worm head visits the same imaginary time with the tail. We thus can measure the static structure factor by simply counting the frequency of extended configurations that contribute to the numerator in Eq .(35) during each worm update.

We have also evaluated the dynamic structure factor at imaginary frequency . It is given by the canonical correlation function as

 C(→q,iω) =1L2β⟨∑j,k∫β0dτ1∫β0dτ2S+j(τ1)S−k(τ2) ×exp{−i[ω(τ2−τ1)+→q⋅(→rk−→rj)]}⟩, (36)

where is the lowest Matsubara frequency in our simulation. The spatial phase factor can be calculated and stored in advance. Then, when the head moves in the worm-update process, the imaginary-time integral is evaluated. In simulations, every time the head reaches a bond operator, a part of the imaginary-time integral is performed. Since the -dependent part of the integrand is simply given by , we can evaluate the part of the integral exactly at each head move. The spatial phase factor is multiplied (if necessary). The matrix elements of the head and tail also need to be considered (they are simply one in the case with ). The contribution to the integral at each head move is summed up until the head returns back to its tail. The average value of the summed integral for each worm insertion will provide the target quantity (36).

The transverse susceptibility (16) is expressed as . The correlation length can be estimated by the second-moment method Cooper et al. (1982); Todo and Kato (2001); the correlation length in the direction, , is expressed as

 ξx=1|δ→qx|√C(→q0,0)C(→q0+δ→qx,0)−1, (37)

where and . Similarly, the correlation length in the imaginary-time direction, , as

 ξτ=1ω√C(→q0,0)C(→q0,iω)−1. (38)

## References

• Sachdev (1999) S. Sachdev, Quantum Phase Transition (Cambridge University Press, Cambridge, 1999).
• Avella and Mancini (2013) A. Avella and F. Mancini, eds., Strongly Correlated Systems: Numerical Methods (Springer-Verlag, Berlin, 2013).
• Suzuki (1976) M. Suzuki, Prog. Theor. Phys. 56, 1454 (1976).
• Vojta (2003) M. Vojta, Rep. Prog. Phys. 66, 2069 (2003).
• Suzuki (1994) M. Suzuki, ed., Quantum Monte Carlo Methods in Condensed Matter Physics (World Scientific, Singapore, 1994).
• Landau and Binder (2005) D. P. Landau and K. Binder, A Guide to Monte Carlo Simulations in Statistical Physics, 2nd ed. (Cambridge University Press, Cambridge, 2005).
• Fisher et al. (1989) M. P. A. Fisher, P. B. Weichman, G. Grinstein,  and D. S. Fisher, Phys. Rev. B 40, 546 (1989).
• Singh and Rokhsar (1992) K. G. Singh and D. S. Rokhsar, Phys. Rev. B 46, 3002 (1992).
• Herbut (1997) I. F. Herbut, Phys. Rev. Lett. 79, 3502 (1997).
• Weichman and Mukhopadhyay (2007) P. B. Weichman and R. Mukhopadhyay, Phys. Rev. Lett. 98, 245701 (2007).
• Sørensen et al. (1992) E. S. Sørensen, M. Wallin, S. M. Girvin,  and A. P. Young, Phys. Rev. Lett. 69, 828 (1992).
• Zhang et al. (1995) S. Zhang, N. Kawashima, J. Carlson,  and J. E. Gubernatis, Phys. Rev. Lett. 74, 1500 (1995).
• Prokofâev and Svistunov (2004) N. Prokofâev and B. Svistunov, Phys. Rev. Lett. 92, 015703 (2004).
• 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).
• Álvarez Zúñiga et al. (2015) J. P. Álvarez Zúñiga, D. J. Luitz, G. Lemarié,  and N. Laflorencie, Phys. Rev. Lett. 114, 155301 (2015).
• Hüvonen et al. (2012) D. Hüvonen, S. Zhao, M. Månsson, T. Yankova, E. Ressouche, C. Niedermayer, M. Laver, S. N. Gvasaliya,  and A. Zheludev, Phys. Rev. B 85, 100410(R) (2012).
• Yu et al. (2012) R. Yu, C. F. Miclea, F. Weickert, R. Movshovich, A. Paduan-Filho, V. S. Zapf,  and T. Roscilde, Phys. Rev. B 86, 134421 (2012).
• Rieger and Young (1994) H. Rieger and A. P. Young, Phys. Rev. Lett. 72, 4141 (1994).
• Pich et al. (1998) C. Pich, A. P. Young, H. Rieger,  and N. Kawashima, Phys. Rev. Lett. 81, 5916 (1998).
• Pollet (2013) L. Pollet, Comptes Rendus Physique 14, 712 (2013).
• Jiang (2011) F.-J. Jiang, Phys. Rev. B 83, 024419 (2011).
• Jiang (2012) F.-J. Jiang, Phys. Rev. B 85, 014414 (2012).
• Machta et al. (1995) J. Machta, Y. S. Choi, A. Lucke, T. Schweizer,  and L. V. Chayes, Phys. Rev. Lett. 75, 2792 (1995).
• Tomita and Okabe (2001) Y. Tomita and Y. Okabe, Phys. Rev. Lett. 86, 572 (2001).
• Wang and Landau (2001) F. Wang and D. P. Landau, Phys. Rev. Lett. 86, 2050 (2001).
• Yasuda and Todo (2013) S. Yasuda and S. Todo, Phys. Rev. E 88, 061301(R) (2013).
• Wenzel and Janke (2009) S. Wenzel and W. Janke, Phys. Rev. B 79, 014410 (2009).
• Wenzel et al. (2008) S. Wenzel, L. Bogacz,  and W. Janke, Phys. Rev. Lett. 101, 127202 (2008).
• Robbins and Monro (1951) H. Robbins and S. Monro, Ann. Math. Stat. 22, 400 (1951).
• Bishop (2006) C. M. Bishop, Pattern Recognition and Machine Learning (Springer, New York, 2006).
• Albert and Gardner (1970) A. E. Albert and L. A. Gardner, Stochastic Approximation and Nonlinear Regression (MIT Press, Cambridge, 1970).
• Lieb et al. (1961) E. Lieb, T. Schultz,  and D. Mattis, Annals of Physics 16, 407 (1961).
• Hen et al. (2010) I. Hen, M. Iskin,  and M. Rigol, Phys. Rev. B 81, 064503 (2010).
• Hen and Rigol (2009) I. Hen and M. Rigol, Phys. Rev. B 80, 134508 (2009).
• Prokof’ev et al. (1998) N. V. Prokof’ev, B. V. Svistunov,  and I. S. Tupitsyn, Sov. Phys. JETP 87, 310 (1998).
• Syljuasen and Sandvik (2002) O. F. Syljuasen and A. W. Sandvik, Phys. Rev. E 66, 046701 (2002).
• Kawashima and Harada (2004) N. Kawashima and K. Harada, J. Phys. Soc. Jpn. 73, 1379 (2004).
• Suwa and Todo (2010) H. Suwa and S. Todo, Phys. Rev. Lett. 105, 120603 (2010).
• Cooper et al. (1982) F. Cooper, B. Freedman,  and D. Preston, Nucl. Phys. B 210[FS6], 210 (1982).
• Todo and Kato (2001) S. Todo and K. Kato, Phys. Rev. Lett. 87, 047203 (2001).
• Campostrini et al. (2001) M. Campostrini, M. Hasenbusch, A. Pelissetto, P. Rossi,  and E. Vicari, Phys. Rev. B 63, 214503 (2001).
• Rieger and Kawashima (1999) H. Rieger and N. Kawashima, Eur. Phys. J. B 9, 233 (1999).
• Fisher (1995) D. S. Fisher, Phys. Rev. B 51, 6411 (1995).
• (45)
• Bauer and et al. (2011) B. Bauer and et al., J. Stat Mech. , P05001 (2011).
• (47)
• (48)
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters