Determinant Monte Carlo algorithms for dynamical quantities in fermionic systems

We introduce and compare three different Monte Carlo determinantal algorithms that allow one to compute dynamical quantities, such as the self-energy, of fermionic systems in their thermodynamic limit. We show that the most efficient approach expresses the sum of a factorial number of one-particle-irreducible diagrams as a recursive sum of determinants with exponential complexity. By comparing results for the two-dimensional Hubbard model with those obtained from state-of-the-art diagrammatic Monte Carlo, we show that we can reach higher perturbation orders and greater accuracy for the same computational effort.

Figure 1

Imaginary part of the Hubbard atom self-energy at order 8 in $U$ as obtained from Dyson's equation (green), the equations of motion approach (orange), and the direct self-energy measurement (blue). We use $\beta =10, U=1, \epsilon =-0.2$. All simulations lasted 120 CPU hours.

Figure 2

Hubbard model self-energy at the first Matsubara frequency ${\tilde{\mathrm{\Sigma }}}_{\vec{k}}\left(i{\omega }_0\right)$ along the $\vec{k}=(0,0)\to (\pi ,0)\to (\pi ,\pi )\to (0,0)$ path at order 3 in $U$, as obtained from Dyson's equation (green), the equations of motion approach (orange), and the direct self-energy measurement (blue). We use a $32\times 32$ lattice with $\beta t=2, U=4t, \mu =0$ and a uniform $\alpha$ shift ${\alpha }_{\uparrow }={\alpha }_{\downarrow }=1.53t$. All simulations lasted 120 CPU hours.

Figure 3

Imaginary part of the Hubbard atom self-energy at order 12 in $U$ as obtained from the equations of motion approach (orange) and the direct self-energy measurement (blue). We use $\beta =10, U=1, \epsilon =-0.2$. All simulations lasted 120 CPU hours.

Figure 4

Variance of the imaginary part of the Hubbard atom self-energy at the first Matsubara frequency. Orange lines with stars is the result of the equations of motion. Blue line with dots corresponds to the direct self-energy measurement. We use $\beta =10, U=1, \epsilon =-0.2$.

Figure 5

Hubbard model self-energy at order 6 in $U$ on a $32\times 32$ lattice with $\beta t=2, U=4t, \mu =0$ and with a uniform $\alpha$ shift ${\alpha }_{\uparrow }={\alpha }_{\downarrow }=1.53t$. Blue symbols are results for the direct self-energy measurement; orange symbols are results from the equations of motion approach. Upper panel: Self-energy at the first Matsubara frequency ${\tilde{\mathrm{\Sigma }}}_{\vec{k}}\left(i{\omega }_0\right)$ along the $\vec{k}=(0,0)\to (\pi ,0)\to (\pi ,\pi )\to (0,0)$ path. Lower panel: Self-energy as a function of $i{\omega }_n$ at $\vec{k}=(\pi ,\pi /2)$. All simulations lasted 120 CPU hours.

Figure 6

Variance of the imaginary part of the Hubbard model self-energy Im ${\tilde{\mathrm{\Sigma }}}_{(\pi ,\pi /2)}\left(i{\omega }_0\right)$. Orange lines with stars is the result equations of motion. Blue line with dots corresponds to the direct self-energy measurement. We use a $32\times 32$ lattice with $\beta t=2, U=4t, \mu =0$ and with a uniform $\alpha$ shift ${\alpha }_{\uparrow }={\alpha }_{\downarrow }=1.53t$.

Figure 7

Hubbard model self-energy at the first Matsubara frequency ${\tilde{\mathrm{\Sigma }}}_{\vec{k}}\left(i{\omega }_0\right)$ along the $\vec{k}=(0,0)\to (\pi ,0)\to (\pi ,\pi )\to (0,0)$ path at order 7 in $U$, as obtained from DiagMC (red) and the direct self-energy measurement (blue). We use a $32\times 32$ lattice with $\beta t=2, U=4t, \mu =0$ and a uniform $\alpha$ shift ${\alpha }_{\uparrow }={\alpha }_{\downarrow }=1.53t$. Simulations lasted 1440 CPU hours for the $\mathrm{\Sigma }\mathrm{Det}$ and 4000 CPU hours for the DiagMC.

Figure 8

Imaginary part of the local lattice self-energy ${\mathrm{\Sigma }}_{\text{loc}}^{\sigma }\left(i{\omega }_n\right)$ as a function of Matsubara frequency, as computed using $k$ orders, with $k=2,...,9$. The red squares are results obtained from DQMC (error bars are smaller than the symbol size). Inset: Zoom on the first Matsubara frequency. It is seen that the results are converged with an error bar smaller than $1\%$. We use a $32\times 32$ lattice with $\beta t=2,U=4t,\mu =0$ and with a uniform $\alpha$ shift ${\alpha }_{\uparrow }={\alpha }_{\downarrow }=1.53t$. The discrete time interval in DQMC is $\mathrm{\Delta }\tau =1/32$.

Figure 9

Benchmark of the contribution to the Matsubara frequency self-energy $\tilde{\mathrm{\Sigma }}\left(i{\omega }_n\right)$ for the Hubbard atom at order 5 in the perturbation series in $U$. Red squares are the analytical solution. Green lines are obtained from a calculation of the Green's function with Eq. (2). Orange line is the result of the equations of motion and lies on top of the blue curve corresponding to the direct self-energy measurement. We use $\beta =10, U=1, \epsilon =-0.2$. All simulations lasted 1200 CPU hours.

Figure 10

Hubbard model self-energy at the first Matsubara frequency ${\tilde{\mathrm{\Sigma }}}_{\vec{k}}\left(i{\omega }_0\right)$ along the $\vec{k}=(0,0)\to (\pi ,0)\to (\pi ,\pi )\to (0,0)$ path at order 4 in $U$, as obtained from DiagMC (red), Dyson's equation (green), the equations of motion approach (orange), and the direct self-energy measurement (blue). We use a $32\times 32$ lattice with $\beta t=2, U=4t, \mu =0$ and a uniform $\alpha$ shift ${\alpha }_{\uparrow }={\alpha }_{\downarrow }=1.53t$. The DiagMC simulation lasted 400 CPU hours, while all other simulations lasted 1440 CPU hours.

Figure 11

Comparison of the time for one Monte Carlo cycle (in microseconds) between the direct accumulation of the self-energy (blue curve with dots) and the computation of the Green's function using CDet (green curve with dots), on a semilog scale, as a function of the perturbation order $n$. Each curve is fitted by its expected high-$n$ behavior: ${\gamma }_{\mathrm{\Sigma }}{n}^2{3}^n$ for the $\mathrm{\Sigma }\mathrm{Det}$ (dotted red line) and ${\gamma }_G{3}^n$ for Dyson (dashed red line), where ${\gamma }_G=0.0464$ and ${\gamma }_{\mathrm{\Sigma }}=0.0012$ are implementation-dependent constants. Inset: Ratio of the time of one MC cycle for $\mathrm{\Sigma }\mathrm{Det}$ ($t_{\mathrm{\Sigma }}$) and for the CDet ($t_G$), as a function of the perturbation order $n$.