Improved Kelbg Potential in Plasma Simulations

• The improved Kelbg potential is a temperature-dependent quantum statistical potential that interpolates between diffraction-dominated and bound-state regimes in plasma simulations.
• It employs a Padé-type parameterization for the fit parameter γ to accurately reproduce numerically exact potentials across atomic numbers (Z=1–54) with errors within 5%.
• Incorporating finite-size and exchange (Pauli) corrections, the potential enables efficient classical and semiclassical MD simulations of warm dense matter and high energy density plasmas.

The improved Kelbg potential is a class of short-range, temperature-dependent quantum statistical potentials designed for semiclassical simulations of plasmas, especially in regimes where quantum diffraction and partial bound-state formation are non-negligible but where full quantum many-body approaches (e.g., path integral Monte Carlo) are computationally prohibitive. These potentials generalize the original Kelbg pseudopotential by introducing fit parameters to interpolate between purely diffractive and bound-state-dominated regimes and correcting for both exchange (Pauli) effects and, in hydrogenic cases, finite electron size. The resulting formalism is applicable across atomic numbers $Z$ and enables efficient, physically faithful classical or semiclassical molecular dynamics (MD) simulations of warm dense matter (WDM), high energy density (HED) plasmas, and related systems, with broad agreement with quantum calculations in the regime of weak to moderate degeneracy, moderate coupling, and sufficiently high ionization (Demyanov et al., 27 Aug 2025, Whitley et al., 23 Jan 2026).

1. Quantum Statistical Derivation and Analytic Formulation

The improved Kelbg potential originates from the diagonal elements of the exact quantum two-body density matrix,

$$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$

normalized by the free-particle density matrix. The quantum statistical effective potential for particle pair $(i, j)$ with charges $q_i$, $q_j$, reduced mass $\mu_{ij}$, and thermal de Broglie length $\lambda_{ij} = \sqrt{\hbar^2 \beta / 2\mu_{ij}}$ is then defined by

$$U_{ij}(r; \beta) = -\frac{1}{\beta} \ln \left[\frac{\rho_2(\mathbf{r}, \mathbf{r}; \beta)}{\rho_F(\mathbf{r}, \mathbf{r}; \beta)}\right].$$

The improved Kelbg analytic ansatz takes the form

$$U_{ij}^{\mathrm{impr}}(r; \beta) = \frac{q_i q_j}{r} \left[1 - e^{-(r/\lambda_{ij})^2} + \sqrt{\pi} \frac{r}{\lambda_{ij} \gamma_{ij}} \left(1-\mathrm{erf}\left(\gamma_{ij} \frac{r}{\lambda_{ij}}\right)\right)\right]^{-1},$$

where $\gamma_{ij}(T, Z)$ is a dimensionless fit parameter that interpolates between the high-temperature (Kelbg) and low-temperature, bound-state-corrected regimes. For repulsive (e–e) pairs, $$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$0 weakens the repulsion at short range, while for attractive (e–ion) pairs, $$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$1 enhances short-range attraction, mimicking K-shell binding (Whitley et al., 23 Jan 2026).

2. Calibration and Universality Across $$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$2 and $$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$3

For $$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$4 (covering hydrogen up to xenon), $$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$5 is determined via nonlinear least-squares fit to numerically exact diagonal pair potentials computed from Slater sums (matrix squaring, Pollock algorithm). The fit variable,

$$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$6

(where $$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$7 is the Hartree energy), enables a Padé-type parameterization: $$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$8 This functional form reproduces the numerically exact $$\rho_2(\mathbf{r}, \mathbf{r}; \beta) = \sum_s e^{-\beta E_s} |\Psi_s(\mathbf{r})|^2,$$9 to within $(i, j)$0 for all $(i, j)$1 and $(i, j)$2 K, outperforming prior attempts (e.g., Filinov et al.’s hydrogen-only fit) at higher $(i, j)$3. In the classical ($(i, j)$4) limit, $(i, j)$5 and the potential reduces to the original Kelbg form (Whitley et al., 23 Jan 2026).

3. Finite-Size and Pauli Corrections

