Velocity of excitations in ordered, disordered and critical antiferromagnets

# Velocity of excitations in ordered, disordered and critical antiferromagnets

Arnab Sen Department of Theoretical Physics, Indian Association for the Cultivation of Science, Jadavpur, Kolkata 700032, India    Hidemaro Suwa Department of Physics, University of Tokyo, Tokyo, Japan Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA    Anders W. Sandvik Department of Physics, Boston University, 590 Commonwealth Avenue, Boston, Massachusetts 02215, USA Department of Physics, National Taiwan University, Taipei 10607, Taiwan
July 12, 2019
###### Abstract

We test three different approaches, based on quantum Monte Carlo simulations, for computing the velocity of triplet excitations in antiferromagnets. We consider the standard one- and two-dimensional Heisenberg models, as well as a bilayer Heisenberg model at its critical point. Computing correlation functions in imaginary time and using their long-time behavior, we extract the lowest excitation energy versus momentum using improved fitting procedures and a generalized moment method. The velocity is then obtained from the dispersion relation. We also exploit winding numbers to define a cubic space-time geometry, where the velocity is obtained as the ratio of the spatial and temporal lengths of the system when all winding number fluctuations are equal. The two methods give consistent results for both ordered and critical systems, but the winding-number estimator is more precise. For the Heisenberg chain, we accurately reproduce the exactly known velocity. For the two-dimensional Heisenberg model, our results are consistent with other recent calculations, but with an improved statistical precision; . We also use the hydrodynamic relation between , the spin stiffness , and the transversal susceptibility , using the smallest non-zero momentum . This method also is well controlled in two dimensions, but the cubic criterion for winding numbers delivers better numerical precision. In one dimension the hydrodynamic relation is affected by logarithmic corrections which make accurate extrapolations difficult. As an application of the winding-number method, for the quantum-critical bilayer model our high-precision determination of the velocity enables us to quantitatively test, at an unprecedented level, field-theoretic predictions for low-temperature scaling forms where enters. We find agreement to within with expansions for the coefficients of the leading susceptibility and specific heat forms; and .

###### pacs:
75.10.Jm, 75.40.Mg, 75.30.Ds, 75.40.Cx

## I Introduction

Ordered quantum antiferromagnets exhibit linear dispersions of their elementary spin-wave (magnon) excitations and the associated velocity is an important parameter characterizing such systems. The prototypical two-dimensional Heisenberg model has an ordered ground state and its low-energy magnon spectrum is well described by spin wave theory.manousakis91 () When corrections are properly taken into account, the velocity and other properties computed within this approximation for the most extreme (and interesting) case of the spin model agrees to within with results of quantum Monte Carlo (QMC) calculations, which can be considered exact to within statistical errors if sufficiently large lattices are used for extrapolations to the thermodynamic limit.runge92 (); wiese94 (); sandvik97 (); sandvik10a (); jiang11a () This good agreement can be traced to the fact that the ground state is strongly ordered, the sublattice magnetization being reduced by quantum fluctuations by only about from the classical value. Upon introducing other interactions which enhance the quantum fluctuations and suppress the order, the quantitative predictive power of spin-wave calculations rapidly deteriorates and the quantum fluctuations have to be treated in more sophisticated ways. singh88 (); millis93 (); chubukov95 (); sandvik94 (); sandvik99a () An extreme case is when a system is driven to criticality. The low-energy critical excitations are still linearly dispersing but the corresponding velocity cannot be reliably calculated in any simple theoretical manner. Lastly, quantum-disordered antiferromagnets also have propagating triplet excitations, which are gapped and often called triplons.

In this paper we will discuss three very different ways to extract the velocity of the elementary excitations of quantum spin models based on ground-state projector and finite-temperature QMC simulations. We consider one-dimensional (1D) and two-dimensional (2D) ordered, disordered and critical quantum antiferromagnets. Using imaginary-time dependent spin correlation functions, the long-time behavior computed at different momenta contain information on the dispersion relation, from which the velocity can be extracted if the limits of the system size going to infinity and the momentum going to zero are treated correctly. One can also in some cases extract the velocity in a simpler, indirect way using winding number fluctuations in the space-time representation of the system sampled in the QMC calculations.kaul08 (); jiang11b () We will develop stable procedure based on these two approaches and compare the results for a number of different cases. In addition, we also test the well known hydrodynamic relationship for an antiferromagnetic state,halperin69 () where is the transversal magnetic susceptibility and the spin stiffness.

We first successfully test the methods on both 1D and 2D systems for which the velocity (of spinons and spin waves, respectively) is previously known, and thereafter study a bilayer system, both in its paramagnetic phase and at its quantum-critical point. In the latter case we subsequently use our high-precision estimate of to investigate in detail, to our knowledge at an unprecedented level of control, the reliability of finite-temperature quantum-critical scaling forms for the magnetic susceptibility and the specific heat.chakravarty89 (); chubukov94 ()

Before turning to the calculations and results, in Sec. I.1 we first provide some more background on the utility of carrying out precise determinations of the velocity in quantum spin systems, focusing on quantum-criticality in dimerized antiferromagnets. In I.2 we provide a brief summary of the different calculations to be presented and in I.3 we outline the organization of the rest of the paper.

### i.1 Effective low-energy descriptions

Quantum field theories are often used to describe universal low-energy properties of quantum magnets.chakravarty89 (); chubukov94 () Numerical techniques, such as QMC simulations, can be used to extract system-dependent parameters appearing in various predicted forms of physical quantities at low energy. Such an approach of combining quantum field theory and numerics has been established over the past several years for certain types of quantum phase transitions in 2D systems.kaul13 (); kim98 () Most well studied are transitions in dimerized models, where for weak inter-dimer coupling there is a tendency to singlet formation on the dimers, leading to a quantum phase transition into a quantum paramagnet at a critical ratio of the inter- and intra-dimer couplings.singh88 (); sandvik94 (); shevchenko00 (); matsumoto01 (); brenig06 (); wang06 (); wenzel09 (); fritz11 () A similar transition (of the same universality class) also take place if the dimers are replaced by some other unit cell of an even number of spins on which the single-cell ground state is a singlet.troyer96 (); albuquerque08 () The many studies of these systems have shown rather convincingly that the phase transition is in the 3D O universality class, in agreement with field-theory descriptions based on the non-linear sigma model.chubukov94 (); chakravarty89 () Other properties associated with quantum-criticality are also well captured by the field theory,chubukov94 () e.g., the uniform magnetic susceptibility is linear in temperature at the critical coupling ratio, , and the specific heat grows quadratically; at low . Perhaps surprisingly, however, the values of the prefactors and have still not been tested quantitatively in a completely unbiased manner.sandvik11 () This may be largely because, as indicated above, they depend on the velocity of the critically damped spin waves (and on no other low-energy parameter), but this parameter has not yet been independently calculated to high precision for model systems (by first-principles methods not depending on other field-theory predictions). As a demonstration of the value of determining to high precision, in this paper we will also provide a test case of a detailed comparison with field theory predictions for one of the prototypical dimerized models; the Heisenberg bilayer.

