Path-integral Monte Carlo estimator for the dipole polarizability of quantum plasma

TL;DR

This work develops a path-integral Monte Carlo estimator for the dipole polarizability of interacting quantum plasmas in the optical, long-wavelength limit, leveraging imaginary-time dipole autocorrelation functions. By using Boltzmann statistics to avoid Fermion sign problems and carefully treating periodic boundaries, the method yields 𝔊̃(τ) and its Matsubara transforms that can be compared to the Drude/Lindhard reference. The results show near-perfect agreement (within ~1%) with the Drude model across metallic densities and moderate temperatures, indicating the absence of strong many-body biases in the tested regime and establishing a robust first-principles baseline for polarizability. The approach holds promise for extending to nonlinear responses, confinement, dispersion forces, and plasmonic materials, and provides a pathway to explore effective damping and higher-order optical properties from quantum simulations.

Abstract

We present a path-integral Monte Carlo estimator for calculating the dipole polarizability of interacting Coulomb plasma in the long-wavelength limit, i.e., the optical region. We present comprehensive details and method validation studies for our approach based on dipole imaginary-time autocorrelation functions. The simulation of thermal equilibrium in imaginary time has exact Coulomb interactions and Boltzmann quantum statistics. For reference, we use the Drude model as the long-wavelength limit of the Lindhard response, presenting its analytical continuation into the imaginary time and Matsubara series. We demonstrate great agreement within $1\%$ between PIMC and the reference model, indicating negligible numerical biases and physical many-body effects in metallic densities. The approach is amenable to nonlinear optical response, quantum confinements, dispersion forces, and to inform applications such as plasmonics and epsilon-near-zero materials.

Figures (6)

• Figure 1: Schematic illustration of the periodic dipole moment estimator and sampling constraints. Blue circles represent Monte Carlo walkers of quantum particles, i.e., positions advancing in discrete steps through imaginary time and forming closed loops $z(0) = z(\beta)$ modulo $L$, where $z$ denotes the real coordinate in one dimension. (left) The mean dipole autocorrelation function is easily measured by summing up individual particles, here represented by $\mu_z(\tau)/qN$ for $N=2$ by the smaller green dots. Their mean squared fluctuation w.r.t. to the mean (green vertical line) induces the relative autocorrelation function $\mathcal{\tilde{G}}(\tau)$ (visualized by the width of a green distribution). (middle) Unfolding of the path. Since the estimator does not treat periodicity with a Fourier phase, paths crossing the simulation box boundary to the nearest image must be unfolded, i.e., moved within $L/2$ of the earlier bead, before measurement. (right) Winding restriction. MC sampling may cause the path to wind across box boundary (branch of red dots), when the thermal wavelength of a multilevel bisection move is a significant fraction of $L$. For distinguishable particles, this is clearly an artifact of the finite cell and also detrimental to the estimation of $\mathcal{\tilde{G}}(\tau)$, and thus, such moves are rejected.
• Figure 2: Comparison overview of the relative dipole correlation function $\mathcal{\tilde{G}}$ from PIMC and the Drude model. In the top, the imaginary-time correlation function $\mathcal{\tilde{G}}(\tau)$ is plotted at two temperatures, $\Theta_1=0.1$ and $\Theta_2=0.2$, resulting in two curves ending at different values for $\beta_1 = 2 \beta_2$. In the bottom, the first few corresponding Matsubara data $\mathcal{\tilde{G}}(i \omega_n)$ are plotted, where the frequencies $\omega_n = 2 \pi n/\hbar \beta$ for $\Theta_1$ are twice those of $\Theta_2$, but the two functions lie on the same curve. In both cases, the data match the analytically continued Drude model within statistical uncertainties.
• Figure 3: Scaling of per-particle kinetic energy $K$ (left) and potential energy $V$ (bottom) with finite $N$ and variable dimensional Ewald cutoff parameter $r_c k_c$. The solid lines correspond to raw PIMC observables and the dotted lines are finite-size corrected values based on Eqs. \ref{['eq:K_corr']} and \ref{['eq:V_corr']}. The case $r_c k_c=0$ means no analytical Ewald summation of any kind, including neutralizing background, and so $V$ is off above the plot.
• Figure 4: Scaling of $\mathcal{\tilde{G}}(\tau)$ (left) and $\mathcal{\tilde{G}}(i \omega_n)$ (right) with finite $N$. The top panels show absolute values of the PIMC data (solid) and the corresponding Drude value (dotted). The bottom panels shows the relative differences $\Delta \mathcal{\tilde{G}}(\tau)$ (solid lines) and $\Delta \mathcal{\tilde{G}}(i \omega_n)$, including statistical uncertainties (dotted lines).
• Figure 5: On the left, $\mathcal{\tilde{G}}(\tau)$ and $\Delta \mathcal{\tilde{G}}(\tau)$ with variable $\Theta$ at $r_s=4$. On the right, $\mathcal{\tilde{G}}(i\omega_n)$ and $\Delta \mathcal{\tilde{G}}(i\omega_n)$ with variable $\Theta$ at $r_s=4$. Solid lines represent the absolute data and dotted lines the statistical uncertainties.