Path integral molecular dynamics for fermions: Alleviating the sign problem with the Bogoliubov inequality Free

Interacting fermions at finite temperature are fundamental to a wide range of chemical and physical phenomena. Examples include ultracold trapped atoms,^1–3 electrons in quantum dots,^4–6 and warm dense matter.⁷ Developing simulation methods that accurately describe correlated fermions at various thermodynamic conditions is therefore very desirable.⁸ The path integral (PI) formulation of quantum mechanics⁹ is a powerful approach for studying many-body systems at finite temperature. It is based on the observation that the partition function of a quantum system is isomorphic to the partition function of a fictitious, extended classical system.¹⁰ Assuming that the quantum particles are distinguishable, each one is represented by a classical “ring polymer” composed of P replicas of the particle (the “beads” of the ring polymer) connected by harmonic springs.¹¹ The frequency of the springs is proportional to the temperature and number of beads. The exact quantum result is obtained in the limit P → ∞. Beads having the same index but representing different particles interact through a scaled interaction potential. The partition function of the fictitious classical system is typically sampled using molecular dynamics (MD)^11,12 or Monte Carlo (MC) methods.¹³ 

For N indistinguishable particles, bosons or fermions, one must consider all N! particle permutations. As a result, the fictitious classical system is no longer composed of a single ring polymer for each particle. Ring polymer configurations in which particles are connected into longer rings must be included,^13,14 as shown in Fig. 1 in the two-particle case. The number of ring polymer configurations scales exponentially with system size, and enumerating them is impractical for more than a few atoms. In Path Integral Monte Carlo (PIMC) simulations, this problem is avoided by using MC moves designed to sample particle permutations, as pioneered by Ceperley and co-workers.^13,15–18 Recently, we showed that the potential and forces required for path integral molecular dynamics (PIMD) simulations for bosons can be evaluated recursively without enumerating or sampling particle permutations.¹⁹ The resulting algorithm scales cubically with system size, which has allowed the first PIMD simulations of large bosonic systems to be performed. In this communication, we show how this method can be extended to obtain thermodynamic properties for fermionic systems. We will focus on the spin-polarized case for clarity.

In comparison to bosons, fermions present an additional challenge since even and odd sequential particle permutations carry a positive and negative sign, respectively. As a result, direct sampling of fermionic partition functions is problematic, either in MC or in MD, since their integrands are not positive definite. One approach to overcome this problem within the framework of PIMC is the fixed node approximation, often referred to as restricted PIMC.¹⁸ Since the exact nodal surface of the fermionic density matrix is usually unknown, an approximate ansatz is employed, such as the free electron density matrix.²⁰ An alternative approach is to simulate the system as if it were composed of bosons and reweight the results according to the sign of each sampled configuration.²¹ However, this procedure leads to the infamous fermionic sign problem—the exponentially slow convergence of expectation values due to the delicate cancellation of positive and negative contributions.²² It was shown by Troyer and Wiese that a general solution to the sign problem—an algorithm that converges in polynomial time—is “almost certainly unattainable.”²³ Even if a complete solution is unlikely, it is highly desirable to develop new methods for simulating fermions, which may provide avenues for alleviating the sign problem.

Recently, several methods have been proposed to alleviate the sign problem within PIMC, such as the multi-level blocking approach,²⁴ configuration PIMC,²⁵ and permutation blocking PIMC.²⁶ In contrast, only a few methods have been proposed to perform fermionic PIMD simulations in general, and even fewer studies have addressed the sign problem in this context. Miura and Okazaki proposed to sample the absolute value of the probability density for fermions.²⁷ The potential is obtained from the absolute value of a determinant, which accounts for all particle permutations. Expectation values for fermions are obtained by reweighting the contribution of different configurations according to the sign of their determinant. However, this method was applied to only three particles, and its sign problem was not investigated in detail. Another limitation of this approach is that the potential has an infinite barrier at the nodes of the fermionic probability density, effectively preventing the system from exploring all regions of phase space. More recently, Runeson et al.²⁸ have shown that by combining enhanced sampling techniques with PIMD simulations of distinguishable particles, the problem of two electrons in a quantum dot can be solved. However, it is not clear how to extend their work to larger systems.

In this communication, we propose a different approach to performing PIMD simulations of fermions. We employ our recently developed method for bosonic PIMD¹⁹ and reweight the results to obtain fermionic expectation values. The method is applicable to any system for which the bosonic simulation can be performed and naturally avoids the problem of the infinite barriers, since the bosonic probability density is nodeless. However, as in PIMC, the method suffers from the sign problem at low temperatures. We propose a simple approach to alleviate the sign problem based on the Bugoliubov inequality for the free energy.²⁹ This allows us to estimate the energies of fermionic systems at lower temperatures than would otherwise have been possible using PIMD. In the following, we first present the method in detail. Then, we benchmark the new approach against recent fermionic PIMC simulations by Dornheim³⁰ for up to N = 7 electrons in a two-dimensional quantum dot (QD). In Sec. IV, we propose to use the Bogoliubov inequality to alleviate the sign problem and apply it to three-electron QDs at low temperatures. Section V presents the conclusions and briefly discusses future directions.

