Permutation blocking path integral Monte Carlo: a highly efficient approach to the simulation of strongly degenerate non-ideal fermions

We consider the canonical ensemble (the particle number N, volume V and inverse temperature $\beta =1/{k}_{{\rm{B}}}T$ are fixed) and write the partition function in coordinate representation as

$$\begin{eqnarray}&&Z=\displaystyle \frac{1}{N!}\displaystyle \sum _{\sigma \in {S}_{N}}{\rm{sgn}}(\sigma )\int {\rm{d}}{\bf{R}}\ \langle {\bf{R}}| {{\rm{e}}}^{-\beta \hat{H}}| {\hat{\pi }}_{\sigma }{\bf{R}}\rangle ,\end{eqnarray} \tag{ 1 }$$

where ${\bf{R}}=\{{{\bf{r}}}_{1},...,{{\bf{r}}}_{N}\}$ contains the coordinates of all particles and ${\hat{\pi }}_{\sigma }$ denotes the exchange operator corresponding to a particular element σ from the permutation group S_N. The Hamiltonian is given by the sum of the kinetic ($\hat{K}$) and potential ($\hat{V}$) energy, $\hat{H}=\hat{K}+\hat{V}$. For the next step, we use the group property of the density operator

$$\begin{eqnarray}&&\hat{\rho }={{\rm{e}}}^{-\beta \hat{H}}=\displaystyle \prod _{\alpha =0}^{P-1}{{\rm{e}}}^{-\epsilon \hat{H}},\end{eqnarray} \tag{ 2 }$$

with $\epsilon =\beta /P$, and insert $P-1$ unities of the form $\hat{1}=\int {\rm{d}}{{\bf{R}}}_{\alpha }\ | {{\bf{R}}}_{\alpha }\rangle \langle {{\bf{R}}}_{\alpha }|$. This gives

$$\begin{eqnarray}&&Z=\int {\rm{d}}{{\bf{R}}}_{0}...{\rm{d}}{{\bf{R}}}_{P-1}\displaystyle \prod _{\alpha =0}^{P-1}\left(\displaystyle \frac{1}{N!}\displaystyle \sum _{\sigma \in {S}_{N}}{\rm{sgn}}(\sigma )\langle {{\bf{R}}}_{\alpha }| {{\rm{e}}}^{-\epsilon \hat{H}}| {\hat{\pi }}_{\sigma }{{\bf{R}}}_{\alpha +1}\rangle \right).\end{eqnarray} \tag{ 3 }$$

Note that we have exploited the permutation operator's idempotency property in equation (3) to introduce antisymmetry on all P imaginary time slices. Following Sakkos et al [38], we introduce the factorization from [37],

$$\begin{eqnarray}&&{{\rm{e}}}^{-\epsilon \hat{H}}\approx {{\rm{e}}}^{-{v}_{1}\epsilon {\hat{W}}_{{a}_{1}}}{{\rm{e}}}^{-{t}_{1}\epsilon \hat{K}}{{\rm{e}}}^{-{v}_{2}\epsilon {\hat{W}}_{1-2{a}_{1}}}{{\rm{e}}}^{-{t}_{1}\epsilon \hat{K}}{{\rm{e}}}^{-{v}_{1}\epsilon {\hat{W}}_{{a}_{1}}}{{\rm{e}}}^{-2{t}_{0}\epsilon \hat{K}},\end{eqnarray} \tag{ 4 }$$

for each of the exponential functions in equation (3). By including double commutator terms of the form

$$\begin{eqnarray}&&[[\hat{V},\hat{K}],\hat{V}]=\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}\displaystyle \sum _{i=1}^{N}{| {{\bf{F}}}_{i}| }^{2},\end{eqnarray} \tag{ 5 }$$

we have to evaluate the total force on each particle, ${{\bf{F}}}_{i}=-{\nabla }_{i}V({\bf{R}})$, and equation (4) is accurate to fourth order in . The explicit form of the modified potential terms $\hat{W}$ is given by

$$\begin{eqnarray}{\hat{W}}_{{a}_{1}} & = & \hat{V}+\displaystyle \frac{{u}_{0}}{{v}_{1}}{a}_{1}{\epsilon }^{2}\left(\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}\displaystyle \sum _{i=1}^{N}{| {{\bf{F}}}_{i}| }^{2}\right)\quad {\rm{and}}\\ {\hat{W}}_{1-{a}_{1}} & = & \hat{V}+\displaystyle \frac{{u}_{0}}{{v}_{2}}(1-{a}_{1}){\epsilon }^{2}\left(\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}\displaystyle \sum _{i=1}^{N}{| {{\bf{F}}}_{i}| }^{2}\right).\end{eqnarray} \tag{ 6 }$$

There are two free parameters in equation (4), namely $0\leqslant {a}_{1}\leqslant 1$, which controls the relative weight of the forces on a particular slice, and $0\leqslant {t}_{0}\leqslant (1-1/\sqrt{3})/2$, which determines the ratio of the, in general, non-equidistant time steps between 'daughter' slices, cf figure 2. All other factors are calculated from these choices:

$$\begin{eqnarray}{u}_{0} & = & \displaystyle \frac{1}{12}\left(1-\displaystyle \frac{1}{1-2{t}_{0}}+\displaystyle \frac{1}{6{(1-2{t}_{0})}^{3}}\right),\\ {v}_{1} & = & \displaystyle \frac{1}{6{(1-2{t}_{0})}^{2}},\\ {v}_{2} & = & 1-2{v}_{1}\quad {\rm{and}}\\ {t}_{1} & = & \displaystyle \frac{1}{2}-{t}_{0}.\end{eqnarray} \tag{ 7 }$$

Zoom In Zoom Out Reset image size

The fourth-order approximation of the imaginary time propagator ${{\rm{e}}}^{-\epsilon \hat{H}}$ is visualized in figure 2. The inverse temperature β has been split into P = 4 intervals of length , which are further divided into three, in general, non-equidistant sub-intervals. Thus, for each main 'bead' ${\tau }_{\alpha }$, there exist two daughter beads, ${\tau }_{\alpha A}$ and ${\tau }_{\alpha B}$.

Let us for a moment ignore the antisymmetry in equation (3) and evaluate the imaginary time propagator in a straightforward way [38]:

$$\begin{eqnarray}\langle {{\bf{R}}}_{\alpha }| {{\rm{e}}}^{-\epsilon \hat{H}}| {{\bf{R}}}_{\alpha +1}\rangle & = & \displaystyle \int {\rm{d}}{{\bf{R}}}_{\alpha A}{\rm{d}}{{\bf{R}}}_{\alpha B}\ \left[{{\rm{e}}}^{-\epsilon {\tilde{V}}_{\alpha }}{{\rm{e}}}^{-{u}_{0}{\epsilon }^{3}\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}{\tilde{F}}_{\alpha }}\right.\\ & & \left.\displaystyle \prod _{i=1}^{N}{\rho }_{\alpha }(i,i){\rho }_{\alpha ,A}(i,i){\rho }_{\alpha B}(i,i)\right],\end{eqnarray} \tag{ 8 }$$

with the definitions of the potential terms

$$\begin{eqnarray}{\tilde{V}}_{\alpha } & = & {v}_{1}V({{\bf{R}}}_{\alpha })+{v}_{2}V({{\bf{R}}}_{\alpha A})+{v}_{1}V({{\bf{R}}}_{\alpha B}),\\ {\tilde{F}}_{\alpha } & = & \displaystyle \sum _{i=1}^{N}\left({a}_{1}{| {{\bf{F}}}_{\alpha ,i}| }^{2}+(1-2{a}_{1}){| {{\bf{F}}}_{\alpha A,i}| }^{2}+{a}_{1}{| {{\bf{F}}}_{\alpha B,i}| }^{2}\right),\end{eqnarray} \tag{ 9 }$$

and the diffusion matrices

$$\begin{eqnarray}{\rho }_{\alpha }(i,j) & = & {\lambda }_{{t}_{1}\epsilon }^{-D}{\rm{exp}}\left(-\displaystyle \frac{\pi }{{\lambda }_{{t}_{1}\epsilon }^{2}}{({{\bf{r}}}_{\alpha ,j}-{{\bf{r}}}_{\alpha A,i})}^{2}\right),\\ {\rho }_{\alpha A}(i,j) & = & {\lambda }_{{t}_{1}\epsilon }^{-D}{\rm{exp}}\left(-\displaystyle \frac{\pi }{{\lambda }_{{t}_{1}\epsilon }^{2}}{({{\bf{r}}}_{\alpha A,j}-{{\bf{r}}}_{\alpha B,i})}^{2}\right),\\ {\rho }_{\alpha B}(i,j) & = & {\lambda }_{2{t}_{0}\epsilon }^{-D}{\rm{exp}}\left(-\displaystyle \frac{\pi }{{\lambda }_{2{t}_{0}\epsilon }^{2}}{({{\bf{r}}}_{\alpha B,j}-{{\bf{r}}}_{\alpha +1,i})}^{2}\right),\end{eqnarray} \tag{ 10 }$$

where ${\lambda }_{\beta }$ denotes the thermal wavelength ${\lambda }_{\beta }^{2}=2\pi {{\rm{\hslash }}}^{2}\beta /m$ and D is the dimensionality of the system. Thus, the matrix elements of equation (10) are equal to the free particle density matrix, ${\rho }_{\alpha }(i,j)={\rho }_{0}({{\bf{r}}}_{\alpha ,j},{{\bf{r}}}_{\alpha A,i},{t}_{1}\epsilon )$. The permutation operator commutes with both $\hat{\rho }$ and $\hat{H}$ and we are, therefore, allowed to artificially introduce the antisymmetrization between all $3P$ slices without changing the result. This transforms equation (8) to

$$\begin{eqnarray}&&\displaystyle \frac{1}{N!}\displaystyle \sum _{\sigma \in {S}_{N}}{\rm{sgn}}(\sigma )\langle {{\bf{R}}}_{\alpha }| {{\rm{e}}}^{-\epsilon \hat{H}}| {\hat{\pi }}_{\sigma }{{\bf{R}}}_{\alpha +1}\rangle \\ &&\quad ={\left(\displaystyle \frac{1}{N!}\right)}^{3}\displaystyle \int {\rm{d}}{{\bf{R}}}_{\alpha A}{\rm{d}}{{\bf{R}}}_{\alpha B}\left[{{\rm{e}}}^{-\epsilon {\tilde{V}}_{\alpha }}{{\rm{e}}}^{-{\epsilon }^{3}{u}_{0}\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}{\tilde{F}}_{\alpha }}{\rm{det}}({\rho }_{\alpha }){\rm{det}}({\rho }_{\alpha A}){\rm{det}}({\rho }_{\alpha B})\right].\end{eqnarray} \tag{ 11 }$$

Finally, this gives the partition function

$$\begin{eqnarray}&&Z=\displaystyle \frac{1}{{(N!)}^{3P}}\int {\rm{d}}{\bf{X}}\displaystyle \prod _{\alpha =0}^{P-1}{{\rm{e}}}^{-\epsilon {\tilde{V}}_{\alpha }}{{\rm{e}}}^{-{\epsilon }^{3}{u}_{0}\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}{\tilde{F}}_{\alpha }}{\rm{det}}({\rho }_{\alpha }){\rm{det}}({\rho }_{\alpha A}){\rm{det}}({\rho }_{\alpha B}),\end{eqnarray} \tag{ 12 }$$

and the integration is carried out over all coordinates on all $3P$ slices:

$$\begin{eqnarray}&&{\rm{d}}{\bf{X}}={\rm{d}}{{\bf{R}}}_{0}...{\rm{d}}{{\bf{R}}}_{P-1}{\rm{d}}{{\bf{R}}}_{0A}...{\rm{d}}{{\bf{R}}}_{P-1A}{\rm{d}}{{\bf{R}}}_{0B}...{\rm{d}}{{\bf{R}}}_{P-1B}.\end{eqnarray} \tag{ 13 }$$

The benefits of the partition function equation (12) are illustrated in the right panel of figure 2 where the beads of two particles are plotted in the τ–x-plane. In the usual PIMC formulation (without the determinants), each of the particles would correspond to a single closed trajectory as visualized by the blue and red connections. To take into account the antisymmetry of fermions, one would also need to sample all configurations with the same positions of the individual beads but different connections between adjacent time slices, which have both positive and negative weights. By indroducing determinants between all slices, we include all $N!$ possible connections between beads on adjacent slices (the green lines) into a single configuration weight and the usual interpretation of mapping a quantum system onto an ensemble of interacting ringpolymers [41] is no longer appropriate. Therefore, a large number of sign changes, due to different permutations, are grouped together resulting in an efficient compensation of many terms (blocking), and the average sign (cf equation (22)) in our simulations is significantly increased [31].

The total energy E follows from the partition function via the familiar relation

$$\begin{eqnarray}&&E=-\displaystyle \frac{1}{Z}\displaystyle \frac{\partial Z}{\partial \beta }.\end{eqnarray} \tag{ 14 }$$

Substituting the expression from equation (12) into (14) and performing a lengthy but straightforward calculation gives the final result for the thermodynamic (TD) estimator

$$\begin{eqnarray}E & = & \displaystyle \frac{3{DN}}{2\epsilon }-\displaystyle \sum _{k=0}^{P-1}\displaystyle \sum _{\kappa =1}^{N}\displaystyle \sum _{\xi =1}^{N}\left(\displaystyle \frac{\pi {\Psi }_{\kappa \xi }^{k}}{\epsilon P{\lambda }_{{t}_{1}\epsilon }^{2}}{({{\bf{r}}}_{k,\kappa }-{{\bf{r}}}_{{kA},\xi })}^{2}\right.\\ & & +\left.\displaystyle \frac{\pi {\Psi }_{\kappa \xi }^{{kA}}}{\epsilon P{\lambda }_{{t}_{1}\epsilon }^{2}}{({{\bf{r}}}_{{kA},\kappa }-{{\bf{r}}}_{{kB},\xi })}^{2}+\displaystyle \frac{\pi {\Psi }_{\kappa \xi }^{{kB}}}{\epsilon P{\lambda }_{2{t}_{0}\epsilon }^{2}}{({{\bf{r}}}_{{kB},\kappa }-{{\bf{r}}}_{k+1,\xi })}^{2}\right)\\ & & +\displaystyle \frac{1}{P}\displaystyle \sum _{k=0}^{P-1}\left({\tilde{V}}_{k}+3{\epsilon }^{2}{u}_{0}\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}{\tilde{F}}_{k}\right),\end{eqnarray} \tag{ 15 }$$

with the definitions

$$\begin{eqnarray}{\Psi }_{\kappa \xi }^{k} & = & {\left({\rho }_{k}^{-1}\right)}_{\kappa \xi }{\left({\rho }_{k}\right)}_{\xi \kappa }\\ {\Psi }_{\kappa \xi }^{{kA}} & = & {\left({\rho }_{{kA}}^{-1}\right)}_{\kappa \xi }{\left({\rho }_{{kA}}\right)}_{\xi \kappa }\\ {\Psi }_{\kappa \xi }^{{kB}} & = & {\left({\rho }_{{kB}}^{-1}\right)}_{\kappa \xi }{\left({\rho }_{{kB}}\right)}_{\xi \kappa }.\end{eqnarray} \tag{ 16 }$$

To split the total energy into a kinetic and a potential part, we evaluate

$$\begin{eqnarray}&&K=\displaystyle \frac{m}{\beta Z}\displaystyle \frac{\partial }{\partial m}Z,\end{eqnarray} \tag{ 17 }$$

and find the TD estimator of the kinetic energy

$$\begin{eqnarray}K & = & \displaystyle \frac{3{ND}}{2\epsilon }-\displaystyle \sum _{k=0}^{P-1}\displaystyle \sum _{\kappa =1}^{N}\displaystyle \sum _{\xi =1}^{N}\left[\displaystyle \frac{\pi {\Psi }_{\kappa \xi }^{k}}{\epsilon P{\lambda }_{{t}_{1}\epsilon }^{2}}{({{\bf{r}}}_{k,\kappa }-{{\bf{r}}}_{{kA},\xi })}^{2}+\displaystyle \frac{\pi {\Psi }_{\kappa \xi }^{{kA}}}{\epsilon P{\lambda }_{{t}_{1}\epsilon }^{2}}{({{\bf{r}}}_{{kA},\kappa }-{{\bf{r}}}_{{kB},\xi })}^{2}\right.\\ & & +\left.\displaystyle \frac{\pi {\Psi }_{\kappa \xi }^{{kB}}}{\epsilon P{\lambda }_{2{t}_{0}\epsilon }^{2}}{({{\bf{r}}}_{{kB},\kappa }-{{\bf{r}}}_{k+1,\xi })}^{2}\right]+\displaystyle \frac{1}{P}\displaystyle \sum _{k=0}^{P-1}\left({\epsilon }^{2}{u}_{0}\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}{\tilde{F}}_{k}\right).\end{eqnarray} \tag{ 18 }$$

Thus, the estimator of the potential energy is given by

$$\begin{eqnarray}&&V=E-K=\displaystyle \frac{1}{P}\displaystyle \sum _{k=0}^{P-1}\left({\tilde{V}}_{k}+2{\epsilon }^{2}{u}_{0}\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}{\tilde{F}}_{k}\right).\end{eqnarray} \tag{ 19 }$$

We notice that the forces contribute to both the kinetic and the potential energy. For completeness, we mention that, for an increasing number of propagators, $P\to \infty$, the first and second terms in equation (15) diverge, which leads to a growing variance and, therefore, statistical uncertainty of both E and K. To avoid this problem, one might derive a virial estimator, e.g. [42], which requires the evaluation of the derivative of the potential terms instead. However, since we are explicitly interested in performing simulations with few propagators to relieve the fermion sign problem, the estimator from equation (15) is sufficient.

To take advantage of the main benefits from the usual continuous space WA, we will temporarily construct artificial trajectories and sample new beads according to standard PIMC techniques, e.g. [43]. The initial situation for our considerations is illustrated in the left panel of figure 4, where a pre-existing trajectory (pink curve) with four missing beads in the middle is shown in the τ–x-plane. We choose the sampling probability to close the configuration as

$$\begin{eqnarray}&&{T}_{{\rm{sample}}}=\displaystyle \frac{{\displaystyle \prod }_{i=0}^{M-1}{\rho }_{0}({{\bf{r}}}_{i},{{\bf{r}}}_{i+1},{\tau }_{i+1}-{\tau }_{i})}{{\rho }_{0}({{\bf{r}}}_{0},{{\bf{r}}}_{M},{\tau }_{M}-{\tau }_{0})},\end{eqnarray} \tag{ 24 }$$

which results in the consecutive generation of $M-1$ new coordinates ${{\bf{r}}}_{i}$, $i\in [1,M-1]$, according to

$$\begin{eqnarray}P({{\bf{r}}}_{i}) & = & \displaystyle \frac{{\rho }_{0}({{\bf{r}}}_{i-1},{{\bf{r}}}_{i},{\tau }_{i}-{\tau }_{i-1}){\rho }_{0}({{\bf{r}}}_{i},{{\bf{r}}}_{M},{\tau }_{M}-{\tau }_{i})}{{\rho }_{0}({{\bf{r}}}_{i-1},{{\bf{r}}}_{M},{\tau }_{M}-{\tau }_{i-1})}\\ & = & {\left(\displaystyle \frac{1}{\sqrt{2\pi {\sigma }_{i}^{2}}}\right)}^{D}{\rm{exp}}\left(-\displaystyle \frac{{({{\bf{r}}}_{i}-{{\boldsymbol{\xi }}}_{i})}^{2}}{2{\sigma }_{i}^{2}}\right),\end{eqnarray} \tag{ 25 }$$

which is a Gaussian (cf the blue curves in figure 4) with the variance

$$\begin{eqnarray}&&{\sigma }_{i}^{2}=\displaystyle \frac{{{\rm{\hslash }}}^{2}}{m}\displaystyle \frac{({\tau }_{i}-{\tau }_{i-1})({\tau }_{M}-{\tau }_{i})}{{\tau }_{M}-{\tau }_{i-1}},\end{eqnarray} \tag{ 26 }$$

around the intersection of the connection between the previous coordinate, ${{\bf{r}}}_{i-1}$, with the end point ${{\bf{r}}}_{M}$ and the time slice ${\tau }_{i}$

$$\begin{eqnarray}&&{{\boldsymbol{\xi }}}_{i}=\displaystyle \frac{{\tau }_{M}-{\tau }_{i}}{{\tau }_{M}-{\tau }_{i-1}}{{\bf{r}}}_{i-1}+\displaystyle \frac{{\tau }_{i}-{\tau }_{i-1}}{{\tau }_{M}-{\tau }_{i-1}}{{\bf{r}}}_{M}.\end{eqnarray} \tag{ 27 }$$

Zoom In Zoom Out Reset image size

In the usual WA-PIMC, the configuration space is defined by the Matsubara Green function (e.g. [44]) which implies that the algorithm does not only allow for the change of the particle number N (grand canonical ensemble) but, in addition, requires the generation of configurations with a single open path, the so-called worm. However, in the PB-PIMC configuration space defined by equation (12), there are no trajectories and, therefore, no direct realization of a worm is possible. Instead, we consider an extended ensemble, which combines closed configurations with a total of $3{NP}$ beads and open configurations, where on some consecutive time slices the number of beads is reduced by one, to $N-1$. Such a configuration is illustrated in the right panel of figure 4. There are two special beads which are denoted as 'head' and 'tail' and the triangles, circles and squares symbolize beads from three different particles. There are eight beads from different particles missing (indicated by the empty symbols at the right boundary) between ${\tau }_{{\rm{head}}}={\tau }_{2A}$ and ${\tau }_{{\rm{tail}}}={\tau }_{1A}$, going forward in imaginary time.

For most slices, the computation of the diffusion matrix allows for no degree of freedom in the extended ensemble. We define the latter in a way, that the head bead does not serve as a starting point for the elements but is treated as if it was missing. This is justified because, otherwise, there does not necessarily exist a large matrix element in this particular row because no artificial connection has been sampled on the next slice. For the configuration from figure 4, the diffusion matrix of the head's time slice is given by

$$\begin{eqnarray}{\rho }_{2A} & = & \left(\begin{array}{rcl}{\rho }_{0}({{\bf{r}}}_{\mathrm{1,2}A},{{\bf{r}}}_{\mathrm{1,2}B},{t}_{1}\epsilon ) & {\rho }_{0}({{\bf{r}}}_{\mathrm{1,2}A},{{\bf{r}}}_{\mathrm{2,2}B},{t}_{1}\epsilon ) & 0\\ 1 & 1 & 1\\ {\rho }_{0}({{\bf{r}}}_{\mathrm{3,2}A},{{\bf{r}}}_{\mathrm{1,2}B},{t}_{1}\epsilon ) & {\rho }_{0}({{\bf{r}}}_{\mathrm{3,2}A},{{\bf{r}}}_{\mathrm{2,2}B},{t}_{1}\epsilon ) & 0\end{array}\right)\\ \Rightarrow {\rm{det}}({\rho }_{2A}) & = & {\rm{det}}\left(\begin{array}{cc}{\rho }_{0}({{\bf{r}}}_{\mathrm{1,2}A},{{\bf{r}}}_{\mathrm{1,2}B},{t}_{1}\epsilon ) & {\rho }_{0}({{\bf{r}}}_{\mathrm{1,2}A},{{\bf{r}}}_{\mathrm{2,2}B},{t}_{1}\epsilon )\\ {\rho }_{0}({{\bf{r}}}_{\mathrm{3,2}A},{{\bf{r}}}_{\mathrm{1,2}B},{t}_{1}\epsilon ) & {\rho }_{0}({{\bf{r}}}_{\mathrm{3,2}A},{{\bf{r}}}_{\mathrm{2,2}B},{t}_{1}\epsilon )\end{array}\right).\end{eqnarray} \tag{ 28 }$$

All diffusion matrices with $N-1$ beads on their slices are computed in the same way. The other degree of freedom for which the extended ensemble allows is the choice whether the tail will be included as the final coordinate in the diffusion matrix or not. Here, it makes sense to allow for this possibility, because there does exist at least a single large element in this particular row anyway. The corresponding matrix for the configuration from figure 4 looks like

$$\begin{eqnarray}{\rho }_{2A}=\left(\begin{array}{ccc}{\rho }_{0}({{\bf{r}}}_{\mathrm{1,1}},{{\bf{r}}}_{\mathrm{1,1}A},{t}_{1}\epsilon ) & {\rho }_{0}({{\bf{r}}}_{\mathrm{1,1}},{{\bf{r}}}_{\mathrm{2,1}A},{t}_{1}\epsilon ) & {\rho }_{0}({{\bf{r}}}_{\mathrm{1,1}},{{\bf{r}}}_{\mathrm{3,1}A},{t}_{1}\epsilon )\\ {\rho }_{0}({{\bf{r}}}_{\mathrm{2,1}},{{\bf{r}}}_{\mathrm{1,1}A},{t}_{1}\epsilon ) & {\rho }_{0}({{\bf{r}}}_{\mathrm{2,1}},{{\bf{r}}}_{\mathrm{2,1}A},{t}_{1}\epsilon ) & {\rho }_{0}({{\bf{r}}}_{\mathrm{2,1}},{{\bf{r}}}_{\mathrm{3,1}A},{t}_{1}\epsilon )\\ 1 & 1 & 1\end{array}\right).\end{eqnarray} \tag{ 29 }$$

