Attenuating the fermion sign problem in path integral Monte Carlo simulations using the Bogoliubov inequality and thermodynamic integration
Tobias Dornheim, Michele Invernizzi, Jan Vorberger, Barak Hirshberg
AAttenuating the fermion sign problem in path integral Monte Carlo simulations usingthe Bogoliubov inequality and thermodynamic integration
Tobias Dornheim, ∗ Michele Invernizzi,
2, 3, 4, † Jan Vorberger, ‡ and Barak Hirshberg
2, 6, § Center for Advanced Systems Understanding (CASUS), D-02826 G¨orlitz, Germany Institute of Computational Sciences, Universit`a della Svizzera italiana, 6900 Lugano, Switzerland National Centre for Computational Design and Discovery of Novel Materials MARVEL,Universit`a della Svizzera italiana, 6900 Lugano, Switzerland Department of Physics, ETH Zurich, 8092 Zurich, Switzerland Helmholtz-Zentrum Dresden-Rossendorf (HZDR), D-01328 Dresden, Germany Department of Chemistry and Applied Biosciences, ETH Zurich, 8092 Zurich, Switzerland
Accurate thermodynamic simulations of correlated fermions using path integral Monte Carlo(PIMC) methods are of paramount importance for many applications such as the description ofultracold atoms, electrons in quantum dots, and warm-dense matter. The main obstacle is thefermion sign problem (FSP), which leads to an exponential increase in computation time bothwith increasing the system-size and with decreasing temperature. Very recently, Hirshberg etal. [
J. Chem. Phys. , 171102 (2020)] have proposed to alleviate the FSP based on the Bo-goliubov inequality. In the present work, we extend this approach by adding a parameter thatcontrols the perturbation, allowing for an extrapolation to the exact result. In this way, we canalso use thermodynamic integration to obtain an improved estimate of the fermionic energy. As atest system, we choose electrons in 2 D and 3 D quantum dots and find in some cases a speed-upexceeding 10 , as compared to standard PIMC, while retaining a relative accuracy of ∼ . I. INTRODUCTION
The accurate estimation of electronic properties is ofparamount importance for many fields such as quantumchemistry, physics, and material science [1]. The mostaccurate results can be obtained using quantum MonteCarlo (QMC) methods which, in principle, allow for a quasi-exact description. Unfortunately, QMC simula-tions of fermions are severely hampered by the notori-ous fermion sign problem (FSP) [2–4], which leads to anexponential increase in the computation time with in-creasing the system-size or decreasing temperature, andhas been shown to be
N P -hard for a specific class ofHamiltonians [3].In the ground-state, the seminal QMC study of theuniform electron gas by Ceperley and Alder [5] has fa-cilitated the success of density functional theory (DFT)regarding the description of real materials [6–8]. Theseresults were obtained on the basis of the fixed-node ap-proximation [9, 10], where the sign problem is avoidedby an a priori decomposition of the wave-function intoa positive and a negative region. Although formally ex-act, the true nodal structure of the wave function is notknown, and one has to rely on approximations. This lim-itation, however, can be alleviated as the ground-stateenergy is variational with respect to the nodes, which canbe exploited for optimization [11–13]. At the same time,there is a broad consensus among the QMC community ∗ [email protected] † [email protected] ‡ [email protected] § [email protected] that the fixed-node approximation has severe limitationsin many cases, and alternative methods [14–17] are de-sirable [18].In addition, a surge of activity has recently emergedin the field of fermionic QMC simulations at finite tem-perature [19–32]. This has been motivated mainly byinterest in warm dense matter (WDM)—an exotic stateat the interface of plasma and solid state physics [33–36]. For example, thermal DFT simulations [37–39] ofWDM require the construction of exchange–correlationfunctionals that explicitly take into account the tem-perature [40, 41], which can be realized on the basis ofQMC data for electrons at these conditions [42, 43]. SeeRef. [35] for a review on recent developments.Other fields for the application of fermionic QMCmethods include dipolar systems such as ultracold atomsor Rydberg dressed states [4, 44], bilayer-systems [45,46], electrons in quantum dots [47–51], and even semi-relativistic quark-gluon plasmas [52, 53]. These sys-tems offer a plethora of interesting effects such as anabnormal superfluid fraction [44, 54], Wigner crystalliza-tion [51, 55], the BCS-BEC transition [56, 57], and col-lective excitations [32, 58–60].Despite this progress, there are still many thermody-namic conditions that are not accessible to QMC meth-ods [29, 61] and their development remains an activetopic of research. In this work, we present an exten-sion of the standard path-integral Monte Carlo (PIMC)method [62] which is motivated by the behavior of thesign for different interaction potentials and is justified bythe well-known Bogoliubov inequality [63]. More specifi-cally, Hirshberg et al. [64] have recently proposed to carryout a path-integral Molecular Dynamics (PIMD) simula-tion [65] of an auxiliary system where the FSP is less a r X i v : . [ phy s i c s . c o m p - ph ] S e p severe and obtain an accurate estimate for the energyof a computationally more challenging system using theBogoliubov inequality. Here, we adapt this idea to thePIMC method. We also extend this approach by addingto the original system a repulsive two-body term thatphenomenologically mimics the effect of the Pauli repul-sion between fermions and allows for a controlled extrap-olation towards the exact result. Moreover, we show thatit is possible to accurately estimate the energy differ-ence between the original and the auxiliary system usingthermodynamic integration [66], which further increasesthe reliability of the method. Our approach results in aspeed-up of up to 10 as compared to standard PIMC,while retaining a relative accuracy of ∼ . D harmonic trap in Sec. III D. Thepaper is concluded by a concise summary and discussion(Sec. IV), where we also indicate possible future direc-tions. II. THEORYA. Path-integral Monte Carlo
The basic idea of the PIMC method [62] is to stochasti-cally sample the thermal density matrix of the canonicalensemble ρ ( R , R (cid:48) , β ) = (cid:104) R | e − β ˆ H | R (cid:48) (cid:105) , (1)where R = ( r , . . . , r N ) T contains the coordinates of all N particles, β = ( k B T ) − is the inverse temperatureand ˆ H denotes the Hamiltonian. The path-integral ex-pression is obtained by performing a Trotter decompo-sition [67], leading to each particle being expressed asan entire path at P discrete positions in imaginary-time τ ∈ [0 , β ]. The collection of the paths of all N particlesis known as a configuration X = ( R , . . . , R P − ) T . Eachconfiguration contributes to the full partition functionaccording to its corresponding weight W ( X ), which is afunction that can be readily evaluated [64], Z = (cid:90) d X W ( X ) . (2) In practice, one uses the metropolis algorithm [68] togenerate a Markov chain of configurations X which aredistributed as P ( X ) = W ( X ) /Z .For indistinguishable particles, one has to explicitlysum over all possible permutations of particle coordi-nates [69]. For bosons, the thermal density matrix is sym-metric under the exchange of particle coordinates, andall terms remain positive. Thus, modern sampling algo-rithms [70, 71] allow for quasi-exact simulations of up to10 particles, which has facilitated profound insights intophenomena such as superfluidity [72–75] and collectiveexcitations [32, 45, 76–79]. Recently, it became possibleto simulate large bosonic systems also using PIMD [65].For fermions, on the other hand, the density matrix isanti-symmetric under particle-exchange, which leads tosign changes in W ( X ) for each pair exchange. Therefore, P = W/Z cannot be interpreted as a probability distri-bution. At this point, one introduces a modified partitionfunction Z (cid:48) = (cid:90) d X | W ( X ) | , (3)where the configurations are generated according to theabsolute value of W ( X ), i.e., P (cid:48) ( X ) = | W ( X ) | /Z (cid:48) . Theexact fermionic expectation value of an observable ˆ A isthen computed as (cid:104) ˆ A (cid:105) = (cid:104) ˆ A ˆ S (cid:105) (cid:48) (cid:104) S (cid:105) (cid:48) , (4)where S ( X ) = W ( X ) / | W ( X ) | denotes the sign associ-ated with a particular configuration, and the denomina-tor of Eq. (4) is the so-called average sign S . S is a measure for the amount of cancellation of posi-tive and negative terms in Z , and exponentially decreasesboth with the system-size N and the inverse temperature β , S = exp ( − βN ( f − f (cid:48) )) , (5)where f and f (cid:48) denote the free energy density of theoriginal and modified system, respectively. Furthermore,the statistical error in estimating the ratio in Eq. (4) isinversely proportional to S ,∆ AA ∼ S √ N MC ∼ exp ( βN ( f − f (cid:48) )) √ N MC . (6)The resulting exponential increase in the Monte Carloerror bar with increasing N or decreasing temperaturecan only be compensated for by increasing the numberof samples N MC . This inevitably becomes unfeasible andone runs into an exponential wall, which is known as thefermion sign problem [2–4]. Methods to overcome theFSP are therefore very desirable. In the following twosections, we describe two approaches for alleviating theFSP in PIMC simulations. B. Extrapolation based on the Bogoliubovinequality
Let ˆ H denote the original Hamiltonian that we wantto simulate using fermionic PIMC,ˆ H = ˆ K + ˆ V ext + ˆ W , (7)with ˆ K , ˆ V ext and ˆ W being the kinetic, external poten-tial, and interaction contribution to the total energy. Wefurther assume that we are interested in the propertiesof this system at relatively low temperature, and thatthe manifestation of quantum degeneracy effects resultsin a low value of the average sign S . In a recent paper,Hirshberg et al. [64] have shown that it is possible toaccurately approximate the energy E ˆ H = (cid:104) ˆ H (cid:105) , by sim-ulating an auxiliary system where ˆ W is replaced by adifferent pair potential ˆ R that more effectively separatesthe particles, ˆ H R = ˆ K + ˆ V ext + ˆ R. (8)This resulted in a substantially less severe manifestationof the sign problem [4], and simulations became feasibleat lower temperatures than for the original system. Then,they used the Bogoliubov inequality [63] F ˆ H − F ˆ H R ≤ (cid:104) ˆ H − ˆ H R (cid:105) ˆ H R , (9)and assumed that the free energy can be approximatedby the energy at low temperatures, to obtain an upperbound on the energy of the original system E ˆ H (cid:46) (cid:104) ˆ H (cid:105) ˆ H R . (10)The subscript indicates that the expectation value of theoriginal Hamiltonian is evaluated in the ensemble of ˆ H R .Since the sign problem is most severe at low temperature[cf. Eq. (6)], this approximation is expected to hold, andthe scheme is highly valuable as R could in principle beoptimized variationally.In the present work, we extend this approach in termsof a coupling parameter η , by re-writing Eq. (8) asˆ H η = ˆ H + η ˆ φ , (11)where ˆ φ is a pair potential that should mimic the effectiverepulsion due to the fermionic degeneracy, as discussed atthe end of this section. Clearly, the PIMC energies com-puted for the Hamiltonian in Eq. (11) are η -dependent,and for any observable ˆ A , it holds that (cid:104) ˆ A (cid:105) = lim η → (cid:104) ˆ A (cid:105) η = lim η → (cid:104) ˆ A ˆ S (cid:105) η (cid:104) ˆ S (cid:105) η . (12)This results in two key advantages: i) The difference be-tween E ˆ H and (cid:104) ˆ H (cid:105) η vanishes as η → η . Due to the added re-pulsion, the simulations converge faster as η is increased.Then, we evaluate E ( η ) ≡ (cid:104) ˆ H (cid:105) η for each one and ex-trapolate it towards η → (cid:104) ˆ A (cid:105) η with respect to η is highly valuable. We show inSec. III B that a simple empirical extrapolation schemeworks well and the results are not very sensitive to therange of η used in the fitting.Throughout this work, we restrict ourselves to two-body correlations and write ˆ φ as the sum over an effectivepair potential, ˆ φ = 12 N (cid:88) k (cid:54) = l Ψ(ˆ r l , ˆ r k ) . (13)Following an observation from Ref. [4] for electrons ina 2 D harmonic confinement, we choose a dipolar short-range repulsion [44],Ψ( r , r ) = 1 | r − r | , (14)which was shown to have a higher average sign in com-parison to Coulomb repulsion between the fermions. Thischoice for electrons in quantum dots guarantees that thefermions feel the Coulomb repulsion at long range, butthat the repulsion due to the dipolar interaction is dom-inant at short range. This additional short range termis what mimics an increased Pauli repulsion and resultsin a larger average sign. We speculate that any short-range potential that is more repulsive than Coulombicinteraction at small separations would be appropriate. C. Thermodynamic integration
As an alternative to direct extrapolation, the free en-ergy difference between the original and auxiliary sys-tem can be evaluated using thermodynamic integration,a widely used free-energy method for atomistic simula-tions [66]. Contrary to the Bogoliubov inequality thatleads to an upper bound [cf. Eq. (9)], using thermody-namic integration we obtain an equality. Given a Hamil-tonian of the form of Eq. (11), the difference in the freeenergy can be estimated as F ˆ H − F ˆ H η = − (cid:90) η d η (cid:48) (cid:104) ˆ φ (cid:105) η (cid:48) . (15)Assuming again that at low temperatures the free en-ergy can be approximated by the energy, we obtain E ˆ H ≈ E ˆ H η − (cid:90) η d η (cid:48) (cid:104) ˆ φ (cid:105) η (cid:48) . (16)The integral in Eq. (16) can be estimated from our PIMCsimulation data at different values of η , whereas the en-tropic contribution, hereafter denoted as ∆ S ( η ), remainsunknown. Still, it is reasonable to assume that the inte-gral term constitutes the dominant contribution at lowtemperatures. This is precisely where the FSP is mostsevere and, consequently, our approach is needed themost. Finally, we notice that also when using thermody-namic integration we need to perform an extrapolationfor η = 0, but this time only of the perturbing potential φ and not of the whole energy as when using the Bogoli-ubov inequality. An extensive discussion of the practicalaspects regarding the application of Eq. (16) is given inin Sec. III C. III. RESULTSA. Model system and speed-up factor
We consider the Hamiltonian of N spin-polarized elec-trons in a harmonic confinement, a commonly employedmodel for quantum dots,ˆ H = − N (cid:88) k =1 ∇ k + 12 N (cid:88) k =1 ˆr k + N (cid:88) k>l λ | ˆr l − ˆr k | , (17)where we assume oscillator units, corresponding to thecharacteristic length l = (cid:112) (cid:126) /m Ω (with Ω being thetrap frequency) and energy scale E = (cid:126) Ω. The firstterm corresponds to the kinetic contribution ˆ K and thelast two terms to the external potential and the Coulombinteraction, ˆ V ext and ˆ W , respectively. The constant λ is the ratio between the screened Coulomb repulsion inthe quantum dot and E [82]. All simulation results inthis work have been obtained for strictly two- and three-dimensional systems.All the PIMC results in this work have been obtainedusing a canonical adaption [83] of the worm algorithmpresented in Refs. [70, 71]. We use a primitive factor-ization of the density matrix (see Refs. [84, 85] for a de-tailed discussion of different factorization schemes) with P ∈ [200 , η = 0,obtained using Eq. 17, as standard PIMC. They will becompared with results obtained with an added dipolarrepulsion, given by Eq. 14, for different values of the cou-pling parameter η . For each of the studied systems, wereport the speed-up factor T ( η ) obtained with respect β =1 E η PIMClinear[0.0:0.5][0.1:0.5][0.2:0.5] β =1.5 E η PIMClinear[0.0:0.5][0.1:0.5][0.2:0.5] β =2 E η PIMClinear[0.1:0.5][0.2:0.5]
FIG. 1. Convergence ( η -dependence) of PIMC results for N = 6 spin-polarized electrons in a 2 D harmonic trap for λ =0 . E ( η ) ≡ (cid:104) ˆ H (cid:105) η for β = 1, β = 1 . β = 2,respectively. The green crosses are PIMC expectation values,the dotted blue and dashed-double-dotted yellow curves fitsaccording to Eq. (19) in different η -ranges, and the dashedblack lines linear fits serving as a guide to the eye. to a standard PIMC calculation. This speed-up factorfollows directly from Eq. (6) and is defined as T ( η ) = (cid:18) S ( η ) S ( η = 0) (cid:19) , (18)where S ( η ) is the average sign obtained with the addi-tional repulsive term in Eq. (11). For example, if S ( η ) isten times larger than S ( η = 0), then we need two ordersof magnitude less Monte-Carlo samples N MC to achievethe same level of statistical uncertainty, i.e., T ( η ) = 100. B. Direct extrapolation scheme
In Fig. 1, we show PIMC results for N = 6 spin-polarized electrons in a 2 D harmonic trap for an inter-mediate value of λ = 0 .
5, and three different tempera-tures, β = 1 (top panel), β = 1 . β = 2 (bottom panel). The green crosses depict thePIMC results for (cid:104) ˆ H (cid:105) η for different values of η . First, wenote that the PIMC data do monotonically converge to-wards the exact η = 0 limit from above, as predicted byEq. (10). For β = 1, the system is substantially out of theground-state, and we find an average sign of S ∼ − for standard PIMC (i.e., η = 0), cf. Fig. 2. Hence, thePIMC simulations can be converged over the entire η -range. The dashed black line corresponds to a linear fitwithin the interval η ∈ [0 , .
5] and has been included asa guide to the eye, although we find it to be surprisinglyaccurate in this case. The dash-dotted red line was ob-tained from a fit within the same interval, but using themodified functional form E ( η ) = a + bη c , (19)which has been found empirically, and with a , b , and c being free parameters. Furthermore, the dotted blueand dash-double-dotted yellow lines have also been ob-tained from fits via Eq. (19), but within the intervals η ∈ [0 . , .
5] and η ∈ [0 . , . η -range. This is an importantintermediate result, as it indicates that Eq. (19) consti-tutes a suitable form to extrapolate results to η = 0 whenPIMC simulations are limited to some finite value of η due to the sign problem.We next examine the center panel, which correspondsto β = 1 .
5, an intermediate temperature. In this case, wefind S ∼ − for standard PIMC, which means that sim-ulations are computationally demanding, but still feasibleover the full η -range. Firstly, we note that the linear fit isnot as good as for β = 1 above, so we have only used datapoints for η ∈ [0 , .
2] in the linear fitting. Secondly, wefind that the blue and red curves are in excellent agree-ment, whereas the yellow curve somewhat deviates in thelimit η →
0. Most probably, this is a consequence of thereduced fit interval of η ∈ [0 . , . η = 0 agrees with the two other curves (with thered curve basically being exact) to a relative accuracy of0 . S ∼ . η = 0 . T ∼ , see Fig. 2.Finally, the bottom panel corresponds to β = 2, where S ∼ − for standard PIMC. Thus, PIMC simulationsare severely hampered by the sign problem and simu-lations are not feasible for η (cid:46) .
1. Still, the dot-ted blue ( η ∈ [0 . , . η ∈ [0 . , . S ( η ) while the bottompanel shows the corresponding speed-up factor T ( η ), asdefined in Eq. (18). S ηβ =1 β =1.5 β =2
0 0.1 0.2 0.3 0.4 0.5 s p ee d - u p T η FIG. 2. Dependence of the average sign S (top panel) andthe speed-up T (bottom panel, cf. Eq. (18)) on the repul-sion strength η for N = 6 spin-polarized electrons in a 2 D harmonic trap for λ = 0 . β = 1, β = 1 .
5, and β = 2,respectively. Note the logarithmic scale of the y -axis in thebottom panel. The green crosses, red circles, and blue diamonds inthe top panel of Fig. 2 show the η -dependence of S for β = 1, β = 1 .
5, and β = 2, respectively, for the same sim-ulations reported in Fig. 1. All three data sets exhibita qualitatively similar progression and monotonically in-crease with η . This growth is more pronounced for thelower temperature, β = 2, where the sign increases bymore than two orders of magnitude at η = 0 . η = 0. Consequently, the corresponding speed-up (Fig. 2, bottom panel) is the highest for β = 2 (dottedblue line), exceeding 10 for η = 0 .
5. From the bottompanel of Fig. 1, it is evident that simulations for η ≥ . ∼ .
1% at this temperature,also provide a speed-up exceeding 10 .For β = 1 . T ∼ for η = 0 . η = 0 was resolved.Finally, the solid green line corresponds to β = 1, wherewe find a speed-up of T ∼ for η (cid:38) .
20 25 30 35 0.5 1 1.5 2 2.5 3 E β PRE 2019 η =0.2 η =0.1 η =0extrapolated
19 20 21 0.5 1 1.5 2 2.5 3 E β -4 -3 -2 -1 S β FIG. 3. Temperature dependence of the energy for N = 6spin-polarized electrons in a 2 D harmonic trap. The top panelshows the total energy E versus β , with the following key:black squares are standard PIMC data from Ref. [4]; greencrosses are standard PIMC from this work; red circles andblue diamonds have been obtained from Eq. (10) for η = 0 . η = 0 .
1, respectively. The center panel shows a mag-nified zoom of the top panel around the lowest temperaturepoints. The bottom panel shows the corresponding data forthe average sign, with the yellow curve depicting an exponen-tial fit according to Eq. (20). Note that the red, blue, andgreen data points have been obtained for the same amountof Monte Carlo samples and, thus, can be directly comparedregarding efficiency. ing the temperature. The results of this analysis areshown in Fig. 3, where the top and bottom panel showthe β -dependence of the total energy E and the averagesign S . The black squares correspond to the standardPIMC data from Ref. [4], and accurate results for E areavailable for β (cid:46) .
3. Data points at lower temperaturespresent very large error bars, and have been omitted forbetter visibility. Looking at S itself, we find a steep decaywhich is of an exponential form for large β , see Eq. (6)and the corresponding analysis in Ref. [4]. The yellowcurve depicts a fit of the form S ( β ) = a S e − βb S , (20)obtained for β ∈ [1 ,
3] and fully confirms this trend.The green crosses in Fig. 3 represent new standard
1 1.5 2 2.5 3 s p ee d - u p T βη =0.2 η =0.1 FIG. 4. Relative speed-up T [cf. Eq. (18)] of our simulationsat η = 0 . η = 0 . η = 0) for the temperature scan fromFig. 3. Note the logarithmic scale of the y -axis. PIMC results for η = 0, but obtained at a substan-tially increased computational cost. Therefore, resultsare available for β (cid:46) .
5, but at β = 2 the statistical un-certainty substantially increases due to the sign problem.We next examine the performance of our new approachbased on the Bogoliubov inequality and the modifiedHamiltonian from Eq. (11). The red circles and bluediamonds represent (cid:104) ˆ H (cid:105) η for η = 0 . η = 0 .
1, re-spectively, and the simulations can be converged downto β = 3. This can be seen particularly well in the cen-ter panel of Fig. 3, showing a magnified segment aroundthe low-temperature points. The yellow diamonds showthe results which have been extrapolated to η = 0 as de-scribed in the discussion of Fig. 1. These results showthat our scheme allows to double the feasible β -range de-spite the exponential wall in compute time given by theFSP.The speed-up T ( η ) [cf. Eq. (18)] is shown in Fig. 4,with the red circles and blue diamonds depicting the β -dependence for T (0 .
2) and T (0 . η . More-over, this increase appears to be of an exponential formfor large β , which helps to explain the remarkable exten-sion of the parameter space that can be covered with thismethod. In particular, we find T ∼ ( T ∼ ) for η = 0 . η = 0 .
1) for the lowest depicted temperature, β = 3.At the same time, it is important to note that thisexponentially growing speed-up is still not sufficient tofully counter the sign problem, since S does still mono-tonically (and, indeed, exponentially) decrease with β for every fixed value of η , cf. the bottom panel of Fig. 3.Therefore, a full solution of the fermion sign problem(which we would define as a simulation scheme withoutan exponential increase in compute time with decreas-ing temperature) would require that the minimum valueof η that is needed for the extrapolation to η → β =3 E η PIMCc- fi tSpline fi t[0.1:0.5]exact β =5 E η PIMCc- fi tSpline fi t[0.2:0.5] Φ η PIMCSpline[0.1:0.5] fi t[0.2:0.5] fi t[0.1:0.5] Φ η PIMCSpline[0.1:0.5] fi t[0.2:0.5] fi t[0.1:0.5] FIG. 5. Results for the thermodynamic integration correction from Sec. II C for N = 4 spin-polarized electrons in a 2 D harmonic trap with λ = 0 . β = 3 (left), β = 5 (right). Top row: energy estimates E . Bottom row: fit of (cid:104) ˆ φ (cid:105) η toevaluate the integral in Eq. (16). The dashed black and dash-dotted yellow lines have been fitted according to Eq. (21), andthe dash-dotted red line corresponds to a spline-fit to η (cid:104) ˆ φ (cid:105) η for η ∈ [0 : 0 .
5] (without the data point at η = 0 . C. Thermodynamic integration scheme
We estimate the energy of the original system by ther-modynamic integration using Eq. (16), as an alternativeto extrapolation. The results are shown in the left col-umn of Fig. 5 for N = 4 spin-polarized electrons in a 2 D quantum dot with λ = 0 . β = 3. The top panelshows data for the total energy E , and the bottom panelthe expectation value φ ( η ) ≡ (cid:104) φ (cid:105) η , the argument of thethermodynamic integration. The green crosses depict thePIMC data for E η = (cid:104) ˆ H (cid:105) η (top) and φ ( η ) (bottom), thatare available at discrete η values. Since the evaluationof Eq. (16) requires the computation of the area under φ ( η (cid:48) ) for η (cid:48) ∈ [0 , η ], two practical obstacles have to beovercome: i) In order to avoid performing many simula-tions, the computation of the integral would benefit froma continuous representation of φ ( η ) by fitting a modestnumber of data points and ii) we need to know φ ( η ) inthe limit of small η , where PIMC simulations might nolonger be feasible due to the sign problem.Overcoming the first problem by itself is relatively easy,and the dotted blue curve in the bottom panel of Fig. 5corresponds to a cubic spline fit to the PIMC data inthe interval η ∈ [0 . , . η → φ is needed for the solution to ii). We find empirically thata suitable choice is φ ( η ) = a φ + b φ η + c φ η / , (21)where a φ , b φ , and c φ are the free parameters. The re- sulting fits are shown in Fig. 5 as the dash-dotted yel-low and dashed black curves, which have been obtainedtaking into account PIMC data for η ∈ [0 . , .
5] and η ∈ [0 . , . φ ( η ) over the entire η -range given as input only fourdata points at η ∈ [0 . , . E η = (cid:104) ˆ H (cid:105) η . For completeness, we also include an ex-trapolation of these data to η = 0 according to Eq. (19,as described in Sec. III B (dash-dotted red line). Thegrey circle corresponds to a standard PIMC simulation( η = 0) that is exact within the given error bars.Using the representation of φ ( η ) according to Eq. (21)to estimate the correction leads to the yellow stars. Asexpected, the entire η -dependence has been removed bythe correction, and the data are in perfect agreementwith the exact standard PIMC results for all depictedvalues of η . In addition, we also show the corrected PIMCdata that have been obtained by using the spline as arepresentation of φ ( η ) instead, see the blue diamonds inthe top panel. It is important to notice that even whenusing such a poor extrapolation scheme, the systematicerror introduced in the energy estimate is only around0 . φ ( η ) will only account for a small contribution to theoverall correction obtained via Eq. (16), while most of itcomes from interpolation of PIMC data. This is not thecase when using instead the direct extrapolation method,where any inaccuracy in the chosen functional form mightmore strongly impact the quality of the final result.Finally, we mention that at the conditions consideredin the left column of Fig. 5, the entropic contributionto Eq. (16) does indeed vanish within the given level ofaccuracy, as expected. This changes only for higher tem-peratures, see the discussion of Fig. 7 below.We have also used this approach to tackle a harderexample, shown in the right column of Fig. 5, where wehave investigated a substantially lower temperature, β =5, for which S (cid:46) − . Therefore, standard PIMC is notavailable in this case, and PIMC simulations are onlyfeasible for η (cid:38) . φ ( η ) and,also in this case, the fits from Eq. (21) are indistinguish-able for the two different intervals of input data, whichsubstantiates the high quality of this representation. Thespline, on the contrary, significantly deviates at low η val-ues.The corresponding energies are shown in the top rightpanel of Fig. 5, and the dash-dotted red line depicts thedirect extrapolation of the PIMC data to η → φ ( η ) only plays aminor role for the overall level of accuracy of the energy.Secondly, the corrected energies do not exhibit any resid-ual dependence on η and fluctuate around the horizontalred line, with an uncertainty level of 0.1%.The increase of the average sign and the correspondingspeed-up for both β = 5 and β = 3 is shown in Fig. 6.For the higher temperature (green crosses), we observean increase in S (top panel) by two order of magnitudebetween η = 0 and η = 0 .
5, which results in a speed-upof up to T ∼ (bottom panel). For β = 5 (red circles),the relative gain in the sign is even larger, leading toa speed-up (dash-dotted red) exceeding T ∼ at thelargest value of η .A realistic application for our method performs simu-lations down to η (cid:38) .
2, as this does suffice for a quasi-exact extrapolation to η →
0. This boundary is markedas the vertical grey line in Fig. 6, and the two horizon-tal arrows point to the respective speed-up on the y -axis.For β = 3, standard PIMC simulations are, in principle,possible, but our scheme results in a speed-up by a factorof T ∼ . For β = 5, standard PIMC is not feasible,and it is only our speed-up by a factor of T ∼ thatmakes it possible to obtain accurate data.For completeness, we also examine the application ofthe correction from Eq. (16) for higher temperatures, -6 -5 -4 -3 -2 -1
0 0.1 0.2 0.3 0.4 0.5 S ηβ =3 β =5
0 0.1 0.2 0.3 0.4 0.5 s p ee d - u p T η FIG. 6. PIMC results for the average sign S (top panel) andthe respective speed-up T (bottom panel) for N = 4 spin-polarized electrons in a 2 D harmonic trap with λ = 0 .
5. Thesolid green and dash-dotted red curves correspond to β = 3and β = 5, respectively. The arrows point to the speed-up for η = 0 . which is shown in Fig. 7 for β = 2 (left column) and β = 1(right column). The sign problem is not severe at theseparameters and we find S ≈ .
02 ( S ≈ .
25) for β = 2( β = 1). Thus, PIMC simulations are computationallyfeasible over the entire η -range for both cases. Still, wefind that the functional form from Eq. (21) allows foran accurate representation of φ ( η ) (see the dashed blackcurve in the bottom row) taking only into account thefour data points at η = 0 . , . , . , . (cid:104) ˆ H (cid:105) η , and the red curves have been obtainedfrom a direct fit to these data according to Eq. (19).Finally, the black squares have been obtained from thethermodynamic integration correction Eq. (16). In con-trast to the previous results in Fig. 7, we find a significantdependence of the corrected data points on η , which isdue to the entropic contribution ∆ S ( η ) to Eq. (16).Interestingly, this function can be perfectly reproducedby a linear fit, ∆ S ( η ) = a S + b S η , (22)and the results are depicted by the dotted blue lines inFig. 7.For β = 2, the entropic contribution is quite small and β =2 E η PIMCc- fi t fi t[0.2,0.5] β =1 E η PIMCc- fi t fi t[0.2,0.5] Φ η PIMC fi t[0.2:0.5] Φ η PIMC fi t[0.2:0.5] FIG. 7. Results for the thermodynamic integration correction from Sec. II C for N = 4 spin-polarized electrons in a 2 D harmonic trap with λ = 0 . β = 2 (left), β = 1 (right). Top row: energy estimates E . Bottom row: fit of (cid:104) ˆ φ (cid:105) η toevaluate the integral in Eq. (16). The dashed black and dash-dotted yellow lines have been fitted according to Eq. (21), andthe dash-dotted red line corresponds to a spline-fit to η (cid:104) ˆ φ (cid:105) η for η ∈ [0 : 0 .
5] (without the data point at η = 0 . does not exceed 0 .
3% even for η = 0 .
5. In contrast, ∆ S attains a maximum value of ∼ .
7% for β = 1. To putthese findings into the proper context, we find it usefulto briefly recall the following points: i) while we cannotdirectly estimate ∆ S ( η ) from our PIMC results, the in-fluence on the correction from Eq. (16) decreases for lowtemperature, when the sign problem is most severe; ii)even when ∆ S ( η ) does have an impact on the correctedenergies, the residual dependence on the parameter η ismuch smaller than the direct dependence of (cid:104) ˆ H (cid:105) η , whichmakes a potential extrapolation to η → S on η , which is an additional advantageover the direct extrapolation of the uncorrected energy,where the functional dependence is more complicated,cf. Sec. III B.We thus conclude that the correction introduced inSec. II C using thermodynamic integration constitutes adistinct improvement over the direct extrapolation ex-plored in Sec. III B.We conclude this section with an application of the cor-rection approach at higher λ and lower temperature. InFig. 8, we show PIMC results for the energy (top panel)and φ (bottom panel) for N = 4 spin-polarized electronsat λ = 2 and β = 10. First, we mention that standardPIMC is not available at these conditions, and the signvanishes within a statistical uncertainty of 10 − . We havechosen this particular set of parameters because it waspreviously studied by Egger et al. [49] using the approx-imate multi-level blocking (MLB) method [50, 86, 87].While being potentially biased [88], such a data pointstill constitutes a valuable reference for the development of a new method.As usual, the PIMC results for (cid:104) ˆ H (cid:105) η are depicted asthe green crosses, and the dash-dotted red curve corre-sponds to a direct extrapolation of these data accord-ing to Eq. (19). The solid grey horizontal line depictsthe MLB value from Ref. [49], and the two light dot-ted grey lines depict the corresponding statistical uncer-tainty, which is given by ∆ E/E ≈ . φ ( η ), shownin the bottom panel of Fig. 8. The dash-double-dottedyellow curve has been obtain using Eq. (21) as a fit func-tion for the PIMC data, which was empirically shown tobe accurate over the entire relevant η -range, see above.Using this representation to evaluate Eq. (16) results inthe yellow stars in the top panel.Firstly, we note that these data do not exhibit anyresidual dependence on η , as the entropic contribution isnegligible at such a low temperature. Moreover, we findthat the yellow stars are in excellent agreement to theresult from the direct extrapolation of the PIMC data,and, thus, also in agreement within the given uncertaintyinterval of the MLB results.Therefore, both the MLB method and our new ap-proach have been successfully validated against eachother for this system.0 E η PIMCMLBc- fi t fi t[0.5:3.0] Φ η FIG. 8. Results for the thermodynamic integration correc-tion from Sec. II C for N = 4 spin-polarized electrons in a 2 D harmonic trap with λ = 2 and β = 10. Top panel: energyestimates E . The solid grey line (dotted grey line) corre-sponds to the MLB result (error bar) from Ref. [49]. Bottompanel: fit of (cid:104) ˆ φ (cid:105) η to evaluate the integral in Eq. (16). Thedash-double-dotted yellow line has been fitted according toEq. (21). D. Electrons in a 3D quantum dot
The final example to be investigated in this work isthe application of our method to a 3 D system. This isshown in Fig. 9 for N = 6 electrons in a 3 D harmonicconfinement for λ = 0 . β = 2. As before, the greencrosses depict the raw PIMC data for (cid:104) ˆ H (cid:105) η and the dash-dotted red-line a direct extrapolation thereof accordingto Eq. (19). At these conditions, we find an average signof S ∼ − [cf. Fig. 10], which means that standardPIMC simulations are only feasible for η (cid:38) .
1. Yet, thefit function from Eq. (19) evidently allows for a controlledextrapolation to η →
0, which is further highlighted bythe dotted light grey horizontal line, depicting an uncer-tainty interval of 0 .
3% around the extrapolated value.We next explore the estimation of the thermodynamicintegration correction introduced in Sec. II C. To thisend, we show φ ( η ) in the bottom panel of Fig. 9 andthe green crosses again show the PIMC data. Thedashed black (dashed-double-dotted yellow) line corre-sponds to a fit using the functional form from Eq. (21)for η ∈ [0 . , .
5] ( η ∈ [0 . , . η , whereassome deviations appear for the extrapolation of η → η ∈ [0 . , . E η PIMCc- fi t fi t[0.1:0.5] Φ η PIMC fi t[0.2:0.5] fi t[0.1:0.5] FIG. 9. Top: Convergence ( η -dependence) of PIMC resultsfor N = 6 spin-polarized electrons in a 3 D harmonic con-finement with λ = 0 . β = 2. The green crosses depictPIMC data for E η = (cid:104) ˆ H (cid:105) η and the dash-dotted red line a fitaccording to Eq. (19). The yellow stars have been obtained us-ing the thermodynamic integration correction from Eq. (16),and the dashed-double-dotted yellow line an extrapolation ofthe entropic contribution, cf. Eq. (22). Bottom: PIMC datafor φ (green crosses) and fits via Eq. (21) (dashed black anddashed-double-dotted yellow) for two different η -ranges. representation of φ ( η ) would result in almost indistin-guishable energy values, which further validates the ob-servation from the previous Sec. III C that the partic-ular extrapolation of φ ( η ) to η → η -range, and the η → S and the corresponding speed-up T [cf. Eq. (18)], shown in Fig. 10. The green crosses showthe PIMC data for S (left y -axis) and exhibit a monotonicincrease similar to the 2 D case. Further, the respectivespeed-up is depicted as the solid black curve (right y -axis)and attains values exceeding T = 10 for η = 0 . IV. SUMMARY AND DISCUSSION
We present a new approach for PIMC simulations offermions at finite temperature, that extends and im-1 -4 -3 -2 -1
0 0.1 0.2 0.3 0.4 0.510 S s p ee d - u p η Signspeed-up
FIG. 10. Dependence of the average sign S (green crosses,left y -axis) and the speed-up T (cf. Eq. (18), solid black line,right y -axis) on the repulsion strength η for N = 6 spin-polarized electrons in a 3 D harmonic trap for λ = 0 . β = 2. proves the idea by Hirshberg et al. [64], that was based onthe Bogoliubov inequality. We add to the Hamiltonian arepulsive term that mimics Pauli repulsion at short rangeand is proportional to a coupling parameter η . This sig-nificantly improves the efficiency of the simulations, byincreasing the average sign S . We then propose two sim-ple post-processing schemes to recover the energy of theoriginal system. The first one is based on the Bogoliubovinequality and consists of extrapolating (cid:104) ˆ H (cid:105) η obtainedfor various values of η to the limit η →
0. The secondone instead is based on thermodynamic integration andalso relies on an extrapolation, but only of the perturba-tion term φ and not of the whole energy. Combined withthe fact that it is based on an exact relation and not aninequality, this makes the second scheme generally morereliable. We believe that having two distinct schemesfor evaluating the energy, starting from the same simula-tions, is a great advantage, and makes the method morerobust against possible poor extrapolation choices. Forall the systems considered here, the two schemes pro-vided compatible energy estimates. Most importantly,they allowed evaluating accurate estimates of energiesfor conditions in which standard PIMC simulations werenot feasible.As a practical application, we have investigated elec-trons both in 2 D and 3 D quantum dots, and our methodworks very well in both of these cases. In particular, wehave found that a direct extrapolation by itself allowsfor a speed-up of up to 10 (cf. Fig. 4), while retaininga relative accuracy of 0 . β -range as compared to standard PIMC.The investigation of the thermodynamic integrationcorrection introduced in Sec. II C has revealed that thisapproach makes the estimation of the energy of the orig-inal system even more reliable. Here, the main challengeis given by the construction of a representation of φ ( η ) with no data points being available below some minimumvalue. As a solution, we have introduced a suitable em-pirical fit function [cf. Eq. (21)] that, remarkably, allowsfor a highly reliable extrapolation using only a few datapoints at large η as input. At low temperature, this rep-resentation allows to accurately estimate the energy ofthe original system of interest. For higher temperatures,there emerges an additional entropic contribution, whichremains a priori unknown. Still, we stress that the latteronly constitutes a fraction of the full difference betweenthe modified and the original system, and, empirically,exhibits a linear dependence on η (cf. Fig. 7).While the present results are certainly encouraging,much can be done to improve the method further. Froma technical point of view, implementing the sampling pro-cedure proposed in Ref. 89 would make the method moreefficient, requiring only a single simulation for estimating E ( η ) and φ ( η ) in the whole range of η values. Further-more, we mention that, despite a speed-up of 10 , thefermion sign problem has not been completely removedand, while our approach significantly extends the rangeof accessible parameters, it still suffers from exponentialincrease in computing time at low enough temperatures.We have also restricted ourselves to the investigation ofenergy values, and the adaption of the method to otherobservables such as pair distributions, structure factors,or the different contributions to E like the kinetic or ex-ternal potential energies is of high interest. Similarly, wehave solely used a dipolar repulsive potential Ψ( r , r )and the optimization of Ψ might substantially improvethe approach.A particularly interesting topic for future research isgiven by the application to other systems, with warm-dense matter [34] being a promising candidate. Finally,we stress that our approach is quite general and can read-ily be adapted to other simulation methods for fermionsboth in the ground state and at finite temperate. ACKNOWLEDGMENTS
This work was partly funded by the Center of Ad-vanced Systems Understanding (CASUS) which is fi-nanced by Germany’s Federal Ministry of Education andResearch (BMBF) and by the Saxon Ministry for Science,Art, and Tourism (SMWK) with tax funds on the basis ofthe budget approved by the Saxon State Parliament. M.I.acknowledges support from the Swiss National ScienceFoundation through the NCCR MARVEL. The PIMCcalculations were carried out at the Norddeutscher Ver-bund f¨ur Hoch- und H¨ochstleistungsrechnen (HLRN) un-der grant shp00015 and on a Bull Cluster at the Centerfor Information Services and High Performace Comput-ing (ZIH) at Technische Universit¨at Dresden. [1] G. Giuliani and G. Vignale,
Quantum Theory of the Elec-tron Liquid (Cambridge University Press, Cambridge, 2008). [2] E. Y. Loh, J. E. Gubernatis, R. T. Scalettar, S. R. White,D. J. Scalapino, and R. L. Sugar, “Sign problem in thenumerical simulation of many-electron systems,” Phys.Rev. B , 9301–9307 (1990).[3] M. Troyer and U. J. Wiese, “Computational complex-ity and fundamental limitations to fermionic quantumMonte Carlo simulations,” Phys. Rev. Lett , 170201(2005).[4] T. Dornheim, “Fermion sign problem in path inte-gral Monte Carlo simulations: Quantum dots, ultracoldatoms, and warm dense matter,” Phys. Rev. E ,023307 (2019).[5] D. M. Ceperley and B. J. Alder, “Ground state of theelectron gas by a stochastic method,” Phys. Rev. Lett. , 566–569 (1980).[6] J. P. Perdew and Alex Zunger, “Self-interaction cor-rection to density-functional approximations for many-electron systems,” Phys. Rev. B , 5048–5079 (1981).[7] John P Perdew, Kieron Burke, and Matthias Ernzer-hof, “Generalized gradient approximation made simple,”Physical Review Letters , 3865 (1996).[8] Kieron Burke, “Perspective on density functional the-ory,” The Journal of Chemical Physics , 150901(2012), https://doi.org/10.1063/1.4704546.[9] D. M. Ceperley, “Fermion nodes,” Journal of StatisticalPhysics , 1237–1267 (1991).[10] James B. Anderson, “Fixed-node quantum monte carlo,”International Reviews in Physical Chemistry , 85–112(1995), https://doi.org/10.1080/01442359509353305.[11] W. M. C. Foulkes, L. Mitas, R. J. Needs, and G. Ra-jagopal, “Quantum monte carlo simulations of solids,”Rev. Mod. Phys. , 33–83 (2001).[12] P. L´opez R´ıos, A. Ma, N. D. Drummond, M. D. Towler,and R. J. Needs, “Inhomogeneous backflow transforma-tions in quantum monte carlo calculations,” Phys. Rev.E , 066701 (2006).[13] R J Needs, M D Towler, N D Drummond, and P L´opezR´ıos, “Continuum variational and diffusion quantummonte carlo calculations,” Journal of Physics: CondensedMatter , 023201 (2009).[14] George H. Booth, Alex J. W. Thom, and Ali Alavi,“Fermion monte carlo without fixed nodes: A game oflife, death, and annihilation in slater determinant space,”The Journal of Chemical Physics , 054106 (2009),https://aip.scitation.org/doi/pdf/10.1063/1.3193710.[15] Michio Honma, Takahiro Mizusaki, and Takaharu Ot-suka, “Diagonalization of hamiltonians for many-bodysystems by auxiliary field quantum monte carlo tech-nique,” Phys. Rev. Lett. , 1284–1287 (1995).[16] Mario Motta and Shiwei Zhang, “Ab initio com-putations of molecular systems by the auxiliary-field quantum monte carlo method,” WIREs Com-putational Molecular Science , e1364 (2018),https://onlinelibrary.wiley.com/doi/pdf/10.1002/wcms.1364.[17] J. P. F. LeBlanc, Andrey E. Antipov, Federico Becca,Ireneusz W. Bulik, Garnet Kin-Lic Chan, Chia-MinChung, Youjin Deng, Michel Ferrero, Thomas M. Hen-derson, Carlos A. Jim´enez-Hoyos, E. Kozik, Xuan-WenLiu, Andrew J. Millis, N. V. Prokof’ev, Mingpu Qin, Gus-tavo E. Scuseria, Hao Shi, B. V. Svistunov, Luca F. Toc-chio, I. S. Tupitsyn, Steven R. White, Shiwei Zhang, Bo-Xiao Zheng, Zhenyue Zhu, and Emanuel Gull (SimonsCollaboration on the Many-Electron Problem), “Solu-tions of the two-dimensional hubbard model: Bench- marks and results from a wide range of numerical al-gorithms,” Phys. Rev. X , 041041 (2015).[18] George H. Booth, Andreas Gr¨uneis, Georg Kresse,and Ali Alavi, “Towards an exact description of elec-tronic wavefunctions in real solids,” Nature , 365–370(2013).[19] Ethan W. Brown, Bryan K. Clark, Jonathan L. DuBois,and David M. Ceperley, “Path-integral monte carlo sim-ulation of the warm dense homogeneous electron gas,”Phys. Rev. Lett. , 146405 (2013).[20] T. Schoof, M. Bonitz, A. Filinov, D. Hochstuhl, andJ.W. Dufty, “Configuration path integral monte carlo,”Contributions to Plasma Physics , 687–697 (2011).[21] Tobias Dornheim, Simon Groth, Alexey Filinov, andMichael Bonitz, “Permutation blocking path integralmonte carlo: a highly efficient approach to the simulationof strongly degenerate non-ideal fermions,” New Journalof Physics , 073017 (2015).[22] Tobias Dornheim, Simon Groth, Fionn D. Malone, TimSchoof, Travis Sjostrom, W. M. C. Foulkes, and MichaelBonitz, “Ab initio quantum monte carlo simulation of thewarm dense electron gas,” Physics of Plasmas , 056303(2017), https://doi.org/10.1063/1.4977920.[23] N. S. Blunt, T. W. Rogers, J. S. Spencer, and W. M. C.Foulkes, “Density-matrix quantum monte carlo method,”Phys. Rev. B , 245124 (2014).[24] Yuan Liu, Minsik Cho, and Brenda Rubenstein, “Abinitio finite temperature auxiliary field quantum montecarlo,” Journal of Chemical Theory and Computation ,4722–4732 (2018).[25] Fionn D. Malone, N. S. Blunt, James J. Shepherd,D. K. K. Lee, J. S. Spencer, and W. M. C. Foulkes, “In-teraction picture density matrix quantum monte carlo,”The Journal of Chemical Physics , 044116 (2015),https://doi.org/10.1063/1.4927434.[26] Fionn D. Malone, N. S. Blunt, Ethan W. Brown, D. K. K.Lee, J. S. Spencer, W. M. C. Foulkes, and James J.Shepherd, “Accurate exchange-correlation energies forthe warm dense electron gas,” Phys. Rev. Lett. ,115701 (2016).[27] Jahan Claes and Bryan K. Clark, “Finite-temperatureproperties of strongly correlated systems via variationalmonte carlo,” Phys. Rev. B , 205109 (2017).[28] Tobias Dornheim, Simon Groth, and Michael Bonitz,“Permutation blocking path integral monte carlo simu-lations of degenerate electrons at finite temperature,”Contributions to Plasma Physics , e201800157 (2019),https://onlinelibrary.wiley.com/doi/pdf/10.1002/ctpp.201800157.[29] A. Yilmaz, K. Hunger, T. Dornheim, S. Groth, andM. Bonitz, “Restricted configuration path integral montecarlo,” (2020), arXiv:2007.12498 [physics.comp-ph].[30] Tobias Dornheim, Jan Vorberger, and Michael Bonitz,“Nonlinear electronic density response in warm densematter,” Phys. Rev. Lett. , 085001 (2020).[31] K. P. Driver, F. Soubiran, and B. Militzer, “Path integralmonte carlo simulations of warm dense aluminum,” Phys.Rev. E , 063207 (2018).[32] T. Dornheim, S. Groth, J. Vorberger, and M. Bonitz,“Ab initio path integral Monte Carlo results for the dy-namic structure factor of correlated electrons: From theelectron liquid to warm dense matter,” Phys. Rev. Lett. , 255001 (2018).[33] F. Graziani, M. P. Desjarlais, R. Redmer, and S. B.Trickey, eds., Frontiers and Challenges in Warm Dense Matter (Springer, International Publishing, 2014).[34] M. Bonitz, T. Dornheim, Zh. A. Moldabekov, S. Zhang,P. Hamann, H. K¨ahlert, A. Filinov, K. Ramakrishna,and J. Vorberger, “Ab initio simulation of warm densematter,” Physics of Plasmas , 042710 (2020).[35] T. Dornheim, S. Groth, and M. Bonitz, “The uniformelectron gas at warm dense matter conditions,” Phys.Reports , 1–86 (2018).[36] V. E. Fortov, “Extreme states of matter on earth and inspace,” Phys.-Usp , 615–647 (2009).[37] Aurora Pribram-Jones, Stefano Pittalis, E. K. U. Gross,and Kieron Burke, “Thermal Density Functional Theoryin Context,” (2014) pp. 25–60.[38] Justin C. Smith, Francisca Sagredo, and Kieron Burke,“Warming Up Density Functional Theory,” in Frontiersof Quantum Chemistry (Springer Singapore, Singapore,2018) pp. 249–271.[39] Yael Cytter, Eran Rabani, Daniel Neuhauser, and RoiBaer, “Stochastic density functional theory at finite tem-peratures,” Phys. Rev. B , 115207 (2018).[40] V. V. Karasiev, L. Calderin, and S. B. Trickey, “Im-portance of finite-temperature exchange correlation forwarm dense matter calculations,” Phys. Rev. E ,063207 (2016).[41] Kushal Ramakrishna, Tobias Dornheim, and JanVorberger, “Influence of finite temperature exchange-correlation effects in hydrogen,” Phys. Rev. B ,195129 (2020).[42] T. Dornheim, S. Groth, T. Sjostrom, F. D. Malone,W. M. C. Foulkes, and M. Bonitz, “Ab initio quantumMonte Carlo simulation of the warm dense electron gas inthe thermodynamic limit,” Phys. Rev. Lett. , 156403(2016).[43] S. Groth, T. Dornheim, T. Sjostrom, F. D. Malone,W. M. C. Foulkes, and M. Bonitz, “Ab initio exchange–correlation free energy of the uniform electron gas atwarm dense matter conditions,” Phys. Rev. Lett. ,135001 (2017).[44] Tobias Dornheim, “Path-integral monte carlo simulationsof quantum dipole systems in traps: Superfluidity, quan-tum statistics, and structural properties,” Phys. Rev. A , 023307 (2020).[45] A. Filinov, “Correlation effects and collective excitationsin bosonic bilayers: Role of quantum statistics, superflu-idity, and the dimerization transition,” Phys. Rev. A ,013603 (2016).[46] J. Schleede, A. Filinov, M. Bonitz, and H. Fehske,“Phase diagram of bilayer electron-hole plasmas,”Contributions to Plasma Physics , 819–826 (2012),https://onlinelibrary.wiley.com/doi/pdf/10.1002/ctpp.201200045.[47] T. Dornheim, H. Thomsen, P. Ludwig, A. Fil-inov, and M. Bonitz, “Analyzing quan-tum correlations made simple,” Contribu-tions to Plasma Physics , 371–379 (2016),https://onlinelibrary.wiley.com/doi/pdf/10.1002/ctpp.201500120.[48] Ilkka Kyl¨anp¨a¨a and Esa R¨as¨anen, “Path integral montecarlo benchmarks for two-dimensional quantum dots,”Phys. Rev. B , 205445 (2017).[49] R. Egger, W. H¨ausler, C. H. Mak, and H. Grabert,“Crossover from fermi liquid to wigner molecule behav-ior in quantum dots,” Phys. Rev. Lett. , 3320–3323(1999).[50] R. EGGER and C. H. MAK, “Multilevel blockingmonte carlo simulations for quantum dots,” International Journal of Modern Physics B , 1416–1425 (2001),https://doi.org/10.1142/S021797920100591X.[51] A. V. Filinov, M. Bonitz, and Yu. E. Lozovik, “Wignercrystallization in mesoscopic 2d electron systems,” Phys.Rev. Lett. , 3851–3854 (2001).[52] V. S. Filinov, Yu. B. Ivanov, V. E. Fortov, M. Bonitz,and P. R. Levashov, “Color path-integral monte-carlosimulations of quark-gluon plasma: Thermodynamic andtransport properties,” Phys. Rev. C , 035207 (2013).[53] V.S. Filinov, M. Bonitz, Y.B. Ivanov, E.-M. Il-genfritz, and V.E. Fortov, “Thermodynamics ofthe quark-gluon plasma at finite chemical poten-tial: Color path integral monte carlo results,” Con-tributions to Plasma Physics , 203–208 (2015),https://onlinelibrary.wiley.com/doi/pdf/10.1002/ctpp.201400056.[54] Yangqian Yan and D. Blume, “Abnormal superfluid frac-tion of harmonically trapped few-fermion systems,” Phys.Rev. Lett. , 235301 (2014).[55] A.V. Filinov, Yu.E. Lozovik, and M. Bonitz, “Path in-tegral simulations of crystallization of quantum confinedelectrons,” physica status solidi (b) , 231–234 (2000).[56] B. Zenker, D. Ihle, F. X. Bronold, and H. Fehske,“Electron-hole pair condensation at the semimetal-semiconductor transition: A bcs-bec crossover scenario,”Phys. Rev. B , 121102 (2012).[57] Y. Ohashi and A. Griffin, “Bcs-bec crossover in a gas offermi atoms with a feshbach resonance,” Phys. Rev. Lett. , 130402 (2002).[58] S. Groth, T. Dornheim, and J. Vorberger, “Ab initiopath integral Monte Carlo approach to the static anddynamic density response of the uniform electron gas,”Phys. Rev. B , 235122 (2019).[59] Tobias Dornheim, Zhandos A Moldabekov, Jan Vor-berger, and Simon Groth, “Ab initio path integral montecarlo simulation of the uniform electron gas in the highenergy density regime,” Plasma Physics and ControlledFusion , 075003 (2020).[60] Paul Hamann, Tobias Dornheim, Jan Vorberger, Zhan-dos A. Moldabekov, and Michael Bonitz, “Dynamicproperties of the warm dense electron gas: an abinitio path integral monte carlo approach,” (2020),arXiv:2007.15471 [physics.comp-ph].[61] V. V. Karasiev, S. B. Trickey, and J. W. Dufty, “Status offree-energy representations for the homogeneous electrongas,” Phys. Rev. B , 195134 (2019).[62] D. M. Ceperley, “Path integrals in the theory of con-densed helium,” Rev. Mod. Phys , 279 (1995).[63] A.F. Verbeure, Many-Body Boson Systems: Half aCentury Later , Theoretical and Mathematical Physics(Springer London, 2010).[64] Barak Hirshberg, Michele Invernizzi, and Michele Par-rinello, “Path integral molecular dynamics for fermions:Alleviating the sign problem with the bogoliubov in-equality,” The Journal of Chemical Physics , 171102(2020), https://doi.org/10.1063/5.0008720.[65] Barak Hirshberg, Valerio Rizzi, and MicheleParrinello, “Path integral molecular dynam-ics for bosons,” Proceedings of the NationalAcademy of Sciences