Addressing general measurements in quantum Monte Carlo

Domain: physics
Published: 2026-01-02
Editorial outcome: minor_revision

Abstract
Quantum Monte Carlo is one of the most promising approaches for dealing with large-scale quantum many-body systems. It has played an extremely important role in understanding strongly correlated physics. However, two fundamental problems, namely the sign problem and general measurement issues, have seriously hampered its scope of application. We propose a universal scheme to tackle the problems of general measurement. The target observables are expressed as the ratio of two types of partition functions hOi=Z=Z, where Z=trðOe(cid:2)βHÞ and Z=trðe(cid:2)βHÞ. These two partition functions can be estimated separately within the reweight-annealing frame, and then be connected by an easily solvable reference point. We have successfully applied this scheme to XXZ model and transverse field Ising model, from 1D to 2D systems, from two-body to multi-body correlations and even non-local disorder operators, and from equal-time to imaginary-time correlations. The reweighting path is not limited to physical parameters, but also works for space and time. Essentially, this scheme solves the long-standing problem of calculating the overlap between different distribution functions in mathematical statistics, which can be widely used in statistical problems, such as quantum many-body computation, big data and machine learning.

Introduction
Quantum Monte Carlo (QMC) is a highly promising numerical method without approximations for large-scale or high-dimensional quantum many-body systems, capable of simulating complex systems with an exponential degree of freedom while maintaining polynomial computation complexity1–25. Despite the maturity of QMC techniques after decades of development15,26–41, there remain two essential challenges that greatly limit the application of QMC. The first is the notorious sign problem34,35,42–58, and the second is the issue of general (off-diagonal) measurements7,10,16,30,59.

In this work, we will focus on the enduring challenge of measuring general (off-diagonal) observables. The target is how to extract as more as information from the QMC samplings. Unlike other numerical methods, QMC cannot directly obtain the wave-function of ground state. Typically, the evaluation of a physical quantity 〈O〉 in QMC is derived as follows: hOi=trðOe(cid:2)βHÞ=Z, where Z=trðe(cid:2)βHÞ is the partition function (PF), β is the inverse temperature and H is the Hamiltonian. For simplicity, we define Z=trðOe(cid:2)βHÞ, hence hOi=Z=Z. In a standard QMC framework, the partition function Z can be generally decomposed into the sum of all the weights, i.e. Z=∑W. If the operator O can be treated as a number O under the configuration of W, which corresponds to a diagonal measurement, the physical quantity can be readily estimated in the form

hOi= Z=Z = (∑ O W) / (∑ W) (1)

In this way, the value O can be directly obtained when we sample the configurations W of the PF, making diagonal measurements straightforward in the QMC framework. In the case of diagonal measurement, it is clear that two PFs, Z=∑ O W and Z=∑ W, share the same set of configurations {W}, but differ in their associated values, with O for Z and 1 for Z. Consequently, sampling the configurations {W} is sufficient to capture the expectation value hOi=Z=Z.

However, the situation would deteriorate significantly during off-diagonal measurements. Off-diagonal operators typically alter the existing configurations {W} of Z=∑ W, resulting in new configurations {W'} for Z=∑ W' that are entirely distinct from the original set {W}. This implies that we are unable to obtain samples {W'} within the framework of conventional QMC methods, which are designed to sample from {W}. As shown in Fig. 1(a), two PFs no longer share the same configurations, making it impossible to simulate their ratio directly as in the diagonal case. Furthermore, it is usually impossible to design updates between {W} and {W'} in QMC algorithms (If you can realize the updates between {W} and {W'}, the ratio Z=Z then can be obtained, such as the QMC algorithm for entanglement entropy60). This represents the fundamental challenge in the off-diagonal measurements.

In some special cases, certain off-diagonal observables can be extracted in ingenious ways. For instance, determinant QMC (DQMC) can efficiently calculate two-body Green’s functions, where higher-order correlators can be derived using Wick’s theorem11,31,35,61. This works particularly well in fermionic systems, where off-diagonal Green’s functions such as hccyi are naturally embedded in the sampling weight. We now fix the measurement problem in other Monte Carlo methods, such as diagrammatic QMC62–64, world-line QMC and stochastic series expansion (SSE) formulations, in which off-diagonal measurements pose a much more severe challenge. Despite numerous efforts in the past, two-body Green’s functions can be obtained only in certain cases within the framework of worm-like QMC9,65–73. The reason is that the configurations in the worm-like update process can be treated as samplings of the two-point Green’s function. However, multi-body Green’s functions remain challenging to be extracted even with this specialized approach and the worm-like algorithm only works in several models. In addition, if the off-diagonal operator to be measured is a part of the Hamiltonian, it can be estimated through the sampling process74. As an instance, 〈Sx〉 can be measured in a transverse field Ising model (TFIM)75,76. Another example is that, in the stochastic series expansion (SSE) method, the energy value can be calculated directly by counting the number of operators in the space-time configurations2,7,30. Although the importance of off-diagonal observables in quantum systems, there is currently no general method for measuring arbitrary operators in QMC, even though a lot of effort has been devoted to it over the past decades.

Recently, a newly proposed method – reweight-annealing (RA)77 has been successfully applied to determine the ratio of two same-type PFs at different parameters. In the reweight-annealing method, as shown in Fig. 1(b) and (c), the PF at the parameter J0 can be estimated using the value of PF at another parameter J by resetting the weights.

Z(J0)/Z(J) = 〈W(J0)/W(J)〉J (2)

This equation states that the ratio Z(J0)=Z(J) is obtained by averaging W(J0)=W(J) over configurations sampled from the ensemble at J (i.e. 〈⋅〉J). However, this equation works well only when the distributions Z(J0) and Z(J) are adjacent. In this context, the importance sampling can be maintained78. Therefore, if the target parameters J and J0 are far away from each other, a series of intermediate parameters {Ji} need to be inserted to split the reweighting process by gradually moving from J to J0. This can be expressed as Z(J0)=Z(J) = Z(J0)=Z(J1) × Z(J1)=Z(J2) × ... Z(Jn)=Z(J). Since the entire process involves annealing from one parameter to another with iterative reweighting, it is dubbed as “reweight-annealing"77. The similar spirit of reweighting also has been developed in the high-energy physics and other fields79–84. Once a reference point Z(J) is known, Z(J0) can be calculated through the ratio. It has been proved that the computation complexity of the RA method is polynomial if the ratio of two closest Z(J) and Z(J0) is fixed in the division strategies77. More specifically, in the QMC simulation, maintaining the ratio of two closest partition functions Z(J0)=Z(J) with O(1) range (i.e. 0.1~10) is essential – not only to preserve the polynomial computation complexity under importance sampling77, but also to ensure that system error remains controllable (system error analysis can be found in Supplementary Note 7). In addition, the number of slices should be increased near the phase transition to maintain the ratio with O(1) range, since the distribution of the partition function changes more rapidly in this regime. Motivated by the reweighting scheme, we propose a novel scheme termed “bipartite reweight-annealing (BRA)" method to address the challenges of general measurements in QMC simulations. We will present several examples to demonstrate its feasibility and versatility.

Results
Bipartite reweight-annealing framework
In fact, we realize that the reweighting scheme is not only limited to the standard PF Z(J) but can be applied to any distribution that varies with the related parameters. In practice, an off-diagonal observable can be treated as the ratio of two types of PFs hOi=Z=Z, where Z=trðOe(cid:2)βHÞ. This insight inspires us to reweight different kinds of PFs (the numerator Z(J) and denominator Z(J)) respectively, as Fig. 1(b) and (c) show. The key idea is that we firstly calculate the ratios Z(J0)=Z(J) and Z(J0)=Z(J), and if we have a reference point Z(J)=Z(J) which is easily solvable (as displayed in Fig. 1), then the target measurement hO(J0)i can be estimated in this approach:

hO(J0)i = (Z(J0)/Z(J)) / (Z(J0)/Z(J)) × (Z(J)/Z(J)) (3)

can be calculated by reweighting. Likewise, maintaining the ratio within an O(1) range is imperative (detail can be found in Supplementary Note 7).

This BRA scheme avoids the intractable problem of calculating the ratio between two entirely different PFs (Fig. 1(a)) by translating it into a solvable framework. It is highly general and can be applied to almost all physical quantities. In the following sections, we will employ this scheme to demonstrate several off-diagonal measurements that previously were rather difficult, even impossible to be calculated in QMC. Moreover, we can change the annealing path to measure the off-diagonal observables on the path of physical parameters and the degree of freedom in both space and time.

Equal-time off-diagonal correlations
As an example, we consider the Hamiltonian of the spin-1/2 XXZ model, which is given by:

H_XXZ = ∑_{i,j} (1/2(S+i S-j + S-i S+j) + ΔSz i Sz j) (4)

where 〈i, j〉 denotes the nearest neighbors, Δ is the parameter that controls the anisotropy. The Hamiltonian can be simulated using the directed loop algorithm of the SSE method (detail can be found in Supplementary Note 1)9,71,85,86. In this model, two-point off-diagonal operators hSxSxi=hSySyi=ðhS+S(cid:2)i+hS(cid:2)S+iÞ=4, which can be measured using worm-like sampling trick10,68,70. However, measuring a general off-diagonal correlation function is significantly more challenging. Fig. 2 shows (detail can be found in Supplementary Note 2). Here, taking correlation of Sx operators as an example, we show how to measure it via varying the physical parameter Δ in our scheme,

hSx i Sx j i Δ = trðSx i Sx j e(cid:2)βHÞ / trðe(cid:2)βHÞ = Z(Δ) / Z(Δ) (5)

where Z(Δ) represents a general partition function with extra off-diagonal operators inserted, distinguished from a normal partition function without these extra off-diagonal operators. The calculation of hSySyi is the same as hSxSxi in this frame, as explained in the Supplement note 1.

Firstly, we consider an obvious reference point of this model: Δ0=1, known as the “Heisenberg point" and possessing SU(2) symmetry. This implies that the Sz i Sz j correlation is accessible through the diagonal measurement, which is the same as hSx i Sx j i. In a standard QMC framework, the focus shifts to measuring the ratio of the partition functions. For convenience, we define that Zr=Z(Δ)=Z(Δ0) and Zr=Z(Δ)=Z(Δ0). Then the Eq. (5) can be rewritten as

hSx i Sx j i Δ = (Zr / Zr) × hSz i Sz j i Δ0=1 (6)

In this way, the correlation of Sx operators can be easily calculated as Fig. 2 shows (detail can be found in Supplementary Note 2). The QMC results are also compared with the exact diagonalization (ED) in order to demonstrate the reliability of this scheme. Figure 2 shows the calculation results from ED and BRA. The subfigures (a) and (b) exhibit two-point correlations and four-point correlations, represented by C(r) = hSx i Sx i+r i and Cijkl = hSx i Sx j Sx k Sx l i in the one-dimensional chain. The results of 2D lattice with L=4, Ly=2, β=8 are shown in the subfigure (c). The black line represents the ED results which match well with the QMC data. We have plotted only a few points on the graph for clarity, while the actual simulation data points of BRA are densely distributed.

One may feel that the SU(2) symmetry at Δ0=1 is a strict condition which is not general for an arbitrary model. Actually, it is convenient to introduce an auxiliary Hamiltonian H0 with friendly symmetry or easily solvable property. As what quantum annealing does87–89, we can set the BRA path as tH+(1−t)H0 and anneal from t=0 to t=1. This approach allows us to obtain the observable of the target Hamiltonian H.

System size and distance growth by annealing
Another choice for the reference point is to measure the observable through the ED method in small size, then anneal the small system to large size.

In this approach, the system size L and distance r between Sx 1 and Sx 1+r can be considered as BRA parameters. For instance, we can choose the hSxSx i for L=4 as a reference point, and then we obtain the hSxSx i for larger system L via adding the remaining sites L−4 to the original chain, as shown in the Fig. 3(a). In this procedure, the interaction J2 is tuned to couple 4 sites with L−4 sites. When we fix the system size and choose the hSx 1 Sx 2 i as a reference point, we can obtain the hSx 1 Sx 1+r i (r>1) via adding some sites to the area between Sx and Sx and removing some sites at the end of this chain. We need to tune the coupling J2 from 0+ to 1 and also gradually adjust the coupling J3 from 1 to 0+, in order to keeping the chain length unchanged as displayed in the Fig. 3(b).

A consistent match between the QMC and ED results for off-diagonal correlations is found in Supplementary Note 5. This method certainly can be extended to the simulation of large systems. As shown in Fig. 3, we obtain the off-diagonal correlation Cxx(r)=hSxSx i using the reweighting method of changing system size (benchmark with the worm-trick method can be found in Supplementary Note 8). In Fig. 3(b), |Cxx(r)| has power-law decay, which reflects the physical feature of Luttinger liquid. As Δ decreases, the power-law parameter becomes smaller, which indicates the SxSx correlation is enhanced and SzSz correlation is weakened. Besides, we utilize the method of annealing Δ, which was introduced in the preceding section, with fixed large-size L=48 to obtain the curves of |Cxx(r)| (dashed line in Fig. 3(b)), which agrees well with the results through annealing L and r.

Separability of measurements
The second scheme, which involves annealing from a small system to larger system, inspires us to explore the separability of the general measurement in a large system. Without loss of generality, we consider a scenario where a large system is composed of two decoupled smaller subsystems as shown in Fig. 4. This approach can be easily extended to systems with multiple parts. In the decoupled case, the density matrix of the total system is the tensor product of the two density matrices, i.e. ρ = ρA ⊗ ρB. Typically, we encounter two kinds of measured operators, OA ⊗ OB and OA + OB, they satisfy

hOA ⊗ OB i A∪B = tr((ρA ⊗ ρB)(OA ⊗ OB)) / tr(ρA ⊗ ρB) = (tr(ρA OA) / tr(ρA)) × (tr(ρB OB) / tr(ρB)) = hOA iA hOB iB (7)

and

hOA + OB i A∪B = tr((ρA ⊗ ρB)(OA ⊗ I + I ⊗ OB)) / tr(ρA ⊗ ρB) = hOA iA + hOB iB (8)

where 〈…〉 A∪B denotes the observable is measured in the total system A∪B and the coupling between A and B is zero. 〈…〉 A(B) denotes the measurement in the subsystem A(B).

Based on the above two equations, we can firstly decompose a large system into several independent parts without coupling and measure the observables of each part via ED. By taking the ED result as a reference point, we then employ QMC to reweight the coupling between each parts from zero to the target value. Consequently, the final observable in the total system can be obtained in this way.

For example, we assign SxSx operator to subsystem A and another SxSx operator to subsystem B. The expectation value hSxSxi and hSxSx i can be obtained via ED since the system size of A or B is small. Subsequently, we incrementally adjust the coupling JAB between A and B to obtain the correlation hSxSxSxSx i. As depicted in Fig. 4, we utilize the above annealing method to obtain the four-point off-diagonal correlation with different system size, and the reference points are obtained with small system size L0=L=2 via ED. The QMC results are in excellent agreement with the pure ED results, which demonstrates the reliability of this method. From a technical standpoint, it is related to the annealing of L and r of the chain in the previous section, but this section more clearly demonstrates the power of this method from small to large sizes. In the next section, we will use this approach to calculate disorder operators in 2D systems.

Disorder operator in 2D TFIM
Here we investigate the off-diagonal measurement for the transverse Ising model (TFIM). The Hamiltonian is given as follows,

H_TFIM = −J ∑_{i,j} σz i σz j − h ∑ i σx i (9)

where σz/x is the Pauli spin-1/2 matrix and 〈i, j〉 means the nearest-neighbor coupling. h > 0 is transverse field term and J > 0 is the ferromagnetic term6. Because the TFIM only preserves Z2 symmetry, we choose the J=0+ and h=1 as a reference point. When J=0, the reference point hσxσxi=1 since all the σx=1. In the simulation, we can choose J → 0+ which makes hσxσxi very close to 1. The BRA formula can be expressed as Z(J)=Zr=Zr×hσxσxi, where Zr=Z(J)=Z(J0=0+) and Zr=Z(J)=Z(J0=0+). If we want to measure the many-body off-diagonal observables, we just need to change the Z(J)=hσxσxi into hσxσx...σxi. For TFIM, the QMC results in small system sizes are also well consistent with the ED90,91, detail can be found in Supplementary Note 6.

We then mainly focus on the disorder operator of 2D TFIM on a square lattice. The disorder operator is a non-local operator which can reveal the high-form symmetry breaking and conformal field theory (CFT) information in quantum many-body systems92–102. For 2D TFIM, we define the disorder operator hXi=h ∏ i∈M σx i to detect the non-local information, where M is a R×R square area in the lattice. Its perimeter is l = 4R and it contains R2 off-diagonal operators. This disorder operator, a multi-body off-diagonal observable, was only well measured in the QMC based on σx basis in the past, which is challenging to obtain directly in the σz basis92. Although the operator σx is contained in the TFIM Hamiltonian and can be measured in the σz basis in principle27,74, it suffers from rather large fluctuations due to the requirement of a product of a series of σx in an area. It requires that the series of σx operators must appear connectedly in the time-space manifold in the σz basis, which is a low-probability event.

This difficulty can be overcome via BRA method. As depicted in Fig. 5, we have successfully obtained the disorder operator with different perimeters l in the paramagnetic (PM) phase and ferromagnetic (FM) phase (benchmark with the directed calculation in σx basis can be found in Supplementary Note 9). Here we set h=1 and the critical point becomes J=0.328573,92. For convenience, we firstly utilize the separability method in the above section to measure the disorder operator at J/h=0.18. Taking it as a reference point, we then obtain the disorder operators for different J/h via annealing along J. In the PM phase, the disorder operator satisfies the perimeter law 〈X〉~e−al, which is consistent with the CFT prediction. In the FM phase, the disorder operator satisfies the area law hXi~e−bl2, which reveals the presence of high-form symmetry92.

Imaginary-time off-diagonal correlations
Our goal becomes to extend our method to imaginary time correlation functions involving off-diagonal operators. Our discussions will concentrate on the framework of path-integral-like QMC. The first way based on the physical parameter reweighting is straightforward, which is similar to the method we have employed in the above sections. By fixing two operators at distinct points in imaginary time τ, we have observed the evolution of the imaginary-time correlation function hSx(τ)Sx(0)i with varying parameter Δ, as depicted in Fig. 6. This is achieved by evaluating the correlation function at several distinct imaginary-time points: τ=0.1, τ=1.0, τ=3, and τ=5. Notably, when β=10, τ=β/2=5 corresponds to the maximum separation in imaginary time. The simulated values, directly comparable as Zr=Zr, demonstrates excellent agreement with the ED results as shown in the subfigure (a) of Fig. 6. For a larger size L=32 with β=64, as Δ is tuned from 1 to 0, the imaginary-time off-diagonal correlation hSx(τ)Sx(0)i gradually becomes larger, which is analogous to the equal-time cases.

Nontrivially, we perform the BRA measurement along the imaginary-time axis, where the distance between two inserted operators increases linearly during the annealing process, as illustrated in Fig. 7(a). For example, we focus on the measurement of 〈O(τ)O(0)〉 (the operators are indicated by gradient of colors) currently inserted at time zero and τ, the corresponding partition function for this configuration is Z(τ). We aim to derive Z(τ0) for the off-diagonal operator at τ0 using the reweighting technique. Different from the above schemes for reweighting in which the old/new weight uses a same configuration, the measured operators O(τ) and O(τ0) represent different configurations here. The solution is to construct an extended ensemble Z(τ)∪Z(τ0), where Z(τ) and Z(τ0) are the measured ensembles containing operators O(τ) and O(τ0), as shown in Fig. 7(a). In this frame, the ratio Z(τ0)=Z(τ) can be estimated by the ratio of sampling numbers N(τ0)=N(τ), where the number N(τ0) or N(τ) denotes how many times the sampling belongs to the ensemble Z(τ) or Z(τ0). The similar spirit has been used to calculate the entanglement entropy in QMC60. More details about this scheme are explained in the Supplementary Note 3.

We present the numerical outcomes for the XXZ chain in Fig. 7(b). Our analysis has focused on the behavior of the weight ratio Zr, across three distinct coupling strengths: Δ=0.1, Δ=0.5, and the Heisenberg condition Δ=1.0. We observe that Zr initiates from the same starting point for all three curves, with the Heisenberg coupling exhibiting a more rapid decay which reflects the energy gaps in related cases. As Δ decreases, the imaginary-time off-diagonal correlation decays more slowly, indicating the Sx imaginary-time correlation is enhanced that is similar to the equal-time case. Moreover, the inset illustrates that, at Δ=0.1, larger β makes the ratio Zr closer to zero via reducing the finite size effect in imaginary-time direction.

Furthermore, we can obtain the spectrum of operators from the momentum imaginary-time correlations via stochastic analytical continuation (SAC)38,103,104. The momentum imaginary-time correlation is defined as G αα(q,τ)= 1/L ∑ i,j e−iq(ri−rj) hsα i(τ)sα j(0)i (α=x,y,z). All the real-space off-diagonal imaginary-time correlation can be captured by the above imaginary-time BRA method, which is used to simulate the excitation spectrum Sαα(q, ω). As shown in Fig. 8, the off-diagonal spectrum has sharper lower boundary with weak continuum, which is different from the diagonal spectrum that has strong spinon continuum on the upper boundary (ED results can be found in Supplementary Note 4). Since Δ=0.1 here is close to zero, the difference of the diagonal and off-diagonal spectra can be understood qualitatively from the limit Δ=0. When Δ=0, the off-diagonal excitation can be solved by the Jordan-Wigner transformation, which is related to a single-mode dispersion of free fermion, thus its excitation is sharp. Meanwhile, the diagonal spectrum Szz corresponds to two fermion operators, which contributes a continuum therefore. The results demonstrate that our BRA method can be successfully applied to extract the off-diagonal spectrum, which also reveals the different excitation modes compared to the diagonal spectrum with the anisotropic phase.

Methods
Our work proposes a general measurement framework based on the SSE-QMC method, the principle of which has been introduced in the main text. The XXZ model and TFIM model involved employ directed-Loop algorithm and cluster SSE algorithm (detail can be found in Supplementary Note 1 and Note 6), respectively, which will be found in detail in the supplementary information.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This work presents a solution to the problem of measuring general offdiagonal operators in quantum Monte Carlo\nsimulations using a bipartite reweighting-annealing technique. Previously inaccessible or very challenging measurements\nfor the XXZ and transverse-field Ising models are evaluated by connecting them via a reweighting chain that passes through\nan easily solvable point. The authors present a multitude of ways for building such chains, including tuning coupling\nparameters, system size and imaginary time distances.\nThe paper is well written, and the results and flexibility of the method are impressive, warranting publication in Nat.\nCommun. after the authors have addressed the points below.\n1. The introduction of the offdiagonal measurement problem references a lot of stochastic series expansion works. Are there\nsimilar issues in determinant QMC (I am not sure - here we can get a lot using just Wick’s theorem) or variational Monte\nCarlo/neural quantum states? Referring to a larger variety of methods here may strengthen the manuscript.\n2. In the second column of page 1: “Furthermore, it is usually impossible to design updates between {W_i} and {W_i'}” Why is\nit usually impossible? Either provide additional reasoning or a reference.\n3. Fig 1. is quite hard to understand. I only understood it after reading and understanding the main text. The meaning of J' in\nthe histograms is not clear and the thick arrow left of “Easily Solvable” suggests something gets easier to solve from right to\nleft.\nSuggestions for improvement:\n- Show less of the section right of the easily solvable point\n- Mark the path of the reweight-annealing chain from Z(J0) to Zbar(J0)\n- Instead of J', maybe J+ΔJ?\nThese are all only suggestions.\n4. In Eq. 2, it would be helpful to specify next to the brackets which distribution the expectation value belongs to.\n5. Page 2, right column, concerning the polynomial complexity at fixed parameter spacing: When scaling the system size,\nthe overlap between two distributions can be interpreted as an extensive free-energy difference. Would it be fair to say that in\nthe general case, at fixed spacing, I need to scale the number of distributions with the system size? If that is the case, it\nshould be mentioned.\nAlso from standard parallel tempering, it is known that the presence of a critical point on the chain can decrease the overlap\nsignificantly. I expect it is the same here, and it would be helpful to mention it so that future users of the method choose\nchains that avoid critical points if possible.\n6. Page 3, beginning of Sec. III. It is mentioned that ⟨S⁺i S¯j⟩ can be measured using worm sampling. It should be clarified\nthat ⟨S^x_i S^x_j⟩ = ⟨S^y_i S^y_j⟩ = (⟨S⁺i S¯j⟩ + ⟨S¯i S⁺j⟩)/2 for the XXZ model, which is easy to get via existing updates.\nIn general, non-vanishing two-site correlation functions should be measurable using the worm trick.\nThe higher-order plaquette correlations, however, are very hard to get with worms, and a clear improvement of the state of\nthe art.\n7. In Fig. 2, it should be possible to understand the meaning of (I)-(IV) without going to the appendix. For the chain, one can\nwrite C_1234, C_1245 etc. For the 2x4, either use a similar notation or move the sketch from the appendix to the figure.\n8. Page 5, right column: “For TFIM, the QMC results in small system sizes are also well consistent with the ED” I saw\ncalculations for this in the appendix, but they are not referenced here. Add a reference to the appendix.", "reviewer_2": "In their manuscript “Addressing general measurements in quantum Monte Carlo”, the authors address the issue of how to\nconstruct an estimator for any general observable in quantum Monte Carlo methods. They start from the standard idea to\nview the expectation value of an observable as the ratio of two different types of partition functions $\\bar{Z}$ and $Z$. These\nare estimated separately within a reweight-annealing scheme. The key point is then to include a reference point where the\nratio of the two functions can easily be calculated via a different method. This allows to properly normalise the observable.\nThe authors claim that their scheme solves the long-standing problem of calculating the overlap between different\ndistribution functions in mathematical statistics.\nI find the manuscript well-written and interesting. The algorithm and results are interesting and merit publication. The idea of\nusing an reweight-annealing scheme is not new. It has for example been used to calculate the entropy, which is typically an\nobservable that is not easily accessible numerically. As far as I know, the attempt to construct a general measurement\nscheme out of this is more original. However, I am very sceptical about their claim that the scheme will work generically to\nestimate any observable in any regime.\nFirst of all, one needs a reference point. This reference point is an easy case (a specific value of the coupling where the\nmodel has a certain symmetry, or a small enough system). That point needs to be connected to the point of interest via a\nsequence of parameters. The fluctuations of the ratio of partition functions will typically be very difficult to control, unless the\nparameters are chosen close enough. In practice, if they are not close enough, one might end up having systematic errors.\nMoreover, these errors will propagate as one goes through the sequence. The authors do not discuss this potential problem.\nHow do the authors exclude any systematic errors of this nature? What criterion is used for choosing the number of\nparameters in the sequence?\nA related issue: going to more difficult (and more interesting) regimes will require taking more points along the path. When\ngoing from weak coupling to strong coupling, for example, one might cross a phase transition. In principle, any phase\ntransition will be a smooth transition since one is considering a finite-size system. But it is clear that one might run into\nproblems when the system size is large enough, due to accumulation and propagation of error.\nThe current manuscript gives interesting examples, but it does not discuss propagation of error. I believe that an imporant\nquestion should be addressed: what ensures that the final error bar is controlled and small? It seems very difficult to me to\nensure this in general, as the authors implicitly claim when they write that “this scheme solves the long-standing problem of\ncalculating the overlap between different distribution functions”.\nI cannot recommend publication in Nature Communications of the manuscript as it is. In my opinion, the authors should\neither include an in-depth analysis of the error. Or remove the claim of solving a long-standing problem and publish in a\nmore specialised journal."}, "author_response": "RESPONSE TO REFEREES\nReferee #1\nThis work presents a solution to the problem of measuring general offdiagonal operators in quantum Monte Carlo\nsimulations using a bipartite reweighting-annealing technique. Previously inaccessible or very challenging measure-\nments for the XXZ and transverse-field Ising models are evaluated by connecting them via a reweighting chain that\npasses through an easily solvable point. The authors present a multitude of ways for building such chains, including\ntuning coupling parameters, system size and imaginary time distances.\nThe paper is well written, and the results and flexibility of the method are impressive, warranting publication in\nNat. Commun. after the authors have addressed the points below.\nWe really appreciate your careful reading of our work and grasping its main points, and we are also very grateful\nfor the kind and detailed comments you provided. We will address them one by one in the following response.\nQ1.1 The introduction of the off-diagonal measurement problem references a lot of stochastic series expansion works.\nAretheresimilarissuesindeterminantQMC(Iamnotsure-herewecangetalotusingjustWick’stheorem)or\nvariational Monte Carlo/neural quantum states? Referring to a larger variety of methods here may strengthen\nthe manuscript.\nResponse: We thank the referees for insightful comments and suggestions.\nAs demonstrated in Ref. [1–4], off-diagonal Green’s functions such as c c† are naturally encoded in the de-\n⟨ i j⟩\nterminant QMC (DQMC) weight and can be directly measured, while higher-order correlation functions follow\nstraightforwardly from Wick’s theorem. In this sense, determinant QMC does not encounter the same intrinsic\ndifficulties with off-diagonal correlation functions that we explicitly discuss in the present work. However, this\nmeasurement scheme of DQMC, being parasitic on the algorithm itself, does not extend to other QMC methods\nor to the evaluation of off-diagonal observables.\nVariational Monte Carlo (VMC) or Diffusion Monte Carlo (DMC) can indeed estimate off-diagonal observables\ndirectly[5,6]. Fordynamicaloff-diagonalquantities,onemayworkwithinavariationalexcited-statesubspace[7]\nor use linear-response formulations to obtain Green’s and spectral information [8]. The method is not restricted\nto fermionic systems, it inherently depends on the specific choice of trial wavefunction. The accuracy of the\nresults therefore directly relies on the quality of the ansatz. We have noted in the manuscript, “Unlike other\nnumerical methods, QMC cannot directly obtain the wave-function of ground state,” in the first column of page\n1, meaning that VMC or DMC, with its reliance on a trial wavefunction, lies outside the scope of our discussion.\nWethanktherefereeforpointingoutthisdistinction,andintherevisedmanuscriptwehaveemphasizedtherole\nofDQMCintheintroduction(Page1,column2–Page2,column1),\"determinantQMC(DQMC)canefficiently\ncalculate two-body Green’s functions, where higher-order correlators can be derived using Wick’s theorem. [1–4]\nThisworksparticularlywellinfermionicsystems,whereoff-diagonalGreen’sfunctionssuchas c c† arenaturally\n⟨ i j⟩\nembedded in the sampling weight. We now turn to diagrammatic Monte Carlo, world-line QMC and stochastic\nseries expansion (SSE) formulations, in which off-diagonal measurements pose a much more severe challenge.\nDespite numerous efforts in the past, two-body Green’s functions can be obtained only in certain cases within\nthe framework of worm-like QMC algorithm [9–18].\"\nQ1.2 In the second column of page 1: “Furthermore, it is usually impossible to design updates between W and W′”\ni i\nWhy is it usually impossible? Either provide additional reasoning or a reference.\nResponse: We sincerely thank the referee for insightful comments and suggestions. If the Hamiltonian does not\n2\ncontaintheoff-diagonaloperators,thegeneralizedpartitionfunctionZ¯ =tr Oe−βH (Oisoff-diagonaloperator)\n{ }\ndoes not naturally appear when sampling the original partition function Z = tr e−βH . Since O changes the\n{ }\nquantum numbers of the configuration (for instance, O = σx flips spins), its contribution has vanishing\ni∈M i\nprobability in the Z-ensemble. This implies that the manifold of Z¯ is statistically disconnected from that of Z,\n(cid:81)\nmakingitextremelydifficulttodesignlocalupdatesthatstochasticallyconnectW andW′. Thereisnoreference\ni i\nsince it hasn’t been realized.\nFor example, if the measured operator O = σx where M is the whole region of the system, the update\ni∈M i\nbetween Z¯ and Z needs to change all the σx into identity operators I which is rather difficult for designing a\ni(cid:81) i\nloop/cluster update in the spacetime manifold.\nQ1.3 Fig 1. is quite hard to understand. I only understood it after reading and understanding the main text. The\nmeaning of J′ in the histograms is not clear and the thick arrow left of “Easily Solvable” suggests something gets\neasier to solve from right to left. Suggestions for improvement:\n• Show less of the section right of the easily solvable point\n• Mark the path of the reweight-annealing chain from Z(J ) to Z¯(J )\n0 0\n• Instead of J′, maybe J +∆J?\nThese are all only suggestions.\nResponse: We appreciate the referee’s feedback on the clarity of Fig. 1. To improve readability, we have revised\nthe figure and its caption. Specifically, we have made the following changes:\n• We have reduced the space on the right-hand side of the “Easily Solvable” label to avoid misinterpretation\nof the arrow direction. Our intention is to emphasize that both Z¯(J) and Z(J) converge through the RA\nprocess at the solvable reference point, rather than suggesting a directional change in difficulty.\n• We have increased the separation of the two distributions at the initial parameter point J and explicitly\n0\nmarked the “Lack of Overlap” to highlight the origin of the difficulty in direct sampling.\n• We have added contour lines to the histograms and modified subplots (b) and (c) so that the distribution\nat J (i.e., the distribution before reweighting) is represented as an unfilled outline curve. This makes the\nfigure cleaner and more intuitive.\n• We have added a description of the RA process, annealing direction and the J′ parameter.\nThese changes emphasize the core idea of our work and render the schematic much clearer. We thank the referee\nagain for this valuable suggestion.\nOriginal images:\nQ1.4 In Eq. 2, it would be helpful to specify next to the brackets which distribution the expectation value belongs to.\nResponse: We thank the referee for the suggestion. We have revised Eq. (2) to explicitly indicate that the\nexpectation value W(J′)/W(J) is taken with respect to the distribution at parameter J. This clarification\nJ\n⟨ ⟩\nmakesthenotationconsistentwithstandardpracticeinreweightingmethods. Wealsorevisedtheaccompanying\ntext to emphasize that the average is performed over the ensemble sampled at J, \"Eq.(2) states that the ratio\nZ(J′)/Z(J)isobtainedbyaveragingW(J′)/W(J)overconfigurationssampledfromtheensembleatJ (i.e. ).\"\nJ\n⟨·⟩\nQ1.5 Page2,rightcolumn,concerningthepolynomialcomplexityatfixedparameterspacing: Whenscalingthesystem\nsize, the overlap between two distributions can be interpreted as an extensive free-energy difference. Would it be\nfair to say that in the general case, at fixed spacing, I need to scale the number of distributions with the system\nsize? If that is the case, it should be mentioned. Also from standard parallel tempering, it is known that the\npresence of a critical point on the chain can decrease the overlap significantly. I expect it is the same here, and\nit would be helpful to mention it so that future users of the method choose chains that avoid critical points if\npossible.\nResponse: Weappreciatethisimportantpoint. Ourmethodrequiresthatthedistributionsbesufficientlyadjacent\nso that the reweighting procedure remains stable. As already stated in our Introduction, “Once a reference point\nZ(J) is known, Z(J′) can be calculated through the ratio. It has been proved that the computation complexity\n4\nof the RA method is polynomial if the ratio of two closest Z(J) and Z(J′) is fixed in the division strategies.” In\npractice, this implies that as the system size grows, the number of intermediate parameters along an annealing\nprocessmustincreasetomaintainsufficientoverlap. Thisstrategyofautomaticallyadjustingthenumberofslices\nto ensure controlled accuracy has been discussed in detail by Ref. [19]. Following their guideline, it is preferable\nto choose the partition-function ratios between neighbors to remain of order (1) (empirically, within about\nO\n0.1–10).\nRegarding critical regions, we agree that the overlap can drop markedly near a phase transition. Since we have\nno more innovative idea than increasing the number of samples, inserting more slices, or bypassing the critical\npoint,acombinationofthesethreestrategiesmaybethemostpracticalchoice. Besides,inthepseudo-automatic\nannealing process of Ref. [19], the required number of non-identity operators provides a way to design the\nparameter sequence such that the ratio between successive partition functions remains approximately constant.\nInthissense,onecanroughlyidentifythelocationofcriticalregionsinthispreprocessing(e.g.,wherethenumber\nof non-identity operators changes discontinuously).\nIn summary, we thank the referee for raising this important point. We have added a more detailed discussion in\nAppendix G–I.\nQ1.6 Page 3, beginning of Sec. III. It is mentioned that S+S− can be measured using worm sampling. It should be\n⟨ i j ⟩\nclarified that SxSx = SySy = ( S+S− + S−S+ )/2 for the XXZ model, which is easy to get via existing\n⟨ i j⟩ ⟨ i j⟩ ⟨ i j ⟩ ⟨ i j ⟩\nupdates.\nResponse: Weappreciatethereferee’scommentandhaverevisedthemanuscripttoaddressthisconcern. Wehave\nrevised the original sentence as \"In this model, two-point off-diagonal operators SxSx = SySy =( S+S− +\n⟨ i j⟩ ⟨ i j⟩ ⟨ i j ⟩\nS−S+ )/2, which is easy to be measured using worm-like sampling trick.\"\n⟨ i j ⟩\nQ1.7 In Fig. 2, it should be possible to understand the meaning of (I)-(IV) without going to the appendix. For the\nchain, one can write C ,C etc. For the 2 4, either use a similar notation or move the sketch from the\n1234 1245\n×\nappendix to the figure.\nResponse: We thank the referee kindly for useful suggestions. In Fig. 2, for the 1D chain, we have revised the\nnotations: the correlations at from r =1 to r =3 are now denoted as C , C , C , and the panels (I)–(III) are\n12 13 14\nnow labeled as C ,C and C . For 4 2 lattice, we have revised (I)-(IV) as C ,C ,C and C . In\n1234 1245 1357 13 15 12 1278\n×\naddition, upon re-examining the data, we found that the error bars for the SSE data points in subfigure (b) were\nplotted incorrectly. This has been corrected, with no visible change in the overall results. We have also updated\nthe figures in the Appendix (Fig. S3 and Fig. S4) accordingly to reflect the same labeling. Thank you again\nfor pointing this out.\n5\n2.0\n1.5\n1.0\n0 1\n∆\nrZ/r¯Z\n2.0\n3 3\n(a) ED (b) ED (c) ED\nr=1 (I) (I)\nr=2 (II) (II) 1.5 r=3 (III) (III) 2 2\n(IV)\n1.0\n1 1\n0 1\n0 1 0 1\n∆\n∆ ∆\nC4(L/2)= S 1 xS 2 xSx L Sx L+1\n2 2\n(cid:10) (cid:11)\n0.0225 8 8\n×\n20 20\n×\n0.0200\n0.0175\n0.0150\n(d)\n0.0125\n0.00 0.25 0.50 0.75 1.00\n∆\nrZ/r¯Z\n3 3\n(a) ED (b) ED (c) ED\nC12 C1234 C13\nC13 C1245 C15\nC14 2 C1357 2 C C 1 1 2 278\n1 1\n0 1 0 1\n∆ ∆\nC4(L/2)= S\n1\nxS\n2\nxS\n1\nx\n+L\nS\n2\nx\n+L 2 2\n(cid:10) (cid:11)\n0.0225 8 8\n×\n20 20\n×\n0.0200\n0.0175\n0.0150\n(d)\n0.0125\n0.00 0.25 0.50 0.75 1.00\n∆\nFIG. 1: Original (left) and revised (right) versions of Fig. 2\nFIG. 2: Original (left) and revised (right) versions of Fig. S3 in Appendix\nFIG. 3: Original (left) and revised (right) versions of Fig. S4 in Appendix\ncalculations for this in the appendix, but they are not referenced here. Add a reference to the appendix.\nResponse: We sincerely thanks the referee for the useful suggestion. We have added the reference to this sentence\nin the main text. Thus, we revise this sentence as \"For TFIM, the QMC results in small system sizes are also well\nconsistent with the ED [87,88], detail can be found in Appendix F [89].\"\nIfindthemanuscriptwell-writtenandinteresting. Thealgorithmandresultsareinterestingandmeritpublication.\nThe idea of using an reweight-annealing scheme is not new. It has for example been used to calculate the entropy,\nwhich is typically an observable that is not easily accessible numerically. As far as I know, the attempt to construct\na general measurement scheme out of this is more original. However, I am very sceptical about their claim that the\nscheme will work generically to estimate any observable in any regime.\nFirst of all, one needs a reference point. This reference point is an easy case (a specific value of the coupling where\nthe model has a certain symmetry, or a small enough system). That point needs to be connected to the point of interestviaasequenceofparameters. Thefluctuationsoftheratioofpartitionfunctionswilltypicallybeverydifficult\nto control, unless the parameters are chosen close enough. In practice, if they are not close enough, one might end up\nhaving systematic errors. Moreover, these errors will propagate as one goes through the sequence. The authors do\nnot discuss this potential problem. How do the authors exclude any systematic errors of this nature? What criterion\nis used for choosing the number of parameters in the sequence?\nArelatedissue: goingtomoredifficult(andmoreinteresting)regimeswillrequiretakingmorepointsalongthepath.\nWhen going from weak coupling to strong coupling, for example, one might cross a phase transition. In principle,\nany phase transition will be a smooth transition since one is considering a finite-size system. But it is clear that one\nmight run into problems when the system size is large enough, due to accumulation and propagation of error.\nThe current manuscript gives interesting examples, but it does not discuss propagation of error. I believe that an\nimportant question should be addressed: what ensures that the final error bar is controlled and small? It seems very\ndifficult to me to ensure this in general, as the authors implicitly claim when they write that “this scheme solves the\nlong-standing problem of calculating the overlap between different distribution functions”.\nWe thank the respected referee for careful reading and helpful suggestion on our manuscript. In the reweight-\nannealing framework, the ratio between two partition function can be as an expectation under the ensemble at J :\n2 = 2 (1)\nZ(J ) W(J )\n1 (cid:28) 1 (cid:29)J1\nthe reweight-annealing (RA) estimator is unbiased, when the number of Monte Carlo samples N is sufficiently large.\nAs pointed out and analyzed by Ref. [19], the accuracy can, in principle, be improved to an arbitrary level simply\nby increasing N. In the QMC simulation, It is essential to maintain the ratio Z(J′)/Z(J) with O(1) range (i.e.\n0.1 10), which ensures the polynomial computation complexity [19] and controllable system error. Besides, the\n∼\nnumber of slices should be increased near the phase transition, since the partition function evolves rapidly in this\nregime.\nHowever, in practice, finite Monte Carlo steps generate the systematic error. Specifically, the weights lying\nbelow the 1/N threshold (illustrated by the red dashed line in Fig. 4) cannot be sampled reliably under N samplings.\nIn the annealing scheme, we assume that two adjacent parameters J and J are sufficiently close so that each step\n1 2\nmaintains equilibrium. Accordingly, the distributions can be approximated as overlapping Gaussians: the blue curve\nrepresents the ensemble before reweighting (J ), while the red curve corresponds to the ensemble after reweighting\n1\n(J ).\nC\ngoal is to estimate the distribution under the target parameter J while remaining within the J ensemble. However,\n2 1\nthe far tail of the blue distribution cannot be statistically inaccessible under sampling. We thus use the blue dashed\nline as an effective boundary. This divides the contribution from J into two parts: the configurations-overlapping\n2\nregion W and the beyond-tail region W . With a finite number of samples N, the part W lying below the 1/N\n1 3 2\nthresholdcanoccasionallybesampled,thoughonlywithpooraccuracy. Thetruesystematicbiasalsooriginatesfrom\n8\nthe entirely unsampled region W . Taken together, W +W represent the source of systematic error.\n3 2 3\nconfiguration\nthgiew\nZ(J )\n1\nN\nW2\nFIG. 4: The distribution of partition functions Z(J ) and Z(J ) with different parameter J and J . N is the Monte Carlo\n1 2 1 2\nstep. The green region including W and W denotes the weights that can be sampled for Z(J ) in the reweight-annealing\n1 2 2\nprocess. The region including W and W denotes the weights that can be sampled in the original directed sampling.\n1 3\nThen,tocontrolthiserrorinpractice,weadjustthestepsizebetweenJ andJ sothattheratio Z(J2) = W(J2)\n1 2 Z(J1) ⟨W(J1)⟩ J1\nremains of (1), typically within 0.1 to 10. This ensures that the inaccessible part W is kept at an (1) (1/N)\nscale, similarly the systematic contribution of W also stays at (1)/N. Therefore, the total systematic error can be\nexpressed as (1)/N 1/N. Meanwhile, statistical fluctuations introduce a statistical error of order 1/√N in QMC\nO ∼\nsimulations. For sufficiently large N, the systematic error 1/N is much smaller than the statistical error 1/√N. This\nexplains why the systematic error can be safely neglected compared to the statistical error.\nWhen sampling over M slices (each with an expectation value x), the systematic error, according to error propaga-\ntion,isgivenby√MxM−1/N scalingwith 1/N,whichremainssmallerthanthetotalstatisticalerror√MxM−1/√N\n∼\nscaling with 1/√N. Hence, the systematic error does not accumulate to a significant or uncontrollable level as M\n∼\nincreases, since it is always dominated by the statistical error. These discussions have been added to the Appendix G\nin the manuscript.\nWe have calculated two examples to demonstrate our analysis about the error. For example, in the XXZ chain\nwithvaryinganisotropicparameter∆,weemployedtheworm-trickmethod[14]tomeasuretheoff-diagonaltwo-point\ncorrelationfunction. Theresultsareconsistentwiththoseobtainedfromthereweight-annealingmethod, asshownin\nFig. 5.\nIn the case of the 2D transverse-field Ising model, measuring multi-point off-diagonal operators such as the Z\n2\ndisorder operator using the worm-trick method in the σz basis is extremely challenging. However, by switching to\ntheσx basisandsimulatingtheHamiltonianviadirectedloopalgorithms, theZ disorderoperator X = σx\n2 ⟨ ⟩ ⟨ i∈M i⟩\nbecomesdiagonalandstraightforwardtomeasure. FurtheralgorithmicdetailscanbefoundinRef.[20]. Asillustrated\n(cid:81)\nin Fig. 6, the disorder operator values obtained from the reweighting-annealing method agree well with those from\ndirectmeasurementsintheσx basiswithinerrorbars. Thesefindingsconfirmthatthesystematicerrordoesnotgrow\nto an uncontrollable level even as the number of slices and the system size increase. These numerical results have\nbeen added to the Appendix H and I in the manuscript.\nWhen going from weak to strong coupling during the reweighting annealing process, the system may cross a phase\ntransition. According to fixing the ratio of two closest partition functions Z(J2) (1), using more slices to sample\nZ(J1) ∼ O\nthe partition function ratio is beneficial in such regions, since the distribution of the partition function changes more\nrapidly near the phase transition. Based on the above error analysis, the systematic error remains smaller than the\n9\nstatistical error and does not accumulate significantly even when crossing the phase transition. Therefore, the final\nerror bar remains controllable and does not accumulate to an unacceptable level during this process.\n100\n10 1\n−\n10 2\n−\n100 101\nr\n)r(xxC\n|\n|\n∆=0.8\n∆=0.6\n∆=0.4\n∆=0.8\n∆=0.6\n∆=0.4\nFIG. 5: Two point off-diagonal correlations for system size L=48 obtained from lattice annealing method (annealing the\nsystem size L and distance r between off-diagonal operators), The dashed lines represent the off-diagonal operators\ncorrelation from worm-trick method.\n1.4\n1.2\n1.0\n0.8\n0.6\n0.4\n0.2\n0.0\n6 14 22 30 38\nl\n)\nX\n(nl\ni\nh\nJ/h=0.18\nJ/h=0.18\nJ/h=0.21\nJ/h=0.21\nJ/h=0.24\nJ/h=0.24\nJ/h=0.27\nJ/h=0.27\nJ/h=0.30\nJ/h=0.30\nFIG. 6: Z disorder operator ⟨X⟩=⟨ (cid:81) σx⟩ obtained from the reweight annealing method (circle) and directed\n2 i∈M i\nmeasurement in the σ basis (dashed line) in the 2D TFIM model.\nx\n[1] F. Assaad and H. Evertz, World-line and determinantal quantum monte carlo methods for spins, phonons and electrons,\nin Computational Many-Particle Physics, edited by H. Fehske, R. Schneider, and A. Weiße (Springer Berlin Heidelberg,\nBerlin, Heidelberg, 2008) pp. 277–356.\n[2] C.-C. Chang, S. Gogolenko, J. Perez, Z. Bai, and R. T. Scalettar, Recent advances in determinant quantum monte carlo,\nPhilosophical Magazine 95, 1260 (2015), https://doi.org/10.1080/14786435.2013.845314.\n10\n[3] Z.-X. Li and H. Yao, Sign-problem-free fermionic quantum monte carlo: Developments and applications, Annual Review\nof Condensed Matter Physics 10, 337 (2019).\n[4] F. Sun and X. Y. Xu, Delay update in determinant quantum monte carlo, Phys. Rev. B 109, 235140 (2024).\n[5] R. J. Needs, M. D. Towler, N. D. Drummond, and P. López Ríos, Continuum variational and diffusion quantum monte\ncarlo calculations, Journal of Physics: Condensed Matter 22, 023201 (2009).\n[6] J. Toulouse, R. Assaraf, and C. J. Umrigar, Chapter fifteen - introduction to the variational and diffusion monte carlo\nmethods, in Electron Correlation in Molecules – ab initio Beyond Gaussian Quantum Chemistry, Advances in Quantum\nChemistry, Vol. 73, edited by P. E. Hoggan and T. Ozdogan (Academic Press, 2016) pp. 285–314.\n[7] M.CharleboisandM.Imada,Single-particlespectralfunctionformulatedandcalculatedbyvariationalmontecarlomethod\nwith application to d-wave superconducting state, Phys. Rev. X 10, 041023 (2020).\n[8] B. Mussard, E. Coccia, R. Assaraf, M. Otten, C. J. Umrigar, and J. Toulouse, Chapter fourteen - time-dependent linear-\nresponse variational monte carlo, in Novel Electronic Structure Theory: General Innovations and Strongly Correlated\nSystems, Advances in Quantum Chemistry, Vol. 76, edited by P. E. Hoggan (Academic Press, 2018) pp. 255–270.\n[9] N.Prokof’evandB.Svistunov,Wormalgorithmforproblemsofquantumandclassicalstatistics,UnderstandingQuantum\nPhase Transitions, Lincoln D. Carr, ed. Taylor & Francis, Boca Raton 910 (2010).\n[10] M. Boninsegni, N. V. Prokof’ev, and B. V. Svistunov, Worm algorithm and diagrammatic monte carlo: A new approach\nto continuous-space path integral monte carlo simulations, Phys. Rev. E 74, 036701 (2006).\n[11] P. Gunacker, M. Wallerberger, E. Gull, A. Hausoel, G. Sangiovanni, and K. Held, Continuous-time quantum monte carlo\nusing worm sampling, Phys. Rev. B 92, 155102 (2015).\n[12] A. Dorneich and M. Troyer, Accessing the dynamics of large many-particle systems using the stochastic series expansion,\nPhys. Rev. E 64, 066701 (2001).\n[13] C.Zhou,Z.Yan,H.-Q.Wu,K.Sun,O.A.Starykh,andZ.Y.Meng,Amplitudemodeinquantummagnetsviadimensional\ncrossover, Physical Review Letters 126, 227201 (2021).\n[14] W. Zhu and W. Guo, Measuring off-diagonal correlation function in stochastic series expansion quantum monte carlo\nsimulation, Journal of Beijing Normal University (Natural Science) 57, 593 (2021).\n[15] O. F. Syljuåsen and A. W. Sandvik, Quantum monte carlo with directed loops, Phys. Rev. E 66, 046701 (2002).\n[16] O. F. Syljuåsen, Directed loop updates for quantum lattice models, Phys. Rev. E 67, 046701 (2003).\n[17] P.Henelius,P.Fröbrich,P.J.Kuntz,C.Timm,andP.J.Jensen,Quantummontecarlosimulationofthinmagneticfilms,\nPhys. Rev. B 66, 094407 (2002).\n[18] C.-J.Huang,L.Liu,Y.Jiang,andY.Deng,Worm-algorithm-typesimulationofthequantumtransverse-fieldisingmodel,\nPhys. Rev. B 102, 094101 (2020).\n[19] Y.-M. Ding, J.-S. Sun, N. Ma, G. Pan, C. Cheng, and Z. Yan, Reweight-annealing method for evaluating the partition\nfunction via quantum monte carlo calculations, Phys. Rev. B 110, 165152 (2024).\n[20] J. Zhao, Z. Yan, M. Cheng, and Z. Y. Meng, Higher-form symmetry breaking at ising transitions, Phys. Rev. Res. 3,\n033024 (2021)."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have successfully addressed all my remarks. I recommend publication.", "reviewer_2": "After reviewing the new manuscript and the reply of the authors, I believe this version is suitable for publication in Nature\nCommunications. The addition of appendices G, H, and I has significantly strengthened the work and its claims."}, "author_response": ""}]

Decision letter