Permutation Blocking Path Integral Monte Carlo approach to the Static Density Response of the Warm Dense Electron Gas

TL;DR

This paper employs an advanced permutation blocking path integral Monte Carlo method to accurately compute the static density response of the warm dense electron gas, highlighting the importance of exact correlation treatments.

Contribution

The study introduces the application of PB-PIMC to the warm dense electron gas, providing highly accurate ab initio results for the static density response and local field correction.

Findings

PB-PIMC yields highly accurate static density response functions.

Comparison shows dielectric approximations are insufficient without exact correlations.

Multiple perturbations analysis informs static response calculations.

Abstract

The static density response of the uniform electron gas is of fundamental importance for numerous applications. Here, we employ the recently developed \textit{ab initio} permutation blocking path integral Monte Carlo (PB-PIMC) technique [T.~Dornheim \textit{et al.}, \textit{New J.~Phys.}~\textbf{17}, 073017 (2015)] to carry out extensive simulations of the harmonically perturbed electron gas at warm dense matter conditions. In particular, we investigate in detail the validity of linear response theory and demonstrate that PB-PIMC allows to obtain highly accurate results for the static density response function and, thus, the static local field correction. A comparison with dielectric approximations to our new \textit{ab initio} data reveals the need for an exact treatment of correlations. Finally, we consider a superposition of multiple perturbations and discuss the implications for the calculation of the static response function.

Figures

Equations

$${\chi }_{\text{RPA}}(\mathbf{q},\omega )=\frac{{\chi }_0(\mathbf{q},\omega )}{1-\frac{4\pi }{q^2}{\chi }_0(\mathbf{q},\omega )},$$

$${\chi }_{\text{LFC}}(\mathbf{q},\omega )=\frac{{\chi }_0(\mathbf{q},\omega )}{1-\frac{4\pi }{q^2}[1-G(\mathbf{q},\omega )]{\chi }_0(\mathbf{q},\omega )}.$$

$$K_{\text{xc}}(\mathbf{q},\omega )=-\frac{4\pi }{q^2}G(\mathbf{q},\omega )$$

$$\xi =\frac{N^{\uparrow }-N^{\downarrow }}{N}=0.$$

$$\hat{H}=-\frac{1}{2}\sum _{i=1}^N{\nabla }_i^2+\frac{1}{2}\sum _{i=1}^N\sum _{j\ne i}^N{\Psi }_{\text{E}}({\mathbf{r}}_i,{\mathbf{r}}_j)+\frac{N}{2}{\xi }_{\text{M}},$$

$$Z=\text{Tr}\hat{\rho },$$

$$\hat{\rho }=e^{-\beta \hat{H}}.$$

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

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

$$Z$$

$$e^{-ϵ\hat{H}}\approx e^{-ϵ\hat{V}}{e}^{-ϵ\hat{K}},$$

$$Z$$

$$e^{-ϵ\hat{H}}$$

$$[[\hat{V},\hat{K}],\hat{V}]$$

$${\mathbf{F}}_i$$

$$Z$$

$$D_{\alpha ,\uparrow }$$

$$D_{\alpha ,\downarrow }$$

$${\rho }_{\alpha ,\uparrow }(i,j)={\lambda }_{t_1ϵ}^{-3}\sum _{\mathbf{n}}{e}^{-\frac{\pi }{{\lambda }_{t_1ϵ}^2}({\mathbf{r}}_{\alpha ,\uparrow ,j}-{\mathbf{r}}_{\alpha A,\uparrow ,i}+\mathbf{n}L{)}^2},$$

$$\hat{H}={\hat{H}}_0+{\hat{H}}_{\text{ext}}(t).$$

$${\hat{H}}_{\text{ext}}(t)=2A\sum _{i=1}^N\text{cos}\left({\mathbf{r}}_i\cdot \mathbf{q}-\Omega t\right),$$

$${\varphi }_{\text{ext}}(\mathbf{r},t)=2A\text{cos}\left(\mathbf{r}\cdot \mathbf{q}-\Omega t\right).$$

$$\tilde{\chi }(\mathbf{q},\tau )=\frac{-i}{\hslash }{⟨\left[\rho (\mathbf{q},\tau ),\rho (-\mathbf{q},0)\right]⟩}_0\Theta (\tau ),$$

$$\chi (\omega ,\mathbf{q})=\underset{\eta \to 0}{\lim }{\int }_{-\infty }^{\infty }\text{d}\tau e^{(i\omega -\eta )\tau }\tilde{\chi }(\mathbf{q},\tau ).$$

$$\underset{\omega \to 0}{\lim }\chi (\omega ,\mathbf{q})=\chi (\mathbf{q}),$$

$${\varphi }_{\text{ext}}(\mathbf{r})=2A\text{cos}(\mathbf{r}\cdot \mathbf{q}),$$

$${\rho }_{\text{ind}}(\mathbf{q})={\rho }_{\text{ext}}(\mathbf{q})\frac{4\pi }{q^2}\chi (\mathbf{q}).$$

$${\rho }_{\text{ext}}(\mathbf{r})$$

$$\Rightarrow {\rho }_{\text{ext}}(\mathbf{q})$$

$${\rho }_{\text{ind}}(\mathbf{q})$$

$$\chi (\mathbf{q})=\frac{1}{A}{⟨{\hat{\rho }}_{\mathbf{q}}⟩}_A.$$

