Path Integral Monte Carlo Simulation of Degenerate Electrons: Permutation-Cycle Properties

TL;DR

This paper analyzes permutation-cycle properties in path integral Monte Carlo simulations of degenerate electrons, focusing on finite-size effects and introducing a correlation function to better understand exchange cycles, with implications for improving fermionic PIMC methods.

Contribution

It provides a detailed analysis of permutation-cycle behaviors in fermionic PIMC, including the development of a correlation function and insights into finite-size effects in electron systems.

Findings

Exchange-cycle frequencies do not follow simple exponential laws due to finite size.

Finite-size effects significantly influence permutation-cycle properties.

The permutation-cycle correlation function offers new insights into cycle joint probabilities.

Abstract

Being motivated by the surge of fermionic quantum Monte Carlo simulations at finite temperature, we present a detailed analysis of the permutation-cycle properties of path integral Monte Carlo (PIMC) simulations of degenerate electrons. Particular emphasis is put onto the uniform electron gas in the warm dense matter regime. We carry out PIMC simulations of up to $N=100$ electrons and investigate exchange-cycle frequencies, which are found not to follow any simple exponential law even in the case of ideal fermions due to the finite size of the simulation box. Moreover, we introduce a permutation-cycle correlation function, which allows us to analyse the joint probability to simultaneously find cycles of different lengths within a single configuration. Again, we find that finite-size effects predominate the observed behaviour. Finally, we briefly consider an inhomogeneous system, namely electrons in a $2D$ harmonic trap. We expect our results to be of interest for the further development of fermionic PIMC methods, in particular to alleviate the notorious fermion sign problem.

Equations

$$Z=\text{Tr}\hat{\rho }=\text{Tr}{e}^{-\beta \hat{H}},$$

$$Z=\int \text{d}\mathbf{R}⟨\mathbf{R}|e^{-\beta \hat{H}}|\mathbf{R}⟩=\int \text{d}\mathbf{R}\rho (\mathbf{R},\mathbf{R},\beta ),$$

$$e^{-\beta (\hat{V}+\hat{K})}\ne e^{-\beta \hat{V}}{e}^{-\beta \hat{K}}.$$

$$e^{-\beta \hat{H}}=\prod _{\alpha =0}^{P-1}{e}^{-ϵ\hat{H}},$$

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

$$Z=\int \text{d}\mathbf{R}\text{d}{\mathbf{R}}_1\dots \text{d}{\mathbf{R}}_{P-1}\prod _{\alpha =0}^{P-1}\rho ({\mathbf{R}}_i,{\mathbf{R}}_{i+1},ϵ),$$

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

$$⟨\hat{A}⟩\approx \frac{1}{N_{\text{MC}}}\sum _{\mathbf{X}\in \{\mathbf{X}{\}}_{\text{MC}}}A(\mathbf{X}),$$

$$\Delta A\sim \frac{1}{\sqrt{N_{\text{MC}}}}.$$

$$Z_{\text{B,F}}=\frac{1}{N!}\sum _{\sigma \in S_N}\int \text{d}$$

$$Z^{\prime }=\int \text{d}\mathbf{X}|W_{\text{F}}(\mathbf{X})|=Z_{\text{B}},$$

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

$${⟨\hat{A}⟩}^{\prime }=\frac{1}{Z_{\text{B}}}\int \text{d}\mathbf{X}{W}_{\text{B}}(\mathbf{X})A(\mathbf{X})={⟨\hat{A}⟩}_{\text{B}}.$$

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

$$\Delta A_{\text{F}}\sim \frac{1}{S\sqrt{N_{\text{MC}}}}\sim \frac{e^{\beta N(f_{\text{F}}-f_{\text{B}})}}{\sqrt{N_{\text{MC}}}},$$

$$P(l)=\frac{1}{lN}⟨\sum _{i=1}^N\delta (i,l)⟩,$$

$$P(l)=\frac{1}{N_{\text{MC}}}\sum _{\mathbf{X}\in \{\mathbf{X}{\}}_{\text{MC}}}\left(\frac{1}{N}\sum _{m\in \{m{\}}_{\mathbf{X}}}{\delta }_{m,l}\right),$$

$$P(l)=\frac{Z_{\text{B},1}(l\beta )Z_{\text{B},N-l}(\beta )}{l{Z}_{\text{B},N(\beta )}},$$

$$P(l,k)=\frac{2}{lkN(N-1)}⟨\sum _{i=1}^N\sum _{j>i}^N\delta (i,l)\delta (j,k)⟩,$$

$$P(l,k)=$$

$$P_{\text{u}}(l,k)=P(l)P(k),$$

$$Z_{\text{B},N}(\beta )=\frac{1}{N}\sum _{\eta =1}^N{Z}_{\text{B},1}(\eta \beta )Z_{\text{B},N-\eta }(\beta ),$$

$$Z_{\text{B},1}(\beta )$$

$$E_x$$

$$\sum _x\to \int \text{d}x\Rightarrow Z_{\text{B},1}(\beta )=\frac{V}{{\lambda }_{\beta }^3}.$$

$$\hat{H}=-\frac{{\hslash }^2}{2m}\sum _{k=1}^N{\nabla }_k^2,$$

$${\Phi }_{\mathbf{k}}(\mathbf{r})=⟨\mathbf{r}|\mathbf{k}⟩=\frac{e^{i\mathbf{r}\cdot \mathbf{k}}}{\sqrt{V}},$$

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