In hydrogenic systems, the original Kelbg plus Pauli (exchange) correction manifests spurious “same-spin” electron cluster formation below $(i, j)$6 kK. This is rectified by a finite-size correction: in the triplet channel, the Coulomb force is regularized by introducing an effective core size $(i, j)$7, implementing a short-range hard-core repulsion proportional to the electron thermal wavelength. For spin-polarization, a symmetrization term,

$(i, j)$8

is added, with "+" for singlet and "−" for triplet (Demyanov et al., 27 Aug 2025).

For $(i, j)$9, the application of established Pauli terms remains valid under moderate quantum degeneracy. Two widely-used analytic forms are the Lado (low-density) and Deutsch (Hansen–MacDonald) corrections, both vanishing as $q_i$0.

4. Treatment of Long-Range Interactions

The improved Kelbg potential's correct handling of long-range Coulomb tails in finite periodic simulation cells is achieved using the angular-averaged Ewald potential (AAEP). The AAEP introduces a spherically-averaged, finite-range correction to the basic Kelbg form,

$q_i$1

where $q_i$2 is derived from the Ewald summation, ensuring proper energy and pressure scaling with particle number and correct asymptotics up to the simulation cell radius (Demyanov et al., 27 Aug 2025).

5. Applications in Molecular Dynamics Simulations

In classical and semiclassical MD, the improved Kelbg potential, with Pauli and finite-size terms where appropriate, enables accurate modeling of thermodynamic and microphysical properties in partially-quantum plasmas.

Simulations of hydrogen plasmas confirm that, for degeneracy parameter $q_i$3 and coupling parameter $q_i$4, internal energy and pressure agree with benchmark path integral Monte Carlo (PIMC) to better than $q_i$5 for $q_i$6 50 kK. Below this threshold, MD overpredicts binding and underestimates energy due to residual inaccuracies in the potential gradient, and fails to reproduce correct molecular fractions (Demyanov et al., 27 Aug 2025). For carbon ($q_i$7) and other higher-$q_i$8 plasmas, MD with the improved potential plus a Pauli term (Lado or Deutsch) recovers DFT+PIMC EOS results within $q_i$9 provided the K-shell is more than 50% ionized and quantum degeneracy is weak (Whitley et al., 23 Jan 2026).

System	Temperature Regime	EOS Agreement
Hydrogen	$q_j$0 kK	$q_j$1 1% (P, E)
Carbon	$q_j$2 K, $q_j$3 g/cc	$q_j$4 10% (P, E)
H, C, $q_j$5	$q_j$6 low, $q_j$7 or strong degeneracy	Breakdown; clustering

Clustering (“faux molecules”), EOS deviation, and simulation divergence demarcate the boundary where fully quantum methods are required.

6. Validity Domain, Limitations, and Physical Implications

The improved Kelbg potential's applicability is robust across $q_j$8 for WDM conditions where:

• The bound K-shell (or more generally, tightly-bound subshells) is largely ionized ($q_j$9).
• The system is not strongly quantum degenerate ($\mu_{ij}$0).
• Three-body quantum effects are negligible (below the Barker line).

A key limitation is the inability to describe molecular formation, strong recombination, and high-density quantum effects. At lower temperatures or higher densities (where bound states become dominant or the quantum degeneracy parameter becomes large), unphysical classical clustering emerges, and the MD approach fails. In these cases, path-integral or explicitly quantum many-body simulation methods are required (Whitley et al., 23 Jan 2026).

7. Implementation and Practical Usage

The closed-form improved Kelbg potential (with its single-parameter Padé fit) is readily implemented in classical and semiclassical MD codes (e.g., ddcMD, LAMMPS, Sarkas). Its universal form for $\mu_{ij}$1 ensures convenience for EOS, transport, and microphysical studies across a spectrum of elements. For regimes where fast computation is essential and bound-state formation is sufficiently suppressed by temperature or ionization, the improved Kelbg potential provides an efficient and accurate alternative to ab initio quantum simulation, capturing key quantum diffraction and exchange features at a fraction of the computational cost (Whitley et al., 23 Jan 2026).

A plausible implication is the accelerated exploration and refinement of EOS tables and transport properties for astrophysical, inertial confinement fusion, and planetary modeling, particularly in parameter regions where quantum effects are neither entirely negligible nor fully dominant.