$${⟨n(\mathbf{r})⟩}_A=n_0+2A\text{cos}\left(\mathbf{q}\cdot \mathbf{r}\right)\chi (\mathbf{q}).$$

Full text

††preprint: APS/123-QED

I Introduction

The uniform electron gas (UEG), which is comprised of Coulomb interacting electrons in a homogeneous neutralizing background, is one of the most seminal model system in quantum many-body physics and chemistry loos . In addition to the UEG’s importance for, e.g., the formulation of Fermi liquid theory quantum_theory ; quantum_theory2 and the quasi-particle picture of collective excitations pines ; pines2 , accurate parametrizations of its ground state properties vwn ; perdew ; pw ; gori ; gori2 based on ab initio quantum Monte Carlo calculations gs1 ; gs2 ; ortiz ; ortiz2 ; spink have been pivotal for the arguably unrivaled success of density functional theory simulations of real materials ks ; dft_review ; dft_burke .

The density response of the UEG to a small external perturbation as described by the density response function is of high importance for many applications quantum_theory . The well-known random phase approximation (RPA) rpa_original provides a qualitative description for weak coupling strength (high density),

[TABLE]

where $\chi_{0}(\mathbf{q},\omega)$ denotes the density response function of the ideal (i.e., non-interacting) system. However, since Eq. (1) does not incorporate correlations beyond the mean field level, RPA breaks down even for moderate coupling. This shortcoming is usually corrected in the form of a local field correction (LFC) $G(\mathbf{q},\omega)$ kugler1 , modifying Eq. (1) to

[TABLE]

Hence, by definition, the exact LFC contains all exchange-correlation effects beyond RPA. Common approximations for $G$ include the approaches by Singwi-Tosi-Land-Sjölander (STLS) stls_original and Vashishta and Singwi (VS) vs_original . It is important to note that the accurate determination of $G(\mathbf{q},\omega)$ is an important end in itself as it can be straightforwardly utilized as input for other calculations. For example, it is directly related to the XC kernel

[TABLE]

of density functional theory in the adiabatic-connection fluctuation-dissipation formulation lu ; patrick ; burke2 . This allows for the construction of a true non-local XC functional, which is a promising approach to go beyond the ubiquitous gradient approximations dft_burke ; pbe and thereby increase the predictive capabilities of DFT. Further applications of the LFCs for current warm dense matter (WDM, see below) research include the calculation of the dynamic structure factor fortmann1 ; fortmann2 ; collision1 ; collision2 as it can be obtained with X-ray Thomson scattering from a variety of systems, energy transfer rates energy_transfer1 ; energy_transfer2 , the electrical and optical conductivity cond1 ; cond2 , and equation of state models of ionized plasmas equation_of1 ; equation_of2 ; equation_of3 . Finally, we mention the construction of effective potentials both for WDM saum1 ; saum2 and beyond pseudo_potential ; gravel .

In the ground state, Moroni et al. moroni2 obtained accurate QMC results for the static response function [i.e., $\omega\to 0$, see Eq. (23)] - and thereby the static LFC - by simulating an electron gas with a weak external harmonic perturbation moroni1 ; bowen1 ; bowen2 ; senatore . This has allowed for a systematic assessment of the accuracy of previous approximations. Further, the ab initio data for the LFC have subsequently been parametrized by Corradini et al. cdop , and the zero temperature limit of the static density response is well understood.

However, recently there has emerged a growing interest in matter under extreme conditions, i.e., at high density and temperature, which occurs in astrophysical objects such as brown dwarfs and planet interiors knudson ; militzer . Furthermore, similar conditions are now routinely realized in experiments with laser excited solids ernst or inertial confinement fusion targets nora ; schmit ; hurricane3 ; kritcher . This ’warm dense matter’ (WDM) regime is characterized by two parameters being of the order of unity wdm_book : (i) the Wigner-Seitz radius $r_{s}=\overline{r}/a_{\text{B}}$ and (ii) the reduced temperature $\theta=k_{\text{B}}T/E_{\text{F}}$, where $\overline{r}$, $a_{\text{B}}$ and $E_{\text{F}}$ denote the mean inter particle distance, Bohr radius and Fermi energy fermi_E , respectively. Naturally, accurate data for the static LFC at such extreme conditions are highly desirable. In fact, in lieu of thermodynamic data often ground state results are used at WDM conditions, which might not be appropriate wdm_book .

Yet, a theoretical description of warm dense electrons is notoriously hard since it must account for the nontrivial interplay of (a) the strong quantum Coulomb collisions, (b) excitation effects due to the high temperature, and (c) quantum degeneracy effects (e.g. fermionic exchange). In particular, conditions (a) and (b) rule out perturbation expansions and ground state methods, respectively, leaving thermodynamic quantum Monte Carlo methods as the most promising option. Unfortunately, QMC simulations of degenerate electrons suffer from the fermion sign problem (FSP) loh ; troyer so that the widespread path integral Monte Carlo (PIMC) approach cep is limited to small system sizes and high temperatures, preventing simulations under WDM conditions dornheim_pop . Despite its remarkable success in the ground state, at finite temperature, the fixed node approximation node ; bcdc (which avoids the FSP) can lead to systematic errors exceeding $10\%$ tim3 . This unsatisfactory situation has sparked remarkable progress in the field of fermionic QMC simulations. In particular, the joint usage of two novel complementary approaches (in combination with an improved finite-size correction dornheim_prl ) has recently allowed to obtain the first complete ab initio description of the warm dense electron gas dornheim_prl ; groth_prl : (i) At high density and weak to moderate coupling, the configuration PIMC (CPIMC) approach tim_cpp15 ; tim1 ; groth , which is formulated in Fock space and can be understood as a Monte Carlo calculation of the (exact) perturbation expansion around the ideal system, is capable to deliver exact results over a broad temperature range. (ii) The permutation blocking PIMC (PB-PIMC) approach dornheim ; dornheim2 ; dornheim3 extends standard PIMC towards higher density and lower temperature and allows for accurate results in large parts of the WDM regime. In this work, we use the latter method to carry out simulations of the harmonically perturbed electron gas under warm dense matter conditions.