The Hamiltonian operator for a system of N identical particles of mass m interacting through a potential $V^$ is

The partition function is obtained by taking the trace of the density operator for a thermal state, Z = Tr{e^−βĤ}, where $β=(kBT)−1$ is the inverse temperature and k_B is the Boltzmann constant. The path integral expression for the partition function³¹ is

where I represents whether the particles are bosons (B), fermions (F), or distinguishable (D) and R_l represents collectively the coordinates of P replicas of particle l, $rl1,…,rlP$⁠. In practice, P is increased until the desired expectation values are converged. Here, we focus on sampling the partition function of Eq. (2) using MD simulations.

For distinguishable particles, the potential is given by

where $ωP=P/βℏ$ and $rlP+1=rl1$⁠. The first and second terms arise from the kinetic and potential energy operators of the quantum system, respectively. Equation (3) is isomorphic to the partition function of N classical ring polymers, each one composed of P beads that are connected by harmonic springs.¹⁰ Beads j of different ring polymers (representing different quantum particles) interact through a scaled interaction potential.

In order to treat bosons and fermions, the trace should be taken over a properly symmetrized or anti-symmetrized basis, respectively. Since the potential energy is invariant under permutations of identical particles, only the term arising from the kinetic energy operator is affected. Each permutation leads to a ring polymer configuration that can be directly inferred from the cycle notation for the permutation, as demonstrated in Fig. 1 for two particles. Since the number of permutations is N!, direct enumeration of all ring-polymer configurations is possible only for very small systems. We have recently shown that the potential for bosons can be evaluated recursively, avoiding the need to enumerate the exponentially large number of permutations.¹⁹ The algorithm scales cubically with system size, which allowed the first applications of PIMD to large bosonic systems. It is also possible to combine PIMC for bosons with PIMD for distinguishable particles in a hybrid algorithm^32–36 to treat large systems. However, our method is simpler as it treats both distinguishable and indistinguishable particles on an equal footing.

Extending our bosonic PIMD approach¹⁹ to fermions is achieved by defining the potential as

and evaluating the argument of the logarithm by using the recurrence relation

In Eq. (5), ξ = 1 for bosons, ξ = −1 for fermions, and we set $WB/F(0)=1$⁠. $EN(k)(RN−k+1,…,RN)$ is the spring energy of k particles connected together in a long ring, given by

In Eq. (6), $rlP+1=rl+11$⁠, except for l = N for which $rNP+1=rN−k+11$⁠. We note that Eqs. (4) and (5) lead to a slightly different expression for the potential for fermions than the one obtained by evaluating a determinant to account for particle permutations.^21,27 However, both methods converge to the same quantum-mechanical expectation values after integration over all particles in the partition function. In addition, Eqs. (4) and (5) reduce to a known recursion relation for the partition function in the case of non-interacting particles.^37,38 Expectation values are evaluated using the regular estimators. For a position dependent operator Ô(r₁, …, r_N), such as the potential energy, the estimator ε_O is

and the expectation value is obtained by

For evaluating the kinetic energy, we use the virial estimator³⁹ (see the supplementary material for more details).

The key difference between fermions and bosons is that, in the case of fermions, odd permutations carry a negative sign. Therefore, $WF(N)$ is no longer positive definite. When it becomes negative, the potential $UF(N)$ is complex and cannot be sampled in standard MD or MC simulations. Attempts to overcome this problem by using the complex Langevin equation instead have had limited success.^40,41 As mentioned in Sec. I, using the absolute value of $WF(N)$ will result in a real potential, but with an infinite barrier when $WF(N)=0$⁠.

Here, we propose instead to simulate the system as if it were composed of bosons, similarly to what is done in PIMC simulations.²³ The fermionic expectation values are then recovered from the bosonic simulations by reweighting,

where $s=WF(N)/WB(N)$ is the signed relative weight for each sampled configuration. The evaluation of $WF(N)$ is done concurrently with the evaluation of the bosonic potential and therefore presents only a modest additional computational cost. The average sign ⟨s⟩_B decreases exponentially with β and the number of particles.²³ As a result, the relative error for the expectation values evaluated using Eq. (9) becomes exponentially hard to converge, which is referred to as the fermionic sign problem. In Sec. III, we show that as long as ⟨s⟩_B is large enough, the simulations of electrons in two-dimensional QDs can be converged and are in very good agreement with PIMC results.³⁰ Then, we propose a procedure to alleviate the sign problem and evaluate the thermal energy of fermionic systems at lower temperatures than would otherwise have been possible using PIMD.

