Reaching the continuum limit in finite-temperature ab initio field-theory computations in many-fermion systems

TL;DR

This paper introduces a low-rank factorization method that significantly reduces computational costs in finite-temperature field-theory calculations for many-fermion systems, enabling more effective continuum limit extrapolations.

Contribution

The authors develop a systematic low-rank factorization approach that lowers the computational scaling from cubic to linear in lattice size, improving finite-temperature many-body simulations.

Findings

Achieved exact calculations of 2D Fermi gas properties across temperature ranges.

Demonstrated the method's effectiveness in normal and superfluid states.

Enabled more feasible continuum limit extrapolations in correlated fermion systems.

Abstract

Finite-temperature, grand-canonical computations based on field theory are widely applied in areas including condensed matter physics, ultracold atomic gas systems, and lattice gauge theory. However, these calculations have computational costs scaling as $N_s^3$ with the size of the lattice or basis set, $N_s$. We report a new approach based on systematically controllable low-rank factorization which reduces the scaling of such computations to $N_s N_e^2$, where $N_e$ is the average number of fermions in the system. In any realistic calculations aiming to describe the continuum limit, $N_s/N_e$ is large and needs to be extrapolated effectively to infinity for convergence. The method thus fundamentally changes the prospect for finite-temperature many-body computations in correlated fermion systems. Its application, in combination with frameworks to control the sign or phase problem as needed, will provide a powerful tool in {\it ab initio} quantum chemistry and correlated electron materials. We demonstrate the method by computing exact properties of the two-dimensional Fermi gas with zero-range attractive interaction, as a function of temperature in both the normal and superfluid states.

Equations

$$\begin{array}{r}\hat{H}=\sum _{\mathbf{k}\sigma }{\epsilon }_{\mathbf{k}}{c}_{\mathbf{k}\sigma }^{+}{c}_{\mathbf{k}\sigma }-\mu \sum _{\mathbf{i},\sigma }{\hat{n}}_{\mathbf{i}\sigma }+U\sum _{\mathbf{i}}{\hat{n}}_{\mathbf{i}\uparrow }{\hat{n}}_{\mathbf{i}\downarrow },\end{array}$$

$$Z\simeq \int P(\mathbf{X})\prod _{\sigma =\uparrow ,\downarrow }\det (1_{N_s}+{\mathbf{B}}_M^{\sigma }{\mathbf{B}}_{M-1}^{\sigma }\cdots {\mathbf{B}}_2^{\sigma }{\mathbf{B}}_1^{\sigma })d\mathbf{X},$$

Full text

Computations are playing an increasingly important role in addressing the fundamental challenges of understanding strong correlations in interacting quantum systems. Understandably, a major part of the development and application efforts in both physics and chemistry have focused on ground-state properties. However, experimental conditions are always at finite temperatures, where often rich and new properties can be revealed Damascelli et al. (2003); Fischer et al. (2007); Mueller (2017); Brown et al. (2019). One example is the rapid development in the area of experiments with ultracold atoms, where temperature plays a crucial role and very precise measurements of properties are often possible with exquisite control over interaction strengths, environments, etc Bloch et al. (2008); Tarruell and Sanchez-Palencia (2018). Accurate computations of thermodynamic properties allow direct comparisons with experiments, but are challenging because of the presence of strong coupling and thermal fluctuations. A second example is in strongly correlated materials, including high-temperature superconductors Bednorz and Müller (1986); Saleem and Tajammul Hussain (2013), where some of the outstanding and most interesting physics questions concern finite-temperature properties Lee et al. (2006).