$$\delta P(l,k)=\frac{P(l,k)-P(l)P(k)}{P(l,k)}.$$

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

Full text

Path-integral Monte-Carlo, degenerate electrons, Quantum-Monte-Carlo, Fermi gas

I Introduction

The well-known path-integral Monte-Carlo (PIMC) method berne_JCP_1982 ; pollock_PRB_1984 ; ceperley_path_1995 constitutes a highly successful tool for the simulation of distinguishable particles (often referred to as boltzmannons) and bosons. Due to its exact nature, PIMC has been pivotal for the understanding of important physical effects such as superfluidity sindzingre_path-integral_1989 ; kwon_local_2006 ; dornheim_superfluidity_2015 and Bose-Einstein condensation pilati_dilute_2010 ; noauthor_path-integral_2016 and, using state-of-the-art Monte-Carlo sampling techniques like the worm algorithm by Boninsegni et al. boninsegni_worm_2006-1 ; boninsegni_worm_2006 , simulations of up to $N\sim 10^{4}$ particles are possible.

In stark contrast, PIMC simulations of fermions like He3 or electrons are severely limited by the notorious fermion sign problem loh_sign ; doi:10.1063/1.4977920 . This is a direct consequence of the anti-symmetry under particle-exchange, leading to a cancellation of positive and negative terms, which can result in an almost vanishing signal-to-noise ratio. Consequently, PIMC simulations of electrons are restricted to relatively high temperature or strong coupling, but break down when quantum-degeneracy effects become important. This is very unfortunate, as fermionic quantum Monte-Carlo simulations at finite temperature are highly needed for the description of, e.g., ultracold atoms RevModPhys.80.885 ; PhysRevLett.74.2288 ; Bernu1992 , lattice models PhysRevB.72.035122 ; 0295-5075-90-1-10004 ; Cheuk1260 , or even exotic quark-gluon plasmas PhysRevC.87.035207 ; doi:10.1002/ctpp.201400056 .

An application of particular interest is so-called warm dense matter—an extreme state occurring in astrophysical objects like giant planet interiors PhysRevLett.120.115703 ; PhysRevLett.120.025701 ; PhysRevB.75.024206 , on the pathway towards inertial confinement fusion PhysRevB.84.224109 ; PhysRevE.90.033111 , or in state-of-the art experiments with, e.g., free-electron lasers or diamond anvil cells, see Ref. falk_2018, for a topical review. More specifically, the warm dense matter regime is characterized by two parameters that are both of the order of one: 1) the density parameter (or coupling parameter, Wigner-Seitz radiusseitz1943physics ) $r_{s}=\overline{r}/a_{\textnormal{B}}$, with $\overline{r}$ and $a_{\textnormal{B}}$ being the average inter-particle distance and Bohr radius, and 2) the degeneracy temperature $\theta=k_{\textnormal{B}}T/E_{\textnormal{F}}$, where $E_{\textnormal{F}}$ denotes the usual Fermi energy giuliani2005quantum . Therefore, the intricate and nontrivial interplay of a) Coulomb coupling, b) thermal excitation, and c) quantum degeneracy effects renders a thorough theoretical description most challenging graziani2014frontiers , leaving ab initio quantum Monte Carlo methods as the option of choice10.1007/978-3-319-04912-0_5 . Consequently, over the last years, there has been a remarkable spark of new developments in the field of thermodynamic quantum Monte Carlo simulations of electrons at finite temperature dornheim_review ; dornheim_NJP2015 ; doi:10.1063/1.4936145 ; PhysRevLett.117.156403 ; PhysRevLett.119.135001 ; PhysRevB.89.245124 ; doi:10.1063/1.4927434 ; PhysRevLett.117.115701 ; PhysRevLett.115.050603 ; doi:10.1021/acs.jctc.8b00569 ; PhysRevB.95.205109 ; PhysRevLett.115.176403 ; PhysRevE.91.033108 ; universe4120133 ; PhysRevE.96.023203 ; doi:10.1063/1.4999907 ; PhysRevLett.121.255001 . Yet, despite this impressive progress, there still does not exist a single method that is capable to provide an accurate description of a correlated Fermi system for all parameters doi:10.1063/1.4977920 , which makes the further development of state-of-the-art QMC simulations indispensable.

In this article, we present a detailed analysis of the application of the path-integral Monte-Carlo method to a system of correlated electrons at finite temperature. More specifically, we investigate the permutation properties of the uniform electron gas in the warm dense matter regime PhysRevLett.110.146405 ; PhysRevB.93.085102 ; PhysRevB.93.205134 ; doi:10.1002/ctpp.201700096 ; dornheim_review , an important model system, which has been fully described only recently dornheim_review ; PhysRevLett.119.135001 ; doi:10.1002/ctpp.201600082 . In fact, the particular structure of the permutation space within a path-integral Monte-Carlo calculation fully determines the degree of severity of the fermion sign problem, and, therefore, whether simulations are feasible PhysRevA.48.4075 ; PhysRevE.80.066702 . In particular, if permutations of particles were independent of each other, the full fermionic configuration space could be significantly simplified as has been recently proposed in Ref. dubois_permutation, , and PIMC simulations were possible at parameters that are currently out of reach for other methods. The main goal of the present work is to investigate correlation effects (which includes, but is not limited to effects caused by the Coulomb coupling between electrons) within the permutation cycles and, in this way, to assess if further simplifications are possible.

