# Development of a new quantum trajectory molecular dynamics framework

Pontus Svensson  Department of Physics, University of Oxford, Parks Road,
Oxford OX1 3PU, UK  Thomas Campbell  Department of Physics, University of
Oxford, Parks Road, Oxford OX1 3PU, UK  Frank Graziani  Lawrence Livermore
National Laboratory, Livermore, CA 94550, USA  Zhandos Moldabekov  Center of
Advanced Systems Understanding (CASUS), Helmholtz-Zentrum Dresden-Rossendorf,
Untermarkt 20, 02826 Görlitz, Germany  Ningyi Lyu  Department of Chemistry,
Yale University, New Haven, CT 06520, USA  Scott Richardson  AWE, Aldermaston,
Reading, Berkshire RG7 4PR, UK  Sam M. Vinko  Department of Physics,
University of Oxford, Parks Road, Oxford OX1 3PU, UK  Central Laser Facility,
STFC Rutherford Appleton Laboratory, Didcot OX11 0QX, UK  Gianluca Gregori
Department of Physics, University of Oxford, Parks Road, Oxford OX1 3PU, UK

###### Abstract

An extension to the wave packet description of quantum plasmas is presented,
where the wave packet can be elongated in arbitrary directions. A generalised
Ewald summation is constructed for the wave packet models accounting for long-
range Coulomb interactions and fermionic effects are approximated by purpose-
built Pauli potentials, self-consistent with the wave packets used. We
demonstrate its numerical implementation with good parallel support and close
to linear scaling in particle number, used for comparisons with the more
common wave packet employing isotropic states. Ground state and thermal
properties are compared between the models with differences occurring
primarily in the electronic subsystem. Especially, the electrical conductivity
of dense hydrogen is investigated where a $15\%$ increase in DC conductivity
can be seen in our wave packet model compared to other models.

##  1 Introduction

The establishment of high power lasers facilities during the last decades has
been instrumental in the achievements towards inertial confinement fusion
(ICF) [1, 2, 3], but also for the creation of high-density and high-
temperature conditions [4] otherwise only found in astrophysical objects [5,
6]. Furthermore, X-ray lasers are now able to reach complementary high-
pressure regions in phase space [7, 8]. One of the exotic states now
accessible is warm dense matter (WDM) which exists in gas giants [9, 10, 11,
12, 13], brown [14] and white dwarf stars [15, 16], the crust of neutron stars
[17, 18], and during the compression of an ICF capsule [19]. Warm dense matter
is a strongly coupled quantum plasma, with ions moving in a partially
degenerate electron fluid with kinetic energy comparable to the ion-ion
interaction energy [20]. Consequently, WDM inherits properties from both
condensed matter systems and classical plasmas, a challenging combination to
model. Various computational techniques are commonly used to describe the
system, yielding similar thermodynamic [21] and acoustic properties [22, 23],
although dynamic properties differ by orders of magnitude [24]. These
uncertainties limit our understanding of, for example, the Jovian interior
[25], or the modelling of ICF implosions [26].

The three main complications in modeling WDM are electron degeneracy, strong
ion correlations and the separation in timescales between the electron and ion
dynamics. A full solution would require a quantum mechanical treatment of the
electrons, resolving electron dynamics while considering phenomena on the ion
time scale. Consequently, explicit models of ionic motion span a wide range of
theories, including classical systems with effective ion-ion interactions [27,
28, 29], classical electrons with effective quantum statistical potentials
(QSP) [30, 31, 32], Bohmian mechanics [33], density functional theory
molecular dynamics (DFT-MD) using both orbital-free [34, 35] and Kohn-Sham
[36, 37, 12] DFT-variants, phenomenological quantum hydrodynamics based on
DFT-functionals [38, 39, 40], time-dependent DFT [41, 42, 43], and quantum
Monte-Carlo and path integral Monte-Carlo [44, 45, 13] approaches. Coarse-
grained models with effective interactions are fundamentally based on
reconstructing some equilibrium property, the choice of which is arbitrary and
limited to a specific thermodynamic condition, whereas experimental
realisations are commonly non-stationary [46, 47, 48, 49, 50, 7]. Furthermore,
models rooted in the Born-Oppenheimer approximation – where the electrons are
treated adiabatically – e.g. DFT-MD cannot capture a dynamic electron
response, believed to be important for the description of dynamic properties
such as some transport coefficients [51], stopping power [39] and energy
transfer between the electronic and ionic subsystems. However, time-dependent
approaches are computationally costly, and are typically limited in terms of
particle numbers and time scales of studied phenomena.

Wave packet molecular dynamics (WPMD) [52, 53] is a family of models in which
the electron dynamics are computed explicitly, whiles simulating hundreds to
thousands of particles over ionic time scales. This is made possible by
restricting the wave function of each electron to a parameterised functional
form. We present an extension to existing wave packet formulations –
applicable to the WDM regime – in which the wave packets can be elongated in
arbitrary directions. The model accounts for the long-range behaviour of
electrostatic interactions and of fermionic properties by effective Pauli
interactions, while implemented within the scalable molecular dynamics
framework LAMMPS [54] to treat systems with thousands of particles.

In the following section the theoretical model is described, after which
section 3 outlines the numerical details and performance of the
implementation. The model is compared with other computational techniques in
section 4, where we apply it to ground state and dynamic properties of a dense
hydrogen plasma. We compute some structural and transport properties, which
are compared with an isotropic wave packet model. We conclude with a summary
of our results.

##  2 Theoretical Description

Originally proposed in the 1970s as an approximate solution to Schrödinger’s
equation [55, 56], wave packet models can systematically be derived from
variations of the action,

$\mathcal{S}=\int\\!dt\,\bra{Q}i\hbar\frac{d}{dt}-\hat{H}\ket{Q},$ (1)

where $\hat{H}$ is the system Hamiltonian and the state,
${\ket{Q}=\ket{Q(Q_{\mu})}}$, is restricted to some manifold, $\mathcal{M}$,
defined by the adopted wave packets and parameterised by its parameters
$Q_{\mu}$. The resulting time evolution reproduces the true quantum dynamics
to the best of its ability being restricted to the manifold, $\mathcal{M}$,
during short time scales of length $\delta t$. Concretely, it can be shown
$\bra{\Delta(t,\delta t)}\ket{\Delta(t,\delta t)}$ is minimised to
$\mathcal{O}(\delta t^{3})$ where $\ket{\Delta(t,\delta t)}=\ket{\Psi(t+\delta
t)}-\ket{Q(t+\delta t)}$ and $\ket{\Psi(t+\delta t)}$ is the true solution to
Schrödinger’s equation starting from $\ket{\Psi(t)}=\ket{Q(t)}$ [52]. The
long-time evolution is constrained by appropriate conservation laws, most
notably energy conservation [57].

In general, the equations of motion are quasi-Hamiltonian

$i\hbar\sum_{\nu}\mathcal{C}_{\mu\nu}\frac{dQ_{\nu}}{dt}=\frac{\partial\mathcal{H}}{\partial
Q^{*}_{\mu}}$ (2)

where

$\mathcal{H}=\frac{\bra{Q}\hat{H}\ket{Q}}{\bra{Q}\ket{Q}}\equiv\langle\hat{H}\rangle\quad\text{and}\quad\mathcal{C}_{\mu\nu}=\frac{\partial^{2}}{\partial
Q_{\mu}^{*}\partial
Q_{\nu}}\ln\Big(\bra{Q(Q_{\mu}^{*})}\ket{Q(Q_{\nu})}\Big{missing}).$ (3)

Fermions are described by states anti-symmetric under exchange and Slater-
determinants have been considered in Refs. [58, 59, 60]. However, this
approach scales unfavourably with particle number $N$. Instead, here we employ
a product state

$\ket{Q}=\ket{q_{1}}\otimes\ket{q_{2}}\otimes\dots\otimes\ket{q_{N}}$ (4)

of single-particle orbitals, $\ket{q_{i}}$, and $\otimes$ is the tensor
product. Exchange effects are approximated by Pauli potentials of the type
first introduced by Klakow et al. [61, 62]. This structure simplifies
$\mathcal{C}_{\mu\nu}$, which becomes block diagonal, and the orbitals
$\ket{q_{i}}$ only couple through the energy $\mathcal{H}$.

###  2.1 Wave Packets

The choice of wave packet shape is central to the model, dictating the states
that can be described. Most commonly, isotropic gaussians are utilised –
primarily motivated by computational ease – yet other variants exist, see Ref.
[53] and references therein. To account for local gradients, an anisotropic
wave packets form is introduced

$\bra{\boldsymbol{x}}\ket{q_{i}}=\left(\left(2\pi\right)^{3}\det\left(\boldsymbol{\Sigma}_{i}\right)\right)^{-1/4}\times\exp\left[-\boldsymbol{\xi}_{i}\left(\frac{1}{4}\boldsymbol{\Sigma}_{i}^{-1}-\frac{i}{\hbar}\boldsymbol{\Pi}_{i}\right)\boldsymbol{\xi}_{i}+\frac{i}{\hbar}\boldsymbol{p}_{i}^{\intercal}\boldsymbol{\xi}_{i}\right],$
(5)