A brief introduction of the UEG model (Sec. II.1) is followed by a comprehensive introduction to fermionic QMC simulations at finite temperature. In particular, we explain how the antisymmetry of the density operator leads to the fermion sign problem in standard PIMC (Sec. II.2.1), and how this is addressed by the idea of permutation blocking (Sec. II.2.2). Further, we give a concise overview of linear response theory and how the static density response can be obtained by simulating the harmonically perturbed system (Sec. II.3). In Sec. III, we show extensive PB-PIMC results to investigate the dependence on the perturbation strength (III.1), the convergence with the number of imaginary time propagators (III.2), and the wave vector dependence (III.3), which also allows to address possible finite-size effects. Finally, in Sec. III.5 we consider the response to a superposition of multiple perturbations with different wave vectors and the resulting implications for the calculation of $\chi$.

II Theory

II.1 Uniform Electron Gas

The uniform electron gas is a model system of $N$ electrons in a positive homogeneous background that ensures charge neutrality. Throughout this work, we assume an unpolarized (paramagnetic) system, i.e., $N^{\uparrow}=N^{\downarrow}=N/2$ [with $\uparrow$ ($\downarrow$) denoting the number of spin-up (-down) electrons] and, thus,

[TABLE]

To alleviate the differences between a finite model system and the thermodynamic limit (finite-size effects), we employ Ewald summation for the repulsive pair interaction. Therefore, the Hamiltonian (in Hartree atomic units) is given by

[TABLE]

where $\Psi_{\text{E}}(\mathbf{r},\mathbf{s})$ and $\xi_{\text{M}}$ denote the Ewald pair potential and the well-known Madelung constant, see, e.g., Ref. fraser .

II.2 Quantum Monte Carlo

II.2.1 Path Integral Monte Carlo

Throughout the entire work, we consider the canonical ensemble where the volume $V=L^{3}$ (with $L$ being the box length), particle number $N$ and inverse temperature $\beta=1/k_{\text{B}}T$ are fixed. To derive the path integral Monte Carlo formalism cep , we consider the partition function

[TABLE]

which is defined as the trace over the canonical density operator $\hat{\rho}$

[TABLE]

Let us temporarily restrict ourselves to distinguishable particles and re-write Eq. (6) in coordinate representation:

[TABLE]

where $\mathbf{R}=\{\mathbf{r}_{1},\dots,\mathbf{r}_{N}\}$ contains the all $3N$ particle coordinates. Since the matrix elements of $\hat{\rho}$ are not readily known, we use the group property

[TABLE]

with $\epsilon=\beta/P$. Furthermore, we insert $P-1$ unity operators of the form $\hat{1}=\int\text{d}\mathbf{R}_{\alpha}\ket{\mathbf{R}_{\alpha}}\bra{\mathbf{R}_{\alpha}}$ into Eq. (8) and obtain

[TABLE]

and the integration is carried out over $P$ sets of particle coordinates, $\text{d}\mathbf{X}=\text{d}\mathbf{R}_{0}\dots\text{d}\mathbf{R}_{P-1}$. We stress that Eq. (II.2.1) is still exact. The main benefit of this re-casting is that the new expression involves $P$ density matrix elements, but at a $P$ times higher temperature. Each of these high temperature factors can now be substituted using some suitable high-$T$ approximation, e.g., the simple primitive factorization

[TABLE]

with $\hat{V}$ and $\hat{K}$ being the operators for the potential and kinetic contribution to the Hamiltonian, respectively, and which becomes exact in the limit $P\to\infty$ trotter . The resulting high-dimensional integral is then evaluated using the Metropolis algorithm metropolis (we employ a simulation scheme based on the worm algorithm bon ; bon2 ).

However, to simulate fermions we must extend the partition function from Eq. (8) by the sum over all particle permutations, which, for an unpolarized system, gives

[TABLE]

