Efficient continuous-time quantum Monte Carlo algorithm for fermionic lattice models

Efficient continuous-time quantum Monte Carlo (CT-QMC) algorithms that do not suffer from time discretization errors have become the state of the art for most discrete quantum models. They have not been widely used yet for fermionic quantum lattice models, such as the Hubbard model, nor other fermionic lattice systems due to a suboptimal scaling of $O\left({\beta }^3\right)$ with inverse temperature $\beta$, compared to the linear scaling of discrete-time algorithms. Here we present a CT-QMC algorithm for fermionic lattice systems that matches the scaling of discrete-time methods but is more efficient and free of time discretization errors. This provides an efficient simulation scheme that is free from the systematic errors opening an avenue to more precise studies of large systems at low and zero temperature.

Figure 1

Sketch of the insertion and removal updates. The five-vertex configuration $c$ can be modified adding a vertex at site $x$ and time $\tau$ with spin $\rho$. This leads to a proposed configuration $c^{\prime }$. The reverse move consists of removing the vertex $(x,\tau ,\rho )$. These two basic moves (insertion and removal) are sufficient to reach any term in the series (3) from any other.

Figure 2

Kinetic energy and interaction energy for a $4\times 4$ Hubbard plaquette at half filling with $U=4t$, computed using the discrete-time BSS algorithms for various values of ${\mathrm{\Delta }}_{\tau }t$ and LCT-AUX. The inset compares the CT-QMC results to BSS results extrapolated to ${\mathrm{\Delta }}_{\tau }\to 0$.

Figure 3

Computational time of the simulations from Fig. 2 as a function of $\beta$. Both the present method and discrete-time algorithm have been run for ${10}^6$ sweeps for thermalization, then ${10}^6$ sweeps while measuring observables. Both codes have been similarly optimized with delayed rank-1 updates. The linear dependence is clear in both pictures. The slope of the CT scheme is the same as the discrete time with ${\mathrm{\Delta }}_{\tau }t=1/8$ because the average number of vertices in a time slice of size $t$ is around 120, i.e., roughly the same as $8V=128$. The CT code extrapolates to a much larger constant for $\beta \to 0$ mostly due to several allocations per update which have not been optimized.