### i.2 Technical and Physics Objectives

The first aim of the present paper is to systematically test three completely different ways of extracting the velocity of elementary excitations based on QMC calculations in the following ways: (i) Using the momentum dependent gaps extracted from imaginary-time dependent spin correlation functions. (ii) By monitoring spatial and temporal winding number fluctuations, which are proportional to the spin stiffness and the uniform susceptibility, respectively, and adjusting the space-time aspect ratio such that these fluctuations are equal. At this special inverse temperature the ratio of the spatial and temporal lengths should equal .kaul08 (); jiang11b () (iii) Exploiting the hydrodynamic relation , where one has pay attention to the fact that when the temperature in a finite system, while the spin stiffness remains nonzero. We discuss an approach to circumvent this problem based on where the momentum is small but non-zero.

The method (i) is in principle very direct, being connected to the fundamental definition of in the dispersion relation of the lowest-energy excitation versus momentum. However, the precise determination of the momentum-dependent gap is in practice complicated by the presence of a continuum of excitations above the lowest energy. In some cases, as we will show, one also has to take great care with the way the thermodynamic limit and zero-momentum limits are taken. The method (ii) is rather simple and the only uncertainty is introduced by a final extrapolation of the finite-size velocity definition to the thermodynamic limit. However, as far as we are aware, the correctness of this approach has not been formally proven, except for the case of a long-range ordered antiferromagnet, and the method has not been widely used.kaul08 (); jiang11b () We here confirm that the method continues to give the correct velocity also when the system is critical, even in the case of the Heisenberg chain, where the low-energy excitations are not even spin-waves but deconfined spin- carrying topological defects (“spinons”). The method (iii) based on generalized hydrodynamics is simple to apply and we will argue that it as well continues to work also in critical systems.

Comparing methods (i)-(iii) for the Heisenberg chain as well as for the well-studied 2D Heisenberg model gives us insights into how to best apply the methods in practice. After these tests, we study the bilayer Heisenberg model at its quantum-critical point and in its paramagnetic regime, using methods (i) and (ii) Also in this case we find good agreement between the two methods at the critical point, but one has to be more careful when defining the velocity based on gaps for finite systems, because of a slower convergence of the gaps to their infinite-size values than in the ordered state. We also explicitely show that the winding number estimator gives an incorrect velocity in the gapped paramagnetic phase.

The second aim of the paper connects to the low-energy filed-theory description discussed in Sec. I.1—to compute a high-precision value for for the quantum-critical Heisenberg bilayer model, and to use this value to reliably test the field-theory predictions for and . While many such tests have been performed in the past, on the bilayersandvik94 (); shevchenko00 (); brenig06 () as well as other troyer96 () quantum-critical 2D spin systems, the past studies were not able to completely quantitatively test the level of agreement with the existing large- field theory predictions, because an independent, unbiased determination of was lacking. This obstacle is overcome by the reliable calculation of in this paper.

### i.3 Outline of the Paper

The rest of the paper is organized as follows: In Sec. II we discuss calculations of imaginary-time correlation functions within a ground-state projector QMC method formulated in the basis of valence bonds. We discuss our methods to extract the lowest momentum-dependent gaps from the long-time behavior of these correlations. In addition to fitting a leading exponential decay with additional corrections to account for higher excitations, we also present a generalization to ground-state projection QMC data of a systematic high-order moment approach recently introduced for use with finite-temperature QMC simulations.SuwaT2014 () We compare our results with exact solutions for small Heisenberg chains as well as with the rigorously known velocity of this model in the thermodynamic limit. In Sec. III we discuss the determination of based on winding numbers and demonstrate the method using the Heisenberg chain. These finite-temperature calculations are carried out with the stochastic series expansion (SSE) QMC method. In Sec. IV we discuss details of the hydrodynamic relationship generalized to finite system size and test it using SSE calculations for the Heisenberg chain. In Sec. V we present further tests of all three methods for the standard 2D Heisenberg model. We also show that the scaling of the triplet gap at momentum is entirely consistent with quantum rotor-states carrying spin . In Sec. VI we discuss results for the critical and disordered bilayer systems, including the comparison of determined using both methods. The analysis of the quantum-critical susceptibility and specific heat computed using large-scale SSE calculations is presented in Sec. VII. We briefly summarize and discuss our results further in Sec. VIII.

## Ii Momentum-dependent spin gaps

Here we discuss how to use zero temperature () projector QMC methods to calculate imaginary-time correlation functions which are then used to extract the triplet gap (with respect to the ground state energy) as a function of momentum . Applying the imaginary-time evolution operator to a trial state leads to the ground state when . Alternatively, one can use a high power of the Hamiltonian and reach the ground state when . In practice QMC simulation methods based on these two operators are very similar, but the exponential form allows more direct access to the standard imaginary time. As we will show, this connection is not needed if only the excitation energies are of interest, and then one can use the slightly faster power method. We will use both approaches here, with two different ways of analyzing the correlation functions.

### ii.1 H-power projection

This QMC approach is based on projection with a sufficiently high power of the hamiltonian on a trial state . The process can be conveniently expressed in the energy eigenbasis of , leading to the following expression for the dependence on :

 (−H)m|ψt⟩=c0(−E0)m[|0⟩+Λ∑n=1cnc0(EnE0)m|n⟩]. (1)

Here , are the energy eigenstates of and is assumed to be negative, with its absolute value being the largest in magnitude of all the energies (which can always be achieved by adding a suitable constant to the hamiltonian). Then, if the expansion coefficient (which in practice is essentially guaranteed for any reasonable choisce of ), we have for large , and the expectation value of an operator at can be written as

 ⟨O⟩=⟨ψt|(−H)mO(−H)m|ψt⟩⟨ψt|(−H)2m|ψt⟩. (2)