with $\sigma^{\uparrow,\downarrow}$ denoting particular elements from the permutation groups $S_{N}^{\uparrow,\downarrow}$, and $\hat{\pi}_{\sigma^{\uparrow,\downarrow}}$ being the corresponding permutation operators. In practice, this leads to the occurrence of so-called exchange cycles within the PIMC simulations, which are paths incorporating more than a single particle, see Fig. 1. The problem is that the sign of each configuration depends on the parity of the permutations involved which can be both positive and negative. At low temperature and high density, permutation cycles with both positive and negative signs appear with a similar frequency and, thus, the signal to noise ratio vanishes. This is the notorious fermion sign problem loh ; troyer , which limits standard PIMC to weak degeneracy where fermionic exchange plays only a minor role and, therefore, precludes its application to warm dense matter dornheim_pop . This is illustrated in Fig. 1, where we show random configurations from standard PIMC simulations of the UEG with $N=19$ spin-polarized electrons at a density parameter $r_{s}=1$ and three different temperatures. Each particle is represented by $P=32$ so-called ’beads’, which are connected by the (red) kinetic density matrix elements and thus form the eponymous paths. At high temperature, $\theta=8$ (top panel), each particle is represented by a distinct, separate path and exchange cycles occur only infrequently. Therefore, the FSP is not severe and PIMC simulations are feasible. At moderate, WDM temperatures ($\theta=1$, center panel), fermionic exchange is influencing the system significantly, and multiple exchange cycles are visible in the screenshot. Since each pair exchange causes a sign change in the Monte Carlo simulation, a standard PIMC simulation is no longer feasible. Finally, at low temperature ($\theta=0.3$, bottom panel) nearly all particles are involved in exchange cycles, and the system is dominated by the antisymmetric nature of the electrons (i.e., Pauli blocking).

II.2.2 Permutation blocking

The fermion sign problem is $NP$-hard troyer and a general solution is, at the time of this writing, not in sight. Therefore, there does not exist a single QMC method that is applicable for all parameters. Nonetheless, it is possible to go beyond standard PIMC by employing the recently introduced permutation blocking PIMC approach dornheim ; dornheim2 . The first key ingredient is the usage of antisymmetric imaginary time propagators, i.e., determinants, which allows for a combination of positive and negative terms into a single configuration weight det1 ; det2 ; det3 . However, while this ’permutation blocking’ can indeed lead to a significant reduction of the fermion sign problem, with an increasing number of propagators $P$ this advantage quickly vanishes. For this reason, as the second key ingredient, we utilize a higher order factorization of the density matrix ho2 ; ho3

[TABLE]

which allows for sufficient accuracy even for a small number of imaginary time slices, for the definitions of the coefficients $v_{1},t_{1},v_{2},a_{1}$ and $t_{0}$, see Refs. dornheim ; dornheim2 . The $\hat{W}$ operators correspond to modified potential terms combining the standard potential contribution $\hat{V}$ with double commutator terms of the form ho3

[TABLE]

where $\mathbf{F}_{i}$ denotes the total force on a particle ’$i$’. Finally, this allows one to obtain the PB-PIMC partition function dornheim3

[TABLE]

with $\tilde{V}_{\alpha}$ and $\tilde{F}_{\alpha}$ containing all contributions of the potential energy and the forces, respectively, and the exchange-diffusion functions

[TABLE]

Here $\rho_{\alpha,\uparrow}$ denotes the diffusion matrix of a single time slice

[TABLE]

with $\lambda_{t_{1}\epsilon}=\sqrt{2\pi\epsilon t_{1}\hbar^{2}/m}$ being the corresponding thermal wavelength. Observe that Eq. (13) implies that there are three imaginary time slices for each propagator $\alpha=0,\dots,P-1$, with $\mathbf{R}_{\alpha}$, $\mathbf{R}_{\alpha A}$ and $\mathbf{R}_{\alpha B}$ denoting the corresponding sets of particle coordinates.

In a nutshell, in the PB-PIMC approach, we do not have to explicitly sample each positive or negative permutation cycle. Instead, we combine configuration weights with different signs in the determinants, which results in an analytical cancellation of terms and, thus, a significantly alleviated sign problem. This is illustrated in Fig. 2, where we show a random configuration from a PB-PIMC simulation of the warm dense UEG with $N=9$ spin-polarized electrons, $r_{s}=1$ and $\theta=1$ for $P=2$. The green, blue and purple beads correspond to the three different kinds of imaginary time slices due to the higher order factorization of the density operator, cf. Eq. (13). In contrast to the standard PIMC configurations from Fig. 1, every bead can be involved in multiple connections here. In fact, each bead is connected to all $N$ beads on the next and previous slices although the weight of the connection exponentially decreases with spatial difference, which is expressed by the different line widths of the (red) connections. Evidently, many beads of the depicted screen shot exhibit multiple visible connections, which means that a significant amount of analytical cancellation is accomplished within the determinants and, unlike standard PIMC, simulations are still feasible dornheim_pop .

This permutation blocking is most effective when $\lambda_{t_{1}\epsilon}$ is comparable (or larger) than the mean inter-particle distance. However, for $P\to\infty$ the beneficial effect vanishes and the original sign problem from standard PIMC is recovered. This plainly illustrates the paramount importance of a sophisticated higher order factorization scheme such as Eq. (13).

II.3 Linear Response Theory

In linear response theory (LRT), we consider the effect of a small external perturbation on the density of the system of interest

[TABLE]

Note that, in general, $\hat{H}_{\text{ext}}(t)$ is time-dependent. Throughout this work, the unperturbed Hamiltonian $\hat{H}_{0}$ corresponds to the UEG as introduced in Eq. (5) and the perturbation is given by a sinusoidal external charge density of wave vector $\mathbf{q}$

[TABLE]

which corresponds to the potential

[TABLE]

The standard definition of the density response function is given by

[TABLE]

where the expectation value is with respect to the unperturbed system. Note that Eq. (21) only depends on the time difference $\tau=t-t^{\prime}$ and, due to the homogeneity of the unperturbed system, $\chi$ only depends on the modulus of the wave vector. The corresponding Fourier transform is given by

[TABLE]

Throughout this work, we restrict ourselves to the static limit that is defined as

[TABLE]