where $\boldsymbol{\xi}_{i}=\boldsymbol{x}-\boldsymbol{r}_{i}$. The wave
packet is parameterised by 18 degrees of freedom, the position
$\boldsymbol{r}_{i}$, momentum $\boldsymbol{p}_{i}$ and two symmetrical
matrices, $\boldsymbol{\Sigma}_{i}$, describing the elongation and
orientation, and $\boldsymbol{\Pi}_{i}$ the associated momentum to
$\boldsymbol{\Sigma}_{i}$. A similar type of wave packet has been treated
previously, describing molecular binding in water molecules [63], and is
generally believed to improve the description of molecular states [53].

The equations of motion for the functional form (5) have a classical-looking
structure, for the “classical” degrees of freedom [63]

$\frac{d\boldsymbol{r}_{i}}{dt}=\frac{\partial\mathcal{H}}{\partial\boldsymbol{p}_{i}},\hskip
14.22636pt\frac{d\boldsymbol{p}_{i}}{dt}=-\frac{\partial\mathcal{H}}{\partial\boldsymbol{r}_{i}}$
(6a) and the “internal” dynamics of the wave packet follow as
$\frac{d}{dt}\boldsymbol{\Sigma}_{i}=\boldsymbol{\tau}\odot\frac{\partial\mathcal{H}}{\partial\boldsymbol{\Pi}_{i}},\hskip
14.22636pt\frac{d}{dt}\boldsymbol{\Pi}_{i}=-\boldsymbol{\tau}\odot\frac{\partial\mathcal{H}}{\partial\boldsymbol{\Sigma}_{i}},$
(6b)

where the matrix $\boldsymbol{\tau}$ (unity on the diagonal and one half
elsewhere) accounts for the symmetric structure of $\boldsymbol{\Sigma}_{i}$
and $\boldsymbol{\Pi}_{i}$ by the elementwise multiplication operator, denoted
by $\odot$. Specifically, we consider a charged system of classical ions, with
position $\boldsymbol{R}_{I}$, momentum $\boldsymbol{P}_{I}$, charge $Z_{I}e$
and mass $M_{I}$, and quantum electrons with position
$\hat{\boldsymbol{x}}_{i}$ and momentum $\hat{\boldsymbol{p}}_{i}$ operators
as well as charge $-e$ and mass $m$. The system is described by the
Hamiltonian

$\hat{H}=\sum_{I}\frac{\boldsymbol{P}_{I}^{2}}{2M_{I}}+\sum_{I<J}\frac{Z_{I}Z_{J}e^{2}}{|\boldsymbol{R}_{I}-\boldsymbol{R}_{J}|}+\sum_{i}\frac{\hat{\boldsymbol{p}}_{i}^{2}}{2m}+\sum_{i<j}\frac{e^{2}}{|\hat{\boldsymbol{x}}_{i}-\hat{\boldsymbol{x}}_{j}|}-\sum_{i}\sum_{I}\frac{Z_{I}e^{2}}{|\hat{\boldsymbol{x}}_{i}-\boldsymbol{R}_{I}|},$
(7)

the state average of which is required for the time evolution. The average
kinetic energy

$\left\langle\frac{\hat{\boldsymbol{p}}_{i}^{2}}{2m}\right\rangle=\frac{\boldsymbol{p}_{i}^{2}}{2m}+\frac{2}{m}\Tr\\{\boldsymbol{\Pi}_{i}\boldsymbol{\Sigma}_{i}\boldsymbol{\Pi}_{i}\Big{\\}}+\frac{\hbar^{2}}{8m}\Tr\\{\boldsymbol{\Sigma}_{i}^{-1}\Big{\\}}$
(8)

includes both a classical contribution and a part internal to the wave packet.
The last term in equation (8) is the so-called shape-kinetic energy [64],
which keeps $\boldsymbol{\Sigma}_{i}$ positive definite and the wave-packet
well defined during the time evolution. The interaction terms have not been
evaluated explicit and the following section is dedicated to the treatment of
these terms.

###  2.2 Generalised Ewald summation

Within molecular dynamics, it is desirable to truncate pair-interactions at
some distance such that the computation formally scales as
$\mathcal{O}\left(N\right)$ [65]. However, in our case the electrostatic
interaction is long-range [66] and it is beneficial to perform the split [67,
68]

$\frac{1}{r}=\frac{\text{erfc}\left(gr\right)}{r}+\frac{\text{erf}\left(gr\right)}{r},$
(9)

chosen so that the first term can be truncated at a distance of order
$g^{-1}$, while the second term is regular as $r\rightarrow 0$ and efficiently
evaluated in Fourier space. The Ewald parameter $g$ is chosen to optimise
performance. Below we present a self-consistent treatment of both terms, as
the long-range part has only been mentioned once for isotropic wave packets
[69] and is commonly neglected.

#### Short-range forces:

In the case of a gaussian interaction kernel, the required state average can
promptly be evaluated [70, 71]. Therefore, we construct a gaussian
decomposition of the interaction kernel,

$\frac{\text{erfc}\left(gr\right)}{r}\simeq\sum_{p}c_{p}e^{-\alpha_{p}r^{2}},$
(10)

where the coefficients $c_{p}$ and $\alpha_{p}$ are fitting parameters. A
robust numerical scheme to perform the decomposition is described in appendix
A, where typically only 5-15 modes are required. By approximating the
potential form, the notion of energy conservation is retained. We note that
this is not the case for methods based on Taylor expansions [63] or on direct
numerical evaluations [72], due to either truncation errors or numerical
noise.

#### Long-range forces:

To limit surface effects in a finite size simulation, the simulation box is
periodically repeated. Periodic images are included in accordance with the
standard treatment of Ewald summation with particles positioned at
$\boldsymbol{r}_{i}+L\boldsymbol{n}$ for all $\boldsymbol{n}\in\mathbb{Z}^{3}$
and $L$ being the length of a cubic simulation cell. The interaction energy is

$\sum_{i<j}\left\langle
V(\boldsymbol{x}_{ij})\right\rangle\rightarrow\frac{1}{2}\sum_{i,j}^{\prime}\sum_{\boldsymbol{n}\in\mathbb{Z}^{3}}\left\langle
V(\boldsymbol{x}_{ij}\shortminus
L\boldsymbol{n})\right\rangle=\underbrace{\frac{1}{2}\sum_{i,j}\sum_{\boldsymbol{n}\in\mathbb{Z}^{3}}\left\langle
V(\boldsymbol{x}_{ij}\shortminus
L\boldsymbol{n})\right\rangle}_{\mathcal{E}_{k}}-\underbrace{\frac{1}{2}\sum_{i=j}\left\langle
V(\boldsymbol{x}_{ij})\right\rangle}_{\mathcal{E}_{s}},$ (11)

where $V$ is the long-range part of the Coulomb interaction in equation (9).
The special case $i=j$ is excluded when $\boldsymbol{n}=0$ (denoted by the
primed sum) resulting in two distinct terms the main contribution
$\mathcal{E}_{k}$ and the self-energy $\mathcal{E}_{s}$. In appendix B, we
evaluate $\mathcal{E}_{k}$ in reciprocal space to be

$\mathcal{E}_{k}=\frac{1}{2L^{3}}\sum_{\boldsymbol{k}\neq
0}\frac{4\pi}{k^{2}}e^{-k^{2}/4g^{2}}|\tilde{\rho}_{\text{uc}}(\boldsymbol{k})|^{2},$
(12)

where $\tilde{\rho}_{\text{uc}}(\boldsymbol{k})$ is the charge density

$\tilde{\rho}_{\text{uc}}(\boldsymbol{k})=e\sum_{i}Z_{i}e^{-i\boldsymbol{k}\cdot\boldsymbol{r}_{i}}e^{-\frac{1}{2}\boldsymbol{k}^{\intercal}\boldsymbol{\Sigma}_{i}\boldsymbol{k}},$
(13)

and ${\boldsymbol{k}=2\pi\boldsymbol{n}/L}$. Equation (12) converges rapidly
due to the exponential factor. The self-energy term, $\mathcal{E}_{s}$, is
independent of $\boldsymbol{r}_{i}$ and does not influence the classical
degrees of freedom directly. Although for distributed particles, it does
depend on $\boldsymbol{\Sigma}_{i}$ and influences the internal dynamics.
Based on the decomposition of the long-range interaction kernel (appendix A)
the self-energy is evaluated in appendix B.

###  2.3 Pauli interactions

The difference in the kinetic energy between a pair-wise anti-symmetrised
state and the product state has often been used as a Pauli potential to
correct for the fermionic structure of electrons [61, 62, 69]. The electron
force field model (eFF) introduced fitting parameters in the Pauli interaction
to achieve stable bounds for elements with $Z\leq 6$ [73] and has been widely
used with minor modifications [74, 75, 76, 77, 78, 22, 23, 79]. Ref. [80]
considered exchange contributions to the interaction terms while including a
correlation potential with a free parameter. This treatment of the Pauli
interaction is extended here to anisotropic gaussian states.

We construct the potential by considering two particles, $i$ and $j$, with the
Hamiltonian

