We present a method for performing path integral molecular dynamics (PIMD) simulations for fermions and address its sign problem. PIMD simulations are widely used for studying many-body quantum systems at thermal equilibrium. However, they assume that the particles are distinguishable and neglect bosonic and fermionic exchange effects. Interacting fermions play a key role in many chemical and physical systems, such as electrons in quantum dots and ultracold trapped atoms. A direct sampling of the fermionic partition function is impossible using PIMD since its integrand is not positive definite. We show that PIMD simulations for fermions are feasible by employing our recently developed method for bosonic PIMD and reweighting the results to obtain fermionic expectation values. The approach is tested against path integral Monte Carlo (PIMC) simulations for up to seven electrons in a two-dimensional quantum dot for a range of interaction strengths. However, like PIMC, the method suffers from the sign problem at low temperatures. We propose a simple approach for alleviating it by simulating an auxiliary system with a larger average sign and obtaining an upper bound to the energy of the original system using the Bogoliubov inequality. This allows fermions to be studied at temperatures lower than would otherwise have been feasible using PIMD, as demonstrated in the case of a three-electron quantum dot. Our results extend the boundaries of PIMD simulations of fermions and will hopefully stimulate the development of new approaches for tackling the sign problem.

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.7 Developing simulation methods that accurately describe correlated fermions at various thermodynamic conditions is therefore very desirable.8 The path integral (PI) formulation of quantum mechanics9 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.10 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.11 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.13 

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.19 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.

FIG. 1.

Ring polymer configurations for two indistinguishable particles. The permutations of two particles are shown in green, and the corresponding cycle notations are shown in red. The ring polymer configurations can be directly inferred from the cycle notation for permutations.

FIG. 1.

Ring polymer configurations for two indistinguishable particles. The permutations of two particles are shown in green, and the corresponding cycle notations are shown in red. The ring polymer configurations can be directly inferred from the cycle notation for permutations.

Close modal

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.18 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.20 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.21 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.22 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.”23 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,24 configuration PIMC,25 and permutation blocking PIMC.26 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.27 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.28 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 PIMD19 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.29 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 Dornheim30 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

Ĥ=12ml=1Np^l2+V^(r1,,rN).
(1)

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 kB is the Boltzmann constant. The path integral expression for the partition function31 is

ZIlimPeβUI(N)dR1dRN,
(2)

where I represents whether the particles are bosons (B), fermions (F), or distinguishable (D) and Rl 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

UD(N)=12mωP2l=1Nj=1P(rlj+1rlj)2+1Pj=1PV(r1j,,rNj),
(3)

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.10 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.19 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 algorithm32–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 approach19 to fermions is achieved by defining the potential as

UB/F(N)=1βlnWB/F(N)+1Pj=1PV(r1j,,rNj)
(4)

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

WB/F(N)=1Nk=1Nξk1eβEN(k)WB/F(Nk).
(5)

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

EN(k)=12mωP2l=Nk+1Nj=1Prlj+1rlj2.
(6)

In Eq. (6), rlP+1=rl+11, except for l = N for which rNP+1=rNk+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 Ô(r1, …, rN), such as the potential energy, the estimator εO is

εO=1Pj=1PO(r1j,..,rNj),
(7)

and the expectation value is obtained by

ÔI=1ZIεOeβUI(N)dR1dRNεOI.
(8)