However, we emphasize that the particular choice of the extended ensemble does not influence the extracted canonical expectation values as long as detailed balance is fulfilled in all updates. We have developed a simulation scheme which consists of four different types of moves that ensure detailed balance and ergodicity. The updates are presented in detail in the appendix.

We start the discussion of the simulation results by investigating the effects of the two free parameters a₁ and t₀ on the convergence of two different observables, namely the energy E and radial density $n(r)$.

In figure 5, results are summarized for N = 4 electrons with $\lambda =1.3$ and $\beta =5$, i.e., moderate coupling and low temperature, and panel (A) shows the convergence of the total energy as a function of the inverse number of propagators which is proportional to the imaginary time step, $\epsilon \propto 1/P$. The red diamonds [(a) ${t}_{0}=0.04$, ${a}_{1}=0.0$] and blue circles [(b) ${t}_{0}=0.13$, ${a}_{1}=0.33$] denote two different combinations of free parameters and exhibit a clearly different convergence behavior towards the exact result known from CPIMC, i.e., the black line. For P = 2, the energy with parameter set (a) is too low by almost one percent. With increasing P, E increases and reaches a maximum around P = 5, until the curves approach the exact energy from above. For parameter set (b), the energy converges monotonically from above and, even for P = 2, the deviation from the CPIMC result is as small as $0.2\%$. The selected energies which are listed in table 1 reveal that the total energy is converged for P = 14 within the statistical uncertainty. For the panels (C) and (D), the energy has been split into a potential (V) and kinetic (K) contribution. For both parameter combinations, V converges monotonically, although from different directions. In addition, parameter set (b) gives a much better result for small P. Panel (D) reveals, that the kinetic energy K is responsible for the non-monotonous convergence of E for parameter set (a), which again delivers worse results for P = 2, as compared to the blue circles. Finally, panel (B) shows the average sign S as a function of $1/P$. Both curves exhibit a similar decrease with an increasing number of propagators, as it is expected. However, parameter set (a) always allows for a better sign than (b). The reason for this behavior is the free parameter t₀, which controls the relative spacing between the three time slices of an imaginary time step . For ${t}_{0}=0.04$, there are a single small and two large steps. The latter allow for more blocking, since the corresponding decay length ${\lambda }_{{t}_{1}\epsilon }$ in the diffusion matrices is large as well. For ${t}_{0}=0.13$, on the other hand, there are three nearly equal steps, each of which with a smaller decay length than the two large ones for parameter set (a). Therefore, less blocking is possible and more determinants with a negative sign appear in the Markov chain.

Zoom In Zoom Out Reset image size

Table 1.  Convergence of the energy for N = 4, $\lambda =1.3$ and $\beta =5.0$ for selected parameter combinations shown in figure 5.

Simulation	E	V	K	S
$P=2$ a	$12.1924(3)$	$9.0283(3)$	$3.1641(3)$	$0.4907(3)$
$P=2$ b	$12.3186(2)$	$9.0927(2)$	$3.2258(2)$	$0.3771(2)$
$P=14$ a	$12.293(4)$	$9.083(1)$	$3.210(4)$	$0.02664(1)$
$P=14$ b	$12.292(2)$	$9.0831(6)$	$3.209(2)$	$0.020600(7)$
CPIMC	$12.293(3)$	—	—	—

^a ${t}_{0}=0.04$, ${a}_{1}=0.0$ ^b ${t}_{0}=0.13$, ${a}_{1}=0.33$

The different convergence behaviors of the two free parameter combinations for small P leads to the question how to choose t₀ and a₁ for optimal results. To provide an answer, we consider the same system as in figure 5, and investigate the accuracy of the total energy as a function of t₀, for a fixed ${a}_{1}=0.33$. The simulation results are shown in the left panel of figure 6 for P = 2 (red squares), P = 3 (blue circles) and P = 4 (green diamonds). All three curves exhibit a similar decay towards the exact value starting from small t₀, followed by a minimum around ${t}_{0}=0.14$ and finally an increasing error for larger values. We note that as few as two propagators allow for an accuracy of $| \Delta E| /E\lt 2\times {10}^{-3}$ for the best choice of the free parameters. Figure 6(B) shows the dependency of the average sign S on t₀. Again, we observe that S decreases with increasing t₀ as explained during the discussion of figure 5. In addition, it is revealed that the combination of P = 4 and ${t}_{0}=0.01$ leads to a larger sign than P = 3 and ${t}_{0}\gt 0.10$. However, the optimum free parameters allow for a higher accuracy even for P = 2, compared to small t₀ with more propagators. Therefore, it turns out to be advantagous to use the fourth order factorization with the two free parameters despite the smaller average sign for the same P compared to the factorization with only a single daughter slice for each propagator, i.e., ${t}_{0}=0.0$.

Zoom In Zoom Out Reset image size

Finally, we mention that the optimal choice of a₁ and t₀ depends on the observable of interest. In figure 7, we investigate the effects of the free parameters on the convergence of the radial density distribution $n(r)$ for the same system as in figures 5 and 6. The left panel shows n as a function of the distance to the center of the trap, r, for four different P and the parameter combination ${a}_{1}=0.33$ and ${t}_{0}=0.13$, which has been proven to allow for nearly optimum energy values at P = 2, cf figure 6. The black curve corresponds to P = 10 and is converged within statistical uncertainty. For P = 2 (red diamonds), there appear significant deviations to the latter, in particular n is too large around the maximum $r\approx 1.25$ and too small at the boundary of the system. The P = 3 results (blue squares) exhibit the same trends although the differences towards the black curve are reduced. Finally, the density for P = 4 (green circles) can hardly be distinguished from the converged data. The right panel compares the density for P = 2 with two different combinations of free parameters. The red diamonds (parameter set (a)) correspond to the curve from the left panel and the green circles (parameter set (b)) to ${a}_{1}=0.0$ and ${t}_{0}=0.04$. The latter parameters clearly allow for a density distribution which is much closer to the exact results than the a₁ and t₀ values which provide the optimal energy.

Zoom In Zoom Out Reset image size

In the last section, we have demonstrated that the optimal choice of the free parameters a₁ and t₀ allows for the calculation of energies with an accuracy of $0.1\%$ with as few as two propagators, even at a relatively low temperature, $\beta =5.0$. However, with decreasing T (i.e., increasing β) the number of required propagators must be increased to keep the commutator error fixed. In figure 8, we investigate the effect of a decreasing temperature on the accuracy provided by a few propagators P for N = 4 electrons at indermediate coupling, $\lambda =1.3$. The left panel shows the total energy E as a function of the inverse temperature β. We compare results for P = 2 (green circles), P = 3 (red diamonds) and P = 4 (blue triangles) to exact results from CPIMC (black stars). At larger temperature, $\beta \leqslant 7.0$, all four datasets nearly coincide and exhibit the expected decrease towards the energy of the ground state. With increasing β, the P = 2 results exhibit an unphysical drop because two propagators are not sufficient and the commutator errors become more significant. The red and blue curves exhibit a qualitatively similar trend, however, the energy drop is weaker and shifted to lower temperature. Even at $\beta =10.0$, which is already very close to the ground state, three propagators allow for an accurate description of the system.

Zoom In Zoom Out Reset image size

In the right panel of figure 8, the average sign S is plotted versus the inverse temperature. At small β, the wavefunctions of the electrons do not overlap and, hence, the system is not degenerate. With decreasing temperature, exchange effects become increasingly important which leads to a decrease of S. However, while for standard PIMC the sign is expected to exponentially decrease with β, S seems to converge for PB-PIMC with P = 3 and P = 4 and exhibits an even slightly non-monotonous behavior for P = 2. The application of antisymmetric propagators leads to a competition with respect to S and β. On the one hand, with increasing inverse temperature off-diagonal matrix elements are increased, which leads to more negative determinants and, therefore, more negative weights in the Markov chain. On the other hand, the thermal wavelengths ${\lambda }_{{t}_{1}\epsilon }$ and ${\lambda }_{2{t}_{0}\epsilon }$ are increasing with β, which makes the blocking of large diagonal and off-diagonal elements more effective. Hence, the sign can even become larger with β once the system has reached the ground state, because the particle distribution remains constant while more elements in the diffusion matrix compensate each other in the determinants.

We conclude that few propagators allow for the calculation of accurate results up to low temperature, $\beta \leqslant 10.0$. For higher β, the system is in its ground state and finite temperature PIMC is no longer the method of choice.

In the previous sections, we have restricted ourselves to the investigation of small systems to illustrate the convergence and sign behavior depending on relevant parameters. In this section, we demonstrate that PB-PIMC allows for the calculation of accurate results at parameters where no other ab initio results have been reported, so far. Figure 9 shows results for N = 8 and N = 20 electrons at $\beta =3.0$ over a wide range of coupling parameters, λ. In panel (A), the average sign S is plotted versus λ for standard PIMC (squares), CPIMC (circles) and the present PB-PIMC (diamonds) with P = 2 and the parameter sets ${t}_{0}=0.14$ and ${a}_{1}=0.33$ (N = 8, blue symbols) and ${t}_{0}=0.10$ and ${a}_{1}=0.33$ (N = 20, red symbols), which are known to allow for accurate energies, cf figure 6. It is well understood that PIMC allows for the simulation of strongly coupled fermions, where exchange effects do not play a dominant role. With decreasing λ, the sign exhibits a sharp drop and the sign problem prevents the simulation within feasible computation time for $\lambda \leqslant 2.0$ and $\lambda \leqslant 5.0$, respectively. Evidently, larger systems lead to a more severe decrease of S at larger coupling strength. CPIMC, on the other hand, can be interpreted as a Monte Carlo simulation on a perturbation expansion around the ideal quantum system, i.e., $\lambda =0.0$. Hence, the method efficiently provides exact results for small coupling, where the system is close to an ideal one. For N = 20 around $\lambda \approx 0.3$, the sign almost instantly drops from $S\approx 0.97$ towards zero, and CPIMC is no longer applicable, without further approximation. This means that, in particular for larger systems, there have only been results for systems that are (a) almost ideal or (b) so strongly coupled that fermions and bosons lead to nearly equal physical properties. The physically particularly interesting regime where Coulomb correlations and Fermi statistics are significant simultaneously, has remained out of reach.

Zoom In Zoom Out Reset image size

However, the average sign from PB-PIMC exhibits a much less severe drop with decreasing λ than standard PIMC and saturates for $\lambda \leqslant 0.7$. For N = 8, the average sign remains above $S=0.08$, which allows for good accuracy with relatively low effort. The small sign, $S\sim {10}^{-3}$, for N = 20 indicates that the simulations are computationally involved but, in contrast to PIMC and CPIMC, still feasible. In panel (B) of figure 9, the total energy E for N = 20 is plotted versus λ over the entire coupling range and the statistical uncertainty from the PB-PIMC results is smaller than the size of the data points. Both, at small and large λ, the P = 2 results are in excellent agreement with the exact energy known from the other methods and, in addition, results are obtained for the particularly interesting transition region (region [II] in figure 1). In panel (C), we show the radial density for N = 20 and low coupling, $\lambda =0.10$, calculated with the parameter set ${t}_{0}=0.04$ and ${a}_{1}=0.0$, which has been proven effective for accurate densities $n(r)$. The PB-PIMC results (red diamonds) are in excellent agreement with the exact CPIMC data (blue squares) over the entire r-range. For completeness, we mention that this combination of parameters allows for an approximately three times as high sign as the choice from panels (A) and (B), which was choosen to allow for a good energy, and the results have been obtained within ${t}_{{\rm{CPU}}}\sim {10}^{3}$ core hours. Panel (D) shows the density of a strongly coupled system, $\lambda =15.0$, and N = 20. Again, the two propagators already provide very good agreement with the exact curve. In figure 1(B), we have shown density profiles for coupling parameters over the entire coupling range. At $\lambda =15$ (red pluses), there are three distinct shells and the physical behavior is dominated by the strong Coulomb repulsion. Decreasing the coupling to $\lambda =5$ (green bars) leads to a reduced extension of the system, and the three shells exhibit a much larger overlap. At indermediate coupling, $\lambda =2$ (blue crosses), both the interaction and fermionic exchange govern the system. The density profile is still significantly more extended than the ideal pendant, but n exhibits modulations instead of a flat curve. Decreasing the repulsion further to $\lambda =0.7$ (pink circles) leads to a further reduction of the extension. However, n does not approach a Gaussian-like profile as for ideal boltzmannons or bosons, but continues to exhibit the density modulations which are characteristic for fermions. For $\lambda =0.1$, the system is almost ideal and the density is completely dominated by the quantum statistics.

Finally, in figure 10 we compare density profiles for N = 20 particles at $\beta =3.0$ with Fermi-, Bose- and Boltzmann statistics. Panel (A) shows results for intermediate coupling, $\lambda =2.0$. The distinguishable boltzmannons (blue diamonds) exhibit a nearly flat profile without any shell structure, i.e., a liquid-like behavior. The bosonic particles (green circles) lead to an even smoother curve, with a slightly reduced extension of the system. For fermions (red squares), on the other hand, the exchange already plays a significant role, as the particles exhibit an additional repulsion due to the Pauli principle, and n decays only at larger r. In addition, the fermionic density profile exhibits distinct modulations. In panel (B), we show a comparison for smaller coupling, $\lambda =0.7$. Again, the boltzmannons and bosons lead to smooth density profiles which are very similar, despite a reduced extension of the Bose-system and an increased density around the center of the trap. The fermions exhibit a different behavior as the system is significantly more extended and the density profile again features distinct modulations.

Zoom In Zoom Out Reset image size

In conclusion, we have presented ab initio results for the energy and the density for up to 20 electrons over the entire coupling range. A comparison with standard PIMC and CPIMC has revealed excellent agreement in both the limits of weak and strong coupling. A more detailed investigation of the transition from the classical to the degenerate regime, including systematic comparisons with bosons and boltzmannons, is beyond the scope of this work and will be published elsewhere.

In the last section, we have shown that the sign problem is more severe for larger systems, cf figure 9(A). Here, we provide a more detailed investigation of the performance of our method in dependence on the particle number. In figure 11, the average sign S is plotted versus N for $\lambda =0.1$ and $\beta =3.0$, i.e., a very degenerate system, with two different combinations of free parameters. It is revealed that S exhibits an exponential decay with the system size and, as usual, the smaller t₀ leads to a more effective blocking. Therefore, the PB-PIMC approach still suffers from the fermion sign problem, and feasible system sizes for 2D quantum dots at weak coupling are limited to $N\leqslant 30$. This is a remarkable result since standard PIMC simulations for $\lambda =0.1$ and $\beta =3.0$ are possible only for $N\leqslant 4$.

Zoom In Zoom Out Reset image size