i.e., the response of the electron gas to a time-independent external perturbation

[TABLE]

and, henceforth, the $\omega$-dependence is simply dropped. More precisely, the physical interpretation of $\chi(\mathbf{q})$ is the description of the density response (i.e., the induced charge density $\rho_{\text{ind}}(\mathbf{q})$) due to the external charge density $\rho_{\text{ext}}(\mathbf{q})$

[TABLE]

The external density follows from the Poisson equation as

[TABLE]

and the induced density is the difference between the perturbed and unperturbed systems:

[TABLE]

where we made use of the fact that $\braket{\hat{\rho}_{\mathbf{q}}}_{0}=0$. Thus, it holds

[TABLE]

In order to obtain the desired static density response function, we carry out multiple QMC simulations for each wave vector $\mathbf{q}=2\pi L^{-1}(a,b,c)^{T}$ (with $a,b,c\in\mathbb{Z}$) for different values of $A$ and compute the expectation value from Eq. (II.3). For sufficiently small $A$, $\braket{\hat{\rho}_{\mathbf{q}}}_{A}$ is linear with respect to $A$ with $\chi(\mathbf{q})$ being the slope.

Another way to obtain the response function from the QMC simulation of the perturbed system is via the perturbed density profile in coordinate space:

[TABLE]

In practice, we compute the lhs. of Eq. (30) using QMC and perform a fit of the rhs. with $\chi(\mathbf{q})$ being the only free parameter. Naturally, in the linear response regime both ways to obtain $\chi(\mathbf{q})$ are equal.

For completeness, we mention that the dynamic response can be obtained in a similar fashion by considering explicitly time dependent perturbations, e.g., using non-equilibrium Green function techniques negf1 ; negf2 for quantum systems or molecular dynamics md1 ; md2 in the classical case.

A second strategy to compute the density response from thermodynamic QMC simulations in LRT is by considering imaginary-time correlation functions (ITCF) of the unperturbed system. In particular, the static response function can be obtained from the fluctuation dissipation theorem bowen1

[TABLE]

as an integral over the imaginary time $\tau$. If one is solely interested in the linear response of the system, invoking Eq. (31) constitutes the superior strategy since all $\mathbf{q}$-vectors can be computed from a single simulation. However, this requires a QMC estimation of the ITCF on a sufficient $\tau$-grid, which is straightforward in standard PIMC where $P>100$ is not an obstacle. For PB-PIMC, simulations are only possible for a small number of imaginary-time propagators (typically $P\lesssim 4$), see Sec. II.2.2, which precludes the evaluation of Eq. (31). Nevertheless, we stress that it is only the permutation blocking idea that allows to carry out simulations at warm dense matter conditions in the first place, since standard PIMC simulations are not feasible due to the FSP. In addition, the application of an external perturbation allows to go beyond LRT and to consider arbitrarily strong perturbation strengths.

III Results

III.1 Dependence on Perturbation strength

Let us start our investigation of the harmonically perturbed electron gas by considering the dependence on the perturbation amplitude $A$. In Fig.  3, we show PB-PIMC results for the density profile along the $x$-direction for $N=54$ unpolarized electrons at $r_{s}=10$ and $\theta=1$ for the perturbation wave vector $\mathbf{q}=2\pi L^{-1}(2,0,0)^{T}$. In the top panel, the depicted $A$-values are relatively small. The black lines correspond to fits according to Eq. (30). Evidently, for $A=0.001$ and $A=0.005$ those curves are in perfect agreement with the QMC results, which indicates that here the linear response theory is accurate. In contrast, for $A=0.01$ significant (although small, $\Delta A/A\sim 1\%$) deviations appear, which are most pronounced around the minima and maxima. In the center panel, we systematically increase $A$ up to a factor two. Clearly, with increasing perturbation amplitude the deviations between the exact QMC results and the cosine-fit predicted by LRT become more severe, as it is expected. Finally, in the bottom panel we show the density profiles for even larger perturbations. Eventually, the external potential becomes the dominating feature, resulting in a strongly inhomogeneous electron gas. For the largest depicted perturbation, $A=0.1$, there appear two distinct shells with a vanishing density in between.

To systematically investigate the effect of the perturbation amplitude on our QMC estimation of the static response function $\chi(\mathbf{q})$, we show results in Fig. 4 for the induced density $\rho_{\text{ind}}(\mathbf{q})$ for the same system and two different wave vectors, $\mathbf{q}=2\pi L^{-1}(q_{x},0,0)^{T}$ with $q_{x}=2$ (top panel) and $q_{x}=1$ (bottom panel). The black squares correspond to the direct QMC results, cf. Eq. (II.3), and the green crosses have been obtained by performing a cosine-fit to the density profiles according to Eq. (30). The red lines depict a linear fit to the black squares for $A<0.01$. First and foremost, we observe a perfect agreement between the direct QMC results and the cosine-fits for small $A$ as predicted by the linear response theory. Even for $A=0.01$, where the cosine-fit exhibits significant deviations to the density profile from QMC, we find perfect agreement between the black and green points and also to the fit. With increasing $A$, however, the assumptions of linear response theory are no longer valid. Interestingly, the $\rho$-values obtained from the cosine-fit exhibit significantly larger deviations to the linear response prediction (red line) than the direct QMC results. For example, at $A=0.05$ the deviation of the green points is twice as large as for the black squares.