The paper is organised as follows: in Sec. II, we introduce the theoretical background of the PIMC method (II.1), how it needs to be adapted to the simulation of identical particles (II.1.1), and, in the case of fermions, how this leads to the fermion sign problem (II.1.2). Furthermore, in Sec. II.2 we introduce several quantities to measure the desired permutation cycle properties and give a few useful formulas for the noninteracting (ideal) case (II.3). In Sec. III, we discuss our simulation results starting with noninteracting fermions (III.1), which, remarkably, exhibit a quite nontrivial behaviour due to the finite size of the simulation box. The bulk of our results is devoted to the warm dense electron gas (III.2), where we investigate different system sizes, coupling strengths, and temperatures. Further, in Sec. III.3 we extend our considerations to electrons in a $2D$ harmonic trap, where the finite-size effects are an inherent feature of the system. Lastly, in Sec. IV we summarise our findings and discuss the implications on the future development of PIMC simulations of electrons at finite temperature.

II Theory

II.1 Path Integral Monte Carlo

Let us consider a system of $N$ distinguishable particles in a three-dimensional box of length $L$ and volume $V=L^{3}$ at an inverse temperature $\beta=1/k_{\textnormal{B}}T$. In thermodynamic equilibrium, such a system is fully described by the canonical partition function

[TABLE]

with the Hamiltonian $\hat{H}=\hat{V}+\hat{K}$ being the sum of a potential and a kinetic contribution, respectively. In coordinate space, Eq. (1) becomes

[TABLE]

with $\mathbf{R}=(\mathbf{r}_{1},\dots,\mathbf{r}_{N})^{T}$ containing the coordinates of all $N$ particles, and $\rho(\mathbf{R},\mathbf{R}^{\prime},\beta)$ being the thermal density matrix. Unfortunately, Eq. (2) is not directly useful since the matrix elements of $\hat{\rho}$ cannot be readily evaluated as $\hat{V}$ and $\hat{K}$ do not commute,

[TABLE]

To overcome this obstacle, we use a semi-group property of the density operator,

[TABLE]

with $\epsilon=\beta/P$, and, according to the well-known Trotter formula trotter_product_1959 ; de_raedt_applications_1983 , the factorization error from Eq. (3) vanishes in the limit of large $P$,

[TABLE]

In a nutshell, we can express the density-matrix elements at some temperature $T$, $\rho(\mathbf{R},\mathbf{R}^{\prime},\beta)$, as a product of $P$ density-matrices, but at a $P$-times higher temperature, $\rho(\mathbf{R},\mathbf{R}^{\prime},\epsilon)$. This is advantageous, since, for sufficiently large $P$, we can evaluate $\rho(\mathbf{R},\mathbf{R}^{\prime},\epsilon)$ using a suitable high-temperature approximation, with $P$ being a convergence parameter within the PIMC formalism. The partition function thus becomes

[TABLE]

with $\mathbf{R}=\mathbf{R}_{0}=\mathbf{R}_{P}$ due to the definition of the trace. Formally, the density matrices in Eq. (6) are equivalent to propagators in the imaginary time$\tau\in[0,\beta]$ by a time step $\epsilon$, and the trace is then interpreted as the sum over all closed paths $\mathbf{X}$.

This is illustrated in the left panel of Fig. 1, where we show a configuration of $N=3$ particles—each being represented by an entire closed path of $P=6$ sets of coordinates—in the $x$-$\tau$-plane. In this way, the complicated quantum mechanical system of interest has effectively been mapped onto a system of classical ring-polymers chandler_exploiting_1981 . The only seeming drawback of this procedure is the drastic increase in the dimensionality in Eq. (6), which makes the application of standard quadrature schemes unfeasible. Fortunately, this obstacle can be overcome by the utilization of stochastic, i.e., Monte-Carlo methods, which are not afflicted by this curse of dimensionality yellow_book ; liu2013monte . In particular, the basic idea of the path-integral Monte-Carlo method berne_JCP_1982 ; pollock_PRB_1984 ; ceperley_path_1995 is to use the Metropolis algorithm metropolis_equation_1953 to randomly generate a Markov chain of $N_{\textnormal{MC}}$ particle configurations $\mathbf{X}$, $\{\mathbf{X}\}_{\textnormal{MC}}$, distributed according to their respective contribution $W(\mathbf{X})$ to the partition function,

[TABLE]

with $\textnormal{d}\mathbf{X}=\textnormal{d}\mathbf{R}_{0}\dots\textnormal{d}\mathbf{R}_{P-1}$. The expectation value of a given observable $\hat{A}$ is then computed as the sum over these elements,

[TABLE]

with $A(\mathbf{X})$ being the so-called estimator. For completeness, we note that for any finite number of Monte-Carlo samples $N_{\textnormal{MC}}$, Eq. (8) is subject to a statistical uncertainty $\Delta A$, which vanishes as

[TABLE]

Hence, the Monte-Carlo error can—at least in principle—be made arbitrarily small and the PIMC approach is quasi-exact.

A more detailed introduction of the PIMC method including a description of sampling-schemes and the computation of observables is beyond the scope of the present work and the interested reader is referred to, e.g., Refs. ceperley_path_1995, ; boninsegni_worm_2006-1, ; dornheim_review, .

II.1.1 Identical Particles: PIMC simulation of bosons and fermions

Let us next extend our consideration to the simulation of $N$ spin-polarized bosons or fermions. The partition function is then given by ceperley_path_1995

[TABLE]

