The Fermion Sign Problem in Path Integral Monte Carlo Simulations: Quantum Dots, Ultracold Atoms, and Warm Dense Matter

TL;DR

This paper discusses the fermion sign problem in path integral Monte Carlo simulations, analyzing its effects across various systems and conditions, and provides extensive data to aid future method development and benchmarking.

Contribution

It offers a detailed analysis of the fermion sign problem's manifestation and provides comprehensive PIMC data for different fermionic systems and regimes.

Findings

FSP severity varies with temperature, system size, and interaction type.

Fermionic expectation values can have non-Gaussian, fat-tailed distributions.

Extensive PIMC data serve as benchmarks for future research.

Abstract

The ab initio thermodynamic simulation of correlated Fermi systems is of central importance for many applications, such as warm dense matter, electrons in quantum dots, and ultracold atoms. Unfortunately, path integral Monte Carlo (PIMC) simulations of fermions are severely restricted by the notorious fermion sign problem (FSP). In this work, we present a hands-on discussion of the FSP and investigate in detail its manifestation with respect to temperature, system size, interaction-strength and -type, and the dimensionality of the system. Moreover, we analyze the probability distribution of fermionic expectation values, which can be non-Gaussian and fat-tailed when the FSP is severe. As a practical application, we consider electrons and dipolar atoms in a harmonic confinement, and the uniform electron gas in the warm dense matter regime. In addition, we provide extensive PIMC data, which can be used as a reference for the development of new methods and as a benchmark for approximations.

Tables

Table 1. Table 1: System-size dependence: Ab initio path integral Monte Carlo results for spin-polarized electrons (Coulomb interaction) at β = 1 𝛽 1 \beta=1 (top half) and β = 0.3 𝛽 0.3 \beta=0.3 (bottom half) in a 2 ​ D 2 𝐷 2D harmonic trap with coupling strength λ = 0.5 𝜆 0.5 \lambda=0.5 . All results have been obtained for P = 200 𝑃 200 P=200 imaginary-time propagators (see Appendix  A for details) and are given in oscillator units (i.e., energies in units of E 0 = ℏ ​ Ω subscript 𝐸 0 Planck-constant-over-2-pi Ω E_{0}=\hbar\Omega ). Parts of the data set are shown in Fig.  7 .

$N$	$S$	$E_{\text{HO}}$	$V$	$K$	$K_{\text{vir}}$	$E$
$3$	$0.4746\left(2\right)$	$4.1585\left(9\right)$	$4.9640\left(8\right)$	$3.755\left(3\right)$	$3.756\left(1\right)$	$8.719\left(3\right)$
$4$	$0.2453\left(1\right)$	$6.0593\left(9\right)$	$7.6279\left(8\right)$	$5.275\left(7\right)$	$5.275\left(1\right)$	$12.903\left(7\right)$
$5$	$0.1085\left(1\right)$	$8.168\left(3\right)$	$10.721\left(2\right)$	$6.94\left(2\right)$	$6.892\left(4\right)$	$17.66\left(2\right)$
$6$	$0.04184\left(8\right)$	$10.483\left(8\right)$	$14.234\left(5\right)$	$8.59\left(5\right)$	$8.61\left(1\right)$	$22.82\left(5\right)$
$7$	$0.01425\left(9\right)$	$12.98\left(3\right)$	$18.13\left(1\right)$	$10.6\left(2\right)$	$10.41\left(4\right)$	$28.7\left(2\right)$
$8$	$0.00428\left(6\right)$	$15.68\left(7\right)$	$22.40\left(3\right)$	$12.1\left(5\right)$	$12.3\left(1\right)$	$34.5\left(5\right)$
$9$	$0.00130\left(8\right)$	$17.9\left(3\right)$	$26.7\left(2\right)$	$13\left(2\right)$	$13.5\left(5\right)$	$40\left(2\right)$
$10$	$0.00030\left(3\right)$	$21.4\left(6\right)$	$31.8\left(3\right)$	$17\left(3\right)$	$16.2\left(10\right)$	$49\left(3\right)$
$4$	$0.8413\left(1\right)$	$14.392\left(5\right)$	$15.592\left(5\right)$	$13.782\left(6\right)$	$13.792\left(9\right)$	$29.374\left(8\right)$
$6$	$0.6580\left(2\right)$	$22.492\left(8\right)$	$25.438\left(8\right)$	$21.01\left(1\right)$	$21.02\left(1\right)$	$46.45\left(1\right)$
$8$	$0.4679\left(2\right)$	$31.118\left(10\right)$	$36.527\left(9\right)$	$28.44\left(2\right)$	$28.41\left(2\right)$	$64.97\left(2\right)$
$9$	$0.3816\left(3\right)$	$35.65\left(2\right)$	$42.55\left(1\right)$	$32.24\left(2\right)$	$32.20\left(2\right)$	$74.78\left(3\right)$
$10$	$0.3051\left(3\right)$	$40.32\left(2\right)$	$48.88\left(2\right)$	$36.04\left(4\right)$	$36.04\left(3\right)$	$84.92\left(4\right)$
$11$	$0.2384\left(3\right)$	$45.06\left(2\right)$	$55.45\left(2\right)$	$39.86\left(5\right)$	$39.87\left(4\right)$	$95.31\left(6\right)$
$12$	$0.1833\left(1\right)$	$49.98\left(1\right)$	$62.35\left(1\right)$	$43.85\left(5\right)$	$43.80\left(2\right)$	$106.20\left(6\right)$
$14$	$0.1021\left(2\right)$	$60.11\left(3\right)$	$76.95\left(3\right)$	$51.7\left(2\right)$	$51.70\left(5\right)$	$128.7\left(2\right)$
$16$	$0.0530\left(2\right)$	$70.75\left(6\right)$	$92.67\left(5\right)$	$59.7\left(3\right)$	$59.79\left(9\right)$	$152.4\left(3\right)$
$18$	$0.0259\left(1\right)$	$81.89\left(8\right)$	$109.47\left(6\right)$	$69.6\left(6\right)$	$68.1\left(1\right)$	$179.1\left(6\right)$
$20$	$0.01164\left(6\right)$	$93.3\left(1\right)$	$127.18\left(8\right)$	$75.8\left(10\right)$	$76.4\left(2\right)$	$203.0\left(10\right)$

Table 2. Table 2: Coupling strength dependence: Ab initio path integral Monte Carlo results for N = 6 𝑁 6 N=6 spin-polarized electrons (top half) and ultracold atoms with dipole interaction (bottom half) in a 2 ​ D 2 𝐷 2D harmonic trap at an inverse temperature β = 1 𝛽 1 \beta=1 . All results have been obtained for P = 200 𝑃 200 P=200 imaginary-time propagators (see Appendix  A for details) and are given in oscillator units (i.e., energies in units of E 0 = ℏ ​ Ω subscript 𝐸 0 Planck-constant-over-2-pi Ω E_{0}=\hbar\Omega ). Parts of the data set are shown in Fig.  9 .