In the bottom panel of Fig. 4, the same information is shown for a smaller wave vector, $q_{x}=1$. Firstly, we observe a significantly smaller density response [cf. Fig. 8]. This, in turn, means that linear response theory is accurate up to much larger $A$-values as the system only weakly reacts to such an external perturbation.

To further illustrate this point, in the top panel of Fig. 5 we show the corresponding average signs from the QMC simulations for both wave vectors investigated in Fig. 4. For small perturbations, $S$ is equal for both $\mathbf{q}$ and approaches the result for the unperturbed system. With increasing $A$, the system becomes more inhomogeneous, i.e., there appear regions of increased (and also decreased) density, see the bottom panel of Fig. 5 where we show the corresponding density profiles for strong perturbations, $A=0.1$. This, in turn, leads to increased fermionic exchange, resulting in a significantly decreased average sign in our PB-PIMC simulations. Since the density response is more pronounced for $q_{x}=2$, here $S$ exhibits a faster decrease in dependence of $A$. We conclude that PB-PIMC (and also standard PIMC) simulations of the inhomogeneous electron gas are significantly more computationally demanding than simulations of the UEG at equal conditions. Nevertheless, this is of no consequence for the determination of the static response function as this is only possible for $A$-values that are sufficiently small for the linear response theory to deliver an accurate description, i.e., systems that are close to the uniform case.

III.2 Convergence with propagators

As discussed in Sec. II.2.2, PB-PIMC crucially relies on the higher order factorization of the density operator, Eq. (13), to allow for sufficient accuracy with only few imaginary time propagators. In the following section, this situation is investigated in detail.

In the top panel of Fig. 6, we plot direct QMC results for the induced density for the unpolarized UEG with $r_{s}=10$, $\theta=1$, and $N=34$ electrons versus the inverse number of propagators $P^{-1}$. The perturbation is given by the wave vector $\mathbf{q}=2\pi L^{-1}(1,0,0)^{T}$ and amplitude $A=0.01$, which is well within the linear response regime. Evidently, only the result for $\rho$ with $P=2$ propagators significantly deviates from the rest and, for the $P=4$ propagators used above, the PB-PIMC results are converged within the statistical uncertainty. The bottom panel shows the corresponding density profiles along the $x$-direction. Here, even the results for only $P=2$ propagators exhibits no significant deviations to the other curves.

As a second example, in Fig. 7 we consider the same system as in Fig. 6, but with $N=54$ electrons and a larger wave vector for the perturbation, $\mathbf{q}=2\pi L^{-1}(5,0,0)^{T}$. In the top panel, we again show direct QMC results for $\rho$ in dependence of the inverse number of propagators. However, in contrast to the data depicted in Fig. 6, here we see significant differences for different $P$. The black line corresponds to a parabolic fit of the form

[TABLE]

which reproduces all QMC results within error bars. Nevertheless, we stress that the functional form in Eq. (32) has been empirically chosen and does merely serve as a guide to the eye since, for large $P$, the propagator error is expected to exhibit a fourth-order decay, see Ref. ho3 for a comprehensive discussion. Evidently, for $P=4$ there occurs a systematic bias of $\Delta\rho/\rho\approx 2\%$ at such a large wave vector. This is reflected in the increasing error bars towards large $\mathbf{q}$ in the wave vector dependence plot, i.e., Fig. 8, and can be understood as follows: The propagator error is a direct consequence of the non-commuting of the kinetic ($\hat{K}$) and potential ($\hat{V}$) contributions of the Hamiltonian. The larger the wave vector $\mathbf{q}$, the faster the spatial variations of the external potential and, because $\hat{K}\propto\nabla^{2}$, the larger the error terms, which involve nested commutators of $\hat{K}$ and $\hat{V}$.

The bottom panel of Fig. 7 shows the corresponding results for the total potential energy, i.e., the sum of the Ewald interaction and the external perturbation. Evidently, no deviations can be resolved within the given statistical uncertainty, even for $P=2$ propagators. This is similar to previous findings for the unperturbed UEG dornheim2 ; dornheim3 and reflects the circumstance that for $V$ the particle interaction dominates. In stark contrast, the induced density $\rho$ is particularly sensitive to the small external perturbation which, as explained above, requires a larger number of propagators to be sufficiently incorporated.

III.3 Wave vector dependence of $\chi(\mathbf{q})$ and finite size effects

Due to the momentum quantization in a finite simulation box, QMC calculations are only possible at an $N$-dependent discrete $\mathbf{q}$-grid. Therefore, the investigation of finite-size effects in the static response function requires us to obtain results over a broad wave vector range, as shown in Fig. 8. The grey and red curves correspond to the predictions due to the random phase approximation (RPA), cf. Eq. (1), and with a LFC from the (finite-$T$) STLS formalism stls ; stls2 , respectively. For small $\mathbf{q}$, both approximations exhibit the same exact parabolic behavior kugler2 . With increasing $\mathbf{q}$, however, there appear significant systematic deviations with a maximum of $\Delta\chi/\chi\sim 50\%$ around $q\approx 0.35$ (i.e., around twice the Fermi vector $k_{\text{F}}=(9\pi/4)^{1/3}/r_{s}$). The symbols correspond to our QMC results obtained according to Eq. (29) and the colors distinguish different particle numbers, in particular $N=54$ (blue crosses), $N=34$ (light blue circles), $N=20$ (yellow squares), $N=14$ (black triangles), and $N=8$ (green diamonds). First and foremost, we note that the main effect of different system size is the $\mathbf{q}$-grid, while the functional form itself is remarkably well converged even for as few as $N=8$ particles, cf. the bottom panel showing a magnified segment. This is similar to the analogous behavior of the static structure factor $S(\mathbf{q})$ of the warm dense UEG found in Refs. dornheim_prl ; dornheim_pop . Evidently, momentum shell effects as observed at $T=0$ in Refs. bowen2 ; moroni2 do not appear above $\theta=0.5$. Secondly, we find that the static local field correction due to the STLS closure relation leads to a significant improvement compared to RPA due to the improved treatment of correlations.