and incorporates a summation over all elements $\sigma$ of the permutation group $S_{N}$, with $\hat{\pi}_{\sigma}$ being the corresponding permutation operator. The plus and minus sign corresponds to bosons and fermions, respectively, which means that the sign of the configuration weight function $W(\mathbf{X})$ alternates with the number of pair permutations $T$ in the latter case. This renders fermionic PIMC simulations somewhat tricky, as we shall see in Sec. II.1.2.

In practice, Eq. (II.1.1) implies that we now also have to stochastically generate configurations with trajectories where the start- and end-points of a given particle are not identical. This leads to so-called exchange-cycles, which are illustrated in Fig. 1 for a system of $N=3$ identical particles. The left panel shows a configuration without any permutations, such as we have already encountered in the discussion of PIMC simulations of distinguishable particles. Consequently, the corresponding fermionic configuration weight $W_{\textnormal{F}}(\mathbf{X})$ is positive. In contrast, the central panel depicts a configuration with a single pair-exchange, leading to one path containing two particles and one path containing only one. In the case of fermions, the sign of the weight function is negative, $W_{\textnormal{F}}(\mathbf{X})\leq 0$. Finally, in the right panel we show a configuration with two pair-permutations, leading to only a single path, which contains all $N=3$ particles with a positive configuration weight.

As a side remark, we mention that the efficient realization of such exchange-cycles was solved only relatively recently by the worm algorithm introduced in Refs. boninsegni_worm_2006-1, ; boninsegni_worm_2006, . All results shown in the present work have been obtained using a canonical adaption of this approach.

In the case of bosons, such as ultracold atoms PhysRevLett.96.105301 ; PhysRevLett.97.045301 ; PhysRevLett.105.070401 or composite particles like indirect excitons PhysRevB.84.075130 , all contributions to $Z$ in Eq. (II.1.1) are strictly positive and PIMC simulations have been successfully used to study important phenomena such as superfluidity sindzingre_path-integral_1989 ; kwon_local_2006 ; dornheim_superfluidity_2015 , Bose-Einstein condensation pilati_dilute_2010 ; noauthor_path-integral_2016 , and collective excitations filinov_collective_2012 ; filinov_correlation_2016 . Unfortunately, for fermions, such as the ubiquitous electrons, this does not hold true, and PIMC simulations are afflicted with the notorious fermion sign problem loh_sign ; doi:10.1063/1.4977920 .

II.1.2 The fermion sign problem

As we have seen in the previous section, a direct utilization of the Metropolis algorithm to generate a Markov chain of fermionic configurations $\mathbf{X}$ distributed according to their respective weight function $W_{\textnormal{F}}(\mathbf{X})$ is not possible, as probabilities must not be negative. To work around this issue, we can generate a Markov chain of configurations $\mathbf{X}$ distributed according to the absolute value of the fermionic weight function

[TABLE]

which, in the case of standard PIMC, are nothing else than the bosonic weights, as $|W_{\textnormal{F}}(\mathbf{X})|=W_{\textnormal{B}}(\mathbf{X})$. The fermionic expectation value of an observable $\hat{A}$ is then computed as

[TABLE]

with $S(\mathbf{X})=W_{\textnormal{F}}(\mathbf{X})/W_{\textnormal{B}}(\mathbf{X})$ being the sign of a particular configuration $\mathbf{X}$, and the notation $\braket{\dots}^{\prime}$ referring to the expectation value with respect to the bosonic weights

[TABLE]

The denominator in Eq. (12) is the so-called average sign

[TABLE]

which is readily identified as the ratio of the fermionic and bosonic partition function (with $f$ denoting the free energy per particle), and constitutes a measure for the amount cancellation of positive and negative contributions within the PIMC simulation doi:10.1002/ctpp.201800157 ; doi:10.1063/1.4977920 ; dornheim_review . In particular, the statistical uncertainty of Eq. (12) is inversely proportional to $S$,

[TABLE]

and thus exponentially increases both with system size $N$ and towards low temperature. This can only be compensated by increasing the number of Monte Carlo samples as $\sim 1/\sqrt{N_{\textnormal{MC}}}$, which quickly becomes unfeasible as one runs into an exponential wall. This is the origin of the notorious fermion sign problem loh_sign , which limits standard PIMC simulations of electrons to relatively high temperature or strong coupling and has been shown to be $NP$-hard for a certain class of Hamiltonians troyer_sign .

II.2 Permutation Cycle Properties

Evidently, the feasibility of a standard PIMC simulation of electrons depends on the probability of pair-exchanges or, more specifically, on the prevalence of different exchange-cycles. This is illustrated in Fig. 2, where we show snapshots from a PIMC simulation of $N=19$ spin-polarized electrons with $P=100$ high-temperature factors at a metallic density, $r_{s}=2$. The top panel corresponds to a relatively high temperature, $\theta=4$. In this case, the extension of the paths of different particles, which is proportional to the thermal wavelength $\lambda_{\beta}=\sqrt{\hbar^{2}2\pi\beta/m}$, is significantly smaller than the average inter-particle distance $\overline{r}$. Consequently, fermionic exchange-effects are not of paramount importance, and exchange-cycles only seldom occur within a PIMC simulation, which results in an average sign of $S\approx 0.7$. In stark contrast, the bottom panel depicts a snapshot for $\theta=0.5$, which falls into the important warm dense matter regime. In this case, $\lambda_{\beta}$ is comparable to $\overline{r}$, and fermionic exchange-effects predominate. In fact, as we shall see in Sec. III, at these conditions even macroscopic exchange-cycles containing all $N$ particles have a significant weight. Accordingly, configurations with positive and negative configuration weights occur with a similar frequency, the average sign vanishes within the given statistical uncertainty [$S=0.0004(3)$], and the resulting cancellation renders standard PIMC simulations unfeasible in this regime.