A common finite-temperature formalism is based on field theory in which the thermodynamic properties are computed as path-integrals in field space. Approximations can be used to perform the path integrals, including the simplest, namely mean-field calculations. A more powerful approach is to evaluate the many-dimensional integration by Monte Carlo methods Blankenbecler et al. (1981). This has become a key technique for many-body finite-temperature computations, widely applied in several fields of physics Scalapino and Sugar (1981); Hirsch (1983, 1985); ALHASSID (2001); Chandrasekharan (2006, 2010); Huffman and Chandrasekharan (2017); ALHASSID ; Alhassid (2016); Alhassid et al. (2018); Jensen et al. (2019); Anderson and Drut (2015). For example, many of the sign-problem-free computations in lattice models in condensed matter, and in Fermi gas and optical lattices of ultracold atoms, have been performed this way. For general Hamiltonians, for instance the doped Hubbard model or realistic electronic Hamiltonians in solids or molecules, a constraint can be applied to control the sign or phase problem. This framework, even with simple trial wave functions to impose an approximate constraint, has been shown to be very accurate and has been applied widely in ground-state calculations Zhang et al. (1995, 1997a, 1997b); Carlson et al. (1999); Zhang and Krakauer (2003); Suewattana et al. (2007); Chang and Zhang (2010); Qin et al. (2016); Vitali et al. (2019). Finite-temperature generalization of the approach has also been developed Zhang (1999, 2000); Liu et al. (2018); He et al. (2019).

These finite temperature calculations all have computational complexity of $\mathcal{O}(N_{s}^{3})$ LOH et al. (2005), as they are formulated in the grand-canonical ensemble to analytically evaluate the fermion trace along each path in auxiliary-field space, leading to determinants with dimension $N_{s}$. In the majority of applications, for example dilute Fermi gas and all ab initio real material simulations, it is necessary to reach the continuum (large lattice or complete basis set) limit in order to obtain realistic results. One must take $N_{s}\rightarrow\infty$, while keeping the average number of fermions $N_{e}$ fixed (at the targeted number of electrons in the molecule or cluster, or in the supercell). In contrast, all ground-state calculations, which are formulated in canonical ensemble and only need to retain occupied orbitals along the path in field space, scale as $\mathcal{O}(N_{s}N_{e}^{2})$ Shi et al. (2015). The ratio $N_{s}/N_{e}$ is often $\mathcal{O}(10$-$100)$ or larger for realistic calculations. The discrepancy of $(N_{s}/N_{e})^{2}$ in computational cost can render calculations inaccessible at low, or even modest, temperatures, when the corresponding ground-state calculations can be performed straightforwardly. This poses a fundamental obstacle for finite-temperature studies in interacting fermion systems.

We address this problem in the present paper, introducing a new method which allows finite-temperature computations to be performed with complexity $\mathcal{O}(N_{s}N_{e}^{2})$, with no loss of accuracy. We demonstrate the approach within the framework of sign-problem-free determinantal quantum Monte Carlo (DQMC) calculations, in the two-dimensional (2D) dilute Fermi gas with contact interaction, which has attracted intense experimental interest. The algorithm yields speedups of several orders of magnitude over the standard approach, and allows computations of exact thermodynamic properties in the continuum. This adds a new dimension in our computational repertoire; for example, it allows direct ab initio computations of the superfluid state in strongly interacting Fermi gases and investigations of the corresponding phase transitions.

For concreteness in describing the algorithm and to specify the computational details for our results on the 2D dilute Fermi gas, we use the attractive Hubbard model on a square lattice:

[TABLE]

where $\varepsilon_{\mathbf{k}}=2t(2-\cos k_{x}-\cos k_{y})$, $\sigma$ ($=\uparrow$ or $\downarrow$) denotes spin, and $\hat{n}_{\mathbf{i}\sigma}=c_{\mathbf{i}\sigma}^{+}c_{\mathbf{i}\sigma}$ is the density operator. The periodic supercell is represented by a lattice of size $N_{s}=L^{2}$, with corresponding momentum $k_{x}$ (and $k_{y}$) defined in units of $2\pi/L$. The energy scale of the system is set by $t$. The chemical potential $\mu/t$ is tuned to target the desired number of particles in the supercell, $N_{e}=N_{\uparrow}+N_{\downarrow}$, defining a lattice density $n=N_{e}/N_{s}$. The interaction strength $U/t$ is uniquely determined by $\log(k_{F}a)$ Shi et al. (2015), i.e., the ratio of the two-particle scattering length to the average interparticle spacing (given by the Fermi wave vector $k_{F}$). We will measure energies in units of the Fermi energy $E_{F}=2\pi nt$ and temperatures in units of $T_{F}\equiv E_{F}/k_{B}$.

The partition function $Z=\text{Tr}(e^{-\beta\hat{H}})$, after imaginary-time discretization $\beta=\Delta\tau M$, Trotter decomposition and Hubbard-Stratonovich (HS) transformation, can be written as He et al. (2019)