$\hat{H}_{2}=\frac{\hat{\boldsymbol{p}}^{2}_{i}+\hat{\boldsymbol{p}}^{2}_{j}}{2m}+\frac{e^{2}}{|\hat{\boldsymbol{x}}_{i}-\hat{\boldsymbol{x}}_{j}|}+V_{\text{bg}}\left(\hat{\boldsymbol{x}}_{i},t\right)+V_{\text{bg}}\left(\hat{\boldsymbol{x}}_{j},t\right)$
(14a) with a background interaction from all other particles in the system
$V_{\text{bg}}\left(\boldsymbol{x},t\right)=\\!\\!\sum_{k\neq i,j}e^{2}\int
d^{3}\boldsymbol{x}_{k}\frac{|\bra{\boldsymbol{x}_{k}}\ket{q_{k}}|^{2}}{|\boldsymbol{x}-\boldsymbol{x}_{k}|}-\sum_{I}\frac{Z_{I}e^{2}}{|\boldsymbol{x}-\boldsymbol{R}_{I}|},$
(14b)

where $k$ ($I$) runs over all other electrons (ions) in the system. The two-
electron system can be characterised based on its spin structure either as a
singlet or a triplet state, requiring either a symmetric or anti-symmetric
spatial state

$\text{Singlet:
}\ket{q_{i}}\otimes\ket{q_{j}}+\ket{q_{j}}\otimes\ket{q_{i}}\qquad\text{Triplet:
}\ket{q_{i}}\otimes\ket{q_{j}}-\ket{q_{j}}\otimes\ket{q_{i}}$ (15)

written here in terms of single particle orbitals $\ket{q_{i}}$.

For equal spin particles only the spatially anti-symmetric state is allowed
and the Pauli potential $\mathcal{V}_{ij}^{P}$ is the difference between
$\langle\hat{H}_{2}\rangle$ for the triplet and the product state
$\ket{q_{i}}\otimes\ket{q_{j}}$. In the case of opposite spin particles, the
spatial state depends on the spin structure, but along the lines of Ref. [80]
a correlation potential is introduced based on the singlet state,
$\rho\mathcal{V}_{ij}^{C}$, multiplied by the parameter $\rho$. Therefore, the
potentials are

$\mathcal{V}^{P/C}_{ij}=\mp\frac{\bra{q_{i}q_{j}}\hat{H}_{2}\ket{q_{j}q_{i}}-\bra{q_{i}q_{j}}\hat{H}_{2}\ket{q_{i}q_{j}}\left|\bra{q_{i}}\ket{q_{j}}\right|^{2}}{1\mp\left|\bra{q_{i}}\ket{q_{j}}\right|^{2}},$
(16)

which can be shown to scale as $\left|\bra{q_{i}}\ket{q_{j}}\right|^{2}$ for
large particle separations. Due to gaussians being localised states, this
gives a short-range interaction between particles $i$ and $j$. The background
term $V_{\text{bg}}$ introduces a long-range dependence in terms of the third
particle and is accounted for by a type of Ewald summation, see appendix C.
The state averages of the interaction terms in equation (16) are evaluated
based on the gaussian mode decomposition.

The correlation potential, the spin interaction between opposite spin
electrons, is constructed based on the same premise as the Pauli potential in
a pair-wise approximation, however, with an additional parameter $\rho$ which
needs to be chosen a-priori. In the case of a ground state helium or molecular
hydrogen, the electronic structure is well described by a single state and
$\rho=1$ is an appropriate choice. For a free electron gas, this would
overestimate the correlation effects as the appropriate two-particle state is
not simply the singlet state, and therefore $\rho<1$ is more suitable.

###  2.4 Confining potentials

It has been well documented that at significantly high temperatures wave
packets tend to expand indefinitely [69, 81, 58, 82, 83, 22] and the wave
packet may extend over multiple ions without the ability to localise on
multiple sites [84]. If a wave packet is spread too large, it effectively
ceases to interact with other particles as the charge density effectively
becomes flat. Multiple approaches to counter this expansion has been proposed,
see Refs. [82, 85, 78]. Currently, we employ an additional potential energy
term on the form

$\mathcal{V}_{\Sigma}=\frac{1}{2}A_{w}\sum_{i=1}^{N}\sum_{\alpha=1}^{3}\left(\sigma_{i,\alpha}-l_{w}\right)^{2}\theta\\!\left(\sigma_{i,\alpha}^{2}-l_{w}^{2}\right)$
(17)

where $\sigma_{i,\alpha}^{2}$ is the $\alpha$th eigenvalue of
$\boldsymbol{\Sigma}_{i}$ and $\theta\\!\left(x\right)$ is the Heaviside step
function. The parameters $l_{w}$ and $A_{w}$ set the width of a free particle
$\sigma_{\text{free}}(l_{w},A_{w})\gtrsim l_{w}$ by balancing the shape-
kinetic energy. This potential is rotationally invariant and acts only on wave
packets with a width larger than $l_{w}$. Furthermore, the confinement reduces
to the commonly used potential based on a harmonic potential centred at the
particle position in the limit of $l_{w}\rightarrow 0$. In this specific
limit, the potential has also been used to address the heat capacity in the
classical limit [86].

##  3 Numerical Realisation

The standard velocity-Verlet integrator almost exclusively used for MD
simulations is not appropriate for our model because of the momentum-dependent
Pauli potentials [52] resulting in a non-separable Hamiltonian. This prevents
a straightforward generalisation of the velocity-Verlet algorithm which is
based on the ability to separate the Hamiltonian into terms where the dynamics
following from each term in isolation can be solved exactly [87]. Explicit
Runge-Kutta methods of orders $2$ and $4$ are employed instead. Furthermore,
this momentum dependence of the potential affects the interpretation of
temperature in the system, further described in appendix D.

The time integrator, the generalised Ewald summation, and the Pauli
interaction, are all natively implemented in LAMMPS [54], a MD framework
written in C++ which utilises MPI to distribute the computation [88]. Figure
1$(b)$ shows the computational time for a varying degree of parallelisation
for a fixed system of 2000 protons and an equal number of electrons. In
particular, a good scaling of the pair interaction and the Ewald summation is
established as the computation is distributed. The synchronisation time, the
time different processes need to wait for each other due to an unbalanced load
caused by statistical fluctuations in the number of particles in the region
assigned to each processor, limits the efficiency of the parallelisation when
only a few particles are assigned to each process. In the future, dynamic load
balancing could potentially address this issue.

Furthermore, figure 1$(a)$ also demonstrates the scaling of computational cost
with particle number for a test system where
$r_{s}=(4\pi/3\;na_{B}^{3})^{-1/3}=2$ where $n$ is the proton number density
and $a_{B}$ is the Bohr radius. Close to linear scaling in particle number is
demonstrated, showing the feasibility of employing this modeling technique for
large systems of particles. The exact exponent varies in the range $1.1-1.3$,
depending on the degree of parallelisation, caused by different limiting
factors in the computation. The synchronisation time is most likely one of
these factors for the high parallelisation case.

Figure 1: Scaling of computational cost for a quasi-neutral system with $N$
protons, density corresponding to $r_{s}=2$ and temperature equal to the Fermi
temperature. $(a)$ Scaling with particle number for different levels of
parallelisation, in the range between $N^{1.1}$ and $N^{1.3}$. The Ewald
parameter $g$ was optimised to within $\pm 0.5a_{B}^{-1}$ for each case. $(b)$
MPI parallelisation of the computation for a system of $N=2000$ ions with a
fixed Ewald parameter. Showing the full computation (circles), pair
interaction (squares), Ewald summation (diamonds), communication between MPI
processes (upwards triangles), and the synchronisation time different
processors need to wait for each other due to an unbalanced load (downwards
triangles), which is compared to the optimal scaling based on the single
thread performance (dashed).

##  4 Test Systems

###  4.1 Ground state properties

Ground state properties of the wave packet models can be obtained by the
introduction of a generalised friction term into the equations of motion (6b).
Some care is needed to guarantee continuous energy loss due to the momentum
dependence of the Pauli potential. This is further described in appendix E.

The ground state of isolated atoms is spherically symmetric and does not
utilise the additional degrees of freedom of the wave packets. One of the
simplest physical systems which naturally breaks this symmetry are diatomic
molecules and in particular diatomic hydrogen ($\text{H}_{2}$), where the
ability for the wave packet to stretch is believed to be crucial for molecular
binding [53]. The ground state energy of $\text{H}_{2}$ within the wave packet
model for both elongated and isotropic gaussians is shown in figure 2 for a
varying nuclear separation $\delta$. The elongated wave-packets demonstrate an
improvement over the isotropic ones when $\delta<2.8a_{B}$, above which the
electron density is localised on each ion and close to spherically symmetric.
Furthermore, the isotropic model transitions to an electron density localised
on each nucleus at a significantly shorter nuclear separation, $\delta\approx
1.8a_{B}$, compared to our model at $\delta\approx 2.8a_{B}$.

Figure 2: Ground state energy of $\text{H}_{2}$-molecule, $E$, as a function
of nuclear separation $\delta$, computed for the full wave packet model
(Elongated), the isotropic one (Isotropic) and the virtually exact result from
Ref. [80] (Exact). The ground state is modeled as a triplet state, with
$\rho=1$ in the correlation potential (16). The equilibrium displacement is
marked for the two wave-packet models (vertical lines) and the energy of two
isolated hydrogen atoms for the wave packets as well as the exact result are
shown with horizontal markers.