In the following section, we will introduce two quantities that allow for a more rigorous characterization of the permutation cycle properties of a PIMC simulation.

The probability to find a permutation cycle of length $l$ can be readily defined as

[TABLE]

with $\delta(i,l)$ vanishing, except when particle $i$ is involved in an exchange-cycle of length $l$. Note that the pre-factor $1/l$ makes us count each cycle of length $l$ only once in the definition of $P(l)$. The PIMC expectation value for Eq. (16) is then expressed as

[TABLE]

with the second sum being carried out over all lengths $m$ in the set of permutation cycles $\{m\}_{\mathbf{X}}$ present in a particular Monte-Carlo configuration $\mathbf{X}$. Note that, even in the case of a fermionic PIMC simulation, Eq. (17) is always computed as a bosonic expectation value, i.e., an observable in the modified configuration space $Z^{\prime}=Z_{\textnormal{B}}$.

For the ideal (i.e., non-interacting) system, Eq. (16) can be computed semi-analytically from different partition functions (see Sec II.3) as krauth_book_2006statistical

[TABLE]

which can be used to check the validity of our implementation.

One of the central questions to be investigated in this work is whether the probability to find a cycle of length $l$ within a PIMC simulation depends on the presence of other permutation cycles or, more specifically, on $\{m\}_{\mathbf{X}}$. To answer this question, we define a permutation-cycle correlation function

[TABLE]

which is evaluated in our PIMC simulations as

[TABLE]

In particular, if the probability for a cycle of length $l$ is independent of the presence of other cycles, Eq. (II.2) should coincide with the uncorrelated quantity

[TABLE]

for all $l$ and $k$.

II.3 Ideal Bose and Fermi gas

To check the validity of our implementation and analyze our PIMC results for correlated electrons, it is useful to consider the ideal system where $P(l)$ can be computed from Eq. (18). The ideal canonical partition function obeys a recursion relation with respect to the system size $N$ of the form krauth_book_2006statistical

[TABLE]

with $Z_{\textnormal{B},0}=1$. The initial condition for Eq. (22) is given by the single particle partition function $Z_{\textnormal{B},1}$, which for a three-dimensional periodic box of length $L$ is given by

[TABLE]

and has to be evaluated at different multiples of the inverse temperature $\beta$. In the limit of a large box (compared to the thermal wavelength, $L\gg\lambda_{\beta}$), the sum in Eq. (23) can be transformed into a continuous integral, and we obtain the simple expression

[TABLE]

In Fig. 3, we show $Z_{\textnormal{B},1}(\eta\beta)$ computed from Eq. (23) for the relevant multiples of $\beta$ for the UEG with $N=100$ (green squares) and $N=19$ (red circles) at $\theta=0.5$. With decreasing temperature (i.e., increasing inverse temperature $\eta\beta$), only the lowest state in Eq. (23) is occupied, and $Z$ converges to one. In addition, the dark grey symbols depict the results for $Z_{\textnormal{B},1}(\eta\beta)$ from the continuous approximation Eq. (24). At large temperature (small $\eta$), replacing the discrete sum by an integral is accurate, but with increasing $\eta$ the thermal wavelength eventually becomes comparable to $L$ and the approximation breaks down. Naturally, this happens at slightly lower temperature for the larger system. Thus, we conclude that a numerical evaluation of the infinite sum, Eq. (23), is essential to accurately compute the permutation-cycle distribution $P(l)$, and to perform a meaningful check of our PIMC data.

III Results

III.1 Ideal Fermi gas

Let us start our investigation of the permutation-cycle properties by considering $N$ ideal spin-polarized fermions in a three-dimensional periodic box of length $L$. The corresponding Hamiltonian is simply given by

[TABLE]

which is readily diagonalised by plane-waves of the form

[TABLE]

with $\mathbf{k}=2\pi L^{-1}(a_{0},a_{1},a_{2})^{T}$ and $a_{i}\in\mathbb{Z}$, leading to the partition function in Eq. (23).

In panel a) of Fig. 4, we show the probability of a single particle to be in a permutation-cycle of length $l$, $P(l)l$, for $N=1000$ (yellow triangles), $N=100$ (green squares), $N=33$ (red circles), and $N=19$ (black crosses) ideal fermions at $\theta=0.5$. For completeness, we mention that the density-parameter has been chosen as $r_{s}=2$ in our PIMC simulations, although the ideal results are independent of $r_{s}$. The symbols correspond to our PIMC data evaluated as a histogram according to Eq. (17), and the continuous blue lines to the semi-analytical result, Eq. (18). First and foremost, we note the perfect agreement between the PIMC data and the exact result over the entire $l$-range, which is a strong verification of our implementation. An additional comparison for different temperatures can be found in Fig. 9. Moreover, the $P(l)l$ data for different $N$ are in agreement with each other for small $l$, but there appear large-sized deviations for $l\gtrsim 4$. These finite-size effects manifest as a steep decrease for large $l$, which is a consequence of the periodicity of the simulation box. More specifically, if we have a permutation-cycle of length $l$ and propose a pair-exchange with a particle that is already involved in the same trajectory, this exchange-cycle is split up into two smaller ones. Therefore, large permutation-cycles are significantly suppressed as compared to the macroscopic system in the thermodynamic limit, $N\to\infty$. Lastly, we mention that $P(l)l$ does not exhibit an exponential decay with $l$ for any depicted system size, even in the case of noninteracting fermions.