[TABLE]

where $\mathbf{X}=\{\mathbf{x}_{1},\mathbf{x}_{2},\cdots,\mathbf{x}_{M}\}$ is a path in auxiliary-field space, and $P(\mathbf{X})$ is a probability density function. The auxiliary fields at each time slice, $\mathbf{x}_{\ell}=\{x_{i,\ell}\}$, can be either continuous (e.g., $P(\mathbf{X})$ is a Gaussian) or discrete fields (as we adopt below for our Fermi gas calculations), in which case the integral becomes a sum. The ”$\simeq$” indicates the Trotter error, which can be extrapolated away with calculations using smaller values of $\Delta\tau$. The one-body propagator $\mathbf{B}_{\ell}^{\sigma}$ depends on the auxiliary-field $\mathbf{x}_{\ell}$, and is an $N_{s}\times N_{s}$ matrix. The paths $\mathbf{X}$ are sampled according to the integrand in Eq. (2), for example by Metropolis Monte Carlo in standard DQMC, or by a branching random walk of fixed length $M$ in constrained-path auxiliary-field quantum Monte Carlo (AFQMC) He et al. (2019).

To describe our algorithm, we introduce the notation $\mathbf{R}=\mathbf{B}^{\sigma}_{\ell}\cdots\mathbf{B}^{\sigma}_{2}\mathbf{B}^{\sigma}_{1}$ and $\mathbf{L}=\mathbf{B}^{\sigma}_{M}\mathbf{B}^{\sigma}_{M-1}\cdots\mathbf{B}^{\sigma}_{\ell+1}$, where $\ell$ is an arbitrary time slice along the path. The key components of DQMC and AFQMC calculations include: propagation; numerical stabilization; evaluations of ratios or the derivative of the integrand with respect to $\mathbf{x}_{\ell}$ during updates; and measurements. These procedures all involve operations of the matrix products $\mathbf{L}$ and $\mathbf{R}$. For example, update and measurements at $\tau=\ell\Delta\tau$ require the equal-time single-particle Green’s function matrix $\mathbf{G}^{\sigma}(\tau,\tau)\equiv\{\langle c_{i\sigma}c^{+}_{j\sigma}\rangle_{\tau}\}=(\mathbf{1}_{N_{s}}+\mathbf{R}\mathbf{L})^{-1}$. Because of the manipulation of these matrices and the computation of the determinant, the overall computational complexity scales as $\mathcal{O}(N_{s}^{3})$ LOH et al. (2005) (with a proportionality factor containing the imaginary-time length $M$). This basic structure is depicted in the top row in Fig. 1(a).

The matrices can be written in factorized forms, for example using the column-pivoted QR algorithm to perform $\mathbf{U}\mathbf{D}\mathbf{V}$ decompositions (see e.g., Ref. He et al., 2019), to obtain: $\mathbf{R}=\mathbf{U}_{R}\mathbf{D}_{R}\mathbf{V}_{R}$ and $\mathbf{L}=\mathbf{V}_{L}\mathbf{D}_{L}\mathbf{U}_{L}$, where $\mathbf{D}_{R},\mathbf{D}_{L}$ are diagonal matrices with positive elements in descending order and $\mathbf{U}_{R},\mathbf{U}_{L},\mathbf{V}_{R},\mathbf{V}_{L}$ are $N_{s}\times N_{s}$ matrices, which are used in standard methods as part of the procedure to keep numerical stability  White et al. (1989); LOH et al. (2005); He et al. (2019). The diagonal elements of $\mathbf{D}_{L}$ and $\mathbf{D}_{R}$ typically span an enormous scale (e.g., $10^{+100}\sim 10^{-100}$), controlled by the independent-particle spectrum, interaction strength, and inverse temperature. We now introduce a threshold $\epsilon$ to perform truncations on them in computing $\mathbf{R}$ and $\mathbf{L}$, since contributions from the elements of $\mathbf{D}_{R}$ and $\mathbf{D}_{L}$ smaller than $\epsilon$ will be bounded and can be made smaller than numerical noise if the truncation threshold $\epsilon$ is sufficiently small. Thus, $\epsilon$ is a parameter that we can tune to control the numerical precision of the calculation. Suppose the number of elements in $\mathbf{D}_{R}$ larger than $\epsilon$ is $m_{\text{R}}$, then $\mathbf{U}_{R}$,$\mathbf{D}_{R}$,$\mathbf{V}_{R}$ are effectively $N_{s}\times m_{\text{R}}$, $m_{\text{R}}\times m_{\text{R}}$, $m_{\text{R}}\times N_{s}$ matrices. Similarly, we have a truncation dimension $m_{\text{L}}$ for $\mathbf{D}_{L}$, and low-rank approximation for $\mathbf{L}$. The truncations are illustrated in Fig. 1. We can then decompose $\mathbf{R}\mathbf{L}=\mathbf{U}\mathbf{D}\mathbf{V}$ and an additional truncation on $\mathbf{D}$ can be performed, with truncation dimension $m_{\text{T}}$, as shown in the bottom row of Fig. 1(a). We should note that, methodologically, the above idea fits into a theme of low-rank factorizations which have found broad applications in very different contexts in physics and chemistry (e.g.,density matrix renormalization group White (1992); *White1993 and tensor hyper contraction of quantum chemical Hamiltonians Hohenstein et al. (2012); *Parrish2012).