$\lambda$	$S$	$E_{\text{HO}}$	$V$	$K$	$K_{\text{vir}}$	$E$
$0$	$0.00258\left(1\right)$	$9.27\left(2\right)$	$9.27\left(2\right)$	$9.3\left(1\right)$	$9.27\left(3\right)$	$18.5\left(1\right)$
$0.01$	$0.00288\left(5\right)$	$9.24\left(8\right)$	$9.32\left(8\right)$	$9.2\left(4\right)$	$9.2\left(1\right)$	$18.5\left(5\right)$
$0.05$	$0.00393\left(5\right)$	$9.42\left(6\right)$	$9.82\left(5\right)$	$9.5\left(3\right)$	$9.22\left(9\right)$	$19.3\left(3\right)$
$0.1$	$0.00570\left(5\right)$	$9.51\left(4\right)$	$10.31\left(3\right)$	$8.7\left(2\right)$	$9.11\left(6\right)$	$19.0\left(2\right)$
$0.3$	$0.01861\left(10\right)$	$9.98\left(2\right)$	$12.31\left(1\right)$	$8.7\left(1\right)$	$8.81\left(2\right)$	$21.0\left(1\right)$
$0.5$	$0.04184\left(8\right)$	$10.483\left(8\right)$	$14.234\left(5\right)$	$8.59\left(5\right)$	$8.61\left(1\right)$	$22.82\left(5\right)$
$1$	$0.1475\left(1\right)$	$11.657\left(3\right)$	$18.644\left(2\right)$	$8.16\left(2\right)$	$8.164\left(5\right)$	$26.80\left(2\right)$
$3$	$0.6717\left(2\right)$	$15.859\left(2\right)$	$32.9762\left(9\right)$	$7.297\left(4\right)$	$7.301\left(2\right)$	$40.273\left(4\right)$
$10$	$0.99069\left(4\right)$	$27.149\left(1\right)$	$67.7235\left(7\right)$	$6.864\left(2\right)$	$6.862\left(2\right)$	$74.587\left(2\right)$
$0.01$	$0.02419\left(4\right)$	$9.327\left(5\right)$	$9.401\left(4\right)$	$9.22\left(3\right)$	$9.22\left(1\right)$	$18.62\left(3\right)$
$0.05$	$0.09311\left(7\right)$	$9.541\left(2\right)$	$9.865\left(2\right)$	$9.06\left(1\right)$	$9.056\left(6\right)$	$18.93\left(1\right)$
$0.1$	$0.16501\left(6\right)$	$9.7708\left(10\right)$	$10.3445\left(8\right)$	$8.922\left(5\right)$	$8.910\left(3\right)$	$19.266\left(5\right)$
$0.3$	$0.3721\left(1\right)$	$10.4931\left(8\right)$	$11.7701\left(7\right)$	$8.583\left(3\right)$	$8.578\left(2\right)$	$20.353\left(3\right)$
$0.5$	$0.5070\left(1\right)$	$11.0464\left(7\right)$	$12.8143\left(6\right)$	$8.394\left(2\right)$	$8.395\left(2\right)$	$21.208\left(2\right)$
$1$	$0.7015\left(2\right)$	$12.0977\left(10\right)$	$14.7323\left(8\right)$	$8.148\left(2\right)$	$8.146\left(3\right)$	$22.880\left(3\right)$
$3$	$0.9228\left(2\right)$	$14.707\left(1\right)$	$19.2833\left(9\right)$	$7.842\left(2\right)$	$7.842\left(3\right)$	$27.125\left(2\right)$
$10$	$0.99404\left(7\right)$	$19.4307\left(10\right)$	$27.2642\left(7\right)$	$7.680\left(2\right)$	$7.680\left(3\right)$	$34.944\left(2\right)$

Table 3. Table 3: Temperature dependence: Ab initio path integral Monte Carlo results for N = 6 𝑁 6 N=6 spin-polarized electrons in a 2 ​ D 2 𝐷 2D (top half) and 3 ​ D 3 𝐷 3D (bottom half) harmonic trap with Coulomb interaction and coupling strength λ = 0.5 𝜆 0.5 \lambda=0.5 . All results have been obtained for P = 200 𝑃 200 P=200 imaginary-time propagators (see Appendix  A for details) and are given in oscillator units (i.e., energies in units of E 0 = ℏ ​ Ω subscript 𝐸 0 Planck-constant-over-2-pi Ω E_{0}=\hbar\Omega ). Parts of the data set are shown in Figs.  5 and 10 .

$\beta$	$S$	$E_{\text{HO}}$	$V$	$K$	$K_{\text{vir}}$	$E$
$0.3$	$0.6580\left(2\right)$	$22.492\left(8\right)$	$25.438\left(8\right)$	$21.01\left(1\right)$	$21.02\left(1\right)$	$46.45\left(1\right)$
$0.5$	$0.3616\left(2\right)$	$15.246\left(4\right)$	$18.610\left(3\right)$	$13.55\left(1\right)$	$13.564\left(6\right)$	$32.16\left(1\right)$
$0.6$	$0.2507\left(2\right)$	$13.552\left(4\right)$	$17.040\left(4\right)$	$11.82\left(1\right)$	$11.807\left(6\right)$	$28.86\left(1\right)$
$0.8$	$0.1079\left(1\right)$	$11.557\left(5\right)$	$15.211\left(4\right)$	$9.77\left(2\right)$	$9.730\left(8\right)$	$24.98\left(3\right)$
$1$	$0.04184\left(8\right)$	$10.483\left(8\right)$	$14.234\left(5\right)$	$8.59\left(5\right)$	$8.61\left(1\right)$	$22.82\left(5\right)$
$1.1$	$0.02526\left(9\right)$	$10.12\left(1\right)$	$13.902\left(8\right)$	$8.19\left(7\right)$	$8.23\left(2\right)$	$22.10\left(7\right)$
$1.3$	$0.00894\left(8\right)$	$9.62\left(3\right)$	$13.44\left(2\right)$	$7.7\left(2\right)$	$7.71\left(5\right)$	$21.2\left(2\right)$
$1.5$	$0.00321\left(10\right)$	$9.15\left(9\right)$	$13.07\left(5\right)$	$7.3\left(5\right)$	$7.2\left(1\right)$	$20.4\left(5\right)$
$1.7$	$0.00114\left(6\right)$	$8.7\left(2\right)$	$12.77\left(8\right)$	$8\left(1\right)$	$6.7\left(3\right)$	$21\left(1\right)$
$2$	$0.00023\left(2\right)$	$8.5\left(3\right)$	$12.5\left(1\right)$	$4\left(1\right)$	$6.5\left(5\right)$	$17\left(1\right)$
$0.3$	$0.9298\left(1\right)$	$31.49\left(1\right)$	$33.68\left(1\right)$	$30.414\left(10\right)$	$30.39\left(2\right)$	$64.09\left(2\right)$
$0.5$	$0.7434\left(2\right)$	$20.172\left(5\right)$	$22.833\left(5\right)$	$18.829\left(8\right)$	$18.841\left(8\right)$	$41.66\left(1\right)$
$0.6$	$0.6163\left(2\right)$	$17.465\left(4\right)$	$20.286\left(4\right)$	$16.055\left(8\right)$	$16.054\left(6\right)$	$36.341\left(9\right)$
$0.8$	$0.3610\left(2\right)$	$14.290\left(4\right)$	$17.336\left(3\right)$	$12.778\left(9\right)$	$12.768\left(6\right)$	$30.114\left(9\right)$
$1$	$0.1712\left(1\right)$	$12.547\left(2\right)$	$15.736\left(2\right)$	$10.956\left(9\right)$	$10.952\left(4\right)$	$26.692\left(9\right)$
$1.1$	$0.1102\left(2\right)$	$11.954\left(6\right)$	$15.196\left(5\right)$	$10.34\left(2\right)$	$10.333\left(10\right)$	$25.54\left(2\right)$
$1.3$	$0.0408\left(2\right)$	$11.11\left(1\right)$	$14.429\left(8\right)$	$9.49\left(6\right)$	$9.44\left(2\right)$	$23.92\left(6\right)$
$1.5$	$0.01388\left(8\right)$	$10.56\left(2\right)$	$13.93\left(1\right)$	$8.70\left(7\right)$	$8.87\left(3\right)$	$22.63\left(7\right)$
$1.7$	$0.00452\left(7\right)$	$10.14\left(5\right)$	$13.56\left(3\right)$	$8.5\left(2\right)$	$8.42\left(8\right)$	$22.0\left(2\right)$
$2$	$0.00077\left(3\right)$	$9.9\left(1\right)$	$13.28\left(6\right)$	$8.8\left(5\right)$	$8.1\left(2\right)$	$22.1\left(5\right)$