To test the new method for fermionic PIMD, we performed simulations of electrons in two-dimensional QDs. These materials have received much attention in recent years due to their interesting electronic, optical, and catalytic properties and their potential use in technological applications.^42–44 Landman and collaborators^4–6 have shown that the two-dimensional QDs can be modeled using the dimensionless Hamiltonian

In Eq. (10), λ = e²/(κl₀ℏω₀) is the Wigner parameter, κ is the relative QD dielectric constant, and ω₀ is the frequency of the two-dimensional trap. We denote the effective mass of the electrons as m, and $l0=ℏ/mω0$ is the characteristic length of the dot. The Wigner parameter represents the ratio between the screened Coulomb repulsion in the QD and the confinement of the harmonic trap. When λ < 1, exchange effects dominate since the Coulomb repulsion is screened, and conversely, when λ > 1, the repulsive interaction is dominant. In the simulations below, the trap is isotropic with a frequency ℏω₀ = 5.1 meV and the effective mass was set to m = 0.07m_e, which are realistic values for an electron in a GaAs QD.⁴ We use a development version of the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)⁴⁵ to perform all simulations with a time step of 1 fs and a massive Nose–Hoover chains thermostat.⁴⁶ The equations of motion were integrated numerically in the normal-mode representation (of the distinguishable free ring polymer configuration), and the inertial masses were scaled by the corresponding eigenvalues.³¹ Expectation values were converged with respect to the number of beads using P = 12 at most.

Figure 2 shows the average thermal energy for a system with βℏω₀ = 1 and λ = 0.5 as a function of system size. Excellent agreement is obtained in comparison to PIMC with a maximum absolute error of 0.6% and a mean absolute error of 0.4%. Figures 3 and 4 present the average thermal energies for N = 6 as a function of the inverse temperature and repulsion strength. For the range of temperatures (0.5 ≤ βℏω₀ ≤ 1; λ = 0.5) and interaction parameters (0 ≤ λ ≤ 4; βℏω₀ = 1) considered here, we obtain very good agreement with PIMC results with mean absolute errors of 0.8% and 0.9%, respectively. The simulations span a realistic range of interaction strength as the value of the interaction parameter in a GaAs QD is λ ≈ 1.55.⁴ 

The lower panels of Figs. 2–4 also present the average sign as a function of system size, temperature, and repulsion strength. They show that the average sign decreases with an increase in β and the number of particles, but that the sign problem is more severe for weakly interacting electrons than for strongly interacting electrons. This is in agreement with the detailed investigation by Dornheim³⁰ and physical intuition. Converging the fermionic expectation values becomes harder with a decrease in average sign. Of the results presented in this section, the simulations for N = 6, βℏω₀ = 1, and λ = 0 had the smallest average sign with ⟨s⟩_B ≈ 0.003. They required ten independent trajectories of 9 · 10⁸ MD steps in order to converge the energy to within 2%. In comparison, the same simulations with λ = 3 (⟨s⟩_B ≈ 0.64) required only five simulations of 75 · 10⁶ MD steps in order to converge the energy to within 0.05%. The error bars on expectation values were evaluated using a weighted average over the independent simulations, and full details are given in the supplementary material. Section IV describes the use of the Bogoliubov free energy principle to obtain expectation values for systems with an even smaller average sign.

Overcoming the sign problem is a formidable challenge. As mentioned above, a general solution is probably unattainable.²³ Nevertheless, developing methods that alleviate the sign problem for a wide class of systems and extend the boundaries of current simulations is important. Ideally, this should be done in a controlled manner. Here, we propose a simple approach that provides accurate estimates of fermionic energies at lower temperatures. We demonstrate its usefulness for the model of electrons in two-dimensional QDs introduced in Sec. III, but the method can be more generally applied. The key idea involves using the Bogoliubov inequality²⁹ to obtain an upper bound on the fermionic energies.

If the system Hamiltonian is denoted by Ĥ and the Hamiltonian of a second, auxiliary system is denoted by Ĥ′, an upper bound on the difference in free energies between the two systems is given by the Bogoliubov inequality,²⁹ 

In the following, we also assume that the energy of the system at low temperatures can be approximated by its free energy,