The energy difference between our model and the exact result is close to
constant for nuclear separation larger than the equilibrium position. In this
respect, the presented wave packet model is as good a descriptor of the
molecular ground state, $\text{H}_{2}$, as the atomic one, $2\times\text{H}$,
suggesting that for further improvements one needs a more versatile
description of the isotropic atomistic limit first. The present model has the
ability to treat each direction separately and what is limiting the agreement
is the restriction on functional form.

###  4.2 Binary collisions

The dynamical properties of the wave packet model have been tested for
electron-ion scattering against a full numerical solution of the Schrödinger
equation realised by the SOFT code [89, 90, 91]. The resulting trajectory data
for impact parameters in the range of $0.6-3.0\,a_{B}$ is shown in figure 3,
alongside the time evolution of the extent of the wave function, which is
compared with the result from both isotropic and elongated wave packets. The
centre of mass trajectories agree well between all three different sets of
simulations over long-length scales. In the full numerical solution, the
electron density can split and partially bind to the ion core, a qualitative
feature the wave packets cannot reproduce due to their limited functional
form. However, this is of minimal importance for the trajectories discussed
here. At these energies, only a minor fraction of the electron wave packet
gets bound so that the centre of mass agrees well with the mode position in
the full numerical simulation. The two wave packet models differ in the
internal degrees of freedom, where the isotropic wave packet does not have the
flexibility of the complete model. The isotropic wave packet model cannot
reproduce the internal dynamics of the SOFT computation as well as the
elongated wave packet model for the larger set of impact parameters. In the
case of smaller impact parameters, the full numerical solution has wave
packets with non-gaussian structure during the close approach between the
electron and the ion, impacting the subsequent evolution of the wave function.

Figure 3: Electron scattering of a fixed proton with different impact
parameters $b$ computed using elongated wave packets, isotropic wave packets
and SOFT. The trajectory of electron mass-centre $(a)$ and the width of the
electron wave function along $y$-direction, $\sqrt{\langle
y^{2}\rangle-\langle y\rangle^{2}}$, for two different impact parameters
$(b)$. The wave functions were initialised as isotropic gaussians with a width
$1\,a_{B}$ and velocity such that the total classical energy of each
trajectory were $\frac{3}{2}k_{B}T$ and $k_{B}T=10\,\text{eV}$.

One of the more important characteristics for determining the dynamic
properties of the system is the energy transfer in particle collisions, which
is shown for the electron-proton scattering in figure 4. The energy transfer
is calculated by treating the proton dynamically in contrast to the case shown
in figure 3. Differences between the wave packet models can in particular be
seen for intermediate impact parameters, between the classical behaviour for
large impact parameters and the symmetrical configuration when $b=0$.
Furthermore, the difference is the most pronounced for smaller wave packets
where the asymmetry on the scale of the wave packet is more pronounced. The
result suggests there might be an appreciable difference in the dynamical
properties of the two models.

Figure 4: Transfer of (classical) kinetic energy $\Delta E/E$ from the
electron to the ion in a single electron-ion scattering event as a function of
impact parameter $b$. Initial velocity set as in figure 3 with a temperature
$(a)$ $k_{B}T=5\,\text{eV}$ and $(b)$ $k_{B}T=10\,\text{eV}$. The energy
transferred differs between the isotropic and elongated models, especially for
intermediate impact parameters and smaller wave packets. The result is shown
for different free particle widths $w$ set by the confining potential in
equation (17) by varying $l_{w}$ with $A_{w}=3\,E_{h}/a_{B}^{2}$.

###  4.3 Transport properties

The previously shown demonstrations of the model were limited to the dynamics
of a few particles, however, one of the strengths of the wave packet models is
in the treatment of large collections of particles. We demonstrate this by
studying a hydrogen plasma with degeneracy parameter
$\theta=k_{B}T/E_{\text{F}}=1$ and a density $r_{s}=2$ under WDM conditions.
The system is modeled by a thousand protons and an equal number of electron
wave packets describing a spin un-polarised electron fluid (with an equal
number of spin-up and spin-down electrons). For this initial test we set
$\rho=1$ in equation (16) and the width of the wave packets are regularised by
the width confinement where $A_{w}=3\,E_{h}/a_{B}^{2}$ and $l_{w}=1\,a_{B}$.
The initial random configurations of particles were allowed to thermalise
under the influence of periodic velocity re-scaling (see appendix D) for
$75\,\text{fs}$, after which data was collected for $225\,\text{fs}$, a
procedure repeated five times for each wave packet model.

Figure 5: $(a)$ Radial distribution function for proton-proton ($p$-$p$),
electron-electron ($e$-$e$) and proton-electron ($p$-$e$) in both wave-packet
models described, elongated (solid) and isotropic (dashed). Qualitatively
similar structures is seen between the two models, where for isotropic wave
packets the electrons has a stronger interaction between themselves as well as
the ions. The typical errors – estimated from the distinct runs – are shown in
the shaded area only noticeable at small $r$. $(b)$ Electron current auto
correlation function for the wave packet model (solid), described well by an
exponential decay shown by the corresponding fit (dashed). The system
considered in both $(a)$ and $(b)$ is the one described in section 44.3.

Figure 5$(a)$ shows the static pair-correlation functions [27] for the two
wave packet models under consideration, where they primarily differ in
electronic structure. The isotropic wave packets are seen to interact more
strongly both with ions and between themselves. A part of this stronger
interaction is during a collision between a wave packet and a proton. The
isotropic wave function expands slower, as the appropriate expansion is
averaged over all directions, and the smaller wave packets have a more
localised charge and a stronger interaction.

One of the properties of interest for WDM systems is the electrical
conductivity $\sigma$. The conductivity relates to the microscopic dynamics in
an atomistic simulation via [92],

$\sigma(\omega)=\frac{1}{3k_{B}\\!T\,V}\int_{0}^{\infty}\left\langle\\!\left\langle\boldsymbol{J}(t)\cdot\boldsymbol{J}(0)\right\rangle\\!\right\rangle
e^{i\omega t}\,dt,$ (18)

where $\left\langle\\!\left\langle\cdot\right\rangle\\!\right\rangle$ is a
thermal average and $\boldsymbol{J}$ is the total current. The current is
dominated by the electron contribution, $\boldsymbol{J}_{e}$,

$\boldsymbol{J}\approx\boldsymbol{J}_{e}=-\frac{e}{m}\sum_{i=1}^{N_{e}}\langle\hat{\boldsymbol{p}}_{i}\rangle,$
(19)

where we sum over all electrons. The current-current correlation function has
roughly an exponentially decaying behaviour and therefore in the vein of Ref.
[93] an exponential form is fitted to reduce the influence of noise. Both the
numerical data and the fit is shown in figure 5$(b)$ with a high level of
agreement. An exponential current-current correlation function results in a
Drude-like conductivity,

$\sigma(\omega)=\frac{\sigma_{0}}{1+\left(\omega/\lambda\right)^{2}}+i\,\frac{\sigma_{0}\omega/\lambda}{1+\left(\omega/\lambda\right)^{2}},$
(20)

where $\sigma_{0}$ and $\lambda$ are related to the amplitude and time
constant of the decaying correlation function. Physically, $\sigma_{0}$ is the
DC conductivity while $\lambda$ corresponds to the inverse of the mean free
scattering time. The wave packet models differ in their prediction of both
constants. For elongated wave-packets, $\sigma_{0}^{\text{elo}}=(34\,300\pm
1\,900)~{}\text{\textOmega}^{-1}\text{cm}^{-1}$ and
$\lambda^{\text{elo}}=(1.48\pm 0.06)~{}\text{fs}^{-1}$, while for the
isotropic model $\sigma_{0}^{\text{iso}}=(29\,800\pm
1\,200)~{}\text{\textOmega}^{-1}\text{cm}^{-1}$ and
$\lambda^{\text{iso}}=(1.73\pm 0.06)~{}\text{fs}^{-1}$. This corresponds
roughly to a $15\%$ increase in DC conductivity and a $15\%$ decrease in the
scattering frequency as we extend the wave packet formulation. For a Drude-
like conductivity, $\sigma_{0}$ and $\lambda$ are strongly correlated, which
can be understood by the constraints set by the fluctuation-dissipation
theorem [94]. Within this wave packet formulation the electron-electron
scattering is explicitly included in the dynamic formulation, not the case in
the commonly used Kubo-Greenwod formulation of conductivity for DFT-MD,
preventing the DFT based technique to achieve the appropriate high temperature
limit [95].

##  5 Conclusion

Non-equilibrium modelling of quantum many-body systems is a formidable task
that further suffers from the large proton-to-electron mass ratio, resulting
in vastly different time scales for the evolution of electron and ion
dynamics. Any dynamical treatment must therefore resolve the electron motion,
while extending over comparatively long time scales to investigate ion
dynamics. Wave packet models address this with an ansatz for the electronic
wave function allowing the time evolution of a large number of particles to be
performed over long time scales, while retaining quantum mechanical properties
dynamically within the model.

We have extended the functional form of the wave packets used for modelling
WDM by allowing the wave packet to be elongated with arbitrary rotation. This
in turn allows for a dynamic response to gradients across the wave packet that
should better represent quantum dynamics. As a consequence of the non-
isotropic states, explicit evaluation of the interaction Hamiltonian has not
been possible, and a generalised Ewald summation has been used to
appropriately evaluate both the short- and long-range effect of the Coulomb
interaction. Furthermore, a decomposition of the short-range interaction
kernel into gaussian modes has been constructed and used for an explicit
evaluation scheme.

Crucially, warm dense matter systems are partially degenerate and the exchange
interaction contributes significantly to their evolution. A scalable
approximation for this in terms of Pauli potentials has been derived as they
fundamentally depend on the wave packets used. The interactions are
implemented along with the complete dynamical description in LAMMPS, with good
parallel support and close to a linear scaling with particle number.

The elongated wave packet model is seen to improve the description of ground
state hydrogen molecules compared to isotropic wave packets where spherical
symmetry is naturally broken. The added degrees of freedom also further
improve the dynamical description when compared to full quantum mechanical
treatments, which we demonstrate for electron-proton scattering. Furthermore,
in such collisions different energy transfers are observed between the models,
illustrating an important variation in the predicted behaviour. Finally, the
model is used to investigate a partially degenerate hydrogen plasma.
Differences between wave packet models are seen both in the electronic
structure and in the dynamics. A property of fundamental interest for this
type of system is the electrical conductivity, which we extracted from our
wave packet models. The electrical conductivity is dominated by electron
motion and the DC conductivity was seen to increase by approximately $15\%$
when extending the functional form of the wave packet.

Founding: This work was financially supported by AWE UK via Oxford Centre for
High Energy Density Science (OxCHEDS) scholarship jointly founded by Oxford
Physics Endowment for Graduates (OxPEG).

Acknowledgement:The authors would like to thank S. Azadi and T. Gawne for
valuable discussions.

##  Appendix A Gaussian decomposition

The decomposition of the Coulomb interaction kernels described in section 22.2
is solely introduced to efficiently evaluate state averages and only needed to
be computed once. The state averages of the interaction kernel $V$ can be
reduced to a volume integration weighted by some gaussian profile. Therefore,
the error in gaussian approximation was quantified by the volume average of
the square errors

$\mathcal{L}=\int_{0}^{r_{\text{max}}}\\!\\!dr\,r^{2}\Big{|}V(r)-\sum_{p}c_{p}e^{-\alpha_{p}r^{2}}\Big{|}^{2}$
(21)

which is finite without a small $r$ cutoff. In the case of the short-range
contribution, $V(r)=\text{erfc}\left(gr\right)/r$, the integrals are
convergent and $r_{\text{max}}=\infty$ was chosen. On the other hand for the
long-range correction, $V(r)=\text{erf}\left(gr\right)/r$ a finite
$r_{\text{max}}$ is needed and has been chosen to $r_{\text{max}}=15a_{B}$. As
this decomposition is only applied for terms which are exponentially
suppressed over large distances, it dose not constitute a significant
limitation.

The decomposition in terms of the coefficients $c_{p}$ and $\alpha_{p}$ is
based on the minimisation of $\mathcal{L}$. For a fixed set of $\alpha_{p}$,
the minimisation is solved by inversion of the linear system of equations
${\partial\mathcal{L}/\partial c_{p}=0}$. However, instead of using a
particular set of $\alpha_{p}$ [71] a general minimisation algorithm was used
to find the optimal $\alpha_{p}$’s, preconditioned with the optimal
$c_{p}=c_{p}(\left\\{\alpha_{p}\right\\})$. For the short-range interaction
the search was limited to ${\alpha_{p}>0}$ where in the long-range case
${r_{\text{max}}^{-2}/2<\alpha_{p}<g^{2}}$ was needed for a stable
decomposition. Figure 6 demonstrates the final decomposition resulting from
this procedure described above, typical for the ones used in the main text.
The singular nature of the short-range interaction is approximated by a series
of modes with successively larger amplitude and shorter range.

Figure 6: The gaussian mode decomposition of the short-range interaction $(a)$
and the long range-interaction $(b)$ with a Ewald parameter $g=0.5$.

##  Appendix B Ewald summation

The main contribution $\mathcal{E}_{k}$ to the Ewald summation (11) is the
interaction of the charge distribution in the simulation box
$\rho_{\text{uc}}$ and the total surrounding charge distribution
$\rho_{\text{tot}}$,

$\mathcal{E}_{k}=\int
d^{3}\boldsymbol{x}_{1}d^{3}\boldsymbol{x}_{2}\frac{\text{erf}\left(g|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}|\right)}{2|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}|}\rho_{\text{uc}}\left(\boldsymbol{x}_{1}\right)\rho_{\text{tot}}\left(\boldsymbol{x}_{2}\right),$
(22)

seen by exchanging the particle summation and state average in equation (11).
The charge distributions are based on the functional form of the wave packet
and the image charges, namely

$\rho_{\text{uc}}\left(\boldsymbol{x}\right)=\sum_{i}eZ_{i}|\bra{\boldsymbol{x}}\ket{q_{i}}|^{2}\qquad\text{and}\qquad\rho_{\text{tot}}\left(\boldsymbol{x}\right)=\sum_{\boldsymbol{n}\in\mathbb{Z}^{3}}\rho_{uc}\left(\boldsymbol{x}-L\boldsymbol{n}\right),$
(23)

the Fourier transform of which has been evaluated in equation (13). With the
use of Plancherel’s theorem and the Fourier convolution theorem, one can
arrive at the expression (12) which converges rapidly in $K$-space due to the
exponential dependence on $\boldsymbol{k}$. The $\boldsymbol{k}=0$ term vanish
due to charge neutrality $\tilde{\rho}_{\text{uc}}(0)=0$.

In any given computation only a finite number of $\boldsymbol{k}$-vector may
be used and the associated error incurred by truncating equation (12) was
described for the classical point particle case by Ref. [68]. It can be shown
that the resulting force from the Ewald summation for distributed particles
are suppressed by factors
$\exp\left(-0.5\boldsymbol{k}^{\intercal}\boldsymbol{\Sigma}_{i}\boldsymbol{k}\right)\times\exp\left(-0.5\boldsymbol{k}^{\intercal}\boldsymbol{\Sigma}_{j}\boldsymbol{k}\right)$
for each $\boldsymbol{k}$-vector and the force converge more rapidly – as the
charge distribution now is more smooth – such that the classical estimate
still can be used as an upper bound. The calculations in the manuscript were
performed such that this error estimate was one thousand of the force two ions
experience an average distance apart.

The self-energy term can be formulated similarly:

$\mathcal{E}_{s}=-\sum_{i}\left(eZ_{i}\right)^{2}\\!\int\\!d^{3}\\!\boldsymbol{x}_{1}|\bra{\boldsymbol{x}_{1}}\ket{q_{i}}|^{2}\int
d^{3}\\!\boldsymbol{x}_{2}\,\frac{\text{erf}\left(g|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}|\right)}{2|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}|}|\bra{\boldsymbol{x}_{2}}\ket{q_{i}}|^{2},$
(24)

where the integration could be performed if the kernel were gaussian. The
integrand is exponentially small if
$\boldsymbol{x}_{1}\not\approx\boldsymbol{x}_{2}$ and a gaussian decomposition
of the long-range interaction kernel is utilised, using
$\text{erf}\left(gr\right)/r=\sum_{p}\tilde{c}_{p}\exp\left(-\tilde{\alpha}_{p}r^{2}\right)$.
The self-energy contribution is

$\mathcal{E}_{s}=-\frac{1}{2}\sum_{i}\left(eZ_{i}\right)^{2}\sum_{p}\frac{\tilde{c}_{p}}{\sqrt{\det\left(\boldsymbol{I}+4\tilde{\alpha}_{p}\boldsymbol{\Sigma}_{i}\right)}},$
(25)

yielding the correct classical limit, $\boldsymbol{\Sigma}_{i}\rightarrow 0$,
as $\sum_{p}\tilde{c}_{p}=2g/\sqrt{\pi}$, if the decomposition retain the
correct value at $r=0$.

##  Appendix C Log-range Pauli interaction

The background terms, $V_{\text{bg}}$, in the Pauli potential (14b) have a
long-range behaviour and by splitting the interaction in accordance with the
Ewald summation discussed in section 22.2 it may be truncated and any residual
long-range term is then accounted. Once again the long-range interaction can
be formulated as the interaction between two distributions; one is the
particle overlap and other is the charge density from the background
particles:

$\mathcal{V}^{P/C}_{\text{long}}=\int
d^{3}\boldsymbol{x}_{1}d^{3}\boldsymbol{x}_{2}\frac{\text{erf}\left(g|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}|\right)}{2|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}|}\mu^{P/C}(\boldsymbol{x}_{1})\rho_{\text{tot}}(\boldsymbol{x}_{2})+\mathcal{V}_{s}^{P/C}$
(26)

where

$\displaystyle\mu^{P/C}(\boldsymbol{x})=\pm\sum_{i<j}\frac{e}{1\mp|\bra{q_{i}}\ket{q_{j}}|^{2}}\Big{[}2$
$\displaystyle\real\left\\{\bra{q_{j}}\ket{\boldsymbol{x}}\bra{\boldsymbol{x}}\ket{q_{j}}\bra{q_{j}}\ket{q_{j}}\right\\}$
(27)
$\displaystyle-\left(|\bra{\boldsymbol{x}}\ket{q_{i}}|^{2}+|\bra{\boldsymbol{x}}\ket{q_{j}}|^{2}\right)|\bra{q_{i}}\ket{q_{j}}|^{2}\Big{]}$

and the sum is restricted to electrons with relevant spins. A self-energy term
$\mathcal{V}_{s}^{P/C}$ is introduced when the background terms do not have
terms where $k=i,j$,

$\displaystyle\mathcal{V}_{s}^{P/C}=\pm\sum_{i<j}\int$ $\displaystyle
d^{3}\boldsymbol{x}_{1}d^{3}\boldsymbol{x}_{2}\frac{\text{erf}\left(g|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}|\right)}{2|\boldsymbol{x}_{1}-\boldsymbol{x}_{2}|}\frac{e}{1\mp|\bra{q_{i}}\ket{q_{j}}|^{2}}\left(|\bra{\boldsymbol{x}}\ket{q_{i}}|^{2}+|\bra{\boldsymbol{x}}\ket{q_{j}}|^{2}\right)$
(28)
$\displaystyle\times\left[2\real\left\\{\bra{q_{j}}\ket{\boldsymbol{x}}\bra{\boldsymbol{x}}\ket{q_{j}}\bra{q_{j}}\ket{q_{j}}\right\\}-\left(|\bra{\boldsymbol{x}}\ket{q_{i}}|^{2}+|\bra{\boldsymbol{x}}\ket{q_{j}}|^{2}\right)|\bra{q_{i}}\ket{q_{j}}|^{2}\right]$

while the above expression has been evaluated in K-space and included in the
force computation, it has only a limited impact on the overall dynamics, and
it could be omitted if computational speed is needed.

##  Appendix D Temperature measurements

The system temperature, $T$, is defined through the equipartition theorem:

$\left\langle\\!\left\langle Q_{\mu}\frac{\partial\mathcal{H}}{\partial
Q_{\nu}}\right\rangle\\!\right\rangle=\delta_{\mu\nu}k_{B}T,$ (29)

where $\delta_{\mu\nu}$ is the Kronecker delta and
$\left\langle\\!\left\langle\cdot\right\rangle\\!\right\rangle$ represent the
thermal average in the microcanonical ensemble commonly referred to as the NVE
ensemble, due to a constant particle number, $N$, volume, $V$, and energy $E$.
This corresponds to a time average in molecular dynamics in accordance with
the ergodic hypothesis [65]. We may define a temperature based on the
classical degrees of freedom $Q_{\mu}=Q_{\nu}=\boldsymbol{p}_{i}$,

$\left\langle\\!\left\langle\frac{\boldsymbol{p}_{i}^{2}}{m}+\boldsymbol{p}_{i}^{\intercal}\frac{\partial\mathcal{V}}{\partial\boldsymbol{p}_{i}}\right\rangle\\!\right\rangle=3k_{B}T,$
(30)

and the internal ones,

$\left\langle\\!\left\langle\frac{4}{m}\Tr\\{\boldsymbol{\Pi}_{i}\boldsymbol{\Sigma}_{i}\boldsymbol{\Pi}_{i}\big{\\}}+\Tr\left\\{\boldsymbol{\Pi}_{i}\left(\boldsymbol{\tau}\odot\frac{\partial\mathcal{V}}{\partial\boldsymbol{\Pi}_{i}}\right)\right\\}\right\rangle\\!\right\rangle=6k_{B}T$
(31)

where $\mathcal{V}$ is the state average of the interaction terms of the
hamiltonian. The above two definitions are combined, averaged over the
particle index $i$ and evaluated instantaneously to monitor the temperature
evolution of the system. Instantaneously, equations (30) and (31) may differ,
nonetheless converges to the same value by time averaging. Previously, Ref.
[78] has incorporated the width kinetic energy in the temperature measurement
for isotropic wave packets. When we employ isotropic wave packets we use
equation (31) accounting for the reduced degrees of freedom.

##  Appendix E Generalised friction

To achieve a continuous energy loss from friction, the equations of motion are
modified accordingly:

$\frac{d\boldsymbol{p}_{i}}{dt}=-\frac{\partial\mathcal{H}}{\partial\boldsymbol{r}_{i}}+\boldsymbol{\Gamma}_{i}^{p},\qquad\text{and}\qquad\frac{d\boldsymbol{\Pi}_{i}}{dt}=-\boldsymbol{\tau}\odot\frac{\partial\mathcal{H}}{\partial\boldsymbol{\Sigma}_{i}}+\boldsymbol{\Gamma}_{i}^{\Pi}.$
(32)

The resulting energy loss

$\frac{d}{dt}\mathcal{H}=-\sum_{i}\frac{\partial\mathcal{H}}{\partial\boldsymbol{p}_{i}^{\intercal}}\;\boldsymbol{\Gamma}_{i}^{p}+\Tr\left\\{\boldsymbol{\tau}\odot\frac{\partial\mathcal{H}}{\partial\Pi_{i}}\boldsymbol{\Gamma}_{i}^{\Pi}\right\\},$
(33)

can be guaranteed to result in a monotonic energy decrease with the choice

$\boldsymbol{\Gamma}_{i}^{p}=m\gamma_{p}\;\frac{\partial\mathcal{H}}{\partial\boldsymbol{p}_{i}}\hskip
14.22636pt\boldsymbol{\Gamma}_{\Pi}=\frac{m}{8}\gamma_{\Pi}\left(\boldsymbol{\Sigma}_{i}^{-1}\boldsymbol{\tau}\odot\frac{\partial\mathcal{H}}{\partial\Pi_{i}}+\boldsymbol{\tau}\odot\frac{\partial\mathcal{H}}{\partial\Pi_{i}}\boldsymbol{\Sigma}_{i}^{-1}\right)$
(34)

which reduces the classical case in the absence of a Pauli potential and
$\boldsymbol{\Sigma}_{i}\boldsymbol{\Pi}_{i}=\boldsymbol{\Pi}_{i}\boldsymbol{\Sigma}_{i}$.

## References

  * [1] John Nuckolls, Lowell Wood, Albert Thiessen and George Zimmerman  “Laser compression of matter to super-high densities: Thermonuclear (CTR) applications”  In _Nature_ 239.5368 Nature Publishing Group, 1972, pp. 139–142
  * [2] R Betti and OA Hurricane  “Inertial-confinement fusion with lasers”  In _Nature Physics_ 12.5 Nature Publishing Group, 2016, pp. 435–448
  * [3] H Abu-Shawareb et al.  “Lawson criterion for ignition exceeded in an inertial fusion experiment”  In _Physical Review Letters_ 129.7 APS, 2022, pp. 075001
  * [4] Raymond F Smith et al.  “Equation of state of iron under core conditions of large rocky exoplanets”  In _Nature Astronomy_ 2.6 Nature Publishing Group, 2018, pp. 452–458
  * [5] Katerina Falk  “Experimental methods for warm dense matter research”  In _High Power Laser Science and Engineering_ 6 Cambridge University Press, 2018
  * [6] Adam Frank et al.  “Exoplanets and High Energy Density Plasma Science” https://baas.aas.org/pub/2020n3i036  In _Bulletin of the AAS_ 51.3, 2019  URL: https://baas.aas.org/pub/2020n3i036
  * [7] Anna Levy et al.  “The creation of large-volume, gradient-free warm dense matter with an x-ray free-electron laser”  In _Physics of Plasmas_ 22.3 AIP Publishing LLC, 2015, pp. 030703
  * [8] LB Fletcher et al.  “Ultrabright X-ray laser scattering for dynamic warm dense matter physics”  In _Nature photonics_ 9.4 Nature Publishing Group, 2015, pp. 274–279
  * [9] Manfred Schlanges, Michael Bonitz and A Tschttschjan  “Plasma Phase Transition in Fluid Hydrogen-Helium Mixtures”  In _Contributions to Plasma Physics_ 35.2 Wiley Online Library, 1995, pp. 109–125
  * [10] V Bezkrovniy et al.  “Monte Carlo results for the hydrogen Hugoniot”  In _Physical Review E_ 70.5 APS, 2004, pp. 057401
  * [11] Jan Vorberger, I Tamblyn, B Militzer and SA Bonev  “Hydrogen-helium mixtures in the interiors of giant planets”  In _Physical Review B_ 75.2 APS, 2007, pp. 024206
  * [12] Winfried Lorenzen, Andreas Becker and Ronald Redmer  “Progress in warm dense matter and planetary physics”  In _Frontiers and Challenges in Warm Dense Matter_ Springer, 2014, pp. 203–234
  * [13] Burkhard Militzer et al.  “First-principles equation of state database for warm dense matter computation”  In _Physical Review E_ 103.1 APS, 2021, pp. 013203
  * [14] Didier Saumon, William B Hubbard, Gilles Chabrier and Hugh M Van Horn  “The role of the molecular-metallic transition of hydrogen in the evolution of Jupiter, Saturn, and brown dwarfs”  In _The Astrophysical Journal_ 391, 1992, pp. 827–831
  * [15] Gilles Chabrier  “Quantum effects in dense Coulumbic matter-Application to the cooling of white dwarfs”  In _The Astrophysical Journal_ 414, 1993, pp. 695–700
  * [16] Gilles Chabrier, Pierre Brassard, Gilles Fontaine and Didier Saumon  “Cooling sequences and color-magnitude diagrams for cool white dwarfs with hydrogen atmospheres”  In _The Astrophysical Journal_ 543.1 IOP Publishing, 2000, pp. 216
  * [17] Paweł Haensel, Aleksander Yu Potekhin and Dmitry G Yakovlev  “Neutron stars 1: Equation of state and structure”  Springer Science & Business Media, 2007
  * [18] J Daligault and S Gupta  “Electron–ion scattering in dense multi-component plasmas: Application to the outer crust of an accreting neutron star”  In _The Astrophysical Journal_ 703.1 IOP Publishing, 2009, pp. 994
  * [19] SX Hu et al.  “A review on ab initio studies of static, transport, and optical properties of polystyrene under extreme conditions for inertial confinement fusion applications”  In _Physics of Plasmas_ 25.5 AIP Publishing LLC, 2018, pp. 056306
  * [20] M Bonitz et al.  “Ab initio simulation of warm dense matter”  In _Physics of Plasmas_ 27.4 AIP Publishing LLC, 2020, pp. 042710
  * [21] JA Gaffney et al.  “A review of equation-of-state models for inertial confinement fusion materials”  In _High Energy Density Physics_ 28 Elsevier, 2018, pp. 7–24
  * [22] Ryan A Davis, WA Angermeier, RKT Hermsmeier and Thomas G White  “Ion modes in dense ionized plasmas through nonadiabatic molecular dynamics”  In _Physical Review Research_ 2.4 APS, 2020, pp. 043139
  * [23] Yunpeng Yao et al.  “Reduced ionic diffusion by the dynamic electron–ion collisions in warm dense hydrogen”  In _Physics of Plasmas_ 28.1 AIP Publishing LLC, 2021, pp. 012704
  * [24] PE Grabowski et al.  “Review of the first charged-particle transport coefficient comparison workshop”  In _High Energy Density Physics_ 37 Elsevier, 2020, pp. 100905
  * [25] Sean M Wahl et al.  “Comparing Jupiter interior structure models to Juno gravity measurements and the role of a dilute core”  In _Geophysical Research Letters_ 44.10 Wiley Online Library, 2017, pp. 4649–4659
  * [26] SX Hu et al.  “Impact of first-principles properties of deuterium–tritium on inertial confinement fusion target designs”  In _Physics of Plasmas_ 22.5 AIP Publishing LLC, 2015, pp. 056304
  * [27] Setsuo Ichimaru  “Strongly coupled plasmas: high-density classical plasmas and degenerate electron liquids”  In _Reviews of Modern Physics_ 54.4 APS, 1982, pp. 1017
  * [28] James P. Mithen, Jérôme Daligault and Gianluca Gregori  “Extent of validity of the hydrodynamic description of ions in dense plasmas”  In _Phys. Rev. E_ 83 American Physical Society, 2011, pp. 015401  DOI: 10.1103/PhysRevE.83.015401
  * [29] Hanno Kählert  “Thermodynamic and transport coefficients from the dynamic structure factor of Yukawa liquids”  In _Physical Review Research_ 2.3 APS, 2020, pp. 033287
  * [30] H Minoo, MM Gombert and C Deutsch  “Temperature-dependent Coulomb interactions in hydrogenic systems”  In _Physical Review A_ 23.2 APS, 1981, pp. 924
  * [31] JP Hansen and IR McDonald  “Thermal relaxation in a strongly coupled two-temperature plasma”  In _Physics Letters A_ 97.1-2 Elsevier, 1983, pp. 42–44
  * [32] JN Glosli et al.  “Molecular dynamics simulations of temperature equilibration in dense hydrogen”  In _Physical Review E_ 78.2 APS, 2008, pp. 025401
  * [33] Brett Larder et al.  “Fast nonadiabatic dynamics of many-body quantum systems”  In _Science advances_ 5.11 American Association for the Advancement of Science, 2019, pp. eaaw1634
  * [34] F Lambert, J Clérouin, S Mazevet and D Gilles  “Properties of Hot Dense Plasmas by Orbital-Free Molecular Dynamics”  In _Contributions to Plasma Physics_ 47.4-5 Wiley Online Library, 2007, pp. 272–280
  * [35] TG White et al.  “Orbital-free density-functional theory simulations of the dynamic structure factor of warm dense aluminum”  In _Physical review letters_ 111.17 APS, 2013, pp. 175002
  * [36] L Collins et al.  “Quantum molecular dynamics simulations of hot, dense hydrogen”  In _Physical Review E_ 52.6 APS, 1995, pp. 6202
  * [37] Hannes R Rüter and Ronald Redmer  “Ab initio simulations for the ion-ion structure factor of warm dense aluminum”  In _Physical review letters_ 112.14 APS, 2014, pp. 145007
  * [38] Zh A Moldabekov, M Bonitz and TS Ramazanov  “Theoretical foundations of quantum hydrodynamics for plasmas”  In _Physics of Plasmas_ 25.3 AIP Publishing LLC, 2018, pp. 031903
  * [39] David Michta  “Quantum Hydrodynamics: Theory and Computation with Applications to Charged Particle Stopping in Warm Dense Matter”, 2020
  * [40] Zhandos Moldabekov et al.  “Towards a quantum fluid theory of correlated many-fermion systems from first principles”  In _SciPost Physics_ 12.2, 2022, pp. 062
  * [41] Andrew David Baczewski et al.  “X-ray Thomson scattering in warm dense matter without the Chihara decomposition”  In _Physical review letters_ 116.11 APS, 2016, pp. 115004
  * [42] Rudolph J Magyar, L Shulenburger and AD Baczewski  “Stopping of Deuterium in Warm Dense Deuterium from Ehrenfest Time-Dependent Density Functional Theory”  In _Contributions to Plasma Physics_ 56.5 Wiley Online Library, 2016, pp. 459–466
  * [43] S Frydrych et al.  “Demonstration of X-ray Thomson scattering as diagnostics for miscibility in warm dense matter”  In _Nature communications_ 11.1 Nature Publishing Group, 2020, pp. 1–7
  * [44] Burkhard Militzer and David M Ceperley  “Path integral Monte Carlo calculation of the deuterium Hugoniot”  In _Physical review letters_ 85.9 APS, 2000, pp. 1890
  * [45] SX Hu, B Militzer, VN Goncharov and S Skupsky  “First-principles equation-of-state table of deuterium for inertial confinement fusion applications”  In _Physical Review B_ 84.22 APS, 2011, pp. 224109
  * [46] DA Chapman and DO Gericke  “Analysis of Thomson scattering from nonequilibrium plasmas”  In _Physical review letters_ 107.16 APS, 2011, pp. 165004
  * [47] U Zastrau et al.  “Resolving ultrafast heating of dense cryogenic hydrogen”  In _Physical review letters_ 112.10 APS, 2014, pp. 105002
  * [48] TG White et al.  “Electron-ion equilibration in ultrafast heated graphite”  In _Physical review letters_ 112.14 APS, 2014, pp. 145005
  * [49] Kai-Uwe Plagemann et al.  “Ab initio calculation of the ion feature in x-ray Thomson scattering”  In _Physical Review E_ 92.1 APS, 2015, pp. 013103
  * [50] Jean Clérouin et al.  “Evidence for out-of-equilibrium states in warm dense matter probed by x-ray Thomson scattering”  In _Physical Review E_ 91.1 APS, 2015, pp. 011101
  * [51] Paul Mabey et al.  “A strong diffusive ion mode in dense ionized matter predicted by Langevin dynamics”  In _Nature communications_ 8.1 Nature Publishing Group, 2017, pp. 1–6
  * [52] Hans Feldmeier and Jürgen Schnack  “Molecular dynamics for fermions”  In _Reviews of Modern Physics_ 72.3 APS, 2000, pp. 655
  * [53] Paul E. Grabowski  “A Review of Wave Packet Molecular Dynamics”  In _Frontiers and Challenges in Warm Dense Matter_ Cham: Springer International Publishing, 2014, pp. 265–282
  * [54] Steve Plimpton  “Fast parallel algorithms for short-range molecular dynamics”  In _Journal of computational physics_ 117.1 Elsevier, 1995, pp. 1–19
  * [55] Eric J Heller  “Time-dependent approach to semiclassical dynamics”  In _The Journal of Chemical Physics_ 62.4 American Institute of Physics, 1975, pp. 1544–1555
  * [56] N Corbin and K Singer  “Semiclassical molecular dynamics of wave packets”  In _Molecular Physics_ 46.3 Taylor & Francis, 1982, pp. 671–677
  * [57] J Broeckhove, L Lathouwers and P Van Leuven  “Time-dependent variational principles and conservation laws in wavepacket dynamics”  In _Journal of Physics A: Mathematical and General_ 22.20 IOP Publishing, 1989, pp. 4395
  * [58] Michael Knaup, PG Reinhard, C Toepffer and G Zwicknagel  “Wave packet molecular dynamics simulations of warm dense hydrogen”  In _Journal of Physics A: Mathematical and General_ 36.22 IOP Publishing, 2003, pp. 6165
  * [59] B Jakob, P-G Reinhard, C Toepffer and G Zwicknagel  “Wave packet simulation of dense hydrogen”  In _Physical Review E_ 76.3 APS, 2007, pp. 036406
  * [60] B Jakob, PG Reinhard, C Toepffer and G Zwicknagel  “Wave packet simulations for the insulator–metal transition in dense hydrogen”  In _Journal of Physics A: Mathematical and Theoretical_ 42.21 IOP Publishing, 2009, pp. 214055
  * [61] D Klakow, Ch Toepffer and P-G Reinhard  “Hydrogen under extreme conditions”  In _Physics Letters A_ 192.1 Elsevier, 1994, pp. 55–59
  * [62] D Klakow, Ch Toepffer and P-G Reinhard  “Semiclassical molecular dynamics for strongly coupled Coulomb systems”  In _The Journal of chemical physics_ 101.12 American Institute of Physics, 1994, pp. 10766–10774
  * [63] Junichi Ono and Koji Ando  “Semiquantal molecular dynamics simulations of hydrogen-bond dynamics in liquid water using multi-dimensional Gaussian wave packets”  In _The Journal of chemical physics_ 137.17 American Institute of Physics, 2012, pp. 174503
  * [64] Robert E Wyatt  “Quantum dynamics with trajectories: introduction to quantum hydrodynamics”  Springer Science & Business Media, 2006
  * [65] Dennis C Rapaport  “The art of molecular dynamics simulation”  Cambridge university press, 2004
  * [66] Simon W Leeuw, John William Perram and Edgar Roderick Smith  “Simulation of electrostatic systems in periodic boundary conditions. I. Lattice sums and dielectric constants”  In _Proceedings of the Royal Society of London. A. Mathematical and Physical Sciences_ 373.1752 The Royal Society London, 1980, pp. 27–56
  * [67] Markus Deserno and Christian Holm  “How to mesh up Ewald sums. I. A theoretical and numerical comparison of various particle mesh routines”  In _The Journal of chemical physics_ 109.18 American Institute of Physics, 1998, pp. 7678–7693
  * [68] Jiri Kolafa and John W Perram  “Cutoff errors in the Ewald summation formulae for point charge systems”  In _Molecular Simulation_ 9.5 Taylor & Francis, 1992, pp. 351–368
  * [69] Michael Knaup, P-G Reinhard and Ch Toepffer  “Wave packet molecular dynamics simulations of hydrogen near the transition to a metallic fluid”  In _Contributions to Plasma Physics_ 39.1-2 Wiley Online Library, 1999, pp. 57–60
  * [70] Pavel A Frantsuzov and Vladimir A Mandelshtam  “Quantum statistical mechanics with Gaussians: Equilibrium properties of van der Waals clusters”  In _The Journal of chemical physics_ 121.19 American Institute of Physics, 2004, pp. 9247–9256
  * [71] Yasuteru Shigeta et al.  “Quantal cumulant mechanics and dynamics for multidimensional quantum many-body clusters”  In _International Journal of Quantum Chemistry_ 113.3 Wiley Online Library, 2013, pp. 348–355
  * [72] Victoria Buch  “Exploration of multidimensional variational Gaussian wave packets as a simulation tool”  In _The Journal of chemical physics_ 117.10 American Institute of Physics, 2002, pp. 4738–4750
  * [73] Julius T Su and William A Goddard III  “Excited electron dynamics modeling of warm dense matter”  In _Physical review letters_ 99.18 APS, 2007, pp. 185003
  * [74] Julius T Su and William A Goddard III  “The dynamics of highly excited electronic systems: Applications of the electron force field”  In _The Journal of chemical physics_ 131.24 American Institute of Physics, 2009, pp. 244501
  * [75] Andres Jaramillo-Botero, Julius Su, An Qi and William A Goddard III  “Large-scale, long-term nonadiabatic electron molecular dynamics for describing material properties and phenomena in extreme environments”  In _Journal of computational chemistry_ 32.3 Wiley Online Library, 2011, pp. 497–512
  * [76] Hyungjun Kim, Julius T Su and William A Goddard  “High-temperature high-pressure phases of lithium from electron force field (eFF) quantum electron dynamics simulations”  In _Proceedings of the National Academy of Sciences_ 108.37 National Acad Sciences, 2011, pp. 15101–15105
  * [77] Hai Xiao, Andrés Jaramillo-Botero, Patrick L Theofanis and William A Goddard III  “Non-adiabatic dynamics modeling framework for materials in extreme conditions”  In _Mechanics of Materials_ 90 Elsevier, 2015, pp. 243–252
  * [78] Qian Ma et al.  “Extremely low electron-ion temperature relaxation rates in warm dense hydrogen: Interplay between quantum electrons and coupled ions”  In _Physical review letters_ 122.1 APS, 2019, pp. 015001
  * [79] Yun Liu et al.  “Molecular dynamics investigation of the stopping power of warm dense hydrogen for electrons”  In _Physical Review E_ 103.6 APS, 2021, pp. 063215
  * [80] William A Angermeier and Thomas G White  “An investigation into the approximations used in wave packet molecular dynamics for the study of warm dense matter”  In _Plasma_ 4.2 Multidisciplinary Digital Publishing Institute, 2021, pp. 294–308
  * [81] Michael Knaup, P-G Reinhard and Christian Toepffer  “Wave Packet Molecular Dynamics Simulations of Deuterium in the Region of Laser Shock-Wave Experiments”  In _Contributions to plasma physics_ 41.2-3 Wiley Online Library, 2001, pp. 159–162
  * [82] W Ebeling et al.  “The method of effective potentials in the quantum-statistical theory of plasmas”  In _Journal of Physics A: Mathematical and General_ 39.17 IOP Publishing, 2006, pp. 4309
  * [83] IV Morozov and IA Valuev  “Localization constraints in Gaussian wave packet molecular dynamics of nonideal plasmas”  In _Journal of Physics A: Mathematical and Theoretical_ 42.21 IOP Publishing, 2009, pp. 214044
  * [84] Paul E Grabowski et al.  “Wave packet spreading and localization in electron-nuclear scattering”  In _Physical Review E_ 87.6 APS, 2013, pp. 063104
  * [85] Ya S Lavrinenko, IV Morozov and IA Valuev  “Reflecting Boundary Conditions for Classical and Quantum Molecular Dynamics Simulations of Nonideal Plasmas”  In _Contributions to Plasma Physics_ 56.5 Wiley Online Library, 2016, pp. 448–458
  * [86] W Ebeling and B Militzer  “Quantum molecular dynamics of partially ionized plasmas”  In _Physics Letters A_ 226.5 Elsevier, 1997, pp. 298–304
  * [87] Benedict Leimkuhler and Sebastian Reich  “Simulating hamiltonian dynamics”  Cambridge university press, 2004
  * [88] Edgar Gabriel et al.  “Open MPI: Goals, Concept, and Design of a Next Generation MPI Implementation”  In _Proceedings, 11th European PVM/MPI Users’ Group Meeting_ , 2004, pp. 97–104
  * [89] MD Feit, JA Fleck Jr and A Steiger  “Solution of the Schrödinger equation by a spectral method”  In _Journal of Computational Physics_ 47.3 Elsevier, 1982, pp. 412–433
  * [90] Andreas Markmann, Frank Graziani and Victor S Batista  “Kepler Predictor–Corrector Algorithm: Scattering Dynamics with One-Over-R Singular Potentials”  In _Journal of Chemical Theory and Computation_ 8.1 ACS Publications, 2012, pp. 24–35
  * [91] Samuel M Greene and Victor S Batista  “Tensor-train split-operator fourier transform (TT-SOFT) method: Multidimensional nonadiabatic quantum dynamics”  In _Journal of chemical theory and computation_ 13.9 ACS Publications, 2017, pp. 4034–4042
  * [92] Priya Vashishta, Rajiv K Kalia, Aiichiro Nakano and José Pedro Rino  “Interaction potentials for alumina and molecular dynamics simulations of amorphous and liquid alumina”  In _Journal of Applied Physics_ 103.8 American Institute of Physics, 2008, pp. 083504
  * [93] JP Mithen, J Daligault and G Gregori  “Molecular Dynamics Simulations for the Shear Viscosity of the One-Component Plasma”  In _Contributions to Plasma Physics_ 52.1 Wiley Online Library, 2012, pp. 58–61
  * [94] Rep Kubo  “The fluctuation-dissipation theorem”  In _Reports on progress in physics_ 29.1 IOP Publishing, 1966, pp. 255
  * [95] Martin French et al.  “Electronic transport coefficients from density functional theory across the plasma plane”  In _Physical Review E_ 105.6 APS, 2022, pp. 065204