Table 4. Table 4: Snapshots from PIMC simulations of N = 6 𝑁 6 N=6 spin-polarized fermions with coupling strength λ = 0.5 𝜆 0.5 \lambda=0.5 and P = 200 𝑃 200 P=200 imaginary-time slices. The left, center, and right column corresponds to the inverse temperature β = 0.3 𝛽 0.3 \beta=0.3 , β = 1 𝛽 1 \beta=1 , and β = 5 𝛽 5 \beta=5 , respectively, and the top, center, and bottom row to a system of ultracold atoms with dipole interaction in 2 ​ D 2 𝐷 2D , electrons with Coulomb interaction in a 2 ​ D 2 𝐷 2D harmonic trap, and electrons in a 3 ​ D 3 𝐷 3D harmonic trap. The corresponding full β 𝛽 \beta -dependence of the average sign S 𝑆 S for all three systems in shown in Fig.  10 .

	$\beta =0.3$	$\beta =1$	$\beta =5$
$2D$ Dipole
$2D$ Coulomb
$3D$ Coulomb

Equations

$$Z=\frac{1}{N!}\sum _{\sigma \in S_N}\text{sgn}(\sigma )\int \text{d}\mathbf{R}⟨\mathbf{R}|e^{-\beta \hat{H}}|{\hat{\pi }}_{\sigma }\mathbf{R}⟩,$$

$$Z=\int \text{d}\mathbf{X}W(\mathbf{X}),$$

$$Z^{\prime }=\int \text{d}\mathbf{X}|W(\mathbf{X})|,$$

$$⟨\hat{A}⟩=\frac{{⟨\hat{A}\hat{S}⟩}^{\prime }}{{⟨\hat{S}⟩}^{\prime }},$$

$$S:={⟨\hat{S}⟩}^{\prime }$$

$$\frac{\Delta A}{A}\sim \frac{1}{S\sqrt{M}}\sim \frac{e^{\beta N(f-f^{\prime })}}{\sqrt{M}}.$$

$$\hat{H}=-\frac{1}{2}\sum _{k=1}^N{\nabla }_k^2+\frac{1}{2}\sum _{k=1}^N{\hat{\mathbf{r}}}_k^2+\sum _{k>l}^N\frac{\lambda }{|{\hat{\mathbf{r}}}_l-{\hat{\mathbf{r}}}_k{|}^{\alpha }},$$

$$\lambda =\{\begin{array}{ll}\frac{e^2}{4\pi {ϵ}_0{l}_0{E}_0},& \text{mbox}if\alpha =1\\ \frac{D}{4\pi l_0^3{E}_0}& \text{mbox}if\alpha =3\end{array},$$

$$\hat{H}=-\frac{1}{2}\sum _{k=1}^N{\nabla }_k^2+\sum _{k>l}^{N}w({\hat{\mathbf{r}}}_l,{\hat{\mathbf{r}}}_k).$$

$$P(A)=\frac{e^{-\frac{(A-\mu {)}^2}{2{\sigma }^2}}}{\sqrt{2\pi {\sigma }^2}},$$

$$\Delta {⟨V⟩}^{\prime }={\left(\frac{1}{M}\sum _{i=1}^M(V_i-{⟨V⟩}^{\prime }{)}^2\right)}^{1/2}.$$

$$S_N(\beta )=a_N{e}^{-b_N\beta },$$

$$\frac{\Delta K}{K}(\beta )=a_6{e}^{b_6\beta }{c}_K,\frac{\Delta V}{V}(\beta )=a_6{e}^{b_6\beta }{c}_V,$$

$$P(l)=\frac{Z_1^{\prime }(l\beta )Z_{N-l}^{\prime }(\beta )}{l{Z}_N^{\prime }(\beta )}.$$

$$S_{\beta }(N)=a_{\beta }{e}^{-b_{\beta }N},$$

$$K=E_{\text{HO}}-\alpha \frac{V-E_{\text{HO}}}{2},$$

$${\hat{\rho }}_{ϵ}=e^{-ϵ(\hat{K}+\hat{V})}=e^{-ϵ\hat{K}}{e}^{-ϵ\hat{V}}+\mathcal{O}\left({ϵ}^2\right),$$

$$e^{-\beta (\hat{K}+\hat{V})}=\underset{P\to \infty }{\lim }{\left(e^{-ϵ\hat{K}}{e}^{-ϵ\hat{V}}\right)}^P.$$

$$A(P=\infty )=A(P)+\frac{{\delta }_A}{P^2}.$$

Full text

I Introduction

The numerical solution of the well-known, but highly complex equations that govern quantum mechanics using modern high-performance computers has emerged as one of the most active and successful fields in theoretical physics and chemistry. A particularly useful approach to accomplish this goal for a correlated quantum system in thermodynamic equilibrium (i.e., at finite temperature) was already outlined by Feynman feynman , who proposed to map the complicated quantum system of interest onto a classical ensemble of interacting ring-polymers chandler via the path-integral formalism kleinert . The basic idea of the path-integral Monte-Carlo (PIMC) method berne3 ; imada ; pollock ; cep ; review is to stochastically evaluate the resulting high-dimensional integrals using the Metropolis algorithm metropolis , which, remarkably, does not suffer from the curse of dimensionality curse that renders standard quadrature methods unfeasible in this case filinov_chapter ; liu_monte_carlo .

Since its first application to He4 in the late sixties pimc_original ; pimc_original2 , PIMC has emerged as one of the most successful tools in statistical physics and has allowed for profound insights into exciting physical phenomena such as superfluidity ceperley_superfluid ; sindzingre ; kwon ; dornheim_superfluid and Bose-Einstein-condensation BEC1 ; BEC2 . Moreover, PIMC provides exact simulations at strong coupling, which makes it possible to study crystallization in real quantum systems jones_crystal ; filinov_PRL ; clark_casula , and the direct access to imaginary-time correlation functions berne1 ; berne2 ; dynamic_folgepaper can be used as input for an analytic continuation jarrell ; schoett , which makes even possible the computation of dynamic properties such as collective excitations dornheim_dynamic ; vitali ; supersolid . In fact, novel Monte-Carlo sampling techniques allow for exact calculations of up to $N\sim 10^{4}$ bosons and boltzmannons (i.e., distinguishable particles obeying Boltzmann statistics) boninsegni1 ; boninsegni2 .

On the other hand, the situation is entirely different in the case of fermions. More specifically, the antisymmetry of the fermionic density matrix under the exchange of particles [cf. Eq. (1)] leads to a near cancellation of positive and negative terms both with decreasing temperature and increasing the system size ceperley_fermions . This issue is commonly known as the fermion sign problem (FSP) loh ; troyer ; lyubartsev ; vozn , and generally prevents fermionic PIMC simulations once quantum degeneracy effects start to get important fsp_note . This is very unfortunate, as correlated Fermi systems offer a wealth of interesting effects such as the BCS-BEC crossover fermi1 ; fermi2 ; fermi3 in ultracold atoms and the formation of Wigner molecules in quantum dots wigner_molecule1 ; wigner_molecule2 ; wigner_molecule3 .

Of particular importance is the so-called warm dense matter (WDM) regime review ; ross ; koenig ; fortov_review , an extreme state of matter with high temperatures ($T\sim 10^{4}-10^{8}$K) and extreme densities ($n\sim 10^{21}-10^{27}$cm*-3*). These conditions occur in astrophysical objects such as giant planet interiors militzer1 ; militzer2 ; knudson and brown dwarfs saumon1 ; saumon2 ; becker , and are expected to play an important role on the pathway towards inertial confinement fusion hu1 ; hu2 . Moreover, WDM is now routinely realized in the lab (see Ref. falk_wdm for a topical review article) and constitutes one of the most active frontiers in plasma science wdm_book . The theoretical description of WDM, however, is notoriously difficult due to the nontrivial interplay of (i) Coulomb coupling, (ii) quantum degeneracy effects, and (iii) thermal excitations. This regime is typically characterized by two parameters, which are both of the order of one: the density parameter $r_{s}=\overline{r}/a_{\textnormal{B}}$ (with $\overline{r}$ and $a_{\textnormal{B}}$ being the mean interparticle distance and first Bohr radius) and the degeneracy parameter torben_eur $\theta=k_{\textnormal{B}}T/E_{\textnormal{F}}$ (with $E_{\textnormal{F}}$ being the usual Fermi energy quantum_theory ). Therefore, both perturbation theory and ground state methods are not applicable, which leaves ab initio PIMC simulations as one of the most promising options brown_chapter .

Consequently, there has been a spark of new developments in the field of fermionic quantum Monte-Carlo simulations at finite temperature over the last years cpimc_original ; brown_ethan ; blunt1 ; schoof_prl ; malone1 ; blunt2 ; malone2 ; dornheim ; dornheim2 ; vladimir_UEG ; groth ; dornheim3 ; dornheim_prl ; dornheim_pop ; groth_prl ; dubois ; dornheim_pre ; groth_jcp ; claes ; dornheim_cpp ; brenda ; universe ; dornheim_neu . Despite this exciting progress, a thorough study of the fermion sign problem itself seems to be missing. In this work, we aim to fill this gap by presenting a detailed practical investigation of the FSP within standard PIMC simulations of i) electrons in quantum dots, ii) ultracold dipolar atoms in a harmonic confinement, and iii) the uniform electron gas at warm dense matter conditions review ; groth_prl . More specifically, we study the manifestation of the FSP regarding different parameters (e.g., system size, coupling strength, etc.) and discuss the probability distribution of Monte-Carlo expectation values, which, in the presence of a sign problem, is not necessarily given by a simple Gaussian. In addition, we provide extensive benchmark data, which will aid the development of new methods and can be used to gauge the accuracy of novel approximations.

