Modeling partially-ionized dense plasma using wavepacket molecular dynamics

TL;DR

This paper extends wave packet molecular dynamics (WPMD) to partially ionized dense plasmas by explicitly including bound-state electrons and a mixed bound-free electronic representation. The authors implement Ehrenfest dynamics for ions, a constrained electronic state with bound orbitals attached to nuclei, and a Gaussian, anisotropic description for free electrons, augmented by Pauli potentials to mimic exchange effects and a harmonic confining potential to regulate electron spread. Self-consistent charge-state distributions are obtained through thermodynamic integration free-energy minimization across ionization states, and hydrogen serves as a benchmark by comparing structural observables to path-integral Monte Carlo data. The results show that including bound electrons reduces sensitivity to the confinement parameter and improves agreement with reference data under partially ionized conditions, highlighting the framework’s potential to extend semi-classical plasma models to more complex species and higher-Z systems.

Abstract

We develop a wave packet molecular dynamics framework for modeling the structural properties of partially-ionized dense plasmas, based on a chemical model that explicitly includes bound state wavefunctions. Using hydrogen as a representative system, we compute self-consistent charge state distributions through free energy minimization, following the approach of Plummer et al. [Phys. Rev. E 111, 015204 (2025)]. This enables a direct comparison of static equilibrium properties with path integral Monte Carlo data, facilitating an evaluation of the model's underlying approximations and its ability to capture the complex interplay between ionization and structure in dense plasma environments.

Figures (8)

• Figure 1: Radial distribution function (RDF) dependence on confinement parameter $\sigma_0$ for the fully ionized wave packet model. The results are compared against reference path integral Monte Carlo data from Ref. Dornheim2024 and the ideal (uncorrelated) limit. Panels (a) and (c) correspond to the proton-proton and radially weighted electron-proton RDFs at the partially ionized condition: $r_s = 3.23$ and $T=55700 \, \text{K}$. Panels (b) and (d) correspond proton-proton and radially weighted electron-proton RDFs at the fully ionized condition: $r_s = 2$ and $T=145000 \, \text{K}$. Hartree atomic units are used in the legend.
• Figure 2: Each step of the ionization calculation for $\sigma_0 = 1.1 \, a_0$. Panel (a) shows the time trace of the scaled potential energy function $\mathcal{U}(\lambda)$ and the unscaled energy $V$ across a single molecular dynamics run for $\bar{z}=1$. The thin black dotted lines indicate the times at which the system coupling parameter is modified, with values given on the top axis.
• Figure 3: Minimized ionization state $\bar{z}_{\text{min}}$ plotted against confinement parameter $\sigma_0$ at $r_s = 3.23$ and $T = 55700 \, \text{K}$. Results are presented alongside path integral Monte Carlo inferred results from Ref Bellenbaum2025 and the one component plasma (OCP) model from Ref. Plummer2025. This data may be found in Table \ref{['tab:ionization_calculations']}.
• Figure 4: Radial distribution functions (RDFs) from the bound-WPMD model for different values of the confinement parameter $\sigma_0$. Self-consistent ionization states are computed via free energy minimization and found in Table \ref{['tab:ionization_calculations']}. (a) proton-proton (pp). (b) electron-proton (ep).
• Figure 5: Radial distribution functions (RDFs) at $r_s = 3.23$ and $T = 55{,}700\ \mathrm{K}$ computed with the bound-WPMD model for $\bar{z} = 0.44$, compared with path-integral Monte Carlo data, fully ionized ($\bar{z} = 1$) results, and the ideal uncorrelated case. (a) Ion-ion (${ii}$), ion--neutral (${in}$), and neutral--neutral (${nn}$) RDFs, combined via Eq. \ref{['eq:RDF_pp_decomp']} to yield the total proton-proton RDF ($pp$). (b) Free-proton ($fp$) and bound-proton (${bp}$) RDFs, combined via Eq. \ref{['eq:RDF_ep_decomp']} to yield the total electron-proton RDF ($ep$).