For evaluating the kinetic energy, we use the virial estimator39 (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.23 The fermionic expectation values are then recovered from the bosonic simulations by reweighting,

ÔF=εOsBsB,
(9)

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 ⟨sB decreases exponentially with β and the number of particles.23 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 ⟨sB is large enough, the simulations of electrons in two-dimensional QDs can be converged and are in very good agreement with PIMC results.30 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 collaborators4–6 have shown that the two-dimensional QDs can be modeled using the dimensionless Hamiltonian

Ĥ=12l=1Nl2+12l=1Nrl2+l,m>lNλ|rlrm|.
(10)

In Eq. (10), λ = e2/(κl0ℏω0) is the Wigner parameter, κ is the relative QD dielectric constant, and ω0 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 ℏω0 = 5.1 meV and the effective mass was set to m = 0.07me, which are realistic values for an electron in a GaAs QD.4 We use a development version of the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS)45 to perform all simulations with a time step of 1 fs and a massive Nose–Hoover chains thermostat.46 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.31 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 βℏω0 = 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 ≤ βℏω0 ≤ 1; λ = 0.5) and interaction parameters (0 ≤ λ ≤ 4; βℏω0 = 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.4 

FIG. 2.

The average thermal energy (upper panel) and sign (lower panel) as a function of the number of interacting fermions N in a two-dimensional QD. The interaction strength is λ = 0.5 and βℏω0 = 1. Blue and orange circles represent the results obtained using PIMD for bosons and fermions, respectively. Green squares represent PIMC results of Dornheim.30 Error bars are smaller than the symbol size.

FIG. 2.

The average thermal energy (upper panel) and sign (lower panel) as a function of the number of interacting fermions N in a two-dimensional QD. The interaction strength is λ = 0.5 and βℏω0 = 1. Blue and orange circles represent the results obtained using PIMD for bosons and fermions, respectively. Green squares represent PIMC results of Dornheim.30 Error bars are smaller than the symbol size.

Close modal
FIG. 3.

The average thermal energy (upper panel) and sign (lower panel) for N = 6 interacting fermions (λ = 0.5) in a two-dimensional QD as a function of inverse temperature. Blue and orange circles represent the results obtained using PIMD for bosons and fermions, respectively. Green squares represent PIMC results of Dornheim.30 Error bars are smaller than the symbol size.

FIG. 3.

The average thermal energy (upper panel) and sign (lower panel) for N = 6 interacting fermions (λ = 0.5) in a two-dimensional QD as a function of inverse temperature. Blue and orange circles represent the results obtained using PIMD for bosons and fermions, respectively. Green squares represent PIMC results of Dornheim.30 Error bars are smaller than the symbol size.

Close modal
FIG. 4.

The average thermal energy (upper panel) and sign (lower panel) at βℏω0 = 1 for N = 6 fermions in a two-dimensional QD as a function of interaction strength λ. Blue and orange circles represent the results obtained using PIMD for bosons and fermions, respectively. Green squares represent PIMC results of Dornheim.30 Error bars are smaller than the symbol size.

FIG. 4.

The average thermal energy (upper panel) and sign (lower panel) at βℏω0 = 1 for N = 6 fermions in a two-dimensional QD as a function of interaction strength λ. Blue and orange circles represent the results obtained using PIMD for bosons and fermions, respectively. Green squares represent PIMC results of Dornheim.30 Error bars are smaller than the symbol size.

Close modal

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 Dornheim30 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, βℏω0 = 1, and λ = 0 had the smallest average sign with ⟨sB ≈ 0.003. They required ten independent trajectories of 9 · 108 MD steps in order to converge the energy to within 2%. In comparison, the same simulations with λ = 3 (⟨sB ≈ 0.64) required only five simulations of 75 · 106 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.23 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 inequality29 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,29 

FĤFĤĤĤĤ.
(11)

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

EĤEĤ+FĤFĤ.
(12)

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.

FIG. 5.

For a fixed number of particles, the average sign depends on the inverse temperature β and the repulsion strength λ. The blue line represents some (arbitrary) critical value of the average sign. When the sign of the system is smaller (red), simulations are too expensive to converge. We propose simulating an auxiliary system that has a larger sign (green) and can be converged in practical simulations. Then, we estimate the energy of the original system, using the Bogoliubov inequality, at lower temperatures than otherwise would have been feasible using PIMD.

FIG. 5.

For a fixed number of particles, the average sign depends on the inverse temperature β and the repulsion strength λ. The blue line represents some (arbitrary) critical value of the average sign. When the sign of the system is smaller (red), simulations are too expensive to converge. We propose simulating an auxiliary system that has a larger sign (green) and can be converged in practical simulations. Then, we estimate the energy of the original system, using the Bogoliubov inequality, at lower temperatures than otherwise would have been feasible using PIMD.

Close modal

We demonstrate this approach for a system of three electrons confined to a two-dimensional QD6 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 βℏω0 = 2.5 where the error is visibly much larger (∼10%) than at higher temperatures.

FIG. 6.

The average thermal energy (upper panel) and sign (lower panel) for three non-interacting electrons in a two-dimensional QD as a function of inverse temperature. Orange circles represent the results obtained by converging Eq. (9) directly. Gray circles represent the results obtained for an auxiliary system, and purple circles show the results obtained using the Bogoliubov inequality. Analytical results are represented by green squares. If error bars are not shown, they are smaller than the symbol size. In the lower panel, red and brown diamonds present the average sign of the non-interacting system and the auxiliary system, respectively.

FIG. 6.

The average thermal energy (upper panel) and sign (lower panel) for three non-interacting electrons in a two-dimensional QD as a function of inverse temperature. Orange circles represent the results obtained by converging Eq. (9) directly. Gray circles represent the results obtained for an auxiliary system, and purple circles show the results obtained using the Bogoliubov inequality. Analytical results are represented by green squares. If error bars are not shown, they are smaller than the symbol size. In the lower panel, red and brown diamonds present the average sign of the non-interacting system and the auxiliary system, respectively.

Close modal

Since it was shown by Dornheim30 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

V(|rlrm|)=gπσ2e(rlrm)2σ2.
(13)

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 = me and the frequency of the trap to ℏω0 = 3 meV, as done by Runeson et al.28 For each temperature, we performed five independent simulations of 75 · 106 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 βℏω0 = 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.

1.
I.
Bloch
, “
Ultracold quantum gases in optical lattices
,”
Nat. Phys.
1
,
23
30
(
2005
).
2.
E.
Haller
,
J.
Hudson
,
A.
Kelly
,
D. A.
Cotta
,
B.
Peaudecerf
,
G. D.
Bruce
, and
S.
Kuhr
, “
Single-atom imaging of fermions in a quantum-gas microscope
,”
Nat. Phys.
11
,
738
742
(
2015
).
3.
U.
Schneider
,
L.
Hackermüller
,
J. P.
Ronzheimer
,
S.
Will
,
S.
Braun
,
T.
Best
,
I.
Bloch
,
E.
Demler
,
S.
Mandt
,
D.
Rasch
, and
A.
Rosch
, “
Fermionic transport and out-of-equilibrium dynamics in a homogeneous Hubbard model with ultracold atoms
,”
Nat. Phys.
8
,
213
218
(
2012
).
4.
C.
Ellenberger
,
T.
Ihn
,
C.
Yannouleas
,
U.
Landman
,
K.
Ensslin
,
D.
Driscoll
, and
A. C.
Gossard
, “
Excitation spectrum of two correlated electrons in a lateral quantum dot with negligible Zeeman splitting
,”
Phys. Rev. Lett.
96
,
126806
(
2006
).
5.
C.
Yannouleas
and
U.
Landman
, “
Symmetry breaking and quantum correlations in finite systems: Studies of quantum dots and ultracold Bose gases and related nuclear and chemical methods
,”
Rep. Prog. Phys.
70
,
2067
2148
(
2007
).
6.
Y.
Li
,
C.
Yannouleas
, and
U.
Landman
, “
Three-electron anisotropic quantum dots in variable magnetic fields: Exact results for excitation spectra, spin structures, and entanglement
,”
Phys. Rev. B
76
,
245310
(
2007
).
7.
Y.
Cytter
,
E.
Rabani
,
D.
Neuhauser
,
M.
Preising
,
R.
Redmer
, and
R.
Baer
, “
Transition to metallization in warm dense helium-hydrogen mixtures using stochastic density functional theory within the Kubo-Greenwood formalism
,”
Phys. Rev. B
100
,
195101
(
2019
).
8.
L. K.
Wagner
and
D. M.
Ceperley
, “
Discovering correlated fermions using quantum Monte Carlo
,”
Rep. Prog. Phys.
79
,
094501
(
2016
).
9.
R. P.
Feynman
and
A. R.
Hibbs
,
Quantum Mechanics and Path Integrals
(
Dover pulications
,
2005
).
10.
D.
Chandler
and
P. G.
Wolynes
, “
Exploiting the isomorphism between quantum theory and classical statistical mechanics of polyatomic fluids
,”
J. Chem. Phys.
74
,
4078
4095
(
1981
).
11.
T. E.
Markland
and
M.
Ceriotti
, “
Nuclear quantum effects enter the mainstream
,”
Nature Reviews Chemistry
2
,
0109
(
2018
).
12.
M.
Parrinello
and
A.
Rahman
, “
Study of an F center in molten KCl
,”
J. Chem. Phys.
80
,
860
867
(
1984
).
13.
D. M.
Ceperley
, “
Path integrals in the theory of condensed helium
,”
Rev. Mod. Phys.
67
,
279
355
(
1995
).
14.
A. P.
Lyubartsev
and
P. N.
Vorontsov-Velyaminov
, “
Path-integral Monte Carlo method in quantum statistics for a system of N identical fermions
,”
Phys. Rev. A
48
,
4075
4083
(
1993
).
15.
E. L.
Pollock
and
D. M.
Ceperley
, “
Simulation of quantum many-body systems by path-integral methods
,”
Phys. Rev. B
30
,
2555
2568
(
1984
).
16.
D. M.
Ceperley
and
E. L.
Pollock
, “
Path-integral computation of the low-temperature properties of liquid 4He
,”
Phys. Rev. Lett.
56
,
351
354
(
1986
).
17.
E. L.
Pollock
and
D. M.
Ceperley
, “
Path-integral computation of superfluid densities
,”
Phys. Rev. B
36
,
8343
8352
(
1987
).
18.
D. M.
Ceperley
, “
Path-integral calculations of normal liquid 3He
,”
Phys. Rev. Lett.
69
,
331
334
(
1992
).
19.
B.
Hirshberg
,
V.
Rizzi
, and
M.
Parrinello
, “
Path integral molecular dynamics for bosons
,”
Proc. Natl. Acad. Sci. U. S. A.
116
,
21445
21449
(
2019
).
20.
E. W.
Brown
,
B. K.
Clark
,
J. L.
DuBois
, and
D. M.
Ceperley
, “
Path-integral Monte Carlo simulation of the warm dense homogeneous electron gas
,”
Phys. Rev. Lett.
110
,
146405
(
2013
).
21.
M.
Takahashi
and
M.
Imada
, “
Monte Carlo calculation of quantum systems
,”
J. Phys. Soc. Jpn.
53
,
963
974
(
1984
).
22.
J. L.
DuBois
,
E. W.
Brown
, and
B. J.
Alder
, “
Overcoming the fermion sign problem in homogeneous systems
,”
In Advances in the Computational Sciences
(
World Scientific
,
2017
), pp.
184
192
.
23.
M.
Troyer
and
U.-J.
Wiese
, “
Computational complexity and fundamental limitations to fermionic quantum Monte Carlo simulations
,”
Phys. Rev. Lett.
94
,
170201
(
2005
).
24.
C. H.
Mak
,
R.
Egger
, and
H.
Weber-Gottschick
, “
Multilevel blocking approach to the fermion sign problem in path-integral Monte Carlo simulations
,”
Phys. Rev. Lett.
81
,
4533
4536
(
1998
).
25.
T.
Schoof
,
M.
Bonitz
,
A.
Filinov
,
D.
Hochstuhl
, and
J. W.
Dufty
, “
Configuration path integral Monte Carlo
,”
Contrib. Plasma Phys.
51
,
687
697
(
2011
).
26.
T.
Dornheim
,
S.
Groth
,
A.
Filinov
, and
M.
Bonitz
, “
Permutation blocking path integral Monte Carlo: A highly efficient approach to the simulation of strongly degenerate non-ideal fermions
,”
New J. Phys.
17
,
073017
(
2015
).
27.
S.
Miura
and
S.
Okazaki
, “
Path integral molecular dynamics for Bose–Einstein and Fermi–Dirac statistics
,”
J. Chem. Phys.
112
,
10116
10124
(
2000
).
28.
J.
Runeson
,
M.
Nava
, and
M.
Parrinello
, “
Quantum symmetry from enhanced sampling methods
,”
Phys. Rev. Lett.
121
,
140602
(
2018
).
29.
R. P.
Feynman
,
Frontiers in Physics
(
Perseus Books
,
1972
), pp.
67
71
.
30.
T.
Dornheim
, “
Fermion sign problem in path integral Monte Carlo simulations: Quantum dots, ultracold atoms, and warm dense matter
,”
Phys. Rev. E
100
,
023307
(
2019
).
31.
M.
Tuckerman
,
Statistical Mechanics: Theory and Molecular Simulation
(
Oxford University Press
,
2010
).
32.
Ł.
Walewski
,
H.
Forbert
, and
D.
Marx
, “
Reactive path integral quantum simulations of molecules solvated in superfluid helium
,”
Comput. Phys. Commun.
185
,
884
899
(
2014
).
33.
Ł.
Walewski
,
H.
Forbert
, and
D.
Marx
, “
Solvation of molecules in superfluid helium enhances the “interaction induced localization” effect
,”
J. Chem. Phys.
140
,
144305
(
2014
).
34.
F.
Uhl
and
D.
Marx
, “
Helium tagging of protonated methane in messenger spectroscopy: Does it interfere with the fluxionality of CH5+?
,”
Angew. Chem., Int. Ed.
57
,
14792
14795
(
2018
).
35.
C.
Schran
,
F.
Uhl
,
J.
Behler
, and
D.
Marx
, “
High-dimensional neural network potentials for solvation: The case of protonated water clusters in helium
,”
J. Chem. Phys.
148
,
102310
(
2018
).
36.
F.
Uhl
and
D.
Marx
, “
Quantum microsolvation of protonated methane with 4He: Large-amplitude motion heavily influences boson
,”
Phys. Rev. Lett.
123
,
123002
(
2019
).
37.
P.
Borrmann
, “
Path integral density functional theory
,” arXiv:9412117 [cond-mat] (
1994
).
38.
P.
Borrmann
,
J.
Harting
,
O.
Mülken
, and
E. R.
Hilf
, “
Calculation of thermodynamic properties of finite Bose-Einstein systems
,”
Phys. Rev. A
60
,
1519
1522
(
1999
).
39.
M. F.
Herman
,
E. J.
Bruskin
, and
B. J.
Berne
, “
On path integral Monte Carlo simulations
,”
J. Chem. Phys.
76
,
5150
5155
(
1982
).
40.
G.
Parisi
, “
On complex probabilities
,”
Phys. Lett. B
131
,
393
395
(
1983
).
41.
G.
Aarts
,
E.
Seiler
, and
I. O.
Stamatescu
, “
Complex Langevin method: When can it be trusted?
,”
Phys. Rev. D
81
,
1
13
(
2010
), arXiv:0912.3360.
42.
X.
Wang
,
G.
Sun
,
N.
Li
, and
P.
Chen
, “
Quantum dots derived from two-dimensional materials and their applications for catalysis and energy
,”
Chem. Soc. Rev.
45
,
2239
2262
(
2016
).
43.
Y.
Xu
,
X.
Wang
,
W. L.
Zhang
,
F.
Lv
, and
S.
Guo
, “
Recent progress in two-dimensional inorganic quantum dots
,”
Chem. Soc. Rev.
47
,
586
625
(
2018
).
44.
A.
Manikandan
,
Y.-Z.
Chen
,
C.-C.
Shen
,
C.-W.
Sher
,
H.-C.
Kuo
, and
Y.-L.
Chueh
, “
A critical review on two-dimensional quantum dots (2D QDs): From synthesis toward applications in energy and optoelectronics
,”
Prog. Quantum Electron.
68
,
100226
(
2019
).
45.
S.
Plimpton
, “
Fast parallel algorithms for short-range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
46.
G. J.
Martyna
,
M. L.
Klein
, and
M.
Tuckerman
, “
Nosé–Hoover chains: The canonical ensemble via continuous dynamics
,”
J. Chem. Phys.
97
,
2635
2643
(
1992
).

Supplementary Material