The paper is organized as follows: In Sec. II, we introduce the required theory, in particular the standard path integral Monte Carlo approach (II.1), followed by the fermion sign problem (II.2) and the considered system types and Hamiltonians (II.3). In Sec. III, we present our simulation results, starting with a detailed discussion of the Monte-Carlo sampling and the probability distribution of fermionic expectation values (III.1). In addition, we study the manifestation of the FSP with respect to temperature (III.2), system-size (III.3), interaction-strength and -type (III.4), and the dimensionality (III.5), all for electrons and ultracold atoms in a harmonic confinement. Lastly, we extend our considerations to the uniform electron gas in the warm dense matter regime, where the FSP exhibits a somewhat different manifestation regarding system size. The paper is concluded by a concise summary and discussion in Sec. IV.

II Theory

II.1 Path Integral Monte Carlo

Throughout this work, we restrict ourselves to the discussion of $N$ spin-polarized fermions in the canonical ensemble, i.e., the inverse temperature $\beta=1/k_{\textnormal{B}}T$, volume $V$ (or trap frequency $\Omega$ in case of a harmonic confinement, see Eq. (7) below) and particle number $N$ are fixed. The central quantity in statistical physics is the partition function, which can be written in coordinate space as

[TABLE]

where $\mathbf{R}=(\mathbf{r}_{1},\dots,\mathbf{r}_{N})^{T}$ contains the coordinates of all particles. Since we are interested in fermions, we have to explicitly evaluate the sum over all possible permutations of particle coordinates $\sigma$, with $S_{N}$ denoting the permutation group and $\hat{\pi}_{\sigma}$ being the corresponding permutation operator. Note that the sign $\textnormal{sgn}(\sigma)$ is positive (negative) for an even (odd) number of pair exchanges. For completeness, we mention that we restrict ourselves to the spin-polarized case (i.e., only one species of fermions, like spin-up electrons) throughout this work, but the generalization to multiple particle species is straightforward and does not affect the manifestation of the FSP. To make the evaluation the matrix elements of the density operator $\hat{\rho}=e^{-\beta\hat{H}}$ in Eq. (1) possible, one typically performs a Trotter decomposition trotter and finds that $Z$ can be expressed as the sum over all closed paths $\mathbf{X}$ in the imaginary time $\tau$. However, since both the derivation and final formulas have already been presented elsewhere cep ; review , they need not be repeated here. For the present purposes, it is fully sufficient to work with the abstract expression

[TABLE]

which can be interpreted as follows: The $PND$-dimensional (with $P$ denoting the number of so-called imaginary time-slices, cf. Fig. 1, and $D$ being the dimensionality of the system) variable $\mathbf{X}$ constitutes a so-called configuration, and each configuration contributes to $Z$ with the appropriate configuration weight $W(\mathbf{X})$, which is a function that can be readily evaluated.

This is illustrated in Fig. 1, where we show example configurations from a PIMC simulation of $N=3$ fermions. First and foremost, we note that each particle is now represented by an entire path in the imaginary time $\tau$, with $P=6$ imaginary time-slices. In panel (a), there is no exchange of particle coordinates and, consequently, the configuration weight $W(\mathbf{X})$ is positive. In contrast, the second depicted configuration contains an exchange-cycle comprised of two fermions. Due to this pair-exchange, the corresponding $W$ is negative.

The basic idea of the path integral Monte Carlo method berne3 ; cep is to generate a Markov chain of path-configurations $\mathbf{X}$ that are distributed according to $P(\mathbf{X})=W(\mathbf{X})/Z$. Although the normalization $Z$ is not known, this can be accomplished efficiently using the celebrated Metropolis algorithm metropolis . Indeed, simulations of up to $N\sim 10^{4}$ bosons and boltzmannons (i.e., distinguishable particles obeying classical Boltzmann statistics clark_casula ; dornheim_analyzing ) are feasible using novel Monte-Carlo sampling techniques boninsegni1 ; boninsegni2 without the introduction of any approximation. Unfortunately, since the weight function $W(\mathbf{X})$ is not strictly positive in the case of fermions, it cannot be interpreted as a probability distribution, which, as we shall see in the next section, is the origin of the infamous fermion sign problem.

For completeness, we mention that it is, at least in principle, possible to recast Eq. (1) into a sum over only positive terms by exploiting the nodal structure of the density matrix ceperley_fermions ; fermion_nodes . However, since the exact nodes are a-priori unknown, this simplification comes at the cost of an uncontrolled approximation node_note .

II.2 The fermion sign problem

The first task at hand is to find a way to generate the paths $\mathbf{X}$ using the Metropolis algorithm, although their weight function is negative. In practice, we switch to the modified partition function

[TABLE]

where the paths $\mathbf{X}$ are now generated according to the absolute value of the weight function. We note that in the case of standard PIMC, as it has been introduced above, Eq. (3) coincides with the (symmetrized) bosonic partion function, which has some interesting implications that are discussed later on. To calculate the fermionic expectation value of an observable $\hat{A}$, we then have to evaluate the ratio

[TABLE]

where the operator $\hat{S}$ measures the sign of the configuration weight, i.e., $S(\mathbf{X})=W(\mathbf{X})/|W(\mathbf{X})|$. The problem with this approach is that both the enumerator and the denominator in Eq. (4) vanish simultaneously both towards low temperature (i.e., large $\beta$) and with increasing system size $N$. This is captured by the average value of $\hat{S}$, which is given by the ratio of the fermionic and bosonic partition function

[TABLE]

and which we will simply refer to as the average sign throughout this work. In fact, Eq. (II.2) constitutes a direct measure for the amount of cancellations within a fermionic PIMC simulation, and exponentially decays both with $N$ and $\beta$ (with $f$ and $f^{\prime}$ being the free energy density of the fermionic and modified system, respectively). This is bad news, because a small sign (typically $S\sim 10^{-3}$) means that simulations are no longer feasible. This can be understood by considering the relative Monte-Carlo error of Eq. (4), which is given by ceperley_fermions

[TABLE]

Evidently, the statistical error exponentially increases with $N$ and $\beta$, which can only be compensated by increasing the number of Monte-Carlo samples as $~{}1/\sqrt{M}$. In practice, one thus quickly runs into an exponential wall, which is nothing else than the fermion sign problem.

