Ab initio path integral Monte Carlo simulations of hydrogen snapshots at warm dense matter conditions

We combine ab initio path integral Monte Carlo (PIMC) simulations with fixed ion configurations from density functional theory molecular dynamics (DFT-MD) simulations to solve the electronic problem for hydrogen under warm dense matter conditions [Böhme et al., Phys. Rev. Lett. 129, 066402 (2022)]. The problem of path collapse due to the Coulomb attraction is avoided by utilizing the pair approximation, which is compared against the simpler Kelbg pair potential. We find very favorable convergence behavior towards the former. Since we do not impose any nodal restrictions, our PIMC simulations are afflicted with the notorious fermion sign problem, which we analyze in detail. While computationally demanding, our results constitute an exact benchmark for other methods and approximations within DFT. Our setup gives us the unique capability to study important properties of warm dense hydrogen such as the electronic static density response and exchange-correlation kernel without any model assumptions, which will be very valuable for a variety of applications such as the interpretation of experiments and the development of new XC functionals.

Figure 1

Illustration of imaginary-time paths in the $x\text{−}\tau$ plane for three particles. The red dots denote the so-called beads, which depict the position of each particle at a particular imaginary-time step. The second and third paths from the left illustrate a so-called exchange cycle, where the trajectories of two particles are involved in a single path. Such permutation cycles [35] are required to correctly take into account the indistinguishable nature of both bosons and fermions in the PIMC formalism. This exchange cycle ensures that PIMC abides by the Pauli exclusion principle and identical bosons remain indistinguishable from one another [54].

Figure 2

Comparison of various effective potentials between a test electron and an ion located at $\mathbf{I}={(L/2,0,0)}^T$. Top panel: diagonal potential, $\mathbf{r}={\mathbf{r}}^{\prime }={(x,0,0)}^T$. Bottom panel: off-diagonal potential, $\mathbf{r}={(x,0,0)}^T, {\mathbf{r}}^{\prime }{(x,{\lambda }_{\varepsilon },0)}^T$.

Figure 3

Convergence of (a) the total energy, (b) potential energy, and (c) kinetic energy, with the inverse number of imaginary time slices $P^{-1}$. The calculations have been carried out for $N=4, \mathrm{\Theta }=1$, and $r_s=2$. The convergence behavior of the pair approximation is clearly favorable compared to the simple Kelbg pair potential.

Figure 4

Induced density convergence as a function of $P$ for a perturbation strength of $A=0.1, N=14, r_s=4$, and $\mathrm{\Theta }=1$. The red curve shows the results using the Kelbg pair potential, and the green curve depicts the PA. Since the lower density of $r_s=4$ causes a stronger localization of the electrons near the protons, the off-diagonal elements in the density matrix of the PA significantly improve the convergence with the number of imaginary time steps.

Figure 5

Density convergence with P for $r_s=4$ and $\mathrm{\Theta }=1$ along the $z$ direction (i.e., integrated over $x$ and $y$). The upper panel depicts PIMC results for the density using the PA, while the lower panel shows the Kelbg density for different numbers of imaginary time steps $P$. For both panels the insets have identical limits to highlight the improved convergence behavior of the PA.

Figure 6

Change in the integrated density in the $z$ direction due to an external perturbation of $A=0.1\mathrm{Ha}$ for different $P$ at $N=14, \mathrm{\Theta }=1$, and $r_s=4$. The deviation is defined as $\mathrm{\Delta }n=n_{pert.}\left(P\right)-n_{unpert.}\left(P\right)$ in units of the unperturbed UEG density $n_0$. In the upper panel the case of the PA is shown, while in the lower panel we depict the Kelbg result for different propagator numbers $P$. The vertical gray dashed lines show the ionic positions.

Figure 7

The average sign $\langle S\rangle$ compared both with the Kelbg potential (red crosses) as well as the PA (green circles) at metallic densities $r_s=2$ and $\mathrm{\Theta }=1$. Both curves exhibit an exponential decay of the average sign with increasing particle number. The Kelbg result has been calculated using $P=600$ propagators, where the PA required only $P=200$ propagators. Furthermore, the same conditions have been studied for $r_s=6$. The magenta crosses shows the Kelbg result at $r_s=6$ again with 600 propagators. For the PA we ran the same calculations but for $P=200$ (purple crosses) and $P=600$ (yellow crosses).

Figure 8

Average sign $\langle S\rangle$ as a function of the degeneracy parameter $\mathrm{\Theta }=1$ for $N=14$ particles at metallic density $r_s=2$. The Kelbg (red) and the PA (green) results show the expected increase of the average sign. The PA used $P=200$, and the Kelbg calculations were run with $P=600$ propagators.

Figure 9

Configurations from pair approximation PIMC simulations of a hydrogen snapshot with $N=20, \theta =1$, and $P=200$ for $r_s=2$ (top) and $r_s=6$ (bottom). The green spheres denote the ion positions, and the blue paths depict the imaginary-time paths of the electrons. At $r_s=6$, the localization of the electron paths around the ions is substantially increased compared to $r_s=2$. This heuristically indicates the emergence of bound states around the ions at lower density. The top panel has been taken from the Supplemental Material of Ref. [33].

Figure 10

Comparing PIMC results for the density along the $z$ direction for the same fixed ion snapshot with $N=14$ at $\theta =1$ for different values of the density parameter $r_s$. The dashed vertical lines show the $z$ coordinates of the ions.

Figure 11

Comparing PIMC results for the density along the $z$ direction for an ion snapshot of $N=14$ at $r_s=2$ for different values of the reduced temperature $\theta$. The dashed vertical lines show the $z$ coordinates of the ions.

Figure 12

Top (bottom): comparison of the density along the $z$ direction of an unperturbed snapshot with $N=14$ unpolarized electrons at $r_s=2$ ($r_s=4$) and $\theta =1$ between thermal DFT within LDA (green crosses) and PIMC with the PA and $P=600$ (red circles).