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^3N_t)$, where $N_x$ ($N_t$) is the number of spatial (temporal) sites.

Figures (4)

• Figure 1: The scaling behavior for both decompositions ("QR") and naïve implementation ("no decomp.") of the evaluation of $\log\det M$. The left panel shows the expected $\mathcal{O}(N_t)$ scaling at fixed system size $N_x$$(=882)$, while the right panel shows the expected asymptotic $\mathcal{O}(N_x^3)$ scaling at fixed $N_t$$(=32)$. The grey dashed lines are to guide the eye. Our decompositions are roughly a factor $\sim 5$ slower compared to the naïve implementation.
• Figure 2: Convergence of leapfrog integrator used in HMC simulations of the Perylene molecule (see e.g. Rodekamp:2024ixu) with onsite coupling $U=2$ and inverse temperature $\beta=90$ (corresponding to room temperature). The change in energy $\Delta H$ of a single HMC trajectory is shown as a function of the inverse number of molecular dynamics $N_{\text{md}}$ steps. Calculations were done with two different numbers of time steps, $N_t=512,1024$. The auxiliary field was initially sampled from a Gaussian distribution, $\phi_{x,t}\sim\mathcal{N}(0,\sqrt{\delta U})$. The expected convergence is $\Delta H\propto N_{\text{md}}^{-2}$.
• Figure 3: The error incurred during integration of equations of motion during a single HMC trajectory, using our decomposition ("QR") and naïve matrix multiplications ("no decomp."). For each data point the trajectory length was the same and the number of leapfrog integration steps was constant with $N_{\text{md}}=50$. The calculation was repeated 400 times with different random seeds and the median was calculated. Error bars represent the median absolute deviation. $N_t$ was chosen such that $\beta/N_t=1/4$. The auxiliary field was initially sampled from a Gaussian distribution, $\phi_{x,t}\sim\mathcal{N}(0,\sqrt{\delta U})$. The expected scaling in the error is $\propto \beta$, and this is shown as the dashed blue line with an arbitrary coefficient.
• Figure 4: Two-body correlator $C(\tau)$ for a spin-singlet $(S=0)$ bright exciton, corresponding to a particle-hole excitation, in a $(10,2)$ chiral carbon nanotube at inverse temperature $\beta=30$. The channel shown has total pseudo-spin $I=0$. The data are obtained from HMC simulations using our decomposition method, which enables stable measurements of the excitonic correlator at large inverse temperature.