We thus conclude that our QMC approach allows, for the first time, to unambiguously assess the accuracy of the multitude of existing and widely used dielectric approximations and, in addition, to provide highly accurate data, which can subsequently be used as input for other theories. However, a comprehensive study over a broad parameter range is beyond the scope of this work and will be provided in a future publication.

III.4 Comparison of PB-PIMC to standard PIMC

As an additional benchmark for the static response obtained with PB-PIMC, in Fig. 9 we show $\chi(\mathbf{q})$ for the unpolarized UEG with $N=8$, $r_{s}=10$, and $\theta=4$. Since for such a temperature fermionic exchange plays only a minor role, in addition to PB-PIMC (green crosses) also standard PIMC (black squares) calculations are feasible. Evidently, both independent data sets are in excellent agreement over the entire $\mathbf{q}$-range, as expected. In addition, we again show results from RPA (grey) and STLS (red) and find qualitatively similar behavior to Fig. 8. However, due to the four times higher temperature correlations play a less important role, which means that (i) RPA and STLS exhibit less deviations towards each other, and (ii) the density response from STLS is in much better agreement with the QMC data. For completeness, we note that a more meaningful assessment of the systemic error due to the STLS approximation requires to eliminate the possibility of finite-size effects in the QMC data (as done in Fig. 8 at lower temperature, $\theta=1$) and, thus, to consider larger particle numbers $N$.

III.5 Multiple $\mathbf{q}$-vectors from a single simulation

When we have to perform at least a single (or even a few for different $A$) QMC simulation for each $\mathbf{q}$-value, the investigation of the wave vector dependence as depicted in Fig. 8 is computationally quite involved. However, by definition in linear response theory the response of a system to multiple perturbations is described by a superposition of the responses to each perturbation. Therefore, it should be possible to obtain the response function for multiple $\mathbf{q}$-values from a single QMC simulation where we apply a superposition of $N_{A}$ harmonic perturbations,

[TABLE]

The induced density is then calculated for each wave vector $\mathbf{q}_{k}$ according to Eq. (II.3). Furthermore, the density profile in coordinate space is given by

[TABLE]

which means that we have to perform a fit where the free parameters are given by the $N_{A}$ values of $\chi(\mathbf{q}_{k})$.

In Fig. 10, we show QMC results for the density profile in $x$-direction for $N=54$, $r_{s}=10$, and $\theta=1$. The green squares have been obtained from a simulation with a superposition of $N_{A}=2$ perturbations with $\mathbf{q}_{1}=2\pi L^{-1}(1,0,0)^{T}$ and $\mathbf{q}_{2}=2\pi L^{-1}(2,0,0)^{T}$ and $A_{1}=A_{2}=0.005$, i.e., an amplitude that is expected to be well within the linear response regime. As a comparison, the yellow and red points correspond to the QMC results with a single perturbation with $q_{x}=1$ (yellow) and $q_{x}=2$ (red). Further, the black crosses have been obtained as a superposition of the latter and are in perfect agreement with the green squares. This is a strong indication that the linear response is still valid for multiple perturbations under the present conditions. In addition, we have fitted the RHS. of Eq. (34) to the green squares and in this way obtained $\chi(\mathbf{q}_{k})$ for both $\mathbf{q}_{k}$-values. This, in turn, allows us to reconstruct the density response of the system to a perturbation with only a single $\mathbf{q}_{k}$-value, i.e., the blue curves. Again, we find excellent agreement to the corresponding QMC simulations.

To further pursue this point, in Fig. 11 we show the induced density matrix for different amplitudes $A$. The green squares and red circles have been obtained from a simulation with two $\mathbf{q_{k}}$-vectors and correspond to the direct QMC estimate and the cosine-fit according to Eq. (34), respectively. The blue crosses have been obtained from the QMC simulation with only a single harmonic perturbation and the red line depicts a linear fit. Evidently, all points are in excellent agreement for all $A$-values both for $q_{x}=1$ (top panel) and $q_{x}=2$ (bottom panel). Therefore, we conclude that it is indeed possible to obtain multiple values of the static density response function $\chi(\mathbf{q})$ simultaneously.

Finally, to investigate the perturbation strength dependence for a QMC simulation with a superposition of multiple $\mathbf{q}$-vectors in more detail, we consider a combination of $N_{A}=3$ perturbations with $\mathbf{q}_{1}=2\pi L^{-1}(1,0,0)^{T}$, $\mathbf{q}_{2}=2\pi L^{-1}(0,2,0)^{T}$, and $\mathbf{q}_{3}=2\pi L^{-1}(0,0,3)^{T}$ and equal amplitude, $A_{1}=A_{2}=A_{3}$, over a broad $A$-range. The results are shown in Fig. 12 where direct QMC results for the induced density matrix are shown both from the simulation with the superposition (black squares) and, as a reference, from a simulation with only a single perturbation (green crosses). As usual, the red line corresponds to a linear fit within the linear response regime. For both $q_{x}=1$ (top panel) and $q_{x}=2$ (bottom panel) we observe that the linear response is accurate up to larger $A$. This is expected, since the more perturbations we apply at the same time, the more inhomogeneous the system becomes and, thus, the stronger the total perturbation will be. Further, we note that this effect is more pronounced for $q_{x}=2$. This is again a consequence of the larger $\chi(\mathbf{q})$-value which implies that the density response is even larger in this case.

