Integral algorithm of exponential observables for interacting fermions in quantum Monte Carlo simulations

Exponential observables, formulated as $ln\langle e^{\widehat{X}}\rangle$ where $\widehat{X}$ is an extensive quantity, play a critical role in the study of quantum many-body systems, examples of which include the free energy and entanglement entropy. Given that $e^X$ becomes exponentially large (or small) in the thermodynamic limit, the accurate computation of the expectation value of this exponential quantity presents a significant challenge. In this paper, we propose a comprehensive algorithm to quantify these observables in interacting fermion systems, utilizing the determinant quantum Monte Carlo method. We have applied this algorithm to the two-dimensional square-lattice half-filled Hubbard model and $\pi$-flux t-V model. In the Hubbard model case at the strong-coupling limit, our method showcases a significant accuracy improvement on free energy compared to conventional methods that are derived from the internal energy, and in the t-V model, we indicate that the free energy offers a precise determination of the second-order phase transition. We also illustrate that this approach delivers highly efficient and precise measurements of the $n\mathrm{th}$ Rényi entanglement entropy. Even more noteworthy is that this improvement comes without incurring increases in computational complexity. This algorithm effectively suppresses exponential fluctuations and can be easily generalized to other models.

Figure 1

(a) Internal energy density $E$, free-energy density $F_2$ derived by solving $E\left(\beta \right)=-\frac{\partial [ln(Z\left)\right]}{\partial \beta }=\frac{\partial \left(\beta F_2\right)}{\partial \beta }$, and free-energy density $F_1$ according to Eq. (6). (b),(d),(f) $\langle ln\left(P_s\right)\rangle$, $\langle ln\left[\mathrm{Tr}\left({\rho }_{A;s_1}{\rho }_{A;s_2}\right)\right]\rangle$ and $\langle ln\left[\mathrm{Tr}\left({\rho }_{A;s_1}{\rho }_{A;s_2}{\rho }_{A;s_3}\right)\right]\rangle$ for different $\beta$ points used in (a),(c),(e) from small to large, corresponding to circles from bottom to top. Black lines show the fitting where the small slopes indicate a small fluctuation for $\langle ln\left(P_s\right)\rangle$ at each $t$ point. (c) The second EE $S_A^{\left(2\right)}$ converging with $\beta$ to 0.89, which is consistent with DMRG at zero temperature. (e) The third EE $S_A^{\left(2\right)}$ converging with $\beta$ to 0.7, which is consistent with DMRG at zero temperature.

Figure 2

(a) Schematic phase diagram of a spinless fermion $\pi$-flux half-filled square-lattice t-V model. The orange point indicates a $N$ = 2 chiral-Ising transition in (2++1)D separating Dirac semimetal (SM) and charge density wave (CDW) phases. The blue line indicates a 2D Ising transition and our simulation goes along the dashed red line at $V=4$ where a finite-temperature phase transition occurs. (b) Directly measured internal energy density $E$ and free-energy density $F$ measured by integral algorithm vs temperature $T$ for system size $L\times 2L$ with $L=6$ ($L=8,10$ have a similar internal energy and free-energy density). (c) First-order (${\partial }_{T}F$) and second-order (${\partial }_T^{2}F$) derivative of free energy for system size $L=6$, where the sharp dip at the green arrow at $T_c\approx 2.1$ in ${\partial }_T^{2}F$ indicates the second-order phase transition temperature. (d) Finite-temperature phase transition shown by the green arrow at $T_c\approx 2.1$ determined by 2D Ising finite-size scaling $[M_2\equiv {\sum }_{i,j}\frac{{\alpha }_i{\alpha }_j}{L^4}\langle (n_i-\frac{1}{2})(n_j-\frac{1}{2})\rangle$, where ${\alpha }_i=\pm 1$ for $i\in A,B$ sublattice and $\eta =\frac{1}{4}$ for 2D Ising university class].