With the new formalism, updating all the auxiliary fields on a single imaginary-time slice costs $\mathcal{O}(N_{s}m_{\text{T}}^{2})$ to calculate the force bias Shi et al. (2015) and $\mathcal{O}(m_{\text{T}}^{3})$ to calculate the ratio of the determinant in Eq. (2), applying $\det(\mathbf{1}_{N_{s}}+\mathbf{U}\mathbf{D}\mathbf{V})=\det(\mathbf{1}_{m_{\text{T}}}+\mathbf{D}\mathbf{V}\mathbf{U})$. Propagation along the path, i.e., moving $\ell$ to the left or right, becomes a propagation on $\mathbf{U}_{R}$ and $\mathbf{U}_{L}$, which can be achieved with $\mathcal{O}(N_{s}m_{\text{R},\text{L}}\log N_{s})$ scaling using FFT and locality of the interaction in Eq. (1). For a single numerical stabilization, the complexity is also lowered to $\mathcal{O}(N_{s}m_{\text{R},\text{L}}^{2})$, since we are dealing with $\mathbf{U}_{R}$ and $\mathbf{U}_{L}$ matrices. For measurements, we compute the Green’s function as $\mathbf{G}^{\sigma}(\tau,\tau)=(\mathbf{1}_{N_{s}}+\mathbf{U}\mathbf{D}\mathbf{V})^{-1}=\mathbf{1}_{N_{s}}-\mathbf{U}\mathbf{D}(\mathbf{1}_{m_{\text{T}}}+\mathbf{V}\mathbf{D}\mathbf{U})^{-1}\mathbf{V}$.

In the calculations, $m_{\text{R}}$, $m_{\text{L}}$ and $m_{\text{T}}$ evolve dynamically, as we carry out the truncations following each numerical stabilization procedure. At the start and end of a sweep, $m_{\text{R}}$ or $m_{\text{L}}$ is close to $N_{s}$, but away from the ends they rapidly decay, as illustrated in Fig. 2(a). The “truncation efficiency”, defined as $\eta_{\text{L},\text{R},\text{T}}=(N_{s}-m_{\text{L},\text{R},\text{T}})/(N_{s}-N_{\sigma})$, remains over 90% except for small portions at the ends. The computational complexity is dominated by the main part of the path in the middle, which scales as $\mathcal{O}(N_{s}m_{\text{T}}^{2})$, and the overall $m_{\text{T}}$ is small across the entire path. In standard DQMC, lowering the temperature makes the calculation progressively more challenging; as $T\rightarrow 0$, the computational efficiency reaches its worst in comparison with the ground-state method. In our new method, $m_{\text{T}}$ approaches $N_{e}$ as $T$ is lowered, and the algorithm restores the same computational complexity as the ground-state projection approach. At higher temperatures $m_{\text{T}}$ increases. In the example in Fig. 2(a) $m_{\text{T}}$ is $\sim 1.3\times$ and $\sim 2.5\times$ the ground-state value at $T/T_{F}=0.055$ and $0.274$, respectively, but even in the latter a speedup of $\sim 30$ is achieved. As mentioned earlier, the closer to the continuum limit, the larger the speedup, since $m_{\text{T}}$ only depends on physical parameters such as $N_{e}$ and $T$, and will change little as $N_{s}$ is increased.

