Stable Determinant Monte Carlo Simulations at Large Inverse Temperature $β$

Abstract: At low temperatures $T$ where $1/T=β\gg1$ the naïve implementation of determinant quantum Monte Carlo (DQMC) methods suffers from loss of precision and numerical instabilities when evaluating the fermion determinant. This instability propagates into the calculation of observables that rely on the evaluation of the inverse of the fermion matrix, or the Greens function. For DQMC methods that rely on the Hamiltonian Monte Carlo (HMC) algorithm, an additional complication comes from evaluating the force terms required for integrating Hamilton's equations of motion, since here loss of precision and numerical instabilities are also prevalent. We show how to address all these issues using various choices of matrix decompositions, allowing us to simulate at $β\gtrsim 90$, which corresponds to room temperature for graphene structures. Furthermore, our implementation has numerical costs that scale similarly to the naïve implementation, namely as $\mathcal{O}(N_x<sup>3N_t)$, where $N_x$ ($N_t$) is the number of spatial (temporal) sites.

• The paper achieves a numerically stable Monte Carlo simulation by using recursive matrix decompositions to prevent precision loss at low temperatures.
• The method maintains core operations—determinant, force, and Green's function calculations—with stable UDT factorizations and without extra asymptotic cost.
• The stabilized framework enables practical low-temperature simulations for quantum many-body systems, facilitating studies of materials like graphene and organic molecules.

Stable Determinant Monte Carlo Simulations at Large Inverse Temperature $\beta$

Introduction and Motivation

Simulation of quantum many-body systems at low temperatures using determinant quantum Monte Carlo (DQMC) approaches is fundamentally limited by numerical instabilities that arise due to the ill-conditioning of the underlying fermionic matrices. This paper provides a robust framework to systematically stabilize DQMC and Hybrid/Hamiltonian Monte Carlo (HMC) calculations at large inverse temperatures ($\beta \gg 1$), which are particularly relevant for modeling systems near and below room temperature, such as graphene and organic molecules. The main technical achievement is the introduction of a maximally stable algorithm for all core operations—determinant, force, and Green's function calculations—without incurring additional asymptotic computational cost.

Numerical Instabilities in Standard DQMC

DQMC methods transform the fermionic path integral into a stochastic sampling problem involving auxiliary fields, where all fermion degrees are integrated analytically, resulting in determinants of large dense matrices. At high $\beta$, the products in these determinants span exponentially separated scales, so naive floating-point arithmetic rapidly loses accuracy, resulting in loss of physical observables, unreliable Metropolis sampling, and ultimately invalid simulation results.

The exponential ill-conditioning can be traced back to expressions involving products such as $\exp(\beta \kappa)$, where the condition number grows as $e^{2\beta K}$ (with $K$ the bandwidth of the tight-binding Hamiltonian). Standard approaches fail for $K\beta \gtrsim 35$, making accurate large-$\beta$ simulations intractable without stabilization.

Stabilized Matrix Decomposition Framework

The authors adapt and generalize recursive matrix decomposition schemes, using QR and SVD factorizations, to separate out scale components at all multiplication steps, preventing catastrophic rounding errors. Each time-slice evolution matrix $M_t$ is decomposed and all subsequent products are maintained in a $UDT$-factorized form (unitary, diagonal, and triangular), ensuring the product remains numerically tractable regardless of $\beta \gg 1$0.

The stabilization is maintained for all core operations:

• $\beta \gg 1$1 calculation: By recursively decomposing each $\beta \gg 1$2 and maintaining the diagonal dominance before the final addition with the identity, the scheme retains stability up to $\beta \gg 1$3 with no extra asymptotic scaling over the naive $\beta \gg 1$4 approach.
• Force calculation (HMC): The authors demonstrate that the standard recursive force evaluation mixes scale orderings and rapidly loses precision (even with SVD/QR), as proven by explicit analysis in minimal examples. Instead, "prefix" and "suffix" products ($\beta \gg 1$5 and $\beta \gg 1$6) and careful cyclic permutations yield stable force terms for the HMC molecular dynamics integration.

  Figure 2: Convergence of the leapfrog integrator for HMC in Perylene at $\beta \gg 1$7 demonstrates stable energy scaling with the number of molecular dynamics steps.