In Fig. 4 b), we show PIMC data for the permutation-cycle correlation function $P(1,l)$ [see Eq. (II.2)] for $N=100$ (green squares) and $N=33$ (red circles), which is a measure for the joint probability to find one single particle and one permutation-cycle of length $l$ within the PIMC simulation at the same time. Again, we observe a perfect agreement of both data sets for small $l$, and significant finite-size effects for $l\gtrsim 5$. The dark grey symbols show the corresponding uncorrelated function $P_{\textnormal{u}}(1,l)$ [see Eq. (21) in Sec. II.2], which, for ideal fermions in the thermodynamic limit, should exactly reproduce $P(1,l)$. While this does indeed (approximately) hold for $N=100$ over the entire depicted $l$-range, there appear sizeable deviations for $N=33$ for $l\gtrsim 15$. In fact, it always holds $P(k,l)=0$ for $k+l>N$, which is not the case for the uncorrelated analogue $P_{u}(k,l)$, and the two quantities deviate even for $l+k<N$, cf. Fig. 5. In panels c) and d) of Fig. 4, we show $P(3,l)$ and $P(10,l)$, which exhibit a similar behavior, although the finite-size effects are more pronounced and start for smaller $l$.

For completeness, in Fig. 5 we show the full $k$- and $l$-dependence of $P(k,l)$ (top) and $P_{\textnormal{u}}(k,l)$ (bottom) for the case of $N=33$. In a nutshell, the joint probability to find a permutation-cycle of length $l$ depends on the composition of a configuration $\mathbf{X}$ even for ideal fermions, which is a direct consequence of the finite number of particles.

III.2 Uniform electron gas

Let us next consider the case of the uniform electron gas dornheim_review ; Gill_UEGREVIEW_2011 , which is governed by the Hamiltonian

[TABLE]

with $W(\mathbf{r},\mathbf{r}^{\prime})$ being the pair-interaction. To speed up our calculations, we have replaced the usual Ewald summation fraser_PRB with the spherically-averaged potential by Yakub and Ronchi Yakub_JCP ; Yakub2005 , where the summations in reciprocal and coordinate space can be carried out analytically. In the moderately coupled warm dense matter regime that is of interest in this work, the effect on the permutation-cycle properties is negligible.

In Fig. 6 (top), we show PIMC results for $P(l)l$ for the spin-polarized UEG in the warm dense matter regime at $r_{s}=2$ and $\theta=0.5$ with $N=100$ (green squares), $N=33$ (red circles), $N=19$ (black crosses), and $N=10$ (blue diamonds). First and foremost, we note that the results are qualitatively very similar to the noninteracing data shown in Fig. 4, including the distinct finite-size effects. Moreover, here, too, we do not find an exponential decay of $P(l)$ with the permutation-cycle length $l$. In the bottom panel, we compare PIMC data for the UEG with $N=33$ electrons at $\theta=0.5$ for $r_{s}=2$ (green squares), $r_{s}=10$ (red circles), $r_{s}=20$ (black crosses), and the noninteracting case (blue diamonds). With increasing $r_{s}$, the system becomes more sparse and, consequently, the coupling strength increases giuliani2005quantum . The resulting repulsion-induced inter-particle separation leads to a suppression of exchange-effects. Thus, we find an increased probability of single-particle trajectories, whereas permutation-cycles with $l\geq 2$ occur less often as compared to the noninteracting case.

Let us next consider the permutation-cycle correlation function $P(k,l)$, which is shown in Fig. 7. In the top panel, we investigate the correlation between single-particles and cycles of length $l$, $P(1,l)$, for $N=33$ and $\theta=0.5$ for $r_{s}=2$ (green squares), $r_{s}=20$ (red circles), and the ideal case (blue diamonds). We note that all three curves exhibit a very similar behaviour. The dark grey symbols correspond to the uncorrelated joint probability $P_{\textnormal{u}}(k,l)$, cf. Eq. (21). Remarkably, we find that independent of the coupling strength, $P_{\textnormal{u}}(k,l)$ and $P(k,l)$ are approximately identical for small $l$, which means that the probabilities to find such permutation-cycles are not correlated. In addition, the deviation between the grey and coloured symbols for larger $l$ is qualitatively identical for all three curves. This is a strong indication that the observed correlations in the exchange-cycles are a direct result of the finite number of electrons, and not due to the Coulomb repulsion, even for the strongly coupled case ($r_{s}=20$). In the bottom panel, we show the same information for $P(6,l)$. While the deviations between $P$ and $P_{\textnormal{u}}$ start to appear even at $l=2$, this difference is again already qualitatively present for the noninteracting case, which means that here, too, Coulomb correlations due not constitute the dominating effect.

Finally, in Fig. 8, we show the full $l$- and $k$-dependence of $P(l,k)$ (left column) and $P_{\textnormal{u}}(l,k)$ (right column) for $r_{s}=2$ (top row) and $r_{s}=20$ (bottom row). Firstly, the grey areas in the top-right corner of the $P$-plots directly follow from the fact that $P(l,k)=0$ for $l+k>N$, see the corresponding discussion in Sec. III.1. Apart from this, the uncorrelated function relatively accurately reproduces the real joint probability from our PIMC simulations for small $k$ and $l$, whereas the deviations become more manifest for larger permutation-lengths due to the finite system size.