For the SU(2) invariant spin models considered in this paper, this form of the expectation value can be evaluated by importance-sampling, using a formulation of the projector QMC method in the non-orthogonal valence bond basis.sandvik05 () An efficient way of sampling the contributions to , very similar to “operator-loop” updates developed within the finite- SSE method,sandvik99 () can be formulated using loop updates in a combined space of spin components () and valence bonds.sandvikevertz10 () The trial state is also expressed using valence bonds, in the form of an amplitude-product state.liang88 () The details of the state (the form of the amplitudes) are not important, as the good convergence to the ground state is observed even if the state is not optimized.sandvikevertz10 ()

We will use this valence-bond variant of the projector QMC method for computing appropriate imaginary-time dependent correlation functions. For technical details on the sampling methods we refer to Ref. sandvik10a, . Below we will focus on the definition of the correlation functions we study and how we process them to extract the velocity.

#### ii.1.1 Imaginary-time correlations

We consider correlation functions of the following form:

 CA(t)=⟨0|A†(−H)tA|0⟩⟨0|(−H)t|0⟩=Λ−1∑n=0∣∣∣EnE0∣∣∣t|⟨n|A|0⟩|2, (3)

where is an integer which can be related to imaginary time sandvik92 () and we will loosely refer to it by this term. More precisely, , where is the system size, is proportional to imaginary time in the sense of the standard Schrödinger time evolution operator , as we will explicitly show below. From , we further define

 QA(t)=CA(t)−|⟨0|A|0⟩|2CA(0)−|⟨0|A|0⟩|2, (4)

and note that and . For large , we have

 QA(t)→(|⟨1|A|0⟩|2⟨0|A†A|0⟩−|⟨0|A|0⟩|2)(1−Δ|E0|)t, (5)

where is the energy gap between the first excited state (connected to the ground state by the operator ) and the ground state. To directly show the relationship between and imaginary time, we can introduce , then write , and in the limit of large have

 QA(τ=t/|E0|)∝(1−ΔN|e0|)N|e0|τ→ e−Δτ, (6)

which is the familiar form of the asymptotic decay of an imaginary-time dependent correlation function. Thus, we have shown that, indeed, . We will not explicitly need to make use of the relationship between and here, however, and we will continue to use as the “time” parameter with the -power method.

We can appropriately choose the operator so that it excites the ground state into a state with desired quantum numbers. The ground state of the unfrustrated hamiltonians considered here are total spin singlets with momentum on a finite lattice with an even number of spins and periodic boundary conditions (in one dimension only when the number of sites is a multiple of four—for other even the momentum is ).sandvik10b () Since we are interested in triplet excitations, we can use the simple Fourier-transformed spin operator,

 A(k)=Sz(k)=1√N∑reik⋅rSz(r), (7)

to create an state with momentum and when acting on the ground state. Thus, the following imaginary-time correlation function

 Ck(t) = ⟨0|A(−k)(−H)tA(k)|0⟩⟨0|(−H)t|0⟩ = ⟨ψt|(−H)2m−p−tA(−k)(−H)tA(−H)p|ψt⟩⟨ψt|(−H)2m|ψt⟩

allows us to directly measure the triplet excitation gap as a function of the momentum. The second line in the above equation explicitly shows the form used in the projector QMC calculations, where both and are assumed to be large enough to achieve projection to the ground state for the system sizes considered.

In practice, we will use projection powers and, to achieve good ground-state convergence, require that and both are larger than , where typically ( or higher in some cases). These choices are motivated by extensive tests indicating that no detectable systematical errors remain. The values of are restricted to be multiples of in our simulations. Since the correlation functions at the different values are measured in the same simulation, the data are correlated and measuring at shorter intervals does not significantly increase the amount of statistical information in the data set.

#### ii.1.2 Extracting the gap

We here use two different ways to extract the lowest triplet gap from the correlation function . Since for , Eq. (4) reduces to

 Qk(t)=Ck(t)Ck(0). (9)

In principle, the gap can be extracted by monitoring the long-time behavior, given by the form (5). A systematic way to extract the gap without performing any curve fits is to consider the ratio of at two different times separated by some interval; e.g., by operations:

 (10)

Note that when is large enough, which follows from the long-time behavior of in Eq. (5). However, in our QMC calculations, it is not always possible to reach perfect convergence of this gap estimate for all , because the relative statistical errors often become too large already for moderately large . This problem is related to the existence of a continuum (for large ) of states above the gap, due to which the pure exponential decay cannot be easily observed in practice.

We thus use another method to estimate the value of the gap based on the entire available set of correlation functions. This scheme is more reliable than the ratio scheme when data (with small relative errors) are not available for large values of imaginary time. It is clear from Eq. (3) that one can define a positive-definite spectral function to fit the normalized imaginary-time correlation as

 Qk(t)=∫∞0dωAk(ω)(1−ω|E0|)t. (11)

This is just the analogue for the -power evolution of the standard form,

 Gk(τ)=∫∞0dωSk(ω)e−ωτ, (12)

relating the imaginary-time Schrödinger evolved correlation function

 Gk(τ)=⟨0|Sz−k(τ)Szk(0)|0⟩, (13)

to the dynamic spin structure factor . In Eq. (11) cannot exceed , following from the fact that we have ensured that is the eigenvalue with the largest magnitude. In practice, as in the standard dynamic spin structure factor in (12), the actual dominant spectral weight will be concentrated only to within a window of order .

For any finite system, is a sum of -functions and this can be replaced by a continuum starting at for a large system (or, in some cases, there is an isolated -function at the gap, following by a second gap and then a continuum). With or computed using QMC calculations, the respective relations (11) or (12) can in principle be inverted using numerical analytic continuation. This procedure is very challenging, however, and it is not easy to extract the gap precisely with conventional methods such as the Maximum Entropy method,jarrell96 () though one can extract the main dominant spectral features (and we note that progress in this regard has been made very recently sandvik15 ()). Here our goal is merely to extract the gap, and instead of trying to reproduce the full spectral function we model the excitations by just a small number of -functions. With the precision of typical QMC data, can be normally fitted very well with just a few -functions (typically to ) over the full range of accessible times . With this procedure, we expect the location of the lowest gap to accurately reproduce , while the higher -functions represent approximately the contributions of the continuum. In this fitting procedure, the extracted is to some extent affected by contributions of the higher states but does not change significantly when increasing the number of -functions. We can therefore quite reliably extract the lowest gap, but not higher ones unless they are separated by significant subsequent gaps (which is not expected in the cases of interest here, except well inside the quantum paramagnetic state of the bilayer model).