In a nutshell, we find that, while it is possible to obtain multiple $\mathbf{q}$-values of the response function within a single QMC simulation, this comes at the cost that the linear response is valid only up to smaller perturbation amplitudes $A$. However, the smaller $A$ the larger the relative statistical uncertainty of the induced density, which means that there is a tradeoff between more Monte Carlo steps for a simulation with multiple $\mathbf{q}$-vectors or multiple QMC simulations with only a single perturbation and fewer MC steps. In practice, applying a superposition of $N_{A}\approx 3$ perturbations is reasonable.

IV Summary and Discussion

In summary, we have carried out extensive permutation blocking PIMC simulations of a harmonically perturbed electron gas to investigate the static density response at warm dense matter conditions. To investigate the dependence of the response on the perturbation strength, we varied the amplitude $A$ over three orders of magnitude. For small $A$, linear response theory is accurate and both ways to obtain the response function $\chi(\mathbf{q})$ [i.e., Eqs. (29) and (30)] give equal results. With increasing $A$, the system becomes strongly inhomogeneous which leads to a significantly increased sign problem due to the regions with increased density. The second important issue investigated in this work is the convergence of the PB-PIMC results for $\chi(\mathbf{q})$ with the number of propagators $P$. For small to medium $\mathbf{q}$, we find that $P=4$ propagators are sufficient at WDM conditions, which agrees with previous findings for the uniform system dornheim2 ; dornheim3 . However, for large $\mathbf{q}$, the external potential exhibits fast spatial variations which lead to increased commutator errors and thus require a larger number of propagators to achieve the same level of accuracy. For the largest considered wave vector, $\mathbf{q}=2\pi L^{-1}(5,0,0)^{T}$, at $\theta=1$, $r_{s}=10$, and $N=54$ we find a propagator error of $\Delta\chi/\chi\sim 2\%$. The main effect of system size on the QMC results for the static response function is given by the different $\mathbf{q}$-grid (which is a consequence of momentum quantization in a finite box), whereas the functional form of $\chi(\mathbf{q})$ is remarkably well converged even for small particle numbers. This is in stark contrast to previous findings at zero temperature bowen2 ; moroni2 and can be ascribed to the absence of momentum shell effects at WDM conditions.

Our first brief comparison of the wave vector dependence of $\chi(\mathbf{q})$ computed from QMC to the approximate results from RPA and STLS for $r_{s}=10$ and $\theta=1$ reveals the stark breakdown of the former when coupling effects are non-negligible. The LFC from the STLS closure relation, on the other hand, constitutes a significant improvement, although there remain significant deviations at intermediate $\mathbf{q}$-values. Finally, we have investigated the possibility to obtain the static response function at multiple wave vectors from a single QMC simulation. As predicted by the linear response theory, we found that the density response of the electron gas to a superposition of $N_{A}$ external harmonic perturbations is given by a linear combination of the responses to each of the perturbations. Unfortunately, however, this means that the linear response is valid only up to smaller perturbation amplitudes $A$ as the system becomes increasingly inhomogeneous for multiple $N_{A}$. Thus, there is a tradeoff between $N_{A}$ and $A$, and applying a superposition of $N_{A}=3$ perturbations is a reasonable strategy.

As mentioned in the introduction, accurate QMC results for the static density response function – and, thus, for the static local field correction – are of high importance for contemporary warm dense matter research. Based on the findings of this work, the construction of a comprehensive set of QMC results for $\chi(\mathbf{q})$ over the entire relevant $r_{s}$-range and temperatures $\theta\geq 0.5$ appears to be within reach. First and foremost, this will allow one to systematically benchmark previous approximate results for the warm dense UEG, such as STLS stls ; stls2 (and ”dynamic STLS” bohm ; arora ), VS stls2 ; stolzmann , or the recent improved LFC by Tanaka tanaka_hnc that is based on the hypernetted chain equation, as well as semi-empirical quantum classical mappings pdw_map ; pdw_param . Furthermore, the construction of an accurate parametrization of $G(\mathbf{q};r_{s},\theta)$ with respect to $r_{s}$ and $\theta$ at WDM conditions gregori ; utsumi ; farid is highly desirable due to its utility for, e.g., new DFT exchange-correlations functionals burke2 ; lu ; patrick , the description of Thomson scattering experiments fortmann1 ; fortmann2 , and the construction of pseudopotentials saum1 ; saum2 ; pseudo_potential . Finally, accurate QMC results for the (weakly and strongly) inhomogeneous electron gas can be used as a highly needed benchmark for different exchange-correlation functionals that are used at WDM conditions groth2 ; gga ; gga2 ; hse ; pbe ; groth_prl ; karasiev2 .

Acknowledgements

This work was supported by the Deutsche Forschungsgemeinschaft via project BO1366-10 and via SFB TR-24 project A9 as well as grant shp00015 for CPU time at the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN).