2026-09-11
Notes on Imaginary-Time Correlations in SSE-QMC: Basic Ideas and Measurement Techniques
This blog post follows Flynn2026ImaginaryTime, a recent systematic guide on measuring imaginary-time correlation functions in Stochastic Series Expansion Quantum Monte Carlo (SSE-QMC).
A quick heads-up: this post is not meant for complete beginners–sorry about that! We won’t be covering basic SSE setup here, and I’ll assume you’re already familiar with the foundational equations. Deriving everything from scratch would take up way too much space and distract from our main focus.
If you have any questions along the way, feel free to contact me. Now enjoy!
We begin with the standard SSE partition function expansion:
As a quick recap, SSE-QMC relies on a direct Taylor expansion of the Boltzmann factor along the imaginary-time direction to simulate quantum partition functions. Because of this, it avoids Trotter decomposition entirely, which is a defining feature and key advantage that sets SSE-QMC apart from other QMCs, such as Path-Integral QMC.
Here, indexes the operator position (or “layer”) along the SSE operator string. Equation gives the exact partition function after carrying out the Taylor expansion, inserting identity operators, and applying fixed-length truncation.
How do we actually extract imaginary-time correlation functions
An imaginary-time correlation function needs to know the exact . The most straightforward, maybe rough, idea is to treat the string index as a proxy for real imaginary time .
Mathematically, this linear relationship is exact in the large cutoff limit where the maximum string cutoff approaches infinity (). We know that, setting is obviously impossible. As long as we choose a sufficiently large cutoff such that (where is the actual number of non-identity operators in the string), in the averaging sense, we obtain .
A more rigorous approach, is respecting the binomial distribution between and .
Once we collect the discrete correlations across all string separations , the exact correlation function is reconstructed by this weighted summation. However, there is a clear practical downside to this approach. Because the binomial sum mixes all discrete string distances , we are forced to measure and accumulate the discrete correlation function across the entire discrete operator space before we can extract the correlation at a specific .
Time slice
The basic idea: first, obtain exact imaginary time, so we slice the time .
Cut into slices, ,
Thus, for each small “Taylor expansion”, every slice boundary has a well-defined imaginary time. The total partition function is,
Well now, the configuration weight is indexed independently by these slices, where each slice contains its own independent number of non-identity operators , and small imaginary-time interval .
Notice: you can not swap and , no matter how you piece or slice, the total PF:
always keeps the sum over configurations on the very outside, which means,
Modified diagonal update
The most direct change brought by time slices is the diagonal update probability. Since this is not the main focus of this section, we quickly summarize it by analogy to the full- case: replacing all variables with their slice-specific counterparts. The diagonal update probabilities are controlled by , becoming:
Next, all strategies revolve around these slice boundaries, as we want to leverage them to obtain exact imaginary-time quantities.
Easy diagonal measurement
Yes, we can measure diagonal imaginary-time correlation functions just like normal diagonal observables. We measure at the slice boundaries, keeping the relative operator distance fixed while sliding the starting point. Measuring the full-space correlation function takes only computational complexity, which is slices, sites.
Easy off-diagonal measurement
Even for those off-diagonal operators, present in the Hamiltonian, for example, in TFIM model, analogous to measuring the energy operator, we can quickly obtain their correlation functions.
then ,
where and are the number of non-identity operators in their respective slices.
well, there may be one subtle point to keep in mind: if we consider an operator sitting right on a slice boundary, which adjacent slice should it actually belong to?
The answer is, you can pick either side freely; whichever side you choose, you just check its corresponding number of non-identity operator in the slice.
Apparently, this method is restricted to slice boundaries. Consider a terrible scenario: when the transverse field is very small, operators are sparse in the configuration, and those falling directly on slice boundaries are rarer still—drastically reducing the number of valid measurement events. This naturally brings us to the core question: How can we average this over the entire slice operator string, freeing our measurements from reliance on slice boundaries?
Improved Estimator, Slice-Level Averaging
We want to increase the number of measurement events and reduce the variance (but this method will introduce systematic bias , which can be a worthwhile in some cases). Still, we want as possible as much string operator to contribute the statistics.
For example, within slice slice, the measurement (imagined insertion) can take place at any of the available operator insertion positions. Summing over all possible positions introduces a normalization factor of . When absorbing the configuration weight of the -th slice, , this extra factor accounts for the virtual insertion of an additional operator, yielding:
Then, the imaginary-time correlation becomes from
to the final form:
Where is the number of operator about in slice k, site . 1
1 Standard energy measurements rely on the expansion:
The feeling of this is to first imagine fixing an insertion position:
Thus, the measured value is actually embodied in a single Hamiltonian operator. Since the expansion only handles the number of non-identity operators as , here must be a specific bond in the Hamiltonian, which can be either a bond or a single ; Now, let us look at the result of averaging over imaginary time:
This also displays the two branching “feelings” from the main text:
- Without averaging: The final statistic is .
- With averaging: The final statistic is , where is the number of occurrences of , which forces a direct dependence on the specific nature of . For example, if is , then simply counts the number of times appears.
However, the method introduces bias.
Consider the actual coordinates of the two measured operators, we can model them as a two-particle () system. Each particle’s equilibrium position be at the midpoint of its respective slice, . Constructing the center-of-mass coordinate and relative coordinate, the fluctuation of their relative distance has a zero mean and non-zero variance .
The bias comes from . If we Taylor expand the imaginary-time correlation around exact with small :
Substituting the relation between and the time-slice length , we obtain,
So, a question regarding this approach is: how do we balance statistical error against systematic bias?
Bias–Variance Tradeoff
(I will omit this part for now.)
Tricky off-diagonal measurement
When off-diagonal operators are not present in the Hamiltonian, things usually get a bit tricky, as it means we genuinely need extra techniques to access them. While countles methods have been developed for this purpose, here we focus specifically on measuring Green’s functions on the fly during the directed loop updates.
Green’s functions:
Here we briefly point out that the mechanism of directed loops can itself be viewed as carrying a pair of charge defects at the head and tail. The tail is generally treated as a fixed starting point , while the head is seen as a defect continuously propagating along the world line. After undergoing various vertex-scattering events, when it eventually returns to the starting point, the defects annihilate and the loop closes.
Thus, as you might imagine, one can randomly choose a starting point , which simultaneously determines the defect type at the tail. As the update progresses, by tracking the position and defect type of the head while continuously accumulating frequencies as it moves, one eventually builds up a distribution akin to a histogram.
Candidate Tails
If we want exact time positions, both the head and tail must simultaneously reside on a certain slice boundary. However, we do not wish to alter the rules of Directed Loop updates, which originate from a randomly chosen leg.
The idea is to shift the position of the tail before it undergoes any scattering. If it can be effectively moved onto a slice boundary, then every location where the head subsequently crosses a boundary during its propagation will contribute a valid count to the Green’s function histogram.
There are several key points worth keeping in mind here:
- The tail can only move vertically along the imaginary-time direction until it encounters another vertex. This constraint arises because we must preserve the charge defect at the tail and avoid any scattering that would alter its charge—such as a transition from to . Consequently, shifting the tail along continuous, unscattered world lines provides the ideal trajectory.
- If multiple candidate positions are available for moving the tail, we can only select one of them to serve as the imagined starting point.
The probability of selecting a starting point depends on the current configuration ; thus, its inverse weight must be applied as a reweighting factor during measurements.
The total probability of selecting a specific “leg + boundary ” pair as the starting point is given by:
Under importance sampling, the update estimator becomes:
If this accumulated value is deemed too large (scaling as ), it can be scaled by , which serves purely as a configuration-independent normalization factor. Ultimately, the resulting histogram is unnormalized and is scaled physical measurements using the boundary condition .
Considering that moving the tail might present multiple candidate boundary positions, I also refer to this approach as the “candidate tail” method.
Several implementation details, such as how to explicitly compensate for the configuration-dependent selection probability, are beyond the scope of this discussion. I may detail them in a dedicated blog post if I have time.