Given a set of -functions at energies with associated amplitudes normalized so that , one can compute the associated time dependent correlation function in analogy with Eq. (11) as

 Qk(t)=n∑i=1Ai(1−ωi|E0|)t. (14)

Now denoting the corresponding QMC-computed function by and their statistical error by , the goodness of the fit is quantified in the standard way by , based on a set of time points :

 (15)

We here use a uniform grid of time points with separation operations, , up to a point where the relative statistical error of exceeds . The choice of cut-off is not very important as, in any case, the noisy data at very large will not affect the fit from the definition of .

For the extracted gap to be reliable, the contribution of the lowest -function to the fit must be significant at the longest times included. To monitor this long-time weight, we compute the relative contribution of the lowest -function, denoted by , at the time included in the fit:

 A1(tmax)=A1(1−ω1/|E0|)tmax∑ni=1Ai(1−ωi/|E0|)tmax. (16)

This quantity should approach for if the lowest -function is at the gap. It is close to in all the fits reported here, indicating stable extraction of .

The statistical error of the extracted gap is estimated using a bootstrap error analysis. With the underlying QMC data for the correlation function saved as bin averages (with typically ), bootstrap averages are constructed by selecting bins at random (i.e., allowing for the same bin to be selected multiple times). The above fitting procedures are then carried out repeatedly for a large number of these samples, and the standard deviation of the estimates is the statistical error of the gap in our procedure.

As already mentioned, the fluctuations of the QMC data at different times are significantly correlated since these are measured in the same simulation. A statistically correct treatment of the data would require the inclusion of the full covariance matrix (instead of just its diagonal elements) in the definition of . However, much of the covariance is already removed when the time-correlations are normalized [by the denominator Eq. (9)], because the errors are correlated primarily by overall fluctuations in the normalization. Based on test cases, including ones reported below, to obtain fully reliable results it is sufficient to use only the diagonal elements of the covariance matrix and define as in (15).

#### ii.1.3 Tests on the Heisenberg chain

We here illustrate the gap extraction method described in the previous subsection using the example of the Heisenberg spin chain with periodic boundary conditions, where spins interact with nearest neighbor exchange constant ;

 H=JL∑i=1Si⋅Sj. (17)

First, we compare the results of our numerics for the lowest triplet gap at for chain sizes and with exact diagonalization (Lanczos-method) results. As can be seen in Fig 1(a), the quantity defined in Eq. (10) indeed converges to the correct gap value in both the cases. Fig 1(b) shows the normalized imaginary time correlation function for the two system sizes, and Fig 1(c) shows the distribution of the gap error obtained using a bootstrap method (i.e., the simulation data stored as bin averages are re-sampled by selecting bins at random, and the fitting procedures are carried out for each such averaged data set) for , using a fit to three -functions.

This analysis show that the gap obtained from the fit agrees statistically with the exact result, which here is within a standard deviation of the distribution obtained in the bootstrapping procedure, and the distribution itself closely matches a normal distribution. Thus, even with only three -functions in the spectrum (which is clearly a much smaller number than what is contained in the full spectrum) we detect no systematic errors introduced by the fitted functional form, supporting our assertion that the simplified description of the spectrum does not significantly affect the location of the lowest -function (the gap). The number of -functions that should be included in a given case depends on the statistical errors of the QMC data and the actual form of the spectral continuum. To determine , we monitor as a function of increasing and stop when no improvements in the fit are observed.

For the Heisenberg chain, the lower edge of the spectrum of triplet excitations is known rigorously based on the Bethe Ansatz solution.cloizeaux62 () For momentum and the spectrum is linear with velocity . For an infinite chain the triplets are degenerate with singlet excitations, due to the fact that the excitations consist of essentially independently propagating deconfined spinons, each carrying spin . The -spinon continuum is maximally broad at and shrinks to zero at . This leads to larger contributions to the time correlations from the continuum close to , which is directly visible in the QMC data for long chains as, e.g., a slower rate of convergence of to a constant. For example, converges much faster to its asymptotic constant value compared to for a chain of size , as shown in Fig. 2. In the inset of the same figure we also show results for , which is the momentum we use to extract as discussed further below. Here the chain is longer, , and does not converge sufficiently before the error bars become too large. However, the method of fitting a simplified spectral function to still works very well and delivers a gap consistent with , but with much better statistical precision. We therefore exclusively uses the fitting method to obtain the results to be discussed next.

In a finite chain, the ground state is non-degenerate and has momentum . However, there is also a quasi-degenerate state with , which is obtained by adding an Umklapp to the true ground state, and this state becomes exactly degenerate with the ground state when takahashi99 () There is also a corresponding finite-size shift in the excitations close to , which is characteristic of the Heisenberg chain but not present in the model in higher dimensions. The lowest triplet gap in the neighborhood of , where we will use only , behaves for large as caux ()

 Δ(π+2π/L)=Δ(π)+c(L)2πL, (18)

where but with a multiplicative logarithmic correction. In Fig 3 we show the behavior of the corresponding velocity estimate,

 c(L)=L2π[Δ(π+2π/L)−Δ(π)], (19)

as a function of the inverse chain length. A smooth monotonic (asymptotically linear) approach to the known velocity, can be observed as . There are no signs of any remaining log corrections as a function of in this estimate. Our results are consistent with a remaining linear finite-size correction.

The velocity of the spinons can also be extracted from the lowest triplet gap at by using the simple estimator as . Results are shown in Fig 4. Also in this case we observe the estimate approaching the correct value of as , but with a non-monotonic behavior with a maximum for before an apparently linear asymptotic approach to the correct value. It should be noted that the spectral weight (before normalizing the correlation function) vanishes as , which implies that the correlation function computed in the QMC simulations becomes very noisy for large system sizes. It is worth noting here that, because of the non-monotonic behavior seen in Fig 4, a calculation using exact (numerical) diagonalization of the Hamiltonian would appear inconsistent with the known velocity, because large enough system sizes required to go beyond the maximum at cannot be reached.

### ii.2 Exponential-form projection

In this section we will discuss a different way of extracting the gap, using a systematic way to analyze moments of the spectrum based on information contained in the imaginary-time correlations. Such a scheme was recently introduced within QMC calculations,SuwaT2014 () and we here present a generalization to projector QMC. In addition, we discuss improvements in the extrapolations required to obtain unbiased results.

To make the connection with the previous version of the method more transparent we will here use the exponential-form projection with continuous imaginary time. The ground state is again projected out of a trial state in the valence-bond basis, but now using instead of . Taylor expanding the exponential (as was also done previously in projector QMC, e.g., in Ref. farhi12, ), the scheme then closely resembles SSE QMC algorithms.sandvik99 () The normalization, which replaces the partition function in methods, is expressed as

 Z′ = ⟨ψt|e−βH|ψt⟩ = ⟨ψt|Tτexp(−∫β0dτH(τ))|ψt⟩ = ⟨ψt|1+∞∑m=1(−1)m∫β0dτ1∫βτ1dτ2⋯ ⋯∫βτm−1dτmm∏i=1H(τi)|ψt⟩

where indicates time ordering and is a Hamiltonian acting at imaginary time . If the Hamiltonian here is actually time dependent, we obtain the non-equilibrium QMC scheme developed in Ref. degrandi11, . Here is time-independent, and we formally use the time-dependent formalism only as a convenient way of accessing imaginary time in the configuration space directly. In principle, we can also use the time dependence corresponding to the interaction representation,beard96 (); prokofev98 (); sandvik97b () where only the - and parts of the interaction appears in the operator product, but here we stay close to the formulation also used in the SSE method and in the -power method discussed above, and expand with the full Hamiltonian. Note that the -power method is exactly recovered if the time integrals are completed, and one can in practice also equivalently use the power method and just generate the ordered time sequences at random,sandvik97b () or generate them in the process of the updates as in the world-line continuous time algorithm.beard96 ()

#### ii.2.1 Generalized Moment Method

Our approach introduced here to extract the gap is a generalization to projection QMC of the moment method proposed recently.SuwaT2014 () The Fourier-transformed dynamical correlation is exploited to derive a sequence of gap estimators. We will first explain how the moment of the dynamical correlation function is related to the gap, and then derive increasingly precise gap estimators. Finally, we will show how these estimates are extrapolated to the limit in which the exact gap is recovered.

We will use the following dynamical correlation function computed with the exponential projector:

 C(τ,τ0) = 1Z′⟨ψt|e−(β−τ−τ0)HA†e−τHAe−τ0H|ψt⟩ (21) = 1Z′∑ℓ,p,qbℓ,p,qe−βEpe−τ(Eℓ−Ep)e−τ0(Eq−Ep) → ∑ℓ≥1bℓe−τΔℓ(β,τ0→∞),

where we use the definitions

 bℓ,p,q = ¯cpcq⟨p|A†|ℓ⟩⟨ℓ|A|q⟩, (22) bℓ ≡ bℓ,0,0=|c0⟨ℓ|A|0⟩|2, (23) Δℓ = Eℓ−E0. (24)

We assume that , , and , all of which apply here.

Let us first consider the moment of the asymptotic correlation function in the limit :

 I∞n = ∫∞0dττnC(τ,τ0) (25) = ∫∞0dττn∑ℓ≥1bℓe−τΔℓ = ∑ℓ≥1bℓΔn+1ℓn! ∼ b1Δn+11n!(n≫1).

Then the lowest gap can, in principle, be obtained using an appropriate ratio, e.g.,

 (n+1)I∞nI∞n+1→Δ1(n→∞). (26)

In practice, however, the projection time cannot be infinite in simulations, and we have to consider carefully the effects of finite . In addition, the range of integration, , is different from (less than) because the imaginary-time correlations are measured from the reference point at the center of the projected trial state [noting that Eq. (II.2) corresponds to projection of the bra and ket state with ], or, alternatively, between two points located symmetrically within the projection range , and one has to stay well away from the boundaries (trial states) for unbiased measurements. Moreover, in many cases the statistical errors grow too large at large times, as mentioned previously. Thus, in practice .

The moments for finite and take the form

 In =∫βint0dττnC(τ,τ0) =1Z′∫βint0dττn∑ℓ,pdℓ,pe−τΔℓ,p =1Z′∑ℓ,pdℓ,pn!Δn+1ℓ,p(1−n∑m=0P(m)), (27)

where

 dℓ,p = e−βintEp∑qbℓ,p,qe−τ0(Eq−Ep), (28) Δℓ,p = Eℓ−Ep, (29)

and

 P(m)=(βintΔℓ,p)me−βintΔℓ,pm! (30)

is a properly normalized Poisson distribution;

 ∞∑m=0P(m)=1. (31)

Owing to the finite integration range, the Poisson term does not vanish. As a consequence, the ratio of the moments for finite does not contain information of the gap in the limit . Instead, we have a completely different limiting behavior,

 InIn+1∼n+2βint(n+1)→1βint(n→∞), (32)

independent of the gap. We can overcome this difficulty and devise a proper gap estimator by using the Fourier transformation. Here we express the moments in a different way using the following expansion:

 I∞2n(2n)! =limβint,τ0→∞(−1)nω2n1n∑k=0xn,k,0R(ωk) (33) I∞2n−1(2n−1)! =limβint,τ0→∞(−1)n−1ω2(n−1)1n∑k=1xn,k,1J(ωk)ωk, (34)

where and are the real and imaginary parts, respectively, of the Fourier-transformed correlation function, i.e.,

 ∫βint0dτeiτωkC(τ,τ0)=R(ωk)+iJ(ωk), (35)

where , and the key coefficients are and

 xn,k,m=1∏m≤j≤n,j≠k(k+j)(k−j). (36)

When deriving these equations we have considered the expansion in on the right-hand side of Eq. (33) and Eq. (34), where the lowest orders then cancel. The coefficients can be solved for by using the inverse of the Vandermonde matrix.Neagoe1996 ()

Combining the results above, we obtain the following improved gap estimator:

 ^Δ(n,βint)=−ω21n∑k=1xn,k,1J(ωk)ωkn∑k=0xn,k,0R(ωk). (37)

Remarkably, this estimator is asymptotically unbiased and the limits are interchangeable (see Appendix for the analytical derivation):

 Δ1 =limn→∞limβint,τ0→∞^Δ(n,βint) =limβint,τ0→∞limn→∞^Δ(n,βint). (38)

Due to the commuting limits, observing convergence of the estimator (37) in large and taken in any convenient fashion will deliver an unbiased result for the gap.

In principle there is also a dependence on the value used in the exponential projection, in addition to the dependence on and . Above we have assumed that is sufficiently large for quantities computed by integration up to have converged to the limit, and in the QMC simulations we also monitor this convergence. Naturally, the value of required grows with .

#### ii.2.2 Performance Check for the Heisenberg Chain

We have checked the validity of the estimator (37) for the Heisenberg chain of 16 spins with periodic boundary conditions. The projection length is set large enough to see the asymptotic behavior of the dynamic correlations on the imaginary-time axis. The correlations were measured from the center of the QMC configuration; that is, from . As in the previous section, the loop algorithm and improved correlation-function estimator sandvik10a () are used. We consider the triplet gap at and the operator in Eq. (21) is then explicitly given by and

 A†eτHA=∑α=x,y,z∑r,r′SαreτHSαr′(−1)r−r′. (39)

The Fourier-transformed correlation functions (35) were directly measured in the simulation, and no discretization error is introduced. The gap estimators (37) for several and were calculated and the errors estimated by the jackknife analysis Berg2004 () (bootstrapping would produce statistically equivalent result).

As shown above, our estimator (37) converges to the exact gap in the limit of infinite and . In practice, if good convergence is observed within statistical errors, a gap value obtained from a sufficiently large and within the converged range can be used as the final estimation, or some appropriate functional form can be used for extrapolation. However, an important issue here is that the statistical errors grow with the two parameters. We therefore inevitably encounter a trade-off problem between systematical and statistical extrapolation errors. Here we will demonstrate that the statistical error can be optimized by proper extrapolations while keeping the estimation unbiased. In the procedure used below, the extrapolation is taken for first and then for , though, as discussed above, other ways to accomplish the limit are also possible.

Examples, based on () Monte Carlo samples, of extrapolations are shown in Fig. 5. The leading finite- correction is linear in [in accord with Eq. (67) and (73) in the Appendix], and we include also a quadratic term to obtain good fits to the data for the full -range. The limit is taken next using the resulting values of the extrapolations, as shown in Fig 6. Here an exponential function is used for the data fit (again, according to results derived in the Appendix) to extrapolate the final result for the gap. Though other ways of extrapolating to are possible, the above protocol is convenient because the extrapolation for , which is taken first, is relatively easier than that for .

As a test of the unbiased nature of the extrapolation scheme, distributions of the relative gap error are shown in Fig. 7. Histograms were collected based on 2048 independent simulations of the Heisenberg chain. Results based on the extrapolation procedure discussed above are compared with those of individual gap estimators for and . The estimator clearly has a non-zero systematic error remaining, which is similar to (but smaller than) the conventional second moment estimator.SuwaT2014 () The estimator has a small enough systematical error (the histogram being centered very close to zero) but has a large statistical error (wide distribution). The extrapolated estimation is unbiased and has a small statistical error.

## Iii The velocity from winding number fluctuations

We here discuss how to use winding numbers to compute the velocity. This method has been known for some time,kaul08 () and was recently applied to the 2D XY model jiang11b () and the 2D Heisenberg model.jiang11a () We begin by briefly recollecting some key aspects about winding numbers in finite-temperature QMC simulations, in particular the SSE method we use for these calculations. We then present results for the Heisenberg chain.

### iii.1 Winding numbers in QMC simulations

QMC simulations at finite temperature are based on some mapping of the partition function of a quantum system in dimensions to an effectively equivalent -dimensional classical system, where the additional dimension of the configurations corresponds to imaginary (Euclidean) time. The effective length of the system in the time dimension is , where (setting ) and the configurations are time-periodic. For a system with conserved particle number, which in the case of spin model corresponds to conserved magnetization along the quantization () axis, imposing periodic boundary in the spatial dimensions leads to another, topological number associated with the configurations—the winding number representing permutations when particles are propagated once or multiple times around the periodic system. The winding numbers were first used in QMC calculations of the superfluid stiffness of interacting bosons.pollock87 ()

The SSE QMC algorithm,sandvik91 (); sandvik92 (); sandvik97 () is based on a Taylor expansion of the imaginary-time evolution operator (the Boltzmann operator) , and similarities with the projector approach discussed in the previous section were already pointed out. Each SSE configuration is associated with some power of the Hamiltonian propagating a basis state (here in the standard computational basis of spin components), and these powers are sampled stochastically to all contributing orders. The trace over all basis state is also sampled. The average expansion power in this procedure is proportional to ; . In simulations, the state propagation is broken up into individual paths corresponding to strings of of the individual local terms of the Hamiltonian, forming successions of evolving basis states, similar to those in path integrals. For a Heisenberg model the Hamiltonian terms are the diagonal operators (in practice with a constant added) and off-diagonal operators, the latter of which transport spin and are associated with currents in the lattice direction corresponding to the site-pair . The winding number in the lattice direction is defined in terms of the currents as

 Wa=1Lan∑p=1Ja(p), (40)

where the index corresponds to the location of the transport “event” in the string of operators and is the length of the lattice in the direction. Defined in this way, the winding numbers are integers. It should be noted that the definition of the winding number is exactly the same in SSE and world-line methods suzuki77 (); hirsch82 () and one can also think of the SSE configurations as consisting of world lines (for up and down spins in the case of quantum spin systems).

The spatial winding number measures the net spin transported around the periodic lattice in the direction in the course of the periodic propagation in imaginary time. Equivalently, this is the number of world lines (up ones minus down ones divided by two) crossing through a plane drawn along the time axis perpendicular to the axis. Since the total magnetization is conserved, one can also think of as a winding number; the net number of world lines crossing a plane drawn at an arbitrary time point perpendicular to the time axis, which is just the magnetization computed in the stored basis state;

 Wτ=Mz=N∑i=1Szi. (41)

The expectation values of the squared winding numbers (i.e., the winding number fluctuations) are related to two important thermodynamic quantities; the spin stiffness,

 (42)

and the uniform magnetic susceptibility,

 (43)

The technicalities of implementing these observables in SSE calculations have been discussed extensively in the literature; see Ref. sandvik10b, for a recent review.

### iii.2 The cubic criterion and the velocity

In the high-temperature limit , the magnetization fluctuations of any system are maximized and therefore according to Eq. (43). For an unfrustrated antiferromagnet the ground state is a singlet, and, on account of the presence of a singlet-triplet finite-size gap, when for any finite system. In contrast to these limits of , for the spatial winding number in direction () we have when on account of there being no quantum fluctuations when the imaginary-time length and there are no contributions from expansion powers (in the case of SSE—in world-line methods there will similarly be no transport events causing shifts of the world lines). In the limit , for a system with long-range order (or a ”quasi-ordered” 1D system with power-law decaying correlations), the stiffness constant converges to a non-zero value for any , and according to Eq. (42) we must then have a divergence . These different behaviors of the spatial and temporal winding numbers versus temperature guarantees that there is a crossing point,

 ⟨W2r(β∗)⟩=⟨W2τ(β∗)⟩, (44)

at some unique value of for given size .

The winding numbers characterize global fluctuations of the system in the different spatial and temporal directions. It is then natural to define a system as having cubic space-time geometry when Eq. (44) holds (with the lattice length the same in all spatial directions). The aspect ratio should then directly correspond to the velocity of the long-wave-length excitations. In some cases this can be shown directly based on low-energy field theory kaul08 (); jiang11a (); jiang11b () but even in the absence of such descriptions the arguments are very general and one can expect the conclusion to always hold for a system with linear dispersion, though we are not aware of any formal proofs in the general case.

### iii.3 Test on the Heisenberg chain

To find the point where the cubic criterion is satisfied, we simulate a system at several values of in the region where based on initial explorations and knowledge of the approximate value of the velocity. We fit a low-order polynomial (typically second- or cubic-order) to the difference and solve the resulting equation for the -value for which the cubic criterion is satisfied. This procedure is illustrated in Fig. 8 for a Heisenberg chain with spins. From this procedure we obtain , which can be extrapolated to . With independent data points, the statistical error of at fixed can be determined by repeating the procedure multiple times with added Gaussian noise of standard deviation equal to the error bars of the data points, whence the standard deviation of the extracted crossing point is the error bar.

Fig. 9 shows results of such a procedure for several chain lengths , graphed versus . We are not aware of any theoretical predictions for the size dependence of this definition of , but the data for the larger systems are well described by a constant (the final infinite-size value of ) plus a term proportional to . Using this fitting form leads to a value of completely consistent with the known value , as is clear from Fig. 9.

## Iv Hydrodynamic Relationship

A well known way to extract the velocity in an antiferromagnet is to use the analogue of a hydrodynamic relationship between the velocity, the spin stiffness (helicity modulus), and the transversal susceptibility (which is the analogue of a mass density):halperin69 ()

 c=√ρsχ⊥. (45)

In a QMC calculation in which the spin-rotation symmetry is not explicitly broken, one cannot compute directly, but one can use the fact that the rotationally-averaged uniform susceptibility computed in the basis, as in Eq. (43), becomes of a transversal (e.g., ) component when in the thermodynamic limit (since the longitudinal component vanishes at ). Thus, one can obtain as by taking the limit before the limit, while taking the limits in the opposite order does not work because then due to the finite-size gap between the and magnetization sectors. In contrast, for the limit has to be taken before , because the system in the thermodynamic limit only has stiffness (Néel order) exactly at . In this case, too, a factor has to be included in the finite-size estimate to account for rotational averaging of the relevant transversal components.

The procedure to obtain in the thermodynamic limit is relatively straight-forward with an extrapolation of using a polynomial in , while the extrapolations requiring first in the calculation is more cumbersome. Results obtained for in this way sandvik10b () have large error bars compared to the result of the winding number method presented above in Sec. V.1.

In order to define a finite-size velocity estimate based on Eq. (45) directly in the limit we here use a modification of the relationship. We use the susceptibility at finite momentum ,

 χ(q)=1N∫β0⟨Mz(−q,τ)Mz(q,0)⟩, (46)

where the magnetization at non-zero momentum is the same as in Eq. (7) but without the normalization by . We can use the smallest momentum for a given system size and approach the limit as . If we here define without the rotational factor discussed above, then the same rotational factor should also not be used in the susceptibility, and we define the velocity for finite size as

 c(L)=√ρs(L)χ(q=2π/L). (47)

We will compute this quantity using SSE simulations at sufficiently large to achieve ground state convergence for each studied.

### iv.1 Test on the Heisenberg chain

The hydrodynamic relationship (45) is normally applied in the magnetically ordered Néel state and one may then question its use in a critical system. In the case of the Heisenberg chain we can rigorously see that it is a valid way to extract . Since there is no long-range order, there is no distinction between longitudinal and transversal modes, but Eq. (45) defined with rotationally averaged quantities should remain valid. According to the exact Bethe Ansatz solution, in the thermodynamic limit we have, in the units with and all relevant constants set to , (derived in Ref. hamer87, ) and (see., e.g., Ref. eggert94, ) and, thus, according to Eq. (47), we obtain the correct result .

For finite system size, it is well known that both and exhibit strong logarithmic corrections to their infinite-size values which can be traced to the presence of marginal operators and it is very difficult to extrapolate them; see, e.g., Ref. byrnes02, for and Ref. eggert94, for (where in the latter case the logarithmic correction in temperature is discussed but one should expect similar corrections for at ). The logarithmic corrections do not cancel in the ratio of the two quantities in Eq. (47), and as a consequence we also find that it is difficult to extrapolate precisely to the thermodynamic limit. Results are shown in Fig. 10. Although it is not possible to fully extrapolate to the exact results based on these data for system size up to , the result for the largest size nevertheless deviates by only from the exact result.

It is interesting to note that also this definition of exhibits a non-monotonic behavior, with a maximum at , the same as in the value obtained previously from the gap at in Fig. 4. In the latter case, the maximum corresponds to an elevated excitation energy at , which one should expect to lead to a reduction in because this quantity is given by a sum rule of . By this sum-rule, in a single-mode approximation a larger would lead to a smaller , and even beyond the single-mode approximation one should expect such an effect because the dominant spectral weight is at the gap. If is not appreciably affected by this finite-size anomaly, then indeed an elevated gap and reduced in Eq. (47) can explain the maximum in Fig. 10.

In light of the logarithmic corrections to defined according to Eq. (45), the apparent complete lack of any challenging corrections or non-monotonic behavior in the definition based on the cubic criterion for winding numbers (Fig. 9) is even more remarkable, since also this estimate involves quantities directly related to the susceptibility (temporal winding number) and spin stiffness (spatial winding number). There are also no apparent logarithmic corrections in the velocities defined based on gaps in Figs. 3 and 4. In Fig. 3, which is based on the triplet gap close to , the logarithmic corrections are avoided by subtracting the gap in Eq. (18), while the gap close to used in Fig. 4 is not expected to be affected by logarithms. The latter gap at scales as for large according to our results in Fig. 4, which can be explained from the known dispersion , if there is a finite-size correction (due to irrelevant fields only).

## V 2D Heisenberg model

The 2D spin- Heisenberg model on the square lattice spontaneously breaks spin rotation symmetry at in the thermodynamic limit. This leads to gapless linearly dispersing Goldstone modes in the vicinity of the momenta and . The ground state in a finite periodic system is a total spin singlet. The long-range antiferromagnetic order in the thermodynamic limit is reflected in the energies of the quantum rotor states rotorrefs () which collapse onto the ground state as much faster than than the spin wave excitations, . The rotor states, thus, become degenerate with the ground state as the system size increases, and have momenta and for even and odd respectively. Combinations of rotor states with up to can then be formed which are ground states with fixed direction of the Néel vector (the staggered magnetization), thus allowing for symmetry breaking in the thermodynamic limit. Here we provide an accurate determination of the spin-wave velocity and also directly investigate the scaling of the rotor state.

### v.1 Velocity from winding numbers

Since we expect the winding numbers with the cubic criterion to provide the best determination of we start with this approach. Results based on the procedures discussed in the preceding section are shown in Fig. 11 versus the inverse system length. Carrying out polynomial fits, we find that no linear and quadratic terms are required. Fourth- and fifth-order polynomials fits excluding these terms and using all the data give almost identical results for the extrapolated velocity, with only the error bar somewhat larger for the fifth-order fit. The figure shows the fifth-order fit and the extrapolated velocity with it is . While we do not know the physical reason for the leading cubic correction, we use it as the simplest empirical description of the data. Including also a quadratic term, it comes out equal to 0 within statistical errors and the extrapolated result does not change appreciably.

The above value of agrees within errors bars with the recent result using the same method by Jiang and Wiese,jiang11a () but our error bar is almost an order of magnitude smaller. We note that in Ref. jiang11a, no systematic extrapolation was carried out to the thermodynamic limit—instead an average was taken of results for system sizes in the range . Looking at the data in Fig. 11 it is clear that, with the small error bars on the SSE data we have achieved here, an extrapolation is necessary to obtain a result with no remaining finite-size effects. To our knowledge, the above result is the most precise spin-wave velocity reported to date for the 2D Heisenberg model. Spin-wave theory with corrections up to order gave ,canali92 () where the uncertainty in the last digit reflects estimated numerical errors from evaluations of challenging integrals. Thus, to this order, the spin-wave result deviates by only from the correct value.

### v.2 Gap scaling of rotor states

The quantum rotor excitation gap can be directly accessed by measuring the lowest triplet gap at . This energy scale is related to the uniform (transverse) magnetic susceptibility as:rotorrefs ()

 E(S,L)=S(S+1)2L2χ⊥. (48)

From the behavior of the lowest triplet gap up to system sizes , as shown in Fig. 12, we indeed observe that at large , but there are also large corrections which we fit with additional higher-order powers of . The extrapolation to infinite size gives . This is consistent with the value of susceptibility obtained using QMC calculations in small external magnetic fields to extract gaps between different spin sectors.Syljausen02 ()

### v.3 Velocity from gaps

The velocity of the spin waves can be estimated by measuring the triplet gap in the vicinity of and . We here choose to measure the triplet gap at [or, equivalently, at , which we use for averaging] since the lowest triplet excitation energy is at and is the closest allowed wave-vector to for a periodic system with linear size . For this model we have carried out only -power QMC simulations. We illustrate the gap extraction with both the ratio method and the simplified spectral function with fitted -functions in Fig. 13, using the largest system sizes considered; . Again, working with the spectral function produces much more stable results, though clearly the correlation ratio also converges to a constant consistent with the same gap.

Since in the thermodynamic limit, the excitation energy of the spin waves equals in the vicinity of , where is the spin-wave velocity, the estimator

 c(L)=L2πΔ(π+2π/L,π) (49)

should converge to as . We graph this quantity in Fig. 14, and it converges to the value of obtained above from the winding-number method when . Note that there is again a non-monotonicity as a function of , similar to the case of the Heisenberg chain in Fig. 4. In the latter case this behavior only was observed in the gap extracted close to , however, while here we have used the spectrum close to .

### v.4 Hydrodynamic relationship

We next consider the definition of by the hydrodynamic relationship, Eq. (47). We again go to sufficiently low temperature in SSE simulations for any remaining finite- corrections to be negligible compared to the statistical errors, which typically meant . Results are graphed in Fig. 15. Here again we expect the finite size corrections in the Néel state to be described asymptotically by a polynomial in , but a rather high order of the polynomial is required to fit the data for the smaller system sizes. The behavior for the largest sizes again indicate (as in the case of the winding-number calculation discussed above) that the leading correction in is cubic. To get a statistically acceptable fit for all data we include terms up to order , which gives , which deviates by approximately error bars from the statistically much more precise value obtained in Sec. V.1 based on the winding-number fluctuations. By comparing Figs. 11 and 15, it is clear that the analysis in the former is more reliable, with a larger number of data points used and an overall much weaker size dependence. The above result from the hydrodynamic relationship is still likely somewhat affected by systematical errors, as the behavior still appears to be flattening out more than the fit suggests, and the fit even including the 6th-order term is only marginally good, with . It would clearly be desirable to go to larger system sizes, but given the much better behavior of the winding number data it is already clear that this is the preferred method.

## Vi Bilayer Heisenberg model

We now consider the bilayer Heisenberg model which is a prototypical system to study quantum phase transitions in 2D.singh88 (); millis93 (); sandvik94 () The Hamiltonian is given by

 H=J1