II.3 System types and Hamiltonians

II.3.1 Harmonic confinement

The most widely used model system that is considered in this work are fermions in a harmonic confinement, which is governed by the Hamiltonian

[TABLE]

where we assume oscillator units, i.e., the characteristic length $l_{0}=\sqrt{\hbar/m\Omega}$ (with $\Omega$ being the trap frequency) and energy scale $E_{0}=\hbar\Omega$. Of particular importance is the exponent $\alpha\in\{1,3\}$, which distinguishes between Coulomb interaction ($\alpha=1$, corresponding to electrons in a quantum dot dornheim ; dornheim_analyzing ; reimann ) and dipole interaction ($\alpha=3$, corresponding to ultracold atoms stuhler ; dynamic_alex1 ; jan_willem ). In addition, the coupling constant $\lambda$ is defined as the ratio of the interaction energy to $E_{0}$,

[TABLE]

with $D$ being the usual dipole-dipole interaction constant, see, e.g., Ref. dynamic_alex1 . Finally, we consider two- and three-dimensional systems in this work, and the dimensionality of the harmonic confinement is always equal to the overall number of dimensions.

II.3.2 Uniform electron gas

The second type of model system that we consider in this work is the uniform electron gas (see Refs. review ; loos for topical review articles), which is defined as an ensemble of $N$ electrons in a periodic box of length $L$ and volume $V=L^{3}$. The corresponding Hamiltonian is given by

[TABLE]

Note that we always assume Hartree atomic units (i.e., energies in Hartree and distances in units of the first Bohr radius $a_{0}$) when discussing the UEG. Let us briefly turn our attention to the pair interaction potential $w(\mathbf{\hat{r}}_{l},\mathbf{\hat{r}}_{k})$ in Eq. (9). In order to mitigate finite-size effects, one typically employs the Ewald summation, which takes into account both the interaction between the two electrons $l$ and $k$ (and the respective positive background) and the infinite array of periodic images fraser . In this work, we use a pre-averaged (i.e., with respect to the orientation of the array of images) Ewald potential introduced by Yakub and Ronchi yakub1 ; yakub2 , where the infinite sums both in real and reciprocal space are evaluated analytically beforehand. This leads to a significant saving of time, while the differences to the real Ewald summation are expected to be small under the conditions considered in this work.

For completeness, we mention that a complete thermodynamic description of the UEG at warm dense matter conditions was achieved only recently groth_prl on the basis of novel configuration PIMC and permutation-blocking PIMC simulation data, see Ref. review for a comprehensive discussion.

III Results

All results in this work have been obtained using a canonical adaption mezza of the worm algorithm boninsegni1 ; boninsegni2 . Further, we use $P=200$ imaginary-time propagators based on the primitive action, see Appendix A for details, and Refs. brualla ; sakkos for an accessible discussion. All fermionic results listed in Tabs. 1, 2, and 3 are conveged with respect to $P$ within the given statistical uncertainty.

III.1 Monte-Carlo sampling and probability distribution of expectation values

Let us start our discussion of the fermion sign problem with an illustration of the sampling of the expectation value of an observable. In Fig. 2, we show PIMC results for a simulation of $N=6$ electrons in a $2D$ harmonic trap [cf. Eq. (7)] at moderate coupling $\lambda=0.5$ and two inverse temperatures, $\beta=0.2$ (red) and $\beta=1.6$ (blue). Panel (a) shows a series of $M=300$ measurements for the signed potential energy $VS$, i.e., the enumerator from Eq. (4). The solid red line corresponds to $\beta=0.2$, which is a comparatively high temperature, where fermionic exchange-effects are not that important. Consequently, the sign stays mostly positive (with $S\approx 0.82$), and sign-changes due to permutation cycles appear as brief negative spikes in the series of measurements. In stark contrast, the blue line corresponds to $\beta=1.6$, and the situation looks completely different: at this low temperature, positive and negative signs appear with a similar frequency and the average sign has decreased to $S\approx 0.002$.

To further illustrate the origin of these cancellations, it is instructive to consider the modified (bosonic) probability distribution, which is used to generate the paths $\mathbf{X}$. To this end, we show in Fig. 2 (b) ($\beta=0.2$) and Fig. 2 (c) ($\beta=1.6$) the radial density distribution $n(r)$ both for Bose (red squares) and Fermi (blue crosses) statistics. At high temperature, the two data sets are very similar and the most significant deviations occur around the center of the trap, where the density is at the maximum. At $\beta=1.6$, on the other hand, the two densities exhibit severe discrepancies over the entire $r$-range. While bosons tend to cluster around the center of the trap, the fermions are pushed outward by the Pauli blocking. Since the paths in our simulations are distributed according to the bosonic density, the difference in the results for fermions at low temperature can only be accomplished by the cancellation and subsequent division by the small value for $S$ [cf. Eq. (4)].

Let us next consider the probability distribution of the expectation values within a fermionic PIMC simulation. According to the central limiting theorem monte_carlo_book , the average value of a Metropolis Monte-Carlo calculation of an expectation value $\braket{\hat{A}}$ with $M$ measurements (and $M$ being large) is normally distributed around the exact value, and the standard deviation decreases as $\sigma_{M}\sim 1/\sqrt{M}$. This is verified in Fig. 3, where we show histograms for the Monte-Carlo average of $S$ (a) and $VS$ (b) for $N_{\textnormal{s}}=600$ independent seeds with $M=5\cdot 10^{6}$ measurements per seed, for a system of $N=6$ noninteracting fermions in a $2D$ harmonic trap at $\beta=1$. The blue bars correspond to our PIMC data, and the solid red curves to Gaussian fits according to

[TABLE]

with $\sigma$ and $\mu$ being the free parameters. Evidently, we do indeed find the expected normal distribution for both cases, which means that the statistical uncertainty for a single seed can be straightforwardly estimated from the Monte-Carlo data via

[TABLE]

For completeness, we note that the error bars given in all tables and figures have been obtained by evaluating Eq. (11) for $N_{\textnormal{s}}$ statistically independent seeds, instead of $M$ measurements from a single seed, see Ref. error_note for details.

Let us next consider the distribution of the fermionic observable $\braket{V}=\braket{VS}^{\prime}/\braket{S}^{\prime}$, which is shown in Fig. 3 (c). Remarkably, the histogram does not exhibit a Gaussian form, and the corresponding normal fit is not in agreement with the PIMC results. More specifically, the simulation results show a distinct tail towards large values of $VS/S$, with the two largest outliers (see the two blue arrows in the plot) being located around $VS/S=24$, i.e., around $16$ standard deviations (assuming the $\sigma$ value from the Gaussian fit) away from the mean.

The reason for this peculiar finding is the nonlinear nature of the fermionic expectation value from Eq. (4). In fact, it can be shown hatano that the probability distribution of the ratio of $\braket{VS}^{\prime}$ and $\braket{S}^{\prime}$ is the superposition of a Lorentzian (also known as Cauchy distribution) and a Gaussian, with the former one being responsible for the tail. Moreover, the sample deviation as defined in Eq. (11) actually diverges, and, therefore, does not constitute a good measure for the real uncertainty in the fermionic expectation value $\braket{V}$. For completeness, we mention that a similar behavior has been found in other fields, most notably financial modelling mandelbrot .

To further illustrate the occurrence of these tail events in our fermionic PIMC simulation, we show the series of the average values for all $N_{\textnormal{s}}=600$ seeds of both the ratio $VS/S$ (blue, left $y$-axis) and the sign $S$ (red, right $y$-axis) in Fig. 3 (d). Let us first consider the blue curve: evidently, the expectation values of most seeds are located somewhere around the mean value, with a few spikes corresponding to the upward outliers. In contrast, the red curve does not exhibit any spikes, and we have already seen that $S_{i}$ follows a normal distribution, see Fig. 3 (a). A comparison of both curves reveals that the spikes in the ratio appear in those seeds with the smallest values of $S$, and the two smallest values, which are responsible for the blue arrows in Fig. 3 (c), are highlighted by red crosses. Indeed, these $S_{i}$ are more than an order of magnitude smaller than the corresponding mean value of the distribution.

Up to this point, one might conclude that fermionic PIMC simulations appear to be doomed as 1) we do not have a good measure for the statistical uncertainty, which would make the Monte-Carlo expectation value an uncontrolled approximation, and 2) the distribution of the ratio $P(VS/S)$ is fat-tailed and outliers exceeding $16$-times the standard deviation (often called black swan events taleb ) do appear with finite probability. However, as we will see next, all is not lost.

In Fig. 4, we show PIMC results for the same conditions as in Fig. 3, but with dipole-interaction and a coupling constant $\lambda=0.1$ (i.e., ultracold atoms). Due to the dipolar repulsion, fermionic exchange is suppressed, and we find an average sign of $S\approx 0.165$, as compared to $S\approx 0.0025$ for the noninteracting case. Panel (a) shows results for $P(S)$, and we again find the expected normal distribution. In contrast to the noninteracting case, this time $P(S)$ has a significantly smaller relative deviation $\sigma/S$, and no expectation values $S_{i}$ with a value that is an order of magnitude smaller than the average appear. Consequently, there are no spikes in the seed-averages of the ratio, and the corresponding distribution $P(VS/S)$ [Fig. 4 (b)] cannot be distinguished from a Gaussian.

In summary, the nonlinear nature of the fermionic expectation value Eq. (4) causes the distribution $P(VS/S)$ to be non-Gaussian, with a fat tail towards larger values. To put it in another way, if the relative uncertainty of the denominator (i.e., $S$) is large, the smallest signs $S_{i}$ lead to spikes in the fermionic observable. In contrast, if the relative error of the sign is small (as in Fig. 4), these spikes do not appear (or are sufficiently unlikely), and the resulting distribution $P(VS/S)$ cannot be distinguished from a normal distribution in practice. Since $S$ itself does obey a normal distribution in any case, this condition can always be checked, and fermionic PIMC results can safely be labelled as quasi-exact after all. Thus, the statistical uncertainty for the fermionic expectation values given in all figures and data tables has been computed assuming the Gaussian form from Eq. (11), which is reliable for all presented cases.

III.2 Temperature dependence

Let us next investigate the manifestiation of the FSP upon decreasing the temperature. To this end, we simulate spin-polarized electrons in a $2D$ harmonic trap at intermediate coupling $\lambda=0.5$. Fig. 5 (a) shows PIMC data for the $\beta$-dependence of the average sign $S$ for $N=6$ (red squares) and $N=9$ (blue crosses). First and foremost, we note that both data sets exhibit a qualitatively similar behavior: for small $\beta$, the system is nearly classical and $S$ is large, whereas it monotonically decreases with increasing $\beta$. In addition, the sign for $N=9$ is always smaller than for $N=6$, as it is expected. To verify the predicted exponential decrease of $S$ with $\beta$ (see Eq. (II.2) in Sec. II.2), we perform fits (starting at $\beta\geq 1$) of the form

[TABLE]

with $a_{N}$ and $b_{N}$ being the free parameters. The results are shown as the dashed lines and are indeed in excellent agreement with the PIMC data for $\beta\geq 1$. Note that at higher temperature, the free energy density $f$ [cf. Eq. (II.2 )] changes significantly with $\beta$, which leads to the deviation from Eq. (12) in this regime.

Panel (b) shows the kinetic energy $K$ for the case of $N=6$ both for fermions (blue crosses) and bosons (red squares). Firstly, we mention that the relative deviations between Bose and Fermi statistics increase towards low temperature, as it is expected. Secondly, the red curve is very smooth over the entire depicted $\beta$-range, and the error bars cannot be seen with the naked eye. In contrast, the fermionic data are accurate for small $\beta$, but eventually the error bars markedly increase when $S$ becomes small.

To check if we are really running into the exponential wall as predicted by Eq. (6), we show the corresponding relative statistical uncertainty of $K$ (blue crosses) and the total potential energy $V$ (i.e., both interaction and external potential, red squares) in Fig. 5 (c). The dashed lines depict exponential fits (for $\beta\geq 1$) of the form

[TABLE]

with $C_{K}$ ($C_{V}$) being the only free parameter, as $a_{6},b_{6}$ have already been determined by a fit to $S$. Evidently, the data and the fit are in excellent agreement both for $V$ and $K$, which (sadly) confirms the severity of the FSP.

Extensive PIMC data for the temperature-dependence of electrons in $2D$ and $3D$ (cf. Sec. III.5) are given in Tab. 3.

Let us conclude this section on the temperature dependence with a brief excursion to the distribution of permutation-cycles. In Fig. 6, we investigate the probability to find a particle involved in a permutation cycle of length $l$, $P(l)l$, see Ref. dornheim_permutation_cycles for a topical introduction and extensive discussion. Panel (a) shows simulation results for $N=6$ noninteracting fermions in a $2D$ harmonic confinement at $\beta=0.3$ (red squares), $\beta=1$ (blue crosses), and $\beta=5$ (green circles). At the highest temperature, the paths resemble classical particles (see also Tab. 4), pair-exchanges are quite improbable and approximately $90\%$ of particles are not involved in any exchange. Therefore, we find an average sign of $S\approx 0.51$. At $\beta=1$, the situation has already drastically changed, and the distribution has become significantly flatter. Due to the resulting cancellation of positive and negative terms, the sign has decreased to $S\approx 0.002$. At the lowest temperature, $\beta=5$, the distribution has become almost completely flat and the sign vanishes within the given statistical uncertainty. In fact, it does hold $P(l)l=1/N$ in the zero temperature limit, which means that PIMC simulations are not possible in the ground state since the sign vanishes krauth_book .

To verify the correctness of our implementation, we compare our PIMC data to the theoretical result for $P(l)$, which can be phrased in terms of the noninteracting partition function at different temperature and system-size as krauth_book ; dornheim_permutation_cycles

[TABLE]

The corresponding dashed lines are in perfect agreement with our PIMC data for all temperatures $\beta$ and cycle-lengths $l$.

In Fig. 6 (b), we show results for $P(l)l$ for the same conditions as in panel (a), but with Coulomb (red) and dipole (blue) interaction and coupling strength $\lambda=0.5$. For $\beta=0.3$, we observe a qualitatively similar behavior as for the noninteracting case shown above. Still, the repulsion between the particles leads to a steeper decay of $P(l)l$ towards large $l$, which is even more pronounced in the case of dipoles. This is a direct consequence of the stronger repulsion at small distances in the latter case, which renders the formation of exchange-cycles within the simulation even more improbable, cf. the discussion of Fig. 9. For $\beta=1$, the distribution is significantly less flat than in the noninteracting case, which is again more pronounced for the dipolar interaction.

III.3 System-size dependence

Another question that is of fundamental importance regarding fermionic PIMC simulations is the manifestation of the FSP with the system size. This topic is investigated in Fig. 7 (a), where we show PIMC results for the average sign $S$ for $N=6$ electrons in a $2D$ harmonic trap with the coupling strength $\lambda=0.5$ and the inverse temperature $\beta=1$ (red squares), and $\beta=0.3$ (blue crosses). Both data sets exhibit a steep decay with increasing $N$, which is significantly more pronounced for the lower temperature, as it is expected. To check the predicted exponential decay with $N$, we perform fits of the form

[TABLE]

with $a_{\beta},b_{\beta}$ being the free parameters. The results for Eq. (15) are shown as the dashed lines, and are in qualitative agreement with the PIMC data. Still, the simulation results appear to exhibit an even faster decay than the exponential function from Eq. (15).

To explain this finding, we plot the radial density $n(r)$ for $\beta=0.3$ and three different particle numbers in Fig. 8. Evidently, the addition of particles leads to an increased density, in particular around the center of the trap. Therefore, the system becomes more quantum degenerate, and the average sign decreases even faster than the exponential fit from Eq. (15). It is important to note that the situation is entirely different for a uniform, periodic system like the UEG (see Sec. III.6), where a change in system size does hardly affect the degree of degeneracy because the density remains constant. Therefore, one does indeed find an exponential decay of $S$ with $N$ in that case, cf. Fig. 12.

Let us conclude this discussion of the system-size dependence of the FSP with the consideration of an observable. To this end, we show the $N$-dependence of the total energy per particle $E/N$ for the case of $\beta=1$ in Fig. 7 (b) both for Fermi (blue crosses) and Bose (red squares) statistics. Evidently, the energy per particle does not remain constant for both cases, but increases, as it is expected. This trend is even more pronounced for the case of fermions, which are subject to the Pauli blocking. Thus, they get pushed away from the center of the trap, where the energy due to the external harmonic potential is large. In addition, we note the increasing error bars in the blue curve, which do not appear for bosons, and are a direct consequence of the corresponding decrease in $S$, cf. Eq. (6).

Extensive PIMC results for the $N$-dependence of spin-polarized electrons are given in Tab. 1.

III.4 Interaction and coupling-strength dependence

A somewhat less well understood question is the dependence of the FSP on the interaction-type and coupling strength. In Sec. III.2, we have already seen that ultracold atoms with dipole interaction [$\alpha=3$, cf. Eq. (7)] exhibit a comparatively less severe sign problem than electrons at the same value of the coupling parameter $\lambda$, cf. Fig. 6. In Fig. 9, we present a more systematic investigation of this issue by performing PIMC simulations of $N=6$ electrons (red squares) and ultracold atoms (blue crosses) at $\beta=1$. Panel (a) shows the $\lambda$-dependence of the average sign $S$ over more than three orders of magnitude in the coupling strength. At $\lambda=10$, the particles are spatially separated by the strong repulsion for both types of interaction and fermionic exchange is suppressed. With decreasing $\lambda$, the particles get increasingly close to each other and the sign decreases for both data sets, although it does so significantly faster in the case of the Coulomb interaction. More specifically, the red curve has already almost attained the noninteracting limit ($\lambda=0$, dash-dotted black line) at $\lambda=0.01$, whereas the corresponding blue data point is still one order of magnitude larger. This is a direct consequence of the comparatively larger repulsion for the dipole-interaction at small distances, as we have already discussed in Sec. III.2.

Let us next consider the corresponding $\lambda$-dependence of the potential energy $V$, which is shown in Fig. 9 (b). The squares and crosses depict data for Coulomb- and dipole-interaction, respectively, and the grey points show the corresponding results for Bose statistics. At strong coupling, quantum statistics are negligible, the grey and colored points are in perfect agreement, and the system resembles a semi-classical Coulomb- or dipole-system. With decreasing $\lambda$, there appears a transition region until eventually both the fermions and the bosons attain their respective noninteracting limit (dash-dotted black lines). Remarkably, this happens much faster for fermions, which are already in good agreement for both types of interaction at $\lambda=0.1$, than for bosons, which still significantly deviate for $\lambda=0.01$.

The reason for this striking difference is illustrated in Fig. 9 (c), where we show the radial density $n(r)$ at $\lambda=0.1$ for Coulomb- (red), dipole- (blue), and no interaction (green) and for both Fermi (solid) and Bose (dashed) statistics. Let us first consider the three fermionic curves, which are in good agreement with each other, as it is by now expected from the observed corresponding agreement in $V$ [cf. Fig. 9 (b)]. In stark contrast, the bosonic curve for the dipole-interaction significantly deviates from the other two, which explains the observed behavior in both $S$ and $V$: for Coulomb-interaction (or the noninteracting case), the paths that are sampled within our PIMC simulation are clustered around the center of the trap. The fermionic density, which remains large for much higher values of $r$, must subsequently be recovered by the cancellation and division by a small average sign $S$ according to Eq. (4). For dipole-interaction, on the other hand, the strong repulsion at small distances has a very similar effect to the Pauli blocking, so that already the bosonic density is very close to its fermionic analogue. Consequently, the bosonic and fermionic configuration spaces and partition functions are almost equal, and the average sign $S\approx 0.17$ is large.

In summary, we have found that the fermion sign problem is much less severe for interaction-types with a strong short-range repulsion. This makes the future systematic study of ultracold fermionic dipolar atoms stuhler (in the trap, in periodic boundary conditions, or in other geometries like bilayers dynamic_alex2 ) a promising project for future research.

Extensive PIMC data for the coupling-strength dependence of both electrons and ultracold atoms are given in Tab. 2.

III.5 Dimensionality versus Interaction-type, virial theorem

The last question to be investigated in this work regarding fermions in a harmonic confinement is the impact of the dimensionality. In Fig. 10, we show the $\beta$-dependence of the average sign $S$ for $N=6$ and $\lambda=0.5$. The red squares, blue crosses, and green circles depict our PIMC results for $2D$ Coulomb, $2D$ dipoles, and $3D$ Coulomb, respectively. The corresponding dashed lines depict exponential fits according to Eq. (12), which are in excellent agreement with the data for all types of systems. As usual, the dipole interaction leads to a significantly less steep decay of $S$ with $\beta$, cf. Sec. III.4. In addition, we find that the Coulomb systems exhibit a very similar behavior of $S$, although the exponential decrease starts at somewhat lower temperatures in $3D$. This is most likely due to the additional degree of freedom in this case, which allows the electrons to avoid each other more effectively.

In Tab. 4, we compare snapshots from our PIMC simulation for all three kinds of system types at three different temperature regimes. For $\beta=0.3$ (left column), all systems exhibit a very similar behavior, with the extension of the paths, which is proportional to the thermal wave length $\lambda_{\beta}=\sqrt{2\pi\hbar^{2}\beta/m}$, being significantly smaller than the average inter-particle distance $\overline{r}$. At $\beta=1$ (center column), the paths are clustered more closely around the center of the trap in all three cases (the scale is equal for all three depicted values of $\beta$), and $\lambda_{\beta}$ is comparable to $\overline{r}$. In the case of dipole interaction (top row), the paths of individual particles are still mostly separated by the strong short-range repulsion (cf. Sec. III.4), and no exchange-cycle is present in the snapshot (the two particles in the front are close, but not connected). For $2D$ Coulomb (center row), on the other hand, fermionic exchange already plays a dominant role, and there appear two permutation-cycles with $N=3$ and $N=2$ particles in it. Going to $3D$ (bottom row), the situation looks qualitatively the same as in $2D$, and permutation-cycles are present, too. At low temperature, $\beta=5$, the thermal wavelength is larger than $\overline{r}$ in all three cases and the system is fully quantum degenerate. Yet, the dipole interaction manages to push the particles away from each other, and the corresponding average sign is several orders of magnitude larger than for Coulomb interaction, cf. Fig. 10. For Coulomb interaction in $2D$ and $3D$, the particles form an entangled knot of paths around the center of the trap, the probability to find an exchange cycle of length $l$ is almost constant (cf. Fig. 6), and the average sign vanishes within the given statistical uncertainty.

Let us conclude this section by investigating the virial theorem greiner_book , which gives a relation between the different contributions to the total energy. For example, it holds

[TABLE]

with $V$ and $E_{\textnormal{HO}}$ being the total potential energy and the energy due to the external potential, respectively. Recall that $\alpha\in\{1,3\}$ distinguishes between Coulomb and dipolar interaction, cf. Eq. (7).

In Fig. 11, we show the relative difference between Eq. (16) and the kinetic energy as evaluated using the standard PIMC thermodynamic estimator (for an extensive discussion on energy estimation in PIMC simulations, see Ref. janke ). Panels (a) and (b) show results for bosons and fermions, and the red squares, blue crosses, and green circles depict data for $2D$ Coulomb, $2D$ dipoles, and $3D$ Coulomb, respectively, with the red and green curves having been shifted for better visibility. Due to the absence of the FSP for Bose statistics, the statistical uncertainty is of the order of $\Delta K/K\sim 10^{-4}$, and the difference between both results for $K$ vanishes within the error bars for all three data sets. For fermions, the error eventually explodes with increasing $\beta$, but Eq. (16) still holds within the given uncertainty. Since Eq. (16) typically exhibits a smaller variance than the thermodynamic estimator for $K$, this route constitutes the method of choice and the results have been included as an extra column in all data tables as $K_{\textnormal{vir}}$.

III.6 The uniform electron gas

Let us conclude this investigation of the fermion sign problem with a study of the uniform electron gas, which is shown in Fig. 12. The top abscissa corresponds to the blue crosses, which depict the system-size dependence for the UEG at metallic density ($r_{s}=2$) in the warm dense matter regime review , $\theta=0.75$. Note that the density is kept constant by increasing the volume $V$ of the simulation cell when adding more electrons. Therefore, increasing $N$ only mitigates finite-size effects, but does not significantly affect the degree of quantum degeneracy. The dashed blue line depicts an exponential fit according to Eq. (15), which is in excellent agreement with our data points even for surprisingly small system size. Thus, the FSP does indeed constitute an exponential wall in terms of particle number $N$ for the UEG as predicted in Sec. II.2, and the situation becomes only worse for the harmonic confinement. The red squares in the same plot show the decrease of $S$ with the inverse temperature for $N=14$ electrons at $r_{s}=2$ (bottom abscissa). Again, we find an exponential decay with $\beta\sim\theta^{-1}$, and simulations become unfeasible for $\theta\lesssim 0.5$ even for such a comparatively small system size (a typical system size for the UEG are $N=33$ electrons brown_ethan ; groth ; dornheim2 ).

Lastly, we show snapshots from our PIMC simulation of the UEG in Fig. 13 for $N=14$ electrons at $r_{s}=2$ and $\theta=4$ (a), $\theta=1$ (b), and $\theta=0.25$ (c). At the highest temperature, the UEG resembles a semi-classical one-component plasma and the average sign $S\approx 0.77$ is large. Panel (b) depicts a configuration from the interesting transition regime, where $\lambda_{\beta}$ becomes comparable to $\overline{r}$ and fermionic exchange-effects are important, but do not yet dominate. At $\theta=0.25$, the system is fully degenerate, the sign vanishes within the given statistical uncertainty, and standard PIMC simulations are unfeasible.

IV Summary and discussion

In summary, we have presented a comprehensive, hands-on discussion of the fermion sign problem in path integral Monte Carlo simulations of degenerate Fermi systems. In particular, we have investigated the manifestation of the FSP regarding different parameters and have found the following: i) our PIMC data for the average sign $S$ are consistent with an exponential decrease in $S$ with increasing the inverse temperature $\beta$ for all considered system- and interaction-types; ii) while we do find an exponential decrease of $S$ with system size for the uniform electron gas, it decreases even faster for the case of the harmonic trap. This is explained by the increase in the radial density distribution $n(r)$ around the center of the trap, which leads to a higher degree of quantum degeneracy; iii) both the coupling strength $\lambda$ and the interaction-type have a large impact on the manifestation of the FSP. Firstly, there is a transition with decreasing $\lambda$ from the strongly coupled, quasi-classical regime (with $S\approx 1$) to the respective noninteracting limit. Secondly, the short-range dipole interaction leads to a significantly less severe FSP compared to the long-range Coulomb repulsion, as the particles are effectively separated from each other within the PIMC simulation, which makes the formation of permutation-cycles less probable; iv) the increase of the dimensionality from $2D$ to $3D$ in the case of electrons in a harmonic confinement leads to a somewhat less severe FSP, although the scaling with $\beta$ is quite similar.

In addition, we have provided a practical example for the Monte-Carlo sampling of a fermionic observable, and have studied the probability distribution $P(VS/S)$ of a fermionic expectation value in the presence of the sign problem. In the case of a severe FSP, when the relative statistical uncertainty of $S$ is large, $P(VS/S)$ is given by a superposition of a Gaussian and a Lorentzian, which leads to a fat tail at large values and a divergence of the variance. For small errors in $S$, on the other hand, the distribution of the fermionic observable cannot be distinguished from a simple Gaussian, and the fermionic PIMC simulation is quasi-exact.

We hope that our results—both regarding the manifestation of the FSP and the extensive data tables—will aid the future development of new simulation approaches for quantum degenerate, correlated Fermi systems. Moreover, the comparatively less severe manifestation of the FSP in the case of dipole interaction makes ab initio PIMC simulations of ultracold dipolar atoms a promising project for future research, which could allow for unprecedented insights into, e.g., the emergence of pairing and fermionic superfluidity for a strongly correlated system.

Acknowledgments

T.D. wishes to thank Simon Groth and Jan Vorberger for their valuable feedback.

This work was partly funded by the Center of Advanced Systems Understanding (CASUS) which is financed by Germany’s Federal Ministry of Education and Research (BMBF) and by the Saxon Ministry for Science and Art (SMWK) with tax funds on the basis of the budget approved by the Saxon State Parliament.

All calculations were carried out on the clusters hypnos and hemera at Helmholtz-Zentrum Dresden-Rossendorf (HZDR), and at the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN) under grant shp00015.

Appendix A Convergence with imaginary-time slices $P$

Since the operators for the kinetic and potential energy, $\hat{K}$ and $\hat{V}$, do not commute, the canonical density matrix within the PIMC formalism is typically decomposed using a suitable factorization scheme, see Refs. brualla ; sakkos for a detailed discussion. In the present work, we restrict ourselves to the primitive factorization

[TABLE]

with $\epsilon=\beta/P$ being the so-called imaginary-time step, which is justified by the Trotter formula trotter

[TABLE]

Therefore, $P$ constitutes a convergence parameter within our simulations, and the factorization error in the expectation value of an observable $A$ due to Eq. (17) scales as brualla

[TABLE]

In the following, we will investigate the convergence with $P$ for a few representative cases.

In Fig. A.1, we show the convergence of the potential energy (a) and kinetic energy (b) with $P$ for $N=6$ spin-polarized electrons with $\beta=1.3$ and $\lambda=0.5$ in a $2D$ harmonic trap, i.e., a data point from Tab. 3. While this parameter combination does not constitute the lowest temperature considered in this work, it is still a good choice for this convergence study. For lower temperatures, the FSP leads to an exponentially increasing statistical uncertainty, and even large factorization errors cannot be resolved. Still, even at $\beta=1.3$ no factorization error can be resolved within the given error bars. For completeness, we mention that the increasing noise in $K$ towards large $P$ is a direct consequence of the utilized thermodynamic estimator, see Ref. janke for an extensive discussion.

A second degree of freedom worth considering is the interaction strength $\lambda$. In particular, one would expect that, for fixed temperature, the factorization error is most pronounced for intermediate coupling, as the system becomes effectively classical or noninteracting in the limits of $\lambda\gg 1$ and $\lambda\to 0$, respectively. To this end, we consider $N=6$ spin-polarized electrons in a $2D$ harmonic trap at $\lambda=3$ and $\beta=1$ (i.e., a parameter set from Tab. 2) in Fig. A.2. The green crosses correspond to the PIMC results, and the dashed red lines to parabolic fits according to Eq. (19) for $20\leq P\leq 1000$. First and foremost, we do find a significant yet small dependence of our PIMC data on $P$, which is fully consistent with the expected factorization error. Moreover, we note that the data points for $P=100,200,500,$ and $1000$ cannot be distinguished within the given error bars, which means that the results for $P=200$ are indeed quasi-exact.

Lastly, we consider the case of dipole interaction (ultracold atoms) in Fig. A.3 for the same parameters as in Fig. A.2. Again, we find good agreement between the PIMC data and Eq. (19), and $P=200$ are converged within the given statistical uncertainty.