• Green's function computation: The methodology extends efficiently to the full fermionic propagator, both diagonal and off-diagonal entries, exploiting the previously computed decomposition products. All time-displaced Green's function components relevant for observable calculations are obtained with numerical stability.

  Figure 1: Integration error in HMC trajectories as a function of $\beta \gg 1$8 for the stabilized vs. naive schemes, highlighting the breakdown of the latter at moderate $\beta \gg 1$9.

Algorithmic Structure and Scaling

All core components, including determinant, force, and Green's function, are evaluated with no increase in leading order cost: $\beta$0 for determinant and diagonal Green's function, $\beta$1 for the full Green's function, with only moderate prefactors (about 5 to 10 times slower than naive multiplication). The essential step is that all matrix multiplications, inversions, and additions maintain the $\beta$2-decomposed structure, which allows correct handling of scale separation, even in the presence of strong quantum fluctuations.

Parallelization is naturally compatible with tree-based scan algorithms; thus, the implementation is efficient on both CPUs and GPUs, maximizing throughput for large-scale simulations.

Strong Results and Numerical Validation

The stabilized approach enables DQMC/HMC computations at $\beta$3, whereas prior approaches break down at $\beta$4. This parameter regime corresponds to room temperature for graphene and related materials, opening access to relevant experimental and technological conditions that were previously inaccessible.

Figure 3: Two-body singlet excitonic correlator in a chiral carbon nanotube at $\beta$5 calculated using stabilized HMC, demonstrating the measurement stability for correlation functions at large $\beta$6.

The authors provide systematic validation: the leapfrog integrator in HMC displays the theoretically expected $\beta$7 scaling; the error in the total energy remains controlled and small even at large $\beta$8—contrasting sharply with the unstable divergence seen with naive matrix multiplication.

Theoretical and Practical Implications

From a theoretical standpoint, the work settles the stabilization problem for a broad class of lattice fermion models, covering spin, charge, and particle-hole bases, with applications ranging from condensed matter (Hubbard, Pariser-Parr-Pople models) to quantum chemistry. The explicit demonstration of when and why the naive recursive force evaluation fails provides new conceptual clarity, enabling further algorithm development for related stochastic simulation frameworks.

On the practical side, the algorithm enables low-temperature simulations critical for interpreting experiment and designing quantum materials. Immediate applications include:

• Simulating Dirac material ground states, high-temperature superconductors, and organic molecular semiconductors at physical temperatures.
• Stable measurement of time-displaced correlators, including the "disconnected" parts that were previously difficult to treat, thereby improving access to dynamical and response observables.

Moreover, all routines are made available through the open-source NSL library, facilitating adoption by the broader computational materials and lattice field theory communities.

Outlook and Future Directions

The work provides a foundation for the next generation of ab initio simulations in strongly correlated electron systems. Combining this stabilization framework with other recent algorithmic advances—Fourier acceleration, radial updates for HMC, advanced analysis of Euclidean correlators—positions the field for realistic simulations of large, complex systems (including away from half-filling and in the presence of electron-phonon or disorder effects). The authors indicate that simulations of increasingly complex molecules (e.g., Perylene, Corannulene) at experimental temperature scales have already become feasible, enabling direct comparison with experiment and new materials predictions.

Conclusion

The algorithm presented achieves numerically stable DQMC and HMC simulations at large $\beta$9, eliminating the longstanding barrier of loss of precision due to scale separation in matrix products. Achieving this without additional leading order computational cost makes physically relevant, low-temperature simulations of correlated electron systems a practical reality. The methodology, publicly implemented, is expected to become the standard for future simulations of two-dimensional materials, organic molecules, and strongly interacting quantum matter.