As shown in the inset in Fig. 2(a), the truncation error is invisible even with an aggressive truncation threshold of $\epsilon=10^{-3}$. (The comparison of different $\epsilon$ vs. $\epsilon=0$, i.e. the standard algorithm, was done on an identical set of random paths. In practical calculations the different choices of $\epsilon$ would cause the calculations to, in sufficiently long runs, cease being correlated in their random number streams, but they will give statistically compatible results.) We tested that these results are not sensitive to the temperature. We have also carried out simulations for different interaction strengths ($U/t=-1$ and $-6$) and obtained similar results. Moreover, these results hold generally for different HS transformations, e.g., total density or spin-$s_{z}$ channels for Eq. (1) Hirsch (1983).

The speedup of our new algorithm over current state-of-the-art is illustrated in Fig. 2(b). The main figure displays the speedup vs. the linear dimension of the lattice for 2D Fermi gas system. The speedups fit well a fourth-order polynomial for $L$ up to $30$, consistent with the earlier conclusion of $(N_{s}/m_{\text{T}})^{2}$ from computational complexity. The speedup is around $50$ at $L=30$, and is extrapolated to $\sim 150$ for $L=45$, the data point for which is obtained for the new algorithm with five days’ computing (on 40 Skylake cores) and which would have required over two years’ computing with the standard algorithm using comparable resources. The inset illustrates the speedup vs. filling by tuning the chemical potential in the context of the Hubbard model. Large speedups are seen at low density, as expected. Even at half-filling, a speedup of $2.65$ is achieved.

We next apply the method to the spin-balanced uniform 2D Fermi gas with a zero-range attractive interaction. Major experimental efforts are being devoted to this and related systems Martiyanov et al. (2010); Ries et al. (2015); Murthy et al. (2015); Boettcher et al. (2016); Fenech et al. (2016); Hueck et al. (2018); Murthy et al. (2018), which offer opportunities for highly controllable and clean experiments in a strongly interacting 2D quantum system. To our knowledge, finite-$T$ ab initio quantum Monte Carlo computations to date have been mostly limited to the normal state Anderson and Drut (2015); Mulkerin et al. (2017). The algorithmic advances presented in this paper allow exact computations to reach much larger lattices and lower $T$’s, which adds a new dimension to our computational and theoretical capabilities, and will help calibrate and guide experiments. Here, as an illustration of the method, we present examples on the equation of state (EOS) and pairing properties.

In Fig. 3, we present the results of EOS at the interaction strength $\log(k_{F}a)=4.346574$, which is on the BCS side. As shown in the main plot, the chemical potential shifted by the bound-state energy, $(\mu+\varepsilon_{B}/2)/E_{F}$ [note $\varepsilon_{B}=4\hbar^{2}/ma^{2}e^{2\gamma}$ with Euler’s constant $\gamma=0.57721$], has perfect linear scaling with $1/L^{2}$ at large $L$ for all the temperatures. After extrapolation to the continuum limit, $(\mu+\varepsilon_{B}/2)/E_{F}$ converges to the the ground-state result when $T\rightarrow 0$, as shown in the inset. Our results indicate that the small discrepancy between experiment Boettcher et al. (2016) and ground-state results Shi et al. (2015) is consistent with finite-$T$ effects, which can now be better understood as experimental resolution improves.

We have also determined the pairing wave function and condensate fraction as a function of temperature. We compute the zero-momentum spin-singlet pairing matrix: $\mathbf{M}_{\mathbf{k}\mathbf{k}^{\prime}}=\langle\Delta_{\mathbf{k}}^{+}\Delta_{\mathbf{k}^{\prime}}\rangle-\delta_{\mathbf{k}\mathbf{k}^{\prime}}\langle c_{\mathbf{k}\uparrow}^{+}c_{\mathbf{k}\uparrow}\rangle\langle c_{-\mathbf{k}\downarrow}^{+}c_{-\mathbf{k\downarrow}}\rangle$ with $\Delta_{\mathbf{k}}^{+}=c_{\mathbf{k}\uparrow}^{+}c_{-\mathbf{k}\downarrow}^{+}$. Its leading eigenvalue divided by $N_{e}/2$ is identified as the condensate fraction. The corresponding eigenvector gives the pairing wavefunction $\phi_{\uparrow\downarrow}(\mathbf{k})$. Results are shown in Fig. 4 for $\log(k_{F}a)=0.50$, in the strongly interacting crossover regime. The mean-field results in the inset illustrate the strong correlation effect. The small difference between the finite-size and thermodynamic limit results also help provide a gauge of the finite-size effects in the many-body results. At high temperatures, the pairing wave function is more extended in momentum space, and the condensate fraction is tiny. As the temperature lowers, the pairing wave function peaks more at the Fermi surface, and the condensate fraction increases rapidly. The Berezinskii-Kosterlitz-Thouless (BKT) transition temperature measured experimentally Ries et al. (2015) is also shown in the inset of Fig. 4. It is clear that the new method now makes possible a quantitative comparison with experiments. A more systematic study of the physics of the 2D Fermi gas and the BKT transition will be published separately.

We have described the algorithm with some specificity in the context of DQMC for a lattice model. As mentioned, the advances are general and can be applied to many forms of finite-temperature field-theory computations. For example, the inclusion of spin-flip terms in the effective single-particle Hamiltonian in Eq. (2), such as spin-orbit coupling, can be accommodated straightforwardly Shi et al. (2016); Rosenberg et al. (2017). In the presence of a sign problem, e.g., repulsive or spin-imbalanced Hubbard-like models, the constrained path AFQMC framework Zhang (1999, 2000); Liu et al. (2018); He et al. (2019) can be used. In that case $\mathbf{L}$ becomes the trial propagator matrix $\mathbf{L}_{T}$, and only ${\mathbf{R}}$ depends on the path, which is active only up to the current $\ell$ and is “grown” stochastically to full length. The low-rank factorizations can be applied in exactly the same way. For electronic Hamiltonians with long-range interactions as in molecules and solids, a phase problem arises and the phaseless approximation Zhang and Krakauer (2003); Liu et al. (2018) is needed. The low-rank factorizations and the rest of the algorithm are identical to the sign problem case. (For a plane-wave basis Suewattana et al. (2007), one can alternate between Fourier and real space Suewattana et al. (2007) similar to the Hubbard model and the computational scaling is similar to what we have discussed. For generic basis sets such as in quantum chemistry, however, the $\mathbf{B}$’s have a more general structure and the propagation scales as $\mathcal{O}(N_{s}^{2}m_{\rm L,R})$, which means the reduction in computational complexity is one power of $N_{s}/N_{e}$.) Beyond fermion systems, the factorization and truncation also apply to bosons and Fermi-Bose mixtures Rubenstein et al. (2012). More broadly, the idea should also be applicable to calculations employing similar formalisms in nuclear shell models ALHASSID ; Alhassid (2016) and lattice QCD Chandrasekharan and Jiang (2003); Chandrasekharan (2006); Chandrasekharan and Mehta (2007); Bazavov et al. (2010); Petreczky (2012).

In summary, we have presented a method to perform finite-temperature field-theoretic calculations of many-fermion systems with computational complexity $\mathcal{O}(N_{s}N_{e}^{2})$, i.e,, linear in lattice (or single-particle basis) size. Such calculations are applied widely, and their fundamental complexity of $\mathcal{O}(N_{s}^{3})$ previously has been a major obstacle for approaching the continuum limit ($N_{s}/N_{e}\rightarrow\infty$). We demonstrated our new method in 2D strongly interacting Fermi gas, where the continuum limit can be reached via numerically exact calculations with large lattice sizes even at low temperatures. The method introduced here can be applied to a variety of strongly interacting fermion systems in multiple fields, including ultracold atomic gases, condensed matter and materials, and quantum chemistry.

We thank J. Carlson, P. Dumitrescu and S. Chandrasekharan for helpful discussions. This work was supported by NSF (Grant No. DMR-1409510). The Flatiron Institute is a division of the Simons Foundation.