Equation (11) is very general, and we can choose the auxiliary Hamiltonian such that it alleviates the sign problem. For example, it was shown in Sec. III that the average sign increases with repulsive interaction strength. Thus, it is possible to converge the simulations for strongly repelling fermions at temperatures lower than for weakly repelling ones. If we choose the former as our auxiliary Hamiltonian, we can obtain an upper bound on the energy of the weakly interacting system at low temperatures using the equations above. This is illustrated schematically in Fig. 5.

We demonstrate this approach for a system of three electrons confined to a two-dimensional QD⁶ in the limit of infinite screening (λ = 0). As shown in Sec. III, this leads to the most severe sign problem. This is also emphasized in the upper panel of Fig. 6. Orange circles show that converging the simulations for the non-interacting system becomes difficult at βℏω₀ = 2.5 where the error is visibly much larger (∼10%) than at higher temperatures.

Since it was shown by Dornheim³⁰ that short-range potentials have a larger average sign than systems with long-range interactions, we choose an auxiliary Hamiltonian in which the electrons interact with a repulsive Gaussian pair potential

In Eq. (13), g and σ represent the interaction strength and range, respectively. In the following simulations, we used g = 16 and σ = 0.5 in harmonic oscillator units. The effective mass was set to m = m_e and the frequency of the trap to ℏω₀ = 3 meV, as done by Runeson et al.²⁸ For each temperature, we performed five independent simulations of 75 · 10⁶ steps each. Expectation values were converged with respect to the number of beads using P = 12 for the non-interacting system and P = 72 for the auxiliary system.

The bottom panel of Fig. 6 provides insight into why this approach is expected to succeed. The red symbols show the average sign of the non-interacting system as a function of temperature. It decreases sharply and is equal to ≈0.008 at βℏω₀ = 2.5, which leads to the large error bars at this temperature in the upper panel. In contrast, the average sign for the auxiliary system decays significantly more slowly, which allows simulations at lower temperatures to be converged more easily.

The thermal energies obtained for the auxiliary system are shown in Fig. 6 (gray circles) and compared to the analytical results (green squares) for the non-interacting system. They are higher than the energies of non-interacting particles due to the added repulsion, as expected. Applying the Bogoliubov inequality (purple circles), we obtain upper bounds on the energies of the non-interacting system that agree with the analytical results within the statistical error. Most importantly, our approach allows an accurate estimate of the energy of the three-electron QD to be obtained at temperatures three times lower than would otherwise be feasible.

The results are not very sensitive to the value of the interaction parameter, as shown in the supplementary material, but it must be large enough to converge the simulations at lower temperatures. We note, however, that this is likely to be dependent on the type of repulsive interaction in the auxiliary system. In this case, the energy of the repulsive system grows slowly with the interaction parameter g (see the supplementary material). In other cases, one should remember that the Bogoliubov inequality provides only an upper bound. Therefore, a good strategy is to choose a repulsion large enough to converge the simulations at low temperatures but as small as possible to obtain an upper bound, which is close to the real expectation value. One could also imagine choosing variationally the auxiliary potential, which is beyond the scope of this work.

In this communication, we present a new method for performing PIMD simulations for fermions. Since direct sampling of fermionic partition functions is not feasible in PIMD (or PIMC), we perform simulations as if the system were composed of bosons and reweight the result to obtain fermionic expectation values using Eq. (9). We apply the new method for up to N = 7 electrons in a two-dimensional QD and obtain excellent agreement with PIMC results. However, like PIMC and other exact methods, our approach suffers from the fermionic sign problem. We show that the average sign decreases, as expected, with an increase in β and system size but also with a decrease in repulsion strength. To alleviate the sign problem for weakly interacting systems, we propose performing simulations of an auxiliary system with a larger average sign (e.g., by adding a fictitious repulsive interaction potential) and correcting its thermal energy using the Bogoliubov inequality. The results obtained in this way provide an upper bound to the energy of the weakly interacting system. We show that using this approach, we are able to solve the problem of a three-electron QD in the limit of infinite screening for which the sign problem is most severe. We obtain accurate energies at temperatures that are three times lower than would have been possible by trying to converge the reweighting scheme directly. Challenges still remain in order to perform such simulations for a large class of interesting systems. For example, at even lower temperatures (or for larger systems), our method too will suffer from the sign problem. However, this work extends the boundaries of PIMD simulations of fermions, and we expect that it will encourage development of new methods to attack this important problem.

See the supplementary material for a discussion of the error analysis, the virial estimator, and the sensitivity of the results to the repulsion strength parameter.

Additional raw data that support the findings of this study are openly available in GitHub (https://github.com/BarakHirshberg).

This research was supported by the European Union (Grant No. ERC-2014-ADG-670227/VARMET) and the NCCR MARVEL, funded by the Swiss National Science Foundation. Calculations were carried out on the Euler cluster at ETH Zurich.