Let us conclude our analysis of the permutation-cycle properties of PIMC simulations of the warm dense electron gas with an investigation of the temperature-dependence, which is shown in Fig. 9.

In panel a), we show PIMC results for $P(l)l$ for a system of $N=33$ ideal fermions at $\theta=0.125$ (green squares), $\theta=0.5$ (red circles), $\theta=1$ (black crosses), and $\theta=4$ (yellow triangles). Again, the continuous blue lines correspond to the semi-analytical solution [see Eq. (18)] and are in perfect agreement with our PIMC data over the entire $l$-range for all four temperatures. At $\theta=0.125$, the probability of a particle to be included in a permutation-cycle of length $l$ is almost independent of $l$, resulting in a nearly flat curve. In fact, in the ground state limit, $\theta\to 0$, this probability becomes exactly constant, $P(l)l=1/N$, which leads to a vanishing average sign $S$ in a PIMC simulation krauth_book_2006statistical . Therefore, the associated statistical uncertainty, Eq. (15), diverges and PIMC simulations are even theoretically impossible. With increasing temperature, $P(l)l$ exhibits an increasingly steep descent, whereas the probability to find single-particle trajectories is larger. This is a direct consequence of the smaller extension of the single-particle wave function, cf. the discussion of Fig. 2 in Sec. II.2. Moreover, we note that, at $\theta=1$ and $\theta=4$, $P(l)l$ does exhibit an approximately exponential decay.

In panel b) we show the same information, but for an interacting electron gas at $r_{s}=2$, i.e., at a metallic density in the warm dense matter regime. First and foremost, we stress the qualitative similarity to the noninteracting case shown in panel a), in particular for the lowest temperature. In contrast, for $\theta=1$ and $\theta=4$, $P(l)$ does exhibit an exponential decay similar to the ideal case, but with a significantly steeper slope. Therefore, somewhat counter-intuitively, Coulomb-correlation effects are more visible in the permutation-cycle probability at larger temperature, where the system is more weakly coupled.

Finally, in the bottom row of Fig. 9 we compare the permutation-cycle correlation function $P(k,l)$ for different temperatures. In panel c), we show PIMC results for $P(1,l)$ of the UEG at $\theta=0.125$ (green squares) and $\theta=1$ (red circles). In addition, the corresponding dark grey and yellow symbols depict results for $P_{\textnormal{u}}(1,l)$ and PIMC results for $P(1,l)$ for the noninteracting system at the same conditions, respectively. At the lower temperature, all three curves approximately coincide for $l\lesssim 30$, i.e., almost over the entire $l$-range. Only for large permutation-lengths, finite-size effects lead to a deviation towards the uncorrelated joint probability $P_{\textnormal{u}}(1,l)$, whereas the data for the UEG and the noninteracing system are still in agreement. Evidently, the latter fact is a direct consequence of $P(l)l$ being almost flat in both cases. In contrast, at $\theta=1$, $P(1,l)$ and $P_{\textnormal{u}}(1,l)$ are still in good agreement, while the ideal results significantly deviate for all $l$.

Lastly, panel d) shows results for $P(6,l)$ and $P_{\textnormal{u}}(6,l)$ for the same conditions, with similar results albeit distinctly larger finite-size effects, as it is expected.

III.2.1 Relative importance of permutation cycle correlations

To further quantify the importance of correlations between permutation cycles, we consider the relative deviation between $P(l,k)$ and $P_{\textnormal{u}}(l,k)$,

[TABLE]

The results for Eq. (28) are shown in Fig. 10 where we plot $\delta P(1,l)$ (top) and $\delta P(3,l)$ (bottom) for $N=33$ particles with $r_{s}=2$ and $\theta=0.5$. The green squares correspond to the case of ideal fermions and the red circles to the uniform electron gas. In both plots, we observe a remarkably small deviation for small permutation cycles, as it is expected. In contrast, for large cycles ($l\sim 20$) there appear significant relative deviations of several hundred percent. This illustrates that the previously observed qualitative agreement between $P(l,k)$ and $P_{\textnormal{u}}(l,k)$ does not mean that correlation effects of the permutation properties within a PIMC simulation are not important.

III.3 $2D$ quantum dot

Let us conclude this paper by extending our considerations to an inhomogeneous system. For this purpose, we simulate electrons in a two-dimensional harmonic trap, which are governed by the Hamiltonian

[TABLE]

where we have used oscillator units (i.e., distances are given in units of $l_{0}=\sqrt{\hbar/(m\omega)}$ and energies in units of $E_{0}=\hbar\omega$), and $\lambda$ denotes the dimensionless coupling constant. Such a system is often used as a simple model for electrons in a quantum dot reimann_RevModPhys2002 ; dornheim_CPP2016 ; kylanpaa_PhysRevB2017 , and constitutes a useful benchmark for the development of quantum Monte Carlo methods egger_PRL1998 ; schoof_CPP2011 ; dornheim_NJP2015 .

In Fig. 11 a), we show radial density profiles for $N=10$ electrons at a moderate temperature and weak coupling, $\beta=0.5$ and $\lambda=0.5$. The green curve corresponds to the bosonic density measured in the modified configuration space $Z^{\prime}=Z_{\textnormal{B}}$ and the red curve to the corresponding fermioninc density profile that has been computed from Eq. (12). Evidently, bosons tend to cluster more around the centre of the trap, whereas the Pauli blocking acts as an effective repulsive force CARSON199699 , which leads to a more pronounced spatial separation and pushes the electrons outwards. We again stress that both data sets have been obtained within the same PIMC simulation and the differences between bosons and fermions are purely a result of the cancellation of positive and negative terms, see Ref. doi:10.1002/ctpp.201800157, for an extensive topical discussion.

In panel b), we show the corresponding permutation-cycle probabilities $P(l)l$ for $\beta=0.5$ (green squares), $\beta=1$ (red circles), and $\beta=3$ (blue diamonds). Note that the average signs are given by $S=0.0615(1)$, $S=0.00022(4)$, and $S=-0.00003(3)$, respectively, which means that PIMC simulations are only feasible at $\beta=0.5$. In contrast, for the lowest considered temperature, the sign vanishes within the given statistical uncertainty. For the largest temperature, we find a smooth decay with $l$, whereas for $\beta=3$ the probability to find a particle in an exchange-cycle of length $l$ is again almost constant, see also the description of Fig. 9 for a more extensive discussion.

Finally, in panels c) and d), we show our PIMC results for the joint probability $P(1,l)$ and $P(3,l)$, and compare them to the corresponding uncorrelated function $P_{\textnormal{u}}(l,k)$ (dark grey symbols). Similar to the previously discussed case of the warm dense electron gas, we find that correlations between exchange-cycles are more pronounced for higher temperature, and increase for larger cycle-lengths $l$. However, in contrast to the UEG, in the case of an inhomogeneous system in an external potential, any finite-size effects have real physical meaning and are not merely an artifact due to the inevitably limited finite simulation box within a PIMC simulation. Therefore, approximation schemes that rely on the absence of correlations between different permutation-cycles within a given configuration as proposed in Ref. dubois_permutation, are not applicable in this case.

IV Summary and discussion

In summary, we have provided a detailed analysis of the permutation-cycle properties of path-integral Monte-Carlo simulations of correlated, quantum degenerate electrons. To verify our implementation, and to further analyse our results for the warm dense electron gas, we have begun our investigations by considering $N$ ideal spin-polarised fermions in periodic box. Even for this most simple case, we find an exponential decay in the permutation-cycle distribution function $P(l)$ only at relatively high temperatures, $\theta\gtrsim 1$. In addition, this quantity exhibits severe finite-size effects, which increase with the cycle length $l$. Similarly, the corresponding analysis of the permutation-cycle correlation-function $P(l,k)$ has revealed that cycles of different lengths are independent only for small $l$, whereas correlations increase for larger $l$, again due to the finite number of particles within the simulation. In contrast, permutation-cycles of noninteracting fermions are completely independent in the thermodynamic limit.

Extending our considerations to the warm dense electron gas, here, too, we have found that exchange-cycles are almost independent for small $l$, and deviations between $P(l,k)$ and $P_{\textnormal{u}}(k,l)=P(k)P(l)$ only start to become important for long cycles. Due to the striking similarity to the previously investigated ideal case, we have concluded that correlation effects of permutation-cycle properties within a PIMC simulation are dominated by finite-size effects and not, as it might have been expected, by the Coulomb interaction. In addition, we have also investigated the temperature-dependence of these properties, which, again, has revealed a remarkable behaviour: at low temperatures, where the system is most strongly correlated, the permutation properties of the electron gas most closely resemble the exact behavior of ideal electrons; at high temperatures, where correlation effects are expected to be less important, the ideal and interacting system exhibit more pronounced deviations. The explanation for this rather peculiar finding is the ground-state limit for the permutation-cycle distribution function, $P(l)l=1/N$, which is equal both in the interacting and noninteracting cases.

In a nutshell, we have found that the permutation-cycle properties of interacting electrons are qualitatively very similar to the noninteracting case. In particular, cycles of different lengths are surprisingly uncorrelated, which can potentially be exploited to significantly simplify the fermionic configuration space dubois_permutation , and to extend fermionic PIMC simulations to lower temperatures and stronger quantum degeneracy. Moreover, the observed correlation effects in $P(l,k)$ for large $l$ and $k$ appear to be an artifact due to the finite number of electrons within a given simulation. Therefore, ignoring them might not necessarily be a shortcoming, but rather a handy way to mitigate finite-size effects PhysRevB.78.125106 ; PhysRevB.94.035126 .

Finally, we have further extended our considerations to a finite system in an external harmonic potential. In this case, however, the observed finite-size and correlation effects in the permutation-cycle properties are not of an artificial nature, but instead represent real physical behaviour. Therefore, neglecting correlations between cycles is not appropriate, even for noninteracting fermions.

We expect our results to be of interest for the further development of PIMC simulations of fermions, which, in turn is of high importance for many fields like the simulation of ultracold atoms RevModPhys.80.885 ; PhysRevLett.74.2288 ; Bernu1992 , quark-gluon plasmas PhysRevC.87.035207 ; doi:10.1002/ctpp.201400056 , and warm dense matter dornheim_review ; 10.1007/978-3-319-04912-0_5 .

Acknowledgments

We thank B. Alder for discussions on the topic presented in Ref. dubois_permutation, . This work has been supported by the Deutsche Forschungsgemeinschaft via project BO1366/13 and by the Norddeutscher Verbund für Hoch- und Höchleistungsrechnen (HLRN) via grant shp00015 for CPU time.