# Theory of Moment Propagation for Quantum Dynamics in Single-Particle
Description

Nicholas Boyer  Christopher Shepard  Ruiyi Zhou  Jianhang Xu  Department of
Chemistry, University of North Carolina at Chapel Hill, Chapel Hill, North
Carolina 27599, USA  Yosuke Kanai<EMAIL_ADDRESS>Department of Chemistry,
University of North Carolina at Chapel Hill, Chapel Hill, North Carolina
27599, USA  Department of Physics and Astronomy, University of North Carolina
at Chapel Hill, Chapel Hill, North Carolina 27599, USA

###### Abstract

We present a novel theoretical formulation for performing quantum dynamics in
terms of moments within the single-particle description. By expressing the
quantum dynamics in terms of increasing orders of moments, instead of single-
particle wave functions as generally done in time-dependent density functional
theory, we describe an approach for reducing the high computational cost of
simulating the quantum dynamics. The equation of motion is given for the
moments by deriving analytical expressions for the first-order and second-
order time derivatives of the moments, and a numerical scheme is developed for
performing quantum dynamics by expanding the moments in the Taylor series as
done in classical molecular dynamics simulation. We propose a few numerical
approaches using this theoretical formalism on a simple one-dimensional model
system, for which an analytically exact solution can be derived. Application
of the approaches to an anharmonic system is also discussed to illustrate
their generality. We also discuss the use of an artificial neural network
model to circumvent the numerical evaluation of the second-order time
derivatives of the moments, as analogously done in the context of classical
molecular dynamics simulations.

##  I Introduction

Over the past few decades, explicit time-dependent simulations have become
increasingly popular for investigating non-equilibrium quantum dynamics of
electrons in chemistry, materials science, and condensed matter physics [1, 2,
3]. In particular, the real-time propagation approach[4, 5, 6] to time-
dependent density functional theory[7] (RT-TDDFT) has significantly advanced
over the last decade [8, 9] for investigating various physical phenomena
including optical excitation[10, 11, 12], electronic stopping[13, 14, 15],
charge transfer dynamics[16, 17, 18], laser-induced water splitting[19, 20,
21], electronic circular dichroism [22, 23, 24, 25], and many other electronic
excitation phenomena[1, 26, 27]. Within the Kohn-Sham (KS) ansatz, the quantum
dynamics of electrons are simulated by propagating a set of single-particle
orbitals using time-dependent KS (TD-KS) equations [8]. The quantum dynamics
is invariant with respect to unitary transformation of the single-particle
orbitals, and this gauge freedom has been exploited for achieving numerical
efficiency of computationally expensive simulations[28] as well as for gaining
physical insights from simulations of complex heterogeneous matter [13]. In
particular, the gauge transformation into spatially-localized maximally
localized Wannier functions (MLWFs)[29, 30, 31] has proved quite useful in RT-
TDDFT simulations not only for efficiently implementing advanced hybrid
exchange-correlation approximations [8] but also in studying topological phase
in the context of the Floquet engineering[32, 33].

In parallel, the emerging field of machine-learning has spurred great interest
in employing artificial neural networks (ANNs) for devising efficient schemes
for performing quantum dynamics simulations[34, 35, 36, 37, 38, 39, 40, 41],
as ANN models have been effectively utilized in the context of classical
molecular dynamics (MD) simulations [42, 43, 44, 45]. Recent developments have
focused on simulating electron dynamics by separating the real and imaginary
parts of the wave function [46, 47, 48, 49]. Secor, et al., for example,
introduced a numerical method to propagate wavepackets using ANNs, and proof-
of-principle demonstrations were presented in their recent work[49].

These recent developments in the fields motivated us to examine how quantum
dynamics based on single-particle wave functions like those of the TD-KS
equation can be reformulated as the dynamics of the moments with increasing
orders. Some earlier works have considered Heisenberg equation of motion in
terms of the position and momentum operators in a related context [50, 51,
52]. In this work, the equation of motion for moment propagation is derived,
and numerical demonstrations are provided for simple proof-of-principle
systems, where analytical solutions can be derived. We also use an anharmonic
potential example to show the generality of the numerical methods based the
new theoretical formulation. We further propose how machine-learning
approaches like ANN can be conveniently used to perform efficient quantum
dynamics simulations with the newly developed moment propagation theory,
analogously to the classical molecular dynamics simulation.

##  II Theoretical Formulation

In this section, we develop the theory of moment propagation. We express
single-particle wave functions in terms of moments of increasing orders, and
we derive the equation of motion such that Taylor series expansion can be used
for propagating the moments in time. This formalism allows us to devise a
numerically convenient method for performing quantum dynamics simulation using
the single-particle description like the TD-KS equation in RT-TDDFT.

###  II.1 Equation of Motion for Moments

We consider the situation in which quantum dynamics of a many-particle system
is described by the propagation of single-particle wave functions via
Schrodinger-like equation, as done through formulating the TD-KS equation. In
order to keep the discussion centered on the theoretical formulation, lengthy
mathematical derivations are given in Appendix. Instead of expressing a
single-particle wave function as $\psi(\mathbf{r},t)$, we explicitly write it
as $\psi(x,y,z,t)$, since the orders of moments are generally not the same in
the three Cartesian coordinates. The central quantity for our theoretical
formulation is the moments of single-particle probability density, given by

$\braket{x^{a}y^{b}z^{c}}(t)=\int\int\int{x^{a}y^{b}z^{c}n(x,y,z,t)dxdydz},$
(1)

where $a,b,c$ are non-negative integers used to denote the a-th, b-th, c-th
moment in $x,y,z$ directions of the Cartesian coordinate space and
$n(x,y,z,t)$ is the particle (probability) density, which is the square
modulus of the single-particle wave function (i.e.
$n(x,y,z,t)=|\psi(x,y,z,t)|^{2}$). Moments are widely used as a statistical
quantity to characterize features of probability functions. In order to
propagate moments, the time derivatives of the moments need to be derived from
time-dependent Schrodinger equation (TD-SE) or alike,

$i\frac{\partial\psi(x,y,z,t)}{\partial
t}=-\frac{1}{2}\nabla^{2}\psi(x,y,z,t)+V(x,y,z,t)\psi(x,y,z,t)$ (2)

where $\psi$ is the time-dependent wave function and $V$ specifies time-
dependent potential. The first-order time derivative of the moments is given
by (see Appendix B)

$\displaystyle\frac{d\braket{x^{a}y^{b}z^{c}}(t)}{dt}=$
$\displaystyle\int\int\int{x^{a}y^{b}z^{c}\dot{n}(x,y,z,t)dxdydz}$ (3)
$\displaystyle=$
$\displaystyle-\frac{i}{2}\int\left[\nabla^{2}\left(x^{a}y^{b}z^{c}\right)n+2\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\nabla\psi\psi^{*}\right]d^{3}r$
$\displaystyle=$
$\displaystyle-i\int\int\int(ax^{a-1}y^{b}z^{c}\partial_{x}\psi\psi^{*}+by^{b-1}x^{a}z^{c}\partial_{y}\psi\psi^{*}+cz^{c-1}x^{a}y^{b}\partial_{z}\psi\psi^{*}$
$\displaystyle+$
$\displaystyle\frac{a(a-1)}{2}x^{a-2}y^{b}z^{c}n(x,y,z,t)+\frac{b(b-1)}{2}y^{b-2}x^{a}z^{c}n(x,y,z,t)$
$\displaystyle+$
$\displaystyle\frac{c(c-1)}{2}z^{c-2}x^{a}y^{b}n(x,y,z,t))dxdydz.$

Here, we used the simplified notation such that
$\partial_{g}\equiv\frac{\partial}{\partial g}$ for brevity. Similarly, the
second-order time derivative of the moments (see Appendix C) is

$\displaystyle\frac{d^{2}\braket{x^{a}y^{b}z^{c}}(t)}{dt^{2}}=$
$\displaystyle\int\operatorname{Re}\left[-\nabla(x^{a}y^{b}z^{c})\cdot\nabla
Vn+\frac{1}{4}\nabla^{4}\left(x^{a}y^{b}z^{c}\right)n-(\nabla\otimes\nabla\left(x^{a}y^{b}z^{c}\right)\cdot(\nabla\otimes\nabla\psi))\psi^{*}\right]d^{3}r$
(4) $\displaystyle=$
$\displaystyle\operatorname{Re}(-\int\int\int(ax^{a-1}y^{b}z^{c}\partial_{x}V(x,y,z,t)n(x,y,z,t)$
$\displaystyle+$ $\displaystyle
bx^{a}y^{b-1}z^{c}\partial_{y}V(x,y,z,t)n(x,y,z,t)+cx^{a}y^{b}z^{c-1}\partial_{z}V(x,y,z,t)n(x,y,z,t))dxdydz$
$\displaystyle-$
$\displaystyle\int\int\int(a(a-1)x^{a-2}y^{b}z^{c}\partial_{x}^{2}\psi\psi^{*}+2bax^{a-1}y^{b-1}z^{c}\partial_{x}\partial_{y}\psi\psi^{*}$
$\displaystyle+$ $\displaystyle
2cax^{a-1}y^{b}z^{c-1}\partial_{x}\partial_{z}\psi\psi^{*}+b(b-1)y^{b-2}x^{a}z^{c}\partial_{y}^{2}\psi\psi^{*}$
$\displaystyle+$ $\displaystyle
2cby^{b-1}x^{a}z^{c-1}\partial_{y}\partial_{z}\psi\psi^{*}+c(c-1)z^{c-2}x^{a}y^{b}\partial_{z}^{2}\psi\psi^{*})dxdydz$
$\displaystyle+$ $\displaystyle
a(a-1)(\frac{(a-2)(a-3)}{4}\braket{x^{a-4}y^{b}z^{c}}+\frac{b(b-1)}{2}\braket{y^{b-2}x^{a-2}z^{c}}$
$\displaystyle+$
$\displaystyle\frac{c(c-1)}{2}\braket{z^{c-2}x^{a-2}y^{b}})+b(b-1)(\frac{(b-2)(b-3)}{4}\braket{y^{b-4}x^{a}z^{c}}$
$\displaystyle+$
$\displaystyle\frac{c(c-1)}{2}\braket{z^{c-2}x^{a}y^{b-2}})+\frac{c(c-1)(c-2)(c-3)}{4}\braket{z^{c-4}x^{a}y^{b}})$

where the shorthand notation
$\partial_{g}\partial_{r}\equiv\frac{\partial^{2}}{\partial g\partial r}$ is
used for brevity. Higher order derivatives can be obtained in principle, but
their utility in practice is limited because of the mathematical complexity of
the resulting expressions. In order to remove explicit reference to the
single-particle wave function in Eqs. 3 and 4, we still need to express the
wave function and its derivatives in terms of the moments.

###  II.2 Single-particle Wave Function and Moments

The expressions for the first-order and second-order time derivatives of the
moments (Eqs. 3 and 4) contain explicit dependence on the wave function. We
show here that the explicit dependence can be removed by expressing the wave
function in the polar form

$\psi(x,y,z,t)=\sqrt{n(x,y,z,t)}e^{i\theta(x,y,z,t)}$ (5)

where $n(x,y,z,t)$ is the particle density and $\theta(x,y,z,t)$ represents
the phase.

####  II.2.1 Particle density using cumulants via Edgeworth series

Instead of using the moments directly, it is convenient to transform them into
the cumulant representation so that the Edgeworth series [53] can be employed
for easily evaluating properties such as the particle density. While the
following formalism is applicable in three-dimensions, we focus here on one-
dimensional case, $x$, for clarity. The cumulants are given in terms of the
moments as [54]

$\kappa_{a}=\braket{x^{a}}-\sum_{i=1}^{a-1}{a-1\choose
i-1}\kappa_{i}\braket{x^{a-i}}$ (6)

where $\kappa_{a}$ is the a-th cumulant and $\braket{x^{a}}$ is the a-th
moment. We express the particle density using cumulants via the Edgeworth
series by expanding in the basis of Cartesian Gaussian functions[55, 56]. We
do so by first generating the n-th order probability density function

$pdf_{n}(x)=z(x)\Big{(}1+\sum^{n-2}_{j=1}(\sum^{\\{P_{j}\\}}_{\\{P_{j_{m}}\\}}\big{(}\prod_{i=1}^{\infty}\frac{1}{k_{i}!}(\frac{\lambda_{i+2}}{(i+2)!})^{k_{i}}\big{)}He_{s}(x))\Big{)}$
(7)

where $\\{P_{j}\\}$ is the set of all positive integer partitions of $j$.
Integer partitions of an integer $j$ are all possible ways to add positive
integers to $j$. For example, if $j=3$, then the integer partitions include
$1+1+1=3$ and $1+2=3$. For each $P_{j}$, these combinations are denoted by
$P_{j_{m}}$ where m is the index of the partition in the set. For each $m$,
Eq. 7 contains several terms. First, $\lambda_{n}$ is given by

$\lambda_{n}=\frac{\kappa_{n}}{(\sqrt{\kappa_{2}})^{n}}$ (8)

where $\kappa_{n}$ is the n-th order cumulant. Secondly, integer $k_{i}$ in
Eq. 7 is the count of $i$ within each integer partition such that
$j=\sum_{i}ik_{i}$. For instance, for the case of the integer partition
$P_{j=3_{m=1}}$:$1+1+1$, we have $k_{1}=3$, $k_{\geq 2}=0$. For the integer
partition $P_{j=3_{m=2}}$:$1+2$, $k_{1}=1$, $k_{2}=1$, and $k_{\geq 3}=0$.
Based on the integer value $k_{i}$, we obtain another integer value $s$ as

$s=j+2\sum_{i=1}^{\infty}k_{i}.\\\ $ (9)

For example, for the integer partition, $P_{j=3_{m=1}}$:$1+1+1(i=1)$, we have
$s=9$. The Hermite polynomial $He_{s}(x)$ in Eq. 7 is then given by

$He_{s}(x)=s!\sum^{\lfloor\frac{s}{2}\rfloor}_{j=0}\frac{(-1)^{j}x^{n-2j}}{2^{j}(n-2j)!j!}$
(10)

where $\lfloor s/2\rfloor$ denotes the floor function of $s/2$[57]. We use
Gaussian functions, $z(x)=\frac{1}{\sqrt{2\pi}}\exp{(-\frac{x^{2}}{2})}$, as
the basis in Eq. 7. The Edgeworth expansion then gives the particle density as

$n(x)=\frac{pdf(\frac{x-\kappa_{1}}{\sqrt{\kappa_{2}}})}{\sqrt{\kappa_{2}}}$
(11)

where the probability density function $pdf(x)$ is given by a n-th order
probability density function $pdf_{n}(x)$ using Eq. 7. This procedure
effectively expresses the particle density, $n(x)$, from the moments
$\braket{x^{a}}$.

####  II.2.2 Phase from the moments

In addition to the particle density $n(x,y,z,t)$ as discussed above, the phase
$\theta(x,y,z,t)$ (which is uniquely determined up to a constant) in Eq. 5 or
rather its spatial derivatives need to be related to the moments[58]. The
derivative of the phase with respect to a particular coordinate $u$ can be
written (see Appendix D) as

$\frac{\partial\theta(x,y,z,t)}{\partial
u}=-n(x,y,z,t)^{-1}\int_{-\infty}^{u}\alpha_{u}(\mathbf{r^{\prime}},t)du^{\prime}$
(12)

where $\alpha_{u}\equiv\frac{i}{2}(\frac{\partial^{2}\psi}{\partial
u^{2}}\psi^{*}-\frac{\partial^{2}\psi^{*}}{\partial u^{2}}\psi)$ is used.
Here, the integration variable $u^{\prime}$ corresponds to the same coordinate
as $u$, and we used $\textbf{r}^{\prime}$ to indicates that the function
depends on all three coordinates $\\{x,y,z\\}$. The phase can be obtained as

$\theta(x,y,z,t)=\phi+\int_{0}^{x}\frac{\partial\theta(u,0,0,t)}{\partial
x}du+\int_{0}^{y}\frac{\partial\theta(x,u,0,t)}{\partial
y}du+\int_{0}^{z}\frac{\partial\theta(x,y,u,t)}{\partial z}du$ (13)

where $\phi$ is an arbitrary constant. It is also useful to note that the time
derivative of the particle density can be expressed in terms of these
functions as $\dot{n}=\alpha_{x}+\alpha_{y}+\alpha_{z}$ (see Appendix E).
Instead of the wave function, $\alpha_{u}$ can be written alternatively using
the moment time derivatives (see Appendix E) as

$\alpha_{u}(x,y,z,t)=\sum_{a,b,c}c_{u}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}\frac{\partial
n}{\partial\braket{x^{a}y^{b}z^{c}}}$ (14)

where $c_{u}^{abc}(t)$ are the coefficients for the moment
$\braket{x^{a}y^{b}z^{c}}(t)$, and these coefficients satisfy
$c_{x}^{abc}(t)+c_{y}^{abc}(t)+c_{z}^{abc}(t)=1$ by construction. The product
of $c_{u}^{abc}(t)$ and the time derivative of the moments in Eq. 14 are
propagated together. For example, one has (see Appendix G)

$\frac{d}{dt}(c_{x}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}(t)}{dt})=-\int\int\int(ax^{a-1}y^{b}z^{c}\partial_{x}V(x,y,z,t)n(x,y,z,t))dxdydz\\\
+\int\int\int(a(a-1)x^{a-2}y^{b}z^{c}(\frac{(\partial_{x}n(x,y,z,t)^{2})}{4n(x,y,z,t)}+\frac{L_{x}(x,y,z,t)^{2}}{n(x,y,z,t)})\\\
+bax^{a-1}y^{b-1}z^{c}(\frac{(\partial_{x}n(x,y,z,t)\partial_{y}n(x,y,z,t))}{4n(x,y,z,t)}+\frac{L_{x}(x,y,z,t)L_{y}(x,y,z,t)}{n(x,y,z,t)})\\\
+cax^{a-1}y^{b}z^{c-1}(\frac{(\partial_{x}n(x,y,z,t)\partial_{z}n(x,y,z,t))}{4n(x,y,z,t)}+\frac{L_{x}(x,y,z,t)L_{z}(x,y,z,t)}{n(x,y,z,t)}))dxdydz\\\
-a(a-1)(\frac{(a-2)(a-3)}{4}\braket{x^{a-4}y^{b}z^{c}}+\frac{b(b-1)}{4}\braket{y^{b-2}x^{a-2}z^{c}}+\frac{c(c-1)}{4}\braket{z^{c-2}x^{a-2}y^{b}})\\\
$ (15)

where we define $L_{x}(x,y,z,t)$ as

$\displaystyle
L_{x}(x,y,z,t)\equiv-\int_{-\infty}^{x}\alpha_{x}(\mathbf{r^{\prime}},t)dx^{\prime}.$
(16)

Note that the time derivative of the coefficients $c_{u}^{abc}(t)$ themselves
are not needed as long as the product of the coefficient and the time
derivative of moments are propagated together as can be done using Eq. 15
since they enter together in Eq. 14 for obtaining $\alpha_{u}(t)$.

The key point of this section is that the single-particle wave function,
expressed as in Eq. 5, can be given in terms of the moments and their first-
order time derivative.

###  II.3 Time Derivatives of Moments

Combining the results of the two preceding Subsections II.2.1 and II.2.2, we
can now write the first-order and the second-order time derivatives of the
moments (Eqs. 3 and 4) without referring to the single-particle wave function.
The first-order time derivative of the moments is

$\frac{d\braket{x^{a}y^{b}z^{c}}(t)}{dt}=\int\int\int{x^{a}y^{b}z^{c}(\alpha_{x}(x,y,z,t)+\alpha_{y}(x,y,z,t)+\alpha_{z}(x,y,z,t))dxdydz}.$
(17)

The second-order time derivative of the moments is (see Appendix F)

$\frac{d^{2}\braket{x^{a}y^{b}z^{c}}(t)}{dt^{2}}=-\int\int\int(ax^{a-1}y^{b}z^{c}\partial_{x}V(x,y,z,t)n(x,y,z,t)\\\
+bx^{a}y^{b-1}z^{c}\partial_{y}V(x,y,z,t)n(x,y,z,t)+cx^{a}y^{b}z^{c-1}\partial_{z}V(x,y,z,t)n(x,y,z,t))dxdydz\\\
+\int\int\int(a(a-1)x^{a-2}y^{b}z^{c}(\frac{(\partial_{x}n(x,y,z,t)^{2})}{4n(x,y,z,t)}+\frac{L_{x}(x,y,z,t)^{2}}{n(x,y,z,t)})\\\
+2bax^{a-1}y^{b-1}z^{c}(\frac{(\partial_{x}n(x,y,z,t)\partial_{y}n(x,y,z,t))}{4n(x,y,z,t)}+\frac{L_{x}(x,y,z,t)L_{y}(x,y,z,t)}{n(x,y,z,t)})\\\
+2cax^{a-1}y^{b}z^{c-1}(\frac{(\partial_{x}n(x,y,z,t)\partial_{z}n(x,y,z,t))}{4n(x,y,z,t)}+\frac{L_{x}(x,y,z,t)L_{z}(x,y,z,t)}{n(x,y,z,t)})\\\
+b(b-1)y^{b-2}x^{a}z^{c}(\frac{(\partial_{y}n(x,y,z,t)^{2})}{4n(x,y,z,t)}+\frac{L_{y}(x,y,z,t)^{2}}{n(x,y,z,t)})\\\
+2cby^{b-1}x^{a}z^{c-1}(\frac{(\partial_{z}n(x,y,z,t)\partial_{y}n(x,y,z,t))}{4n(x,y,z,t)}+\frac{L_{z}(x,y,z,t)L_{y}(x,y,z,t)}{n(x,y,z,t)})\\\
+c(c-1)z^{c-2}x^{a}y^{b}(\frac{(\partial_{z}n(x,y,z,t)^{2})}{4n(x,y,z,t)}+\frac{L_{z}(x,y,z,t)^{2}}{n(x,y,z,t)}))dxdydz\\\
-a(a-1)(\frac{(a-2)(a-3)}{4}\braket{x^{a-4}y^{b}z^{c}}+\frac{b(b-1)}{2}\braket{y^{b-2}x^{a-2}z^{c}}+\frac{c(c-1)}{2}\braket{z^{c-2}x^{a-2}y^{b}})\\\
-b(b-1)(\frac{(b-2)(b-3)}{4}\braket{y^{b-4}x^{a}z^{c}}+\frac{c(c-1)}{2}\braket{z^{c-2}x^{a}y^{b-2}})\\\
-\frac{c(c-1)(c-2)(c-3)}{4}\braket{z^{c-4}x^{a}y^{b}}.\\\ $ (18)

In practical calculations, many terms in the above equation vanish for the
moments of low orders. These analytical expressions will be used later for a
numerical demonstration of this moment propagation theory in Sections III.1
and III.2.

###  II.4 Method of Electron Potential Energy Surface

The above theoretical formulation for quantum dynamics in terms of the moments
is exact in principle, but the numerical evaluation of the time derivatives of
the moments is highly complicated. When/if the principle of the energy
conservation holds for individual single-particle wave functions (which is
generally not the case for TD-KS orbitals in RT-TDDFT simulation), one can
derive an alternative simpler expression. Although it is not universally
applicable, such an idea of designing an effective potential energy for
quantum-mechanical particles (in a many-particle system) has been explored in
the literature [59, 60]. The energy of a quantum mechanical particle,
$\epsilon=\braket{\hat{H}}$, can be generally expressed in terms of the
moments such that

$\ \epsilon(\mathbf{x},\dot{\mathbf{x}})=\epsilon(\\{x_{i},\dot{x}_{i}\\})$
(19)

where we used the short-hand notation $x_{i}\equiv\braket{x^{i}}$ and the
integer index $i$ refers to the order of the moment. In principle, the integer
$i$ goes from 1 to $\infty$, but the dimension of $\mathbf{x}$ would be finite
in practice. Eq. 19 generally holds despite its unknown analytical form
because the single-particle wave function can be obtained from the moments and
their first-order time derivatives as shown in Section II.2. Assuming that
this single-particle energy is conserved (i.e.
$\dot{\epsilon}(\mathbf{x},\dot{\mathbf{x}})=0$), we here examine how the
second-order time derivative of the moments can be derived, using a one-
dimensional case for simplicity. Using the chain rule on the variables in
$\epsilon$ yields

$\sum_{k}(\frac{\partial\epsilon(\mathbf{x},\dot{\mathbf{x}})}{\partial
x_{k}}\dot{x}_{k}+\frac{\partial\epsilon(\mathbf{x},\dot{\mathbf{x}})}{\partial\dot{x}_{k}}\ddot{x}_{k})=0.$
(20)

In certain cases, like the harmonic oscillator model we discuss later, the
second-order time derivative for a particular order $k$, $\ddot{x}_{k}$,
depends only on the terms for which the energy derivative is taken with
respect to $x_{k}$ and $\dot{x}_{k}$. Then, the second-order time derivative
can be solved for each order $k$ using

$\frac{\partial\epsilon(\mathbf{x},\dot{\mathbf{x}})}{\partial
x_{k}}\dot{x}_{k}+\frac{\partial\epsilon(\mathbf{x},\dot{\mathbf{x}})}{\partial\dot{x}_{k}}\ddot{x}_{k}=0.$
(21)

Thus, the expression for the second-order time derivative is

$\ddot{x}_{k}\equiv\frac{d^{2}\braket{x^{k}}}{dt^{2}}=-\frac{\partial\epsilon(\mathbf{x},\dot{\mathbf{x}})}{\partial
x_{k}}\left(\frac{\partial\epsilon(\mathbf{x},\dot{\mathbf{x}})}{\partial\dot{x}_{k}}\right)^{-1}\dot{x}_{k}.$
(22)

According to this Electron Potential Energy Surface (EPES) expression, the
second-order time derivatives can be also evaluated from the energy $\epsilon$
by expressing it in terms of the moments and their first-order time
derivatives. When the single-particle energy does not change in time, this
EPES approach (Eq. 22) is numerically more convenient compared to the more
complicated expression, given by Eq. 18. Eq. 22 shows some similarity to the
expression in classical mechanics, where the second-order time derivative of
the position (i.e. acceleration) is proportional to the spatial derivative of
the potential energy. However, here, quantum effects are accounted for by an
extra term, which would be the unity in classical mechanics analogue.

###  II.5 Numerical Propagation of Moments

Formulating quantum dynamics as the time propagation of the moments offers a
convenient molecular dynamics (MD)-like scheme without explicitly integrating
the TD-SE. The moment propagation becomes particularly attractive when the
single-particle wave functions are spatially localized, as for those in the
MLWF or simply the Wannier gauge in RT-TDDFT [29, 8]. In principle, infinite
orders of the moments are necessary if one wishes to reproduce the exact
quantum dynamics by propagating the moments. In practice, however, essential
aspects of quantum dynamics relevant to calculating many physical properties
can be reproduced with the moments of low orders. For example, the frequency-
dependent dielectric function can be obtained from the dynamics of the first-
order moment of single-particle wave functions in RT-TDDFT[8], and we
anticipate that only a few orders of the moments need to be propagated for
accurately modeling the dynamics of the first-order moment.

Here, we propose a practical method for performing a numerical simulation
using the moment propagation theory developed above. We restrict our
discussion here to a one-dimensional case, for simplicity. As done in
classical MD simulation, the moments at a future time (i.e. $+\Delta t$ ) are
Taylor expanded to the 2nd order,

$\braket{x^{a}}(t+\Delta t)=\braket{x^{a}}(t)+\Delta
t\frac{d\braket{x^{a}}(t)}{dt}+\frac{1}{2}\Delta
t^{2}\frac{d^{2}\braket{x^{a}}(t)}{dt^{2}}+O(\Delta t^{3})$ (23)

where the time-step $\Delta t$ can be taken arbitrarily small so that
$O(\Delta t^{3})$ is negligible. In order to employ a numerical integrator for
using Eq. 23, the time derivatives of the moments need to be evaluated. For
many numerical integrators like Verlet algorithm for classical MD simulation,
an analytical expression is necessary only for the second-order time
derivatives. When the first-order time derivatives are also needed as in the
velocity Verlet algorithm, they can be calculated numerically from the second-
order time derivatives, etc. In classical mechanics, the second-order time
derivative is proportional to the gradient of the potential energy via
Newton’s second law and it is a function of the position of atoms. For the
quantum dynamics, the moment propagation theory gives the necessary expression
for the second-order time derivatives of the moments in terms of the moments
and their first-order time derivatives. The second-order time derivatives are
evaluated using Eq. 18. Alternatively, the EPES expression (Eq. 22) can be
used if the single-particle energy is conserved in the dynamics.

###  II.6 Analytical Solution for Simple Harmonic System

In order to demonstrate the moment propagation theory and examine its
feasibility using numerical methods in the following section, we use a simple
one-dimensional harmonic potential, $V(x)=x^{2}$, for which we can derive the
exact analytical solution. We consider a time-dependent Hamiltonian such that
simple homogeneous electric field is applied to the system. While this might
be considered a oversimplified model, note that the dynamics of individual
single-particle wave functions, localized on chemical bonds, are often
confined in such a localized potential well in the MLWF gauge[29]. We restrict
ourselves to the situation in which the quantum particle is initially in the
equilibrium ground state, and thus only the second-order Edgeworth (E2)
expansion (i.e. $n$=2 for Eq. 11) is needed for deriving the analytical
equation of motion [61, 55, 56]. The particle density is given by

$n(x,t)=\frac{1}{\sqrt{2\pi(\braket{x^{2}}(t)-(\braket{x}(t))^{2})}}\exp{(-\frac{1}{2}\frac{(x-\braket{x}(t))^{2}}{\braket{x^{2}}(t)-(\braket{x}(t))^{2}})}.$
(24)

As a commonly encountered situation in which the light is treated as a
classical electromagnetic field, a spatially homogeneous electric field is
used as the time-dependent external potential that perturbs the quantum system
such that $V(x,t)=x^{2}+c(t)x$ in the length gauge. One can find an exact
analytical solution for the second-order time derivatives using the above
described theory of moment propagation, and the second-order time derivative
of the first-order moment (see Appendix I) is

$\frac{d^{2}\braket{x}(t)}{dt^{2}}=-2\braket{x}(t)-c(t)$ (25)

and, for the second-order moment, we have

$\frac{d^{2}\braket{x^{2}}(t)}{dt^{2}}=-4\braket{x^{2}}(t)-2c(t)\braket{x}(t)+2((\frac{d\braket{x}(t)}{dt})^{2}+\frac{1+(\frac{dS(t)}{dt})^{2}}{4S(t)})$
(26)

where $S(t)=\braket{x^{2}}(t)-\left(\braket{x}(t)\right)^{2}$ is essentially
the spread or the variance of the corresponding single-particle wave function.
For the harmonic oscillator with a homogeneous external electric field, only
the two lowest orders of the moments are necessary to describe the quantum
dynamics exactly. The energy of the quantum particle, which is equivalent to
$\epsilon$ for this one-particle system, can be written (see Appendix I) as

$E(t)=\epsilon(t)=\braket{x^{2}}(t)+c(t)\braket{x}(t)+\frac{1}{2}((\frac{d\braket{x}(t)}{dt})^{2}+\frac{1+(\frac{dS(t)}{dt})^{2}}{4S(t)}).$
(27)

In the context of EPES method discussed in Section II.4, one could
alternatively use this analytical energy expression (Eq. 27) to derive the
second-order time derivatives of the moments as shown in Appendix I. As
expected, the EPES method gives the exact same analytical solutions for the
second-order time derivatives of the moments (i.e. Eq. 25 and 26).

##  III Numerical Demonstration

###  III.1 Harmonic Oscillator Potential

A simple one-dimensional system is used here to demonstrate the numerical
methods using the moment propagation theory described in the above sections.
As a proof-of-principle demonstration, we consider $V(x,t)=x^{2}+c(t)x$. This
represents a simple harmonic oscillator with time-dependent external electric
field applied to the system using the length gauge in Hamiltonian. We employ
spatially homogeneous electric field such that

$c(t)=\begin{cases}0.3\text{ a.u.}&0\text{ a.u.}<t<0.45\text{ a.u.}\\\ 0\text{
a.u.}&\text{else.}\end{cases}$ (28)

The electric field is applied initially to perturb the system for the first
449 time-steps (0.45 a.u.), and the system is propagated for a total of 12001
time steps (12 a.u.). For this harmonic potential with the time-dependent
electric field, the analytical solution can be derived as discussed in
Subsection II.6 as the reference for examining the numerical methods.

For the moment propagation theory (MPT)-based methods, 32 real-space grid
points are used as the basis over the spatial range of 7 a.u., and Eq. 23 was
used for the explicit time integration with a time-step of $\Delta t=$ 0.001
a.u. using the Newmark-beta method, with parameters $\gamma=\frac{1}{2}$ and
$\beta=\frac{1}{8}$ [62]. For the MPT methods, we show the simulation results
using the second-order Edgeworth approximation (MPT-E2) and the fourth-order
Edgeworth approximation (MPT-E4). Here, the order, $n=2,4$, in the Edgeworth
approximation (i.e. E2, E4), corresponds the $n$ in Eq. 7. The second-order
time derivatives are calculated explicitly using Eq. 18. Additionally, we
performed the MPT simulation by employing the EPES approach, in which the the
second-order time derivatives are calculated from the energy according to Eq.
22. The energy here is expressed as (see Appendix H)

$E=\int V(x)n(x)dx+\frac{1}{2}\int(\frac{(\frac{\partial n(x)}{\partial
x}^{2})}{4n(x)}+\frac{L(x)^{2}}{n(x)})dx$ (29)

where $L(x)$ is given by Eq. 16. In order to evaluate the particle density
here, the second-order Edgeworth approximation is employed, and we refer to
this approach as MPT-EPES-E2. For comparison, we also perform the standard
propagation of the wave function by directly integrating the time-dependent
Schrodinger equation, which we refer to as wave function theory (WFT). For the
WFT method, the wave function is represented using the real space grids with
100 evenly spaced points over the same 7 a.u. range as for the MPT methods.
The Crank-Nicholson method was used to numerically integrate the wave function
with a time-step of $\Delta t=$ 0.001 a.u.[63]. In the limit of infinitely
dense grid points and an infinitely small time-step, all the methods converge
to the exact analytical result for this particular harmonic potential case.

(a)

(b)

Figure 1:  (a) Total energy as a function of time for the exactly-solvable
one-dimensional (1D) harmonic potential model, $V(x,t)=x^{2}+c(t)x$. The
homogeneous electric field is applied for t=0-0.45 a.u., and the reference
energy (E=0.0 a.u.) is set to that of the initial ground state at t=0 a.u..
The three MPT methods and the standard wave function theory (WFT) method are
compared to the exact analytical solution. (b) The error on the total energy
with respect to the exact analytical solution. See the main text for more
details.

(a)

(b)

Figure 2:  Comparison of the three MPT methods and the standard wave function
theory (WFT) method to the analytical solution for (a) the first-order moment
as a function of time and (b) its error with respect to the exact analytical
solution. See also Figure 1 for additional details.

Figure 1 shows the total energy, $E$, as a function of time, and the energy is
expected to remain constant after the perturbing electric field is turned off.
Having the exact analytical solution for this system, we know that the energy
change must be $\Delta E=8.771\times 10^{-3}$ a.u. in principle if all
numerical parameters are fully converged and computations are performed
exactly. With the real-space basis of 32 grid points, the MPT methods (MPT-E2
and MPT-E4) yield $\Delta E$ values of $8.772\times 10^{-3}$ a.u. and
$8.770\times 10^{-3}$ a.u., respectively. The MPT method with the EPES
approach (MPT-EPES-E2) gives a slightly smaller value of $\Delta E$ of
$8.766\times 10^{-3}$ a.u.. Compared to the MPT-E2 method, the MPT-EPES-E2
method gives a larger error because of the numerical evaluation of the second-
order time derivative in Eq. 22 through the partial derivative of the energy
with respect to the moments and their time derivatives. With a larger real-
space basis of 100 grid points, the standard approach of propagating the wave
function (WFT) gives a similar value of $\Delta E$=$8.756\times 10^{-3}$ a.u..
While the WFT method gives the largest error for this particular harmonic
oscillator potential, all the numerical methods show sufficient accuracy in
terms of the relative error on the energy. Table 1 also summarizes the energy
conservation for all methods, and all numerical methods show that the energy
drift/deviation is at least several orders of magnitude smaller than the
energy change due to the external electric field.

Method |  WFT |  MPT-E2 |  MPT-E4 |  MPT-EPES-E2
---|---|---|---|---
Energy Change, $\Delta E$ (a.u.) |  $8.756\times 10^{-3}$ |  $8.772\times 10^{-3}$ |  $8.770\times 10^{-3}$ |  $8.766\times 10^{-3}$
Energy Drift (a.u.) per 1 a.u. time |  $1.150\times 10^{-12}$ |  $1.240\times 10^{-7}$ |  $2.073\times 10^{-8}$ |  $8.955\times 10^{-8}$
Oscillation Frequency (a.u.) |  $1.41296$ |  $1.41421$ |  $1.41420$ |  $1.41409$
Table 1:  Energy change, energy drift per 1 a.u. time, and the oscillation
frequency of the first-order moment for the simulation shown in Figure 1. The
exact analytical value of $\Delta E$ is $8.771\times 10^{-3}$ a.u., and the
exact analytical value for the oscillation frequency is $\sqrt{2}=1.41421$.

The first-order moments are of particular importance in quantum dynamics
because they can be directly related to physical properties such as the
frequency-dependent polarizability/conductivity, from which, for instance, the
optical absorption spectrum can be calculated in electronic structure theory
[64, 65]. As expected, the first-order moment of this quantum particle
oscillates similarly to classical particle for this simple harmonic potential.
Figure 2 shows the comparison for the time evolution of the first-order moment
among the newly-developed MPT methods and the standard WFT approach. Although
the basis set size necessary for convergence is system-dependent in general,
the MPT methods exhibit a smaller numerical error with only 32 grid points
while the WFT simulation here used 100 grid points. Note that the differences
among the numerical methods are rather negligible (see Fig. 2(a) ), and they
all show only small errors with respect to the exact analytical result, as
seen in Figure 2(b). The oscillation frequency of the first-order moment,
which would be related to the polarizability/conductivity in electronic
structure calculations, does not deviate appreciatively from the exact
analytical value of $\sqrt{2}$ for any of the numerical methods as shown in
Table 1.

Figure 3 shows the time evolution of the second-order moment according to the
different numerical methods, along with the errors with respect to the exact
analytical solution. The second-order moment is related to the spatial
variance (i.e. the spread) of the quantum particle. In this particular case,
all the MPT methods show much smaller errors than the WFT method although all
the errors are sufficiently small (see Fig. 3(b)). Comparing to the exact
analytical solution, these results numerically confirm the validity of the
newly-developed MPT approach and the practical feasibility of employing
numerical methods based on this new theoretical formulation. As demonstrated
here, the MPT numerical methods perform as well as the WFT method, if not
better, for this simple case with the analytical solution available.

(a)

(b)

Figure 3:  Comparison of the three MPT methods and the standard wave function
theory (WFT) method to the analytical solution for (a) the second-order moment
as a function of time and (b) its error with respect to the exact analytical
solution. See also Figure 1 for additional details.

###  III.2 Anharmonic Morse Potential

Although an exact analytical solution is not available for this anharmonic
potential, we also examined the performance of the numerical MPT method for
the Morse potential,

$V(x,t)=a(1-e^{-\sqrt{\frac{1}{a}}x})^{2}+c(t)x$ (30)

where $a$ determines the anharmonicity of this Morse potential, and we set
$a=10$ a.u. (see Fig. 4). The same external electric field, given by Eq. 28,
is applied to the system through the $c(t)x$ term in the length gauge. All
computational parameters are kept the same as the harmonic oscillator
potential case discussed above, except for the grid size of 6 a.u.. Figure 4
shows the total energy as a function of time for the MPT method with
increasingly higher orders of moments, along with the standard WFT method for
comparison. Although an exact analytical solution is not available for this
model, the energy must be conserved in a numerically exact propagation. Figure
4(a) shows better energy conservation during the dynamics with less
fluctuations for the MPT method with increasingly higher orders of the
moments. In this particular case, the WFT method shows slightly better energy
conservation than the MPT method even when up to the 4th-order moments are
included (i.e. MPT-E4). Convergence of the MPT method and the standard WFT
method can be seen in Figure 4(b), as the Edgeworth expansion becomes more
complete by including higher orders of the moments.

Figure 5 shows the time evolution of the first-order and the second-order
moments in the quantum dynamics. The MPT dynamics becomes increasingly closer
to the conventional WFT method as the higher orders of the moments are
included in the MPT method. This is particularly apparent for the second-order
moment (see Fig. 5(b)) as the MPT-E4 result closely resembles the WFT method
result while significant deviations can be seen for the MPT-E2 result.

(a)

(b)

Figure 4:  (a) Total energy as a function of time for the anharmonic 1D model
with Morse potential (Eq. 30). The inset shows the Morse potential energy
function $V(x)$. The homogeneous electric field is applied for t=0-0.45 a.u.
The MPT method with an increasing order of the moments is compared to the
standard wave function theory (WFT) method. (b) The differences with respect
to the WFT method result are shown.

(a)

(b)

Figure 5:  (a) first-order moment and (b) second-order moment as a function of
time for the anharmonic 1D model (Eq. 30). The MPT method with an increasing
order of the moments is compared to the standard wave function theory (WFT)
method. See also Figure 4 for additional details.

Although the moment propagation theory (MPT) is an exact reformulation of
quantum dynamics in terms of the moments instead of the single-particle wave
function, its usefulness might not be obvious in the above proof-of-principle
demonstrations. A key advantage of this novel MPT is that machine-learning
techniques can be readily used for modeling the second-order time derivatives
of the moments instead of explicitly calculating them. For classical MD
simulations, machine-learning techniques, particularly modern neural network
models, are widely used for evaluation of the potential energy and the atomic
forces, and a similar approach can be conveniently adapted for quantum
dynamics simulation using the MPT. This idea is briefly demonstrated in the
following section.

##  IV Machine-Learning Time Derivatives with Artificial Neural Network

Simulating quantum dynamics is highly nontrivial and computationally
expensive, and this often limits our ability to study long-time dynamics of
quantum systems and/or statistical properties that derive from an ensemble of
quantum dynamics trajectories. There have been a number of works to employ
modern machine-learning approaches like artificial neural network (ANN) models
for efficiently performing such quantum dynamics simulations instead of
numerically integrating the time-dependent Schrodinger equation [48, 49, 66,
35] In this section, we demonstrate how the newly-developed MPT offers an
convenient avenue for developing an ANN model for quantum dynamics simulation.
For classical MD simulations, the use of ANN has become widespread in recent
years [67, 68, 69]. In particular, machine-learning techniques are used to
construct an accurate potential energy function (i.e. force field) as an ANN
through training with first-principles quantum-mechanical calculations [70,
71, 72, 73, 69, 74]. This allows for an efficient computation of the accurate
force on atoms, which is proportional to the second-order time derivative of
the classical particle positions. We here demonstrate a machine-learning of
the second-order time derivative of the moments as an ANN model for performing
the quantum dynamics simulation within the MPT framework. In particular, we
examine the construction of ANN models for the second-order time derivatives
of the moments using $\braket{x}$, $\braket{x^{2}}$,
$(\frac{\partial\braket{x}}{\partial t})^{2}$, and $(\frac{\partial
S}{\partial t})^{2}$ as the input descriptors for the harmonic potential
model. The choice of these descriptors are motivated by the analytical
solution. A simple ANN model was used here, consisting of two hidden layers
with eight Scaled Exponential Linear Units (SELU) activation nodes[75]. The
input descriptors are multiplied by weights and they acquire an added bias
when moving to the next layer. The activation function, SELU, acts on each
input for each node and its output is fed into the next layer. There are many
different approaches for developing ANNs, but we chose this simple proof-of-
concept architecture [67]. The mean squared error (MSE) was used for the loss
function, and the ANN training used the Adam optimizer[75]. Each ANN model was
trained for at least a thousand epochs until the MSE reaches a minimum. The
training data was generated from propagating the wave function according to
the TD-SE with a time-step of $\Delta t=0.01$ a.u.. For training each ANN
model, a total of 9498 data points were used. ANN model then replaces the
numerical evaluation of the second-order time derivative of the moments using
Eq. 18 in the quantum dynamics simulation based on the MPT. Figure 6 shows the
first-order and second-order moments from the MPT simulation using the ANN
model (MPT-ANN) for the harmonic potential, compared against its own training
data as a proof-of-principle demonstration. The MPT-ANN simulation
successfully reproduces the time evolution of the moments to a great accuracy.
For the second-order moment, a slight deviation is observed for longer times
while the oscillation frequency remains accurately reproduced. Improving the
machine-learning techniques for quantum dynamics simulation with the MPT will
be a topic of a separate study in the future. The computational procedure/cost
is numerically comparable to that of performing classical MD simulation using
the ANN model for the force evaluation.

We envision that machine-learning techniques like ANN models can be used to
construct an accurate model of the second-order time derivative of the moments
in general. For the harmonic potential model, the exact analytical expression
can be derived and thus it was possible to deduce what physically-meaningful
input descriptors were in terms of the moments for the ANN construction.
However, in realistic first-principles electronic structure theory
descriptions such as when electrons are represented by time-dependent
maximally-localized Wannier functions in RT-TDDFT simulation [29], the
potential that individual electrons are subjected to have no analytical
solutions in general although confining in many cases. We briefly examine here
how the quantum dynamics might be sensitive to the input descriptors in the
ANN model using this exactly-solvable harmonic potential. Table 2 shows the
oscillation frequency and the root mean squared error (RMSE) of the first-
order moment as well as the RMSE of the second-order moment using different
input descriptors. Even with only $\braket{x}$ and $\braket{x^{2}}$ as the
input descriptors for the ANN model, the RMSE remains quite small for both the
first-order and second-order moments in the dynamics as shown in Figure 7.
Additionally, the oscillation frequencies from the MPT simulation using the
ANN models with different input descriptors are essentially the same as the
frequency obtained from the reference WFT simulation as seen in Table 2. In
summary, we showed here how the newly-developed MPT provides a convenient
theoretical framework for accelerating the quantum dynamics simulation through
machine-learning the second-order time derivatives using a simple ANN model.

(a)

(b)

Figure 6:  The use of ANN model for the calculation of the second-order time
derivative of the moments in the MPT method. Comparison of (a) the first-order
moment and (b) the second-order moment in the MPT simulation based on the ANN
model and explicit WFT simulation. Based on the analytical solution for this
exactly-solvable harmonic potential model, the input descriptors for the ANN
model consist of $\braket{x}$, $\braket{x^{2}}$,
$(\frac{\partial\braket{x}}{\partial t})^{2}$, and $(\frac{\partial
S}{\partial t})^{2}$. The WFT simulation was used also to provide the training
data for the ANN model.

(a)

(b)

Figure 7:  The use of ANN model for the calculation of the second-order time derivative of the moments in the MPT method. Additional details are same as in Figure 6 but only $\braket{x}$ and $\braket{x^{2}}$ are used as the input descriptors for the ANN model. Comparison of (a) the first-order moment and (b) the second-order moment in the MPT simulation based on the ANN model and explicit WFT simulation.  Inputs |  $\braket{x}$ Frequency |  $\braket{x}$ RMSE |  $\braket{x^{2}}$ RMSE
---|---|---|---
$\braket{x}$, $\braket{x^{2}}$, $(\frac{\partial\braket{x}}{\partial t})^{2}$, and $(\frac{\partial S}{\partial t})^{2}$ |  $1.40243$ |  $3.41483\times 10^{-4}$ |  $1.87218\times 10^{-4}$
$\braket{x}$, $\braket{x^{2}}$, and $(\frac{\partial S}{\partial t})^{2}$ |  $1.40260$ |  $2.86096\times 10^{-4}$ |  $2.67096\times 10^{-4}$
$\braket{x}$, $\braket{x^{2}}$, and $(\frac{\partial\braket{x}}{\partial t})^{2}$ |  $1.40266$ |  $1.54922\times 10^{-4}$ |  $3.39435\times 10^{-4}$
$\braket{x}$ and $\braket{x^{2}}$ |  $1.40249$ |  $2.36722\times 10^{-4}$ |  $3.67558\times 10^{-4}$
Table 2:  Root-mean-square error (RMSE) for the first-order and the second-
order moments, together with the oscillation frequency of the first-order
moment in the MPT simulations with different input descriptors in the ANN
model.

##  V Conclusion

We presented a novel theoretical formulation of the quantum dynamics in terms
of the moments within the single-particle description. The aim of propagating
the moments instead of single-particle wave functions is to reduce the high
computational cost by expressing the quantum dynamics in terms of increasing
orders of the moments. By doing so, the equation of motion for the quantum
dynamics resembles that of classical molecular dynamics, and therefore modern
machine-learning techniques might be employed more effectively. We derived
analytical expressions for the first-order and second-order time derivatives
of the moments, and they are used to develop a numerical method for performing
the quantum dynamics by expanding the moments in the Taylor series as done
similarly for atom positions in classical molecular dynamics simulation. We
demonstrated a few numerical schemes using the newly-developed moment
propagation theory (MPT) on the harmonic oscillator potential, for which an
exact analytical solution can be derived. We also applied the MPT method to an
anharmonic Morse potential and compared it to the standard approach of
propagating the wave function according to the time-dependent Schrodinger
equation.

The new moment propagation theory (MPT) offers a particularly compact
representation suitable for developing machine-learned models for performing
efficient quantum dynamics simulations. We demonstrated the use of a simple
artificial neural network (ANN) model for machine-learning the second-order
time derivatives of the moments so that their explicit evaluation can be
circumvented in numerical simulations. Quantum dynamics simulation of
electrons in molecules and condensed matter systems is where this new approach
can find its immediate application. The new MPT framework is particularly
powerful in practice when only low orders of the moments are necessary. In the
context of electronic structure theory, real-time time-dependent density
functional theory (RT-TDDFT) simulation can be cast in the Wannier gauge for
simulating electron dynamics[8, 29]. The maximally-localized Wannier functions
(MLWFs) are spatially localized single-particle orbitals, and the MPT is
likely to offer a convenient theoretical framework for developing machine-
learning models for RT-TDDFT simulations performed in the Wannier gauge. As
for the scaling of the MPT method with respect to the order of moments, the
computational cost per each MLWF scales as
$\frac{1}{6}n^{3}+n^{2}+\frac{11}{6}n$, where $n$ is the highest order of
moments included. The computational scales linearly with the number of MLWFs
in a system. Having the formal theoretical formulation established in this
work, future work will explore the use of this new MPT in the context of
first-principles electronic structure theory for studying real chemical
systems. In particular, the use of ANNs and other related machine-learning
methods for simulating electron dynamics will be examined in the context of
RT-TDDFT simulation.

###### Acknowledgements.

N.B. was supported by the Chapel Family summer research fellowship. This work
was supported by the National Science Foundation, under No. CHE-1954894.

## Author Declarations

### Conflict of Interest

The authors have no conflicts to disclose.

### Author Contributions

N.B. led the work and performed all the calculations. N.B. and Y.K. conceived
of the presented idea. All authors discussed the results and contributed to
the final manuscript.

### Data Availability

The data that support the findings of this study are available from the
corresponding author upon reasonable request.

##  Appendix A Notations

We used the following notations to define several gradient operations in
Appendix,

$\nabla=(\partial_{x},\partial_{y},\partial_{z}),$ (A1)

where $\partial_{x}\equiv\frac{\partial}{\partial x}$. Accordingly,

$\displaystyle\nabla^{2}$ $\displaystyle=\nabla\cdot\nabla,$ (A2)
$\displaystyle\nabla^{3}$ $\displaystyle=\nabla\nabla\cdot\nabla,$
$\displaystyle\nabla^{4}$
$\displaystyle=(\nabla\cdot\nabla)(\nabla\cdot\nabla).$

$\nabla^{2}$ is the Laplacian operator, and therefore we have

$\nabla^{2}(\nabla A\cdot\nabla B)=\nabla^{3}A\cdot\nabla
B+2(\nabla\otimes\nabla A)\cdot(\nabla\otimes\nabla B)+\nabla
A\cdot\nabla^{3}B,$ (A3)

where tensor product of $\nabla$ is

$\nabla\otimes\nabla=(\partial_{xx},\partial_{xy},\partial_{xz},\partial_{xy},\partial_{yy},\partial_{zy},\partial_{xz},\partial_{zy},\partial_{zz}).$
(A4)

##  Appendix B Derivation of Eq. 3

For a single-particle system in three dimensions, the derivative with respect
to time for the moment is defined as:

$\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}\equiv\int
x^{a}y^{b}z^{c}\frac{dn(x,y,z,t)}{dt}d^{3}r.$ (B1)

Here we used $\int d^{3}r$ to indicate the spacial integral $\int\int\int
dxdydz$ for brevity. Using
$n(x,y,z,t)=\psi^{*}(x,y,z,t)\psi(x,y,z,t)=|\psi(x,y,z,t)|^{2}$ and the TD-SE
(Eq. 2), one has

$\displaystyle\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}=$ $\displaystyle\int
x^{a}y^{b}z^{c}\left[\left(\frac{i}{2}\nabla^{2}\psi-
iV\psi\right)\psi^{*}+\psi\left(-\frac{i}{2}\nabla^{2}\psi^{*}+iV\psi^{*}\right)\right]d^{3}r$
(B2) $\displaystyle=$ $\displaystyle\frac{i}{2}\int
x^{a}y^{b}z^{c}\left(\nabla^{2}\psi\psi^{*}-\psi\nabla^{2}\psi^{*}\right)d^{3}r$
$\displaystyle=$
$\displaystyle-\frac{i}{2}\int\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\left(\nabla\psi\psi^{*}-\psi\nabla\psi^{*}\right)d^{3}r$
$\displaystyle=$
$\displaystyle-\frac{i}{2}\int\left[\nabla^{2}\left(x^{a}y^{b}z^{c}\right)n+2\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\nabla\psi\psi^{*}\right]d^{3}r,$

where the explicit spatial and time dependence of the wave function
$\psi(x,y,z,t)$, time-dependent potential $V(x,y,z,t)$, and density
$n(x,y,z,t)$ are omitted for brevity. The explicit expression is:

$\displaystyle\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}=$
$\displaystyle-i\int\int\int(ax^{a-1}y^{b}z^{c}\partial_{x}\psi\psi^{*}+by^{b-1}x^{a}z^{c}\partial_{y}\psi\psi^{*}+cz^{c-1}x^{a}y^{b}\partial_{z}\psi\psi^{*}$
(B3) $\displaystyle+$
$\displaystyle\frac{a(a-1)}{2}x^{a-2}y^{b}z^{c}n(x,y,z,t)+\frac{b(b-1)}{2}y^{b-2}x^{a}z^{c}n(x,y,z,t)$
$\displaystyle+$
$\displaystyle\frac{c(c-1)}{2}z^{c-2}x^{a}y^{b}n(x,y,z,t))dxdydz.$

Using
$\nabla\psi\psi^{*}-\psi\nabla\psi^{*}=2i\operatorname{Im}(\nabla\psi\psi^{*})$,
Eq. B2 can be also written concisely as

$\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}=\operatorname{Re}\left[-i\int\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\nabla\psi\psi^{*}d^{3}r\right].$
(B4)

##  Appendix C Derivation of Eq. 4

Given Eq. B4, the second-order time derivative of the moment is:

$\frac{d^{2}\braket{x^{a}y^{b}z^{c}}}{dt^{2}}=\operatorname{Re}\left[-i\frac{d}{dt}\int\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\nabla\psi\psi^{*}d^{3}r\right].$
(C1)

For the rest of the derivations in this section, the real part is implied.
Analogously to the previous section, $n(x,y,z,t)=|\psi(x,y,z,t)|^{2}$ and TD-
SE (Eq.2) gives,

$\displaystyle\frac{d^{2}\braket{x^{a}y^{b}z^{c}}}{dt^{2}}=$
$\displaystyle-i\int\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\left(\nabla\frac{\partial\psi}{\partial
t}\psi^{*}+\nabla\psi\frac{\partial\psi^{*}}{\partial t}\right)d^{3}r$ (C2)
$\displaystyle=$
$\displaystyle-i\int\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\left[\nabla\left(\frac{i}{2}\nabla^{2}\psi-
iV\psi\right)\psi^{*}+\nabla\psi\left(-\frac{i}{2}\nabla^{2}\psi^{*}+iV\psi^{*}\right)\right]d^{3}r$
$\displaystyle=$
$\displaystyle\int\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\left[-\left(\nabla
V\right)n+\frac{1}{2}\left(\nabla^{3}\psi\psi^{*}-\nabla\psi\nabla^{2}\psi^{*}\right)\right]d^{3}r.$

Expanding Eq. C2, we have

$\frac{d^{2}\braket{x^{a}y^{b}z^{c}}}{dt^{2}}=-\int\nabla(x^{a}y^{b}z^{c})\cdot\nabla
Vnd^{3}r+\int\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\frac{1}{2}(\nabla^{3}\psi\psi^{*}-\nabla\psi\nabla^{2}\psi^{*})d^{3}r.$
(C3)

We now expand the right-hand side and apply the integration by parts,

$\displaystyle\frac{d^{2}<x^{a}y^{b}z^{c}>}{dt^{2}}=$
$\displaystyle-\int\nabla(x^{a}y^{b}z^{c})\cdot\nabla
Vnd^{3}r+\int\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\frac{1}{2}(\nabla^{3}\psi\psi^{*})d^{3}r$
(C4)
$\displaystyle-\frac{1}{2}\int\nabla^{2}(\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\nabla\psi)\psi^{*}d^{3}r$
$\displaystyle=$ $\displaystyle-\int\nabla(x^{a}y^{b}z^{c})\cdot\nabla
Vnd^{3}r-\frac{1}{2}\int(\nabla^{3}\left(x^{a}y^{b}z^{c}\right)\cdot\nabla\psi$
$\displaystyle+2\nabla\otimes\nabla\left(x^{a}y^{b}z^{c}\right)\cdot(\nabla\otimes\nabla\psi))\psi^{*}d^{3}r.$

To make further progress on this equation to the working form, we used the
real part of Eq. B2,

$\operatorname{Re}\left(\int\nabla\left(x^{a}y^{b}z^{c}\right)\cdot\nabla\psi\psi^{*}d^{3}r\right)=-\frac{1}{2}\int\nabla^{2}\left(x^{a}y^{b}z^{c}\right)nd^{3}r.$
(C5)

Then, Eq. C4 is expressed as

$\frac{d^{2}\braket{x^{a}y^{b}z^{c}}}{dt^{2}}=\int\operatorname{Re}\left[-\nabla(x^{a}y^{b}z^{c})\cdot\nabla
Vn+\frac{1}{4}\nabla^{4}\left(x^{a}y^{b}z^{c}\right)n-(\nabla\otimes\nabla\left(x^{a}y^{b}z^{c}\right)\cdot(\nabla\otimes\nabla\psi))\psi^{*}\right]d^{3}r,$
(C6)

or completely expanded as

$\displaystyle\frac{d^{2}\braket{x^{a}y^{b}z^{c}}}{dt^{2}}=$
$\displaystyle\operatorname{Re}(-\int\int\int(ax^{a-1}y^{b}z^{c}\partial_{x}V(x,y,z,t)n(x,y,z,t)+bx^{a}y^{b-1}z^{c}\partial_{y}V(x,y,z,t)n(x,y,z,t)$
(C7) $\displaystyle+$ $\displaystyle
cx^{a}y^{b}z^{c-1}\partial_{z}V(x,y,z,t)n(x,y,z,t))d^{3}r$ $\displaystyle-$
$\displaystyle\int\int\int(a(a-1)x^{a-2}y^{b}z^{c}\partial_{x}^{2}\psi\psi^{*}+2bax^{a-1}y^{b-1}z^{c}\partial_{x}\partial_{y}\psi\psi^{*}$
$\displaystyle+$ $\displaystyle
2cax^{a-1}y^{b}z^{c-1}\partial_{x}\partial_{z}\psi\psi^{*}+b(b-1)y^{b-2}x^{a}z^{c}\partial_{y}^{2}\psi\psi^{*}$
$\displaystyle+$ $\displaystyle
2cby^{b-1}x^{a}z^{c-1}\partial_{y}\partial_{z}\psi\psi^{*}+c(c-1)z^{c-2}x^{a}y^{b}\partial_{z}^{2}\psi\psi^{*})d^{3}r$
$\displaystyle+$
$\displaystyle\frac{1}{4}\left(a(a-1)((a-2)(a-3)\braket{x^{a-4}y^{b}z^{c}}+2b(b-1)\braket{y^{b-2}x^{a-2}z^{c}}\right.$
$\displaystyle+$ $\displaystyle 2c(c-1)\braket{z^{c-2}x^{a-2}y^{b}})$
$\displaystyle+$ $\displaystyle
b(b-1)((b-2)(b-3)\braket{y^{b-4}x^{a}z^{c}}+2c(c-1)\braket{z^{c-2}x^{a}y^{b-2}})$
$\displaystyle+$
$\displaystyle\left.c(c-1)(c-2)(c-3)\braket{z^{c-4}x^{a}y^{b}}\right)).$

##  Appendix D Derivation of Eq. 12

The phase of the wave function in the polar form can be expressed as
$\theta=atan2(\operatorname{Im}\psi,\operatorname{Re}\psi)$, and the real and
imaginary parts of the wave function are
$\operatorname{Re}\psi=\sqrt{n}\cos{\theta}$ and
$\operatorname{Im}\psi=\sqrt{n}\sin{\theta}$. Using the definition (from the
main text)
$\alpha_{u}\equiv\frac{i}{2}(\partial_{u}^{2}\psi\psi^{*}-\psi\partial_{u}^{2}\psi^{*})$,
let us write

$L_{u}(x,y,z,t)=-\int_{-\infty}^{u}\alpha_{u}(\mathbf{r}^{\prime},t)du^{\prime}=\frac{-i}{2}(\frac{\partial\psi}{\partial
u}\psi^{*}-\frac{\partial\psi^{*}}{\partial u}\psi)$ (D1)

where $u$ can be $x$, $y$, or $z$ while $\textbf{r}^{\prime}$ denotes the
dependence of the function on all three coordinates for brevity. We can
rearrange this expression by writing out the real and imaginary parts of the
right-hand side

$\displaystyle L_{u}$ $\displaystyle=\frac{-i}{2}(\frac{\partial\psi}{\partial
u}\psi^{*}-\frac{\partial\psi^{*}}{\partial u}\psi)$ (D2)
$\displaystyle=\frac{-i}{2}(i\operatorname{Im}(\frac{\partial\psi}{\partial
u}\psi^{*})+\operatorname{Re}(\frac{\partial\psi}{\partial
u}\psi^{*})-i\operatorname{Im}(\frac{\partial\psi^{*}}{\partial
u}\psi))-\operatorname{Re}(\frac{\partial\psi^{*}}{\partial u}\psi)))$
$\displaystyle=\frac{-i}{2}(2i\operatorname{Im}(\frac{\partial\psi}{\partial
u}\psi^{*}))$ $\displaystyle=\frac{\partial\operatorname{Im}\psi}{\partial
u}\operatorname{Re}\psi-\frac{\partial\operatorname{Re}\psi}{\partial
u}\operatorname{Im}\psi.$

The derivative of the phase can be expressed as

$\frac{\partial\theta}{\partial
u}=-\frac{\operatorname{Im}\psi\frac{\partial\operatorname{Re}\psi}{\partial
u}}{(\operatorname{Im}\psi)^{2}+(\operatorname{Re}\psi)^{2}}+\frac{\operatorname{Re}\psi\frac{\partial\operatorname{Im}\psi}{\partial
u}}{(\operatorname{Im}\psi)^{2}+(\operatorname{Re}\psi)^{2}}.$ (D3)

Using $n=(\operatorname{Im}\psi)^{2}+(\operatorname{Re}\psi)^{2}$ and Eq. D2,
we have $\frac{\partial\theta}{\partial u}=\frac{L_{u}}{n}$ and thus

$\frac{\partial\theta(x,y,z,t)}{\partial
u}=\frac{L_{u}(x,y,z,t)}{n(x,y,z,t)}=-(n(x,y,z,t))^{-1}\int_{-\infty}^{u}\alpha_{u}(\mathbf{r}^{\prime},t)du^{\prime},\\\
$ (D4)

where $\alpha_{u}\equiv\frac{i}{2}(\frac{\partial^{2}\psi}{\partial
u^{2}}\psi^{*}-\frac{\partial^{2}\psi^{*}}{\partial u^{2}}\psi)$. For each
coordinate variable $u$, we have

$\displaystyle\frac{\partial\theta(x,y,z,t)}{\partial
x}=\frac{L_{x}(x,y,z,t)}{n(x,y,z,t)}=-(n(x,y,z,t))^{-1}\int_{-\infty}^{x}\alpha_{x}(x^{\prime},y,z,t)dx^{\prime}$
(D5) $\displaystyle\frac{\partial\theta(x,y,z,t)}{\partial
y}=\frac{L_{y}(x,y,z,t)}{n(x,y,z,t)}=-(n(x,y,z,t))^{-1}\int_{-\infty}^{y}\alpha_{y}(x,y^{\prime},z,t)dy^{\prime}$
$\displaystyle\frac{\partial\theta(x,y,z,t)}{\partial
z}=\frac{L_{z}(x,y,z,t)}{n(x,y,z,t)}=-(n(x,y,z,t))^{-1}\int_{-\infty}^{z}\alpha_{z}(x,y,z^{\prime},t)dz^{\prime}.$

##  Appendix E Derivation of Eq. 14

We start by expressing the wave function in polar form as
$\psi(x,y,z,t)=\sqrt{n(x,y,z,t)}e^{i\theta(x,y,z,t)}$. While $n(x,y,z,t)$ can
be obtained from the Edgeworth series of the moments straightforwardly, the
expression for the phase in terms of the moments is required. Using the TD-SE,
the first-order time-derivative of the particle density $n$ can be expressed
as

$\displaystyle\frac{dn}{dt}$ $\displaystyle=(\frac{i}{2}\nabla^{2}\psi-
iV\psi)\psi^{*}+\psi(-\frac{i}{2}\nabla^{2}\psi^{*}+iV\psi^{*})$ (E1)
$\displaystyle=\frac{i}{2}(\partial_{x}^{2}\psi\psi^{*}-\psi\partial_{x}^{2}\psi^{*}+\partial_{y}^{2}\psi\psi^{*}-\psi\partial_{y}^{2}\psi^{*}+\partial_{z}^{2}\psi\psi^{*}-\psi\partial_{z}^{2}\psi^{*})$
$\displaystyle=\alpha_{x}+\alpha_{y}+\alpha_{z},$

where
$\alpha_{u}\equiv\frac{i}{2}(\partial_{u}^{2}\psi\psi^{*}-\psi\partial_{u}^{2}\psi^{*})$.
Inserting this expression into Eq. B1, we obtain

$\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}=\int
x^{a}y^{b}z^{c}(\alpha_{x}+\alpha_{y}+\alpha_{z})d^{3}r.\\\ $ (E2)

By defining

$c_{u}^{abc}(t)\equiv\left(\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}\right)^{-1}\int
x^{a}y^{b}z^{c}(\alpha_{u}(x,y,z,t))d^{3}r,\\\ $ (E3)

and Eq. E2 can be expressed as

$\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}=c_{x}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}+c_{y}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}+c_{z}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}.$
(E4)

The time-derivative of the particle density can be also rearranged to read

$\frac{dn}{dt}=\sum_{a,b,c}\frac{\partial
n}{\partial\braket{x^{a}y^{b}z^{c}}}\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}.$
(E5)

By substituting the left-hand side of this equation with Eq. E1 and the time-
derivative on the right-hand side with Eq. E4, Eq. E5 reads

$\alpha_{x}+\alpha_{y}+\alpha_{z}=\sum_{a,b,c}\frac{\partial
n}{\partial\braket{x^{a}y^{b}z^{c}}}\\{c_{x}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}+c_{y}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}+c_{z}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}\\}.$
(E6)

Thus, for the spatial derivative in each Cartesian direction, we have Eq. 14

$\alpha_{u}=\sum_{a,b,c}c_{u}^{abc}(t)\frac{\partial
n}{\partial\braket{x^{a}y^{b}z^{c}}}\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}.$
(E7)

##  Appendix F Derivation of Eq. 18

For brevity, let us derive the equation for the one-dimensional case (using
$x$) in detail. The first step in the derivation is to show

$Re(\psi^{*}(x,t)\frac{\partial^{2}\psi(x,t)}{\partial
x^{2}})=(\frac{1}{2}\frac{\partial^{2}n(x,t)}{\partial
x^{2}}-\frac{(\frac{\partial n(x,t)}{\partial
x})^{2}}{4n(x,t)}-\frac{L(x,t)^{2}}{n(x,t)}),$ (F1)

where the right-hand side is expressed in terms of the particle density and
the moments without explicit dependence on the wave function. Let us first
write the left-hand side of Eq. F1 as

$\displaystyle LHS=$
$\displaystyle\operatorname{Re}(\operatorname{Re}\psi(x,t)\frac{\partial^{2}\operatorname{Re}\psi(x,t)}{\partial
x^{2}}+\operatorname{Im}\psi(x,t)\frac{\partial^{2}\operatorname{Im}\psi(x,t)}{\partial
x^{2}}$ (F2)
$\displaystyle+i\operatorname{Re}\psi(x,t)\frac{\partial^{2}\operatorname{Im}\psi(x,t)}{\partial
x^{2}}-i\operatorname{Im}\psi(x,t)\frac{\partial^{2}\operatorname{Re}\psi(x,t)}{\partial
x^{2}})$ $\displaystyle=$
$\displaystyle\operatorname{Re}\psi(x,t)\frac{\partial^{2}\operatorname{Re}\psi(x,t)}{\partial
x^{2}}+\operatorname{Im}\psi(x,t)\frac{\partial^{2}\operatorname{Im}\psi(x,t)}{\partial
x^{2}}.$

To prove Eq. F1, we also relate the first and second derivatives of $n(x,t)$
to the wave function,

$\frac{\partial n(x)}{\partial
x}=2\operatorname{Re}\psi(x)\frac{\partial\operatorname{Re}\psi(x)}{\partial
x}+2\operatorname{Im}\psi(x)\frac{\partial\operatorname{Im}\psi(x)}{\partial
x},$ (F3)

$\frac{\partial^{2}n(x)}{\partial
x^{2}}=2\operatorname{Re}\psi(x)\frac{\partial^{2}\operatorname{Re}\psi(x)}{\partial
x^{2}}+2(\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})^{2}+2\operatorname{Im}\psi(x)\frac{\partial^{2}\operatorname{Im}\psi(x)}{\partial
x^{2}}+2(\frac{\partial\operatorname{Im}\psi(x)}{\partial x})^{2}.\\\ $ (F4)

Using Eq. D3 as well as the above Eqs.F3/F4, the right-hand side of Eq. F1 is

$\displaystyle RHS$
$\displaystyle=\operatorname{Re}\psi(x)\frac{\partial^{2}\operatorname{Re}\psi(x)}{\partial
x^{2}}+(\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})^{2}+\operatorname{Im}\psi(x)\frac{\partial^{2}\operatorname{Im}\psi(x)}{\partial
x^{2}}+(\frac{\partial\operatorname{Im}\psi(x)}{\partial x})^{2})$ (F5)
$\displaystyle-\frac{(\operatorname{Re}\psi(x)\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})^{2}+2(\operatorname{Re}\psi(x)\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})(\operatorname{Im}\psi(x)\frac{\partial\operatorname{Im}\psi(x)}{\partial
x})+(\operatorname{Im}\psi(x)\frac{\partial\operatorname{Im}\psi(x)}{\partial
x})^{2}}{n(x)}$
$\displaystyle-\frac{(\operatorname{Re}\psi(x)\frac{\partial\operatorname{Im}\psi(x)}{\partial
x})^{2}-2(\operatorname{Re}\psi(x)\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})(\operatorname{Im}\psi(x)\frac{\partial\operatorname{Im}\psi(x)}{\partial
x})+(\operatorname{Im}\psi(x)\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})^{2}}{n(x)}$
$\displaystyle=\operatorname{Re}\psi(x)\frac{\partial^{2}\operatorname{Re}\psi(x)}{\partial
x^{2}}+(\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})^{2}+\operatorname{Im}\psi(x)\frac{\partial^{2}\operatorname{Im}\psi(x)}{\partial
x^{2}}+(\frac{\partial\operatorname{Im}\psi(x)}{\partial x})^{2})$
$\displaystyle-\frac{(\operatorname{Re}\psi(x)\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})^{2}+(\operatorname{Im}\psi(x)\frac{\partial\operatorname{Im}\psi(x)}{\partial
x})^{2}+(\operatorname{Re}\psi(x)\frac{\partial\operatorname{Im}\psi(x)}{\partial
x})^{2}+(\operatorname{Im}\psi(x)\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})^{2}}{n(x)}$
$\displaystyle=\operatorname{Re}\psi(x)\frac{\partial^{2}\operatorname{Re}\psi(x)}{\partial
x^{2}}+(\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})^{2}+\operatorname{Im}\psi(x)\frac{\partial^{2}\operatorname{Im}\psi(x)}{\partial
x^{2}}+(\frac{\partial\operatorname{Im}\psi(x)}{\partial x})^{2})$
$\displaystyle-\frac{n(x)((\frac{\partial\operatorname{Re}\psi(x)}{\partial
x})^{2}+(\frac{\partial\operatorname{Im}\psi(x)}{\partial x})^{2})}{n(x)}$
$\displaystyle=\operatorname{Re}\psi(x)\frac{\partial^{2}\operatorname{Re}\psi(x)}{\partial
x^{2}}+\operatorname{Im}\psi(x)\frac{\partial^{2}\operatorname{Im}\psi(x)}{\partial
x^{2}}.$

This is indeed the left-hand side of Eq. F1 (see Eq.F2), and thus Eq. F1 is
proved.

From Eq. C7, the second-order time-derivative of the moments in a one-
dimensional system can be written as

$\displaystyle\frac{d^{2}\braket{x^{a}}}{dt^{2}}=$
$\displaystyle-\int\int\int(ax^{a-1}\partial_{x}V(x,t)n(x,t))d^{3}r$ (F6)
$\displaystyle-$ $\displaystyle\int\int\int
a(a-1)x^{a-2}\operatorname{Re}(\partial_{x}^{2}\psi\psi^{*})d^{3}r$
$\displaystyle+$ $\displaystyle\frac{1}{4}a(a-1)((a-2)(a-3)\braket{x^{a-4}}.$

By substituting the integrand in the second term with Eq. F1, Eq. F6 can be
expressed without explicit reference to the wave function as

$\displaystyle\frac{d^{2}\braket{x^{a}}}{dt^{2}}=-$ $\displaystyle
a\int{\frac{\partial v(x)}{\partial
x}n(x,t)x^{a-1}dx}+a(a-1)\int{(\frac{(\frac{\partial n(x,t)}{\partial
x})^{2}}{4n(x,t)}+\frac{L(x,t)^{2}}{n(x,t)})x^{a-2}dx}$ (F7) $\displaystyle-$
$\displaystyle\frac{a(a-1)(a-2)(a-3)}{4}\braket{x^{a-4}}.$

Extending this expression to the general case (three-dimensional system) is
tedious but straightforward,

$\frac{d^{2}\braket{x^{a}y^{b}z^{c}}(t)}{dt^{2}}=-\int\int\int(ax^{a-1}y^{b}z^{c}\frac{\partial
V(x,y,z,t)}{\partial x}n(x,y,z,t)\\\ +bx^{a}y^{b-1}z^{c}\frac{\partial
V(x,y,z,t)}{\partial y}n(x,y,z,t)+cx^{a}y^{b}z^{c-1}\frac{\partial
V(x,y,z,t)}{\partial z}n(x,y,z,t))d^{3}r\\\
+\int\int\int(a(a-1)x^{a-2}y^{b}z^{c}(\frac{(\frac{\partial
n(x,y,z,t)}{\partial
x}^{2})}{4n(x,y,z,t)}+\frac{L_{x}(x,y,z,t)^{2}}{n(x,y,z,t)})\\\
+2bax^{a-1}y^{b-1}z^{c}(\frac{(\frac{\partial n(x,y,z,t)}{\partial
x}\frac{\partial n(x,y,z,t)}{\partial
y})}{4n(x,y,z,t)}+\frac{L_{x}(x,y,z,t)L_{y}(x,y,z,t)}{n(x,y,z,t)})\\\
+2cax^{a-1}y^{b}z^{c-1}(\frac{(\frac{\partial n(x,y,z,t)}{\partial
x}\frac{\partial n(x,y,z,t)}{\partial
z})}{4n(x,y,z,t)}+\frac{L_{x}(x,y,z,t)L_{z}(x,y,z,t)}{n(x,y,z,t)})\\\
+b(b-1)y^{b-2}x^{a}z^{c}(\frac{(\frac{\partial n(x,y,z,t)}{\partial
y}^{2})}{4n(x,y,z,t)}+\frac{L_{y}(x,y,z,t)^{2}}{n(x,y,z,t)})\\\
+2cby^{b-1}x^{a}z^{c-1}(\frac{(\frac{\partial n(x,y,z,t)}{\partial
z}\frac{\partial n(x,y,z,t)}{\partial
y})}{4n(x,y,z,t)}+\frac{L_{z}(x,y,z,t)L_{y}(x,y,z,t)}{n(x,y,z,t)})\\\
+c(c-1)z^{c-2}x^{a}y^{b}(\frac{(\frac{\partial n(x,y,z,t)}{\partial
z}^{2})}{4n(x,y,z,t)}+\frac{L_{z}(x,y,z,t)^{2}}{n(x,y,z,t)}))d^{3}r\\\
-a(a-1)(\frac{(a-2)(a-3)}{4}\braket{x^{a-4}y^{b}z^{c}}+\frac{b(b-1)}{2}\braket{y^{b-2}x^{a-2}z^{c}}+\frac{c(c-1)}{2}\braket{z^{c-2}x^{a-2}y^{b}})\\\
-b(b-1)(\frac{(b-2)(b-3)}{4}\braket{y^{b-4}x^{a}z^{c}}+\frac{c(c-1)}{2}\braket{z^{c-2}x^{a}y^{b-2}})\\\
-\frac{c(c-1)(c-2)(c-3)}{4}\braket{z^{c-4}x^{a}y^{b}},\\\ $ (F8)

which is Eq. 18.

##  Appendix G Derivation of Eq. 15

Eq. 15 gives the equation of motion for the product term
$c_{u}^{abc}(t)\frac{\partial\braket{x^{a}y^{b}z^{c}}(t)}{\partial t}$
expressed in terms of the particle density and the moments without explicit
dependence on the wave function. We start by writing Eq. E3 as

$\displaystyle c_{u}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}}{dt}$
$\displaystyle=\int\int\int x^{a}y^{b}z^{c}(\alpha_{u}(x,y,z,t))d^{3}r,$ (G1)
$\displaystyle=\int\int\int
x^{a}y^{b}z^{c}(\frac{i}{2}(\partial_{u}^{2}\psi\psi^{*}-\psi\partial_{u}^{2}\psi^{*}))d^{3}r$
$\displaystyle=\operatorname{Re}\left[-i\int\int\int\partial_{u}(x^{a}y^{b}z^{c})\partial_{u}\psi\psi^{*}d^{3}r\right].$

For the rest of the derivations here in this subsection, the real part is
implied. Analogously to the previous sections, using
$n(\mathbf{r},t)=|\psi(\mathbf{r},t)|^{2}$ and TD-SE (Eq.2) yields the time-
derivative as

$\displaystyle\frac{d}{dt}(c_{u}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}(t)}{dt})=$
$\displaystyle-i\int\partial_{u}\left(x^{a}y^{b}z^{c}\right)\left(\partial_{u}\frac{\partial\psi}{\partial
t}\psi^{*}+\partial_{u}\psi\frac{\partial\psi^{*}}{\partial t}\right)d^{3}r$
(G2) $\displaystyle=$
$\displaystyle-i\int\partial_{u}\left(x^{a}y^{b}z^{c}\right)\left[\partial_{u}\left(\frac{i}{2}\nabla^{2}\psi-
iV\psi\right)\psi^{*}+\partial_{u}\psi\left(-\frac{i}{2}\nabla^{2}\psi^{*}+iV\psi^{*}\right)\right]d^{3}r$
$\displaystyle=$
$\displaystyle\int\partial_{u}\left(x^{a}y^{b}z^{c}\right)\left[-\left(\partial_{u}V\right)n+\frac{1}{2}\left(\nabla^{2}\partial_{u}\psi\psi^{*}-\partial_{u}\psi\nabla^{2}\psi^{*}\right)\right]d^{3}r$
$\displaystyle=$
$\displaystyle-\int\partial_{u}(x^{a}y^{b}z^{c})\partial_{u}Vnd^{3}r+\int\partial_{u}\left(x^{a}y^{b}z^{c}\right)\frac{1}{2}(\nabla^{2}\partial_{u}\psi\psi^{*}-\partial_{u}\psi\nabla^{2}\psi^{*})d^{3}r.$

Applying the integration by parts to the right-hand side, we have

$\displaystyle\frac{d}{dt}(c_{u}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}(t)}{dt})=$
$\displaystyle-\int\partial_{u}(x^{a}y^{b}z^{c})\partial_{u}Vnd^{3}r+\int\partial_{u}\left(x^{a}y^{b}z^{c}\right)\frac{1}{2}(\nabla^{2}\partial_{u}\psi\psi^{*})d^{3}r$
(G3)
$\displaystyle-\frac{1}{2}\int\nabla^{2}(\partial_{u}\left(x^{a}y^{b}z^{c}\right)\partial_{u}\psi)\psi^{*}d^{3}r$
$\displaystyle=$
$\displaystyle-\int\partial_{u}(x^{a}y^{b}z^{c})\partial_{u}Vnd^{3}r-\frac{1}{2}\int(\nabla^{2}\partial_{u}\left(x^{a}y^{b}z^{c}\right)\partial_{u}\psi$
$\displaystyle+2\nabla\partial_{u}\left(x^{a}y^{b}z^{c}\right)\cdot(\nabla\partial_{u}\psi))\psi^{*}d^{3}r.$

To make further progress on this equation above to the working form, let us
express write the following with the integration by parts

$\displaystyle\left(\int\nabla^{2}\partial_{u}\left(x^{a}y^{b}z^{c}\right)\partial_{u}\psi\psi^{*}d^{3}r\right)$
$\displaystyle=-\int\nabla^{2}\partial_{u}^{2}\left(x^{a}y^{b}z^{c}\right)nd^{3}r-\left(\int\nabla^{2}\partial_{u}\left(x^{a}y^{b}z^{c}\right)\psi\partial_{u}\psi^{*}d^{3}r\right)$
(G4)
$\displaystyle\left(\int\nabla^{2}\partial_{u}\left(x^{a}y^{b}z^{c}\right)(\partial_{u}\psi\psi^{*}+\psi\partial_{u}\psi^{*})d^{3}r\right)$
$\displaystyle=-\int\nabla^{2}\partial_{u}^{2}\left(x^{a}y^{b}z^{c}\right)nd^{3}r$
$\displaystyle\left(\int\nabla^{2}\partial_{u}\left(x^{a}y^{b}z^{c}\right)2\operatorname{Re}(\partial_{u}\psi\psi^{*})d^{3}r\right)$
$\displaystyle=-\int\nabla^{2}\partial_{u}^{2}\left(x^{a}y^{b}z^{c}\right)nd^{3}r$
$\displaystyle\operatorname{Re}\left(\int\nabla^{2}\partial_{u}\left(x^{a}y^{b}z^{c}\right)(\partial_{u}\psi\psi^{*})d^{3}r\right)$
$\displaystyle=-\frac{1}{2}\int\nabla^{2}\partial_{u}^{2}\left(x^{a}y^{b}z^{c}\right)nd^{3}r,$

where we used
$(\partial_{u}\psi\psi^{*}+\psi\partial_{u}\psi^{*})=2\operatorname{Re}(\partial_{u}\psi\psi^{*})$.
Using the final expression Eq. G4 for the second term of the RHS in Eq. G3,
Eq. G3 is expressed as

$\displaystyle\frac{d}{dt}(c_{u}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}(t)}{dt})=$
$\displaystyle\int\operatorname{Re}(-\partial_{u}(x^{a}y^{b}z^{c})\partial_{u}Vn+\frac{1}{4}\nabla^{2}\partial_{u}^{2}(x^{a}y^{b}z^{c})n$
(G5)
$\displaystyle-(\nabla\partial_{u}(x^{a}y^{b}z^{c})\cdot(\nabla\partial_{u}\psi))\psi^{*})d^{3}r.$

Using the same approach used to derive Eq. F1 in Sec.F, Eq. G3 in terms of the
density and moments is

$\frac{d}{dt}(c_{u}^{abc}(t)\frac{d\braket{x^{a}y^{b}z^{c}}(t)}{dt})=-\int\int\int(\partial_{u}(x^{a}y^{b}z^{c})\frac{\partial
V(x,y,z,t)}{\partial u}n(x,y,z,t))d^{3}r\\\
+\int\int\int(\partial_{u}\partial_{x}(x^{a}y^{b}z^{c})(\frac{(\frac{\partial
n(x,y,z,t)}{\partial u}\frac{\partial n(x,y,z,t)}{\partial
x})}{4n(x,y,z,t)}+\frac{L_{u}(x,y,z,t)L_{x}(x,y,z,t)}{n(x,y,z,t)})\\\
+\partial_{u}\partial_{y}(x^{a}y^{b}z^{c})(\frac{(\frac{\partial
n(x,y,z,t)}{\partial u}\frac{\partial n(x,y,z,t)}{\partial
y})}{4n(x,y,z,t)}+\frac{L_{u}(x,y,z,t)L_{y}(x,y,z,t)}{n(x,y,z,t)})\\\
+\partial_{u}\partial_{z}(x^{a}y^{b}z^{c})(\frac{(\frac{\partial
n(x,y,z,t)}{\partial u}\frac{\partial n(x,y,z,t)}{\partial
z})}{4n(x,y,z,t)}+\frac{L_{u}(x,y,z,t)L_{z}(x,y,z,t)}{n(x,y,z,t)}))d^{3}r\\\
-\frac{1}{4}\braket{\nabla^{2}\partial_{u}^{2}\left(x^{a}y^{b}z^{c}\right)}.\\\
$ (G6)

##  Appendix H Derivation of Eq. 29

For the one-dimensional system, we have

$\epsilon=\braket{\psi}{\hat{H}}{\psi}=\int-\frac{1}{2}\psi^{*}(x)\frac{\partial^{2}\psi(x)}{\partial
x^{2}}+\psi^{*}(x)V(x)\psi(x)dx.$ (H1)

Using Eq. F1, the wave function dependence can be recast in terms of the
probability density and the moments as

$\epsilon=\int V(x)n(x)dx+\frac{1}{2}\int(\frac{(\frac{\partial n(x)}{\partial
x}^{2})}{4n(x)}+\frac{L(x)^{2}}{n(x)})dx.$ (H2)

##  Appendix I Analytical solution for simple harmonic potential (Eqs. 25, 26,
and 27)

The harmonic potential case, $V(x,t)=x^{2}+c(t)x$, can be solved analytically
within the framework of the moment propagation theory (MPT). For this system,
the second-order Edgeworth basis (E2) is sufficient[61], and it requires only
the first-order and the second-order moments to be propagated.

###  I1  Analytical Solution

The MPT equation of motion for a 1D system is:

$\displaystyle\frac{d^{2}\braket{x^{a}}}{dt^{2}}=$
$\displaystyle-a\int{\frac{\partial V(x)}{\partial
x}n(x,t)x^{a-1}dx}+a(a-1)\int{(\frac{(\frac{\partial n(x,t)}{\partial
x})^{2}}{4n(x,t)}+\frac{L(x,t)^{2}}{n(x,t)})x^{a-2}dx}$ (I1)
$\displaystyle-\frac{a(a-1)(a-2)(a-3)}{4}\braket{x^{a-4}}.$

For the first ($a$=1) and second ($a$=2) order moments, the second time
derivative can be simplified as

$\frac{d^{2}\braket{x}}{dt^{2}}=-\int{\frac{\partial V(x)}{\partial
x}n(x,t)dx},$ (I2)

$\frac{d^{2}\braket{x^{2}}}{dt^{2}}=-2\int{\frac{\partial V(x)}{\partial
x}n(x,t)xdx}+2\int{(\frac{(\frac{\partial n(x,t)}{\partial
x})^{2}}{4n(x,t)}+\frac{L(x,t)^{2}}{n(x,t)})dx}.$ (I3)

With the second-order Edgeworth (E2) basis, the particle density has a
Gaussian form

$n(x,t)=\frac{1}{\sqrt{2\pi(\braket{x^{2}}(t)-(\braket{x}(t))^{2})}}\exp{(-\frac{1}{2}\frac{(x-\braket{x}(t))^{2}}{\braket{x^{2}}(t)-(\braket{x}(t))^{2}})}.$
(I4)

The spatial derivative is

$\frac{\partial n(x,t)}{\partial
x}=-n(x,t)(\frac{(x-\braket{x}(t))}{\braket{x^{2}}(t)-(\braket{x}(t))^{2}}),$
(I5)

and the time derivative is

$\displaystyle\frac{\partial n(x,t)}{\partial t}=$
$\displaystyle\frac{-\frac{\partial\braket{x^{2}}(t)}{\partial
t}+2\braket{x}(t)\frac{\partial\braket{x}(t)}{\partial
t}}{2\sqrt{2\pi(\braket{x^{2}}(t)-(\braket{x}(t))^{2})^{3}}}\exp{(\frac{-1}{2}\frac{(x-\braket{x}(t))^{2}}{\braket{x^{2}}(t)-(\braket{x}(t))^{2}})}$
(I6)
$\displaystyle-\frac{1}{2}n(x,t)(\frac{2(x-\braket{x}(t))(-\frac{\partial\braket{x}(t)}{\partial
t})}{(\braket{x^{2}}(t)-(\braket{x}(t))^{2})}-\frac{(\frac{\partial\braket{x^{2}}(t)}{\partial
t}-2\braket{x}(t)\frac{\partial\braket{x}(t)}{\partial
t})(x-\braket{x}(t))^{2}}{(\braket{x^{2}}(t)-(\braket{x}(t))^{2})^{2}}).$

Rewriting the time derivative as a polynomial, we have

$\frac{\partial n(x,t)}{\partial
t}=k_{1}n(x,t)+k_{2}xn(x,t)+k_{3}x^{2}n(x,t),$ (I7)

where $\\{k_{i}\\}$ are given by

$k_{1}=\frac{-\frac{\partial\braket{x^{2}}(t)}{\partial
t}+2\braket{x}(t)\frac{\partial\braket{x}(t)}{\partial
t}}{2(\braket{x^{2}}(t)-(\braket{x}(t))^{2})}+\frac{\braket{x}(t)(-\frac{\partial\braket{x}(t)}{\partial
t})}{(\braket{x^{2}}(t)-(\braket{x}(t))^{2})}\\\
+\frac{1}{2}\frac{(\frac{\partial\braket{x^{2}}(t)}{\partial
t}-2\braket{x}(t)\frac{\partial\braket{x}(t)}{\partial
t})(\braket{x}(t))^{2}}{(\braket{x^{2}}(t)-(\braket{x}(t))^{2})^{2}},\\\
k_{2}x=\frac{x(\frac{\partial\braket{x}(t)}{\partial
t})}{(\braket{x^{2}}(t)-(\braket{x}(t))^{2})}+\frac{1}{2}\frac{(\frac{\partial\braket{x^{2}}(t)}{\partial
t}-2\braket{x}(t)\frac{\partial\braket{x}(t)}{\partial
t})(-2x\braket{x}(t))}{(\braket{x^{2}}(t)-(\braket{x}(t))^{2})^{2}},\\\
k_{3}x^{2}=\frac{1}{2}\frac{(\frac{\partial\braket{x^{2}}(t)}{\partial
t}-2\braket{x}(t)\frac{\partial\braket{x}(t)}{\partial
t})(x^{2})}{(\braket{x^{2}}(t)-(\braket{x}(t))^{2})^{2}}.\\\ $ (I8)

Let us define $S$ as the variance (spread) as

$S\equiv\braket{x^{2}}-\braket{x}^{2}.$ (I9)

In order to simplify the time derivative using the variance, we use
$\braket{x}=0$ for this harmonic potential case without loss of generality due
to the translational invariance. Eq. I7 is then

$\displaystyle\frac{\partial n(x,t)}{\partial
t}=k_{1}n(x,t)+k_{2}xn(x,t)+k_{3}x^{2}n(x,t),$ (I10) $\displaystyle
k_{1}=\frac{-\dot{S}}{2S},k_{2}=\frac{\frac{\partial\braket{x}(t)}{\partial
t}}{S},k_{3}=\frac{\dot{S}}{2S^{2}}.$

Integrating the time derivative over the position x using the definition
$L(x,t)\equiv-\int_{-\infty}^{x}\alpha_{x}(x^{\prime},t)dx^{\prime}$, we
obtain

$L(x,t)=-(k_{1}+Sk_{3})\Phi(x,t)+Sk_{2}n(x,t)+Sxk_{3}n(x,t),$ (I11)

where $\Phi(x,t)\equiv\int_{-\infty}^{x}n(x^{\prime},t)dx^{\prime}$. Because
$(k_{1}+Sk_{3})=\frac{-\dot{S}}{2S}+\frac{\dot{S}}{2S}=0$, we have

$\displaystyle L(x,t)$ $\displaystyle=a_{1}n(x,t)+a_{2}xn(x,t),$ (I12)
$\displaystyle a_{1}$ $\displaystyle=\frac{\partial\braket{x}(t)}{\partial
t},a_{2}=\frac{\dot{S}}{2S}.$

Substituting Eqs. I12 and I5 into the equation of motion (EOM) for the first-
order and second-order moments (Eqs. I2 and I3), we have

$\displaystyle\frac{d^{2}\braket{x}}{dt^{2}}$
$\displaystyle=-\int{\frac{\partial(x^{2}+c(t)x)}{\partial
x}n(x,t)dx}=-2\braket{x}(t)-c(t),$ (I13)

$\displaystyle\frac{d^{2}\braket{x^{2}}}{dt^{2}}=$
$\displaystyle-2\int{\frac{\partial(x^{2}+c(t)x)}{\partial x}n(x,t)xdx}$ (I14)
$\displaystyle+2\int{(\frac{1}{4S^{2}}n(x,t)x^{2}+a_{1}^{2}n(x,t)+2a_{1}a_{2}xn(x,t)+a_{2}^{2}x^{2}n(x,t))dx}$
$\displaystyle=$
$\displaystyle-2\int(2x+c(t))n(x,t)x+2((\frac{\partial\braket{x}(t)}{\partial
t})^{2}n(x,t)+2\frac{\partial\braket{x}(t)}{\partial
t}\frac{\dot{S}}{2S}xn(x,t)$
$\displaystyle+((\frac{\dot{S}}{2S})^{2}x^{2}+\frac{1}{4S^{2}})n(x,t))dx.$
$\displaystyle=-4\braket{x^{2}}(t)-2c(t)\braket{x}(t)+2((\frac{\partial\braket{x}(t)}{\partial
t})^{2}+\frac{\dot{S}^{2}}{4S}+\frac{1}{4S})).$

Succinctly put, the second-order time derivatives of the moments read

$\frac{d^{2}\braket{x}}{dt^{2}}=-2\braket{x}-c,$ (I15)

$\frac{d^{2}\braket{x^{2}}}{dt^{2}}=-4\braket{x^{2}}-2c\braket{x}+2((\frac{(\frac{\partial
S}{\partial t})^{2}}{4S})+\frac{1}{4S}+(\frac{\partial\braket{x}(t)}{\partial
t})^{2}).$ (I16)

The single-particle energy is $\epsilon=\braket{\hat{H}}$, and it is
equivalent to the system’s total energy for this one-dimensional harmonic
oscillator model. The expression for the energy is obtained by using Eqs. I12
and I5 in Eq. H2, resulting in

$E=\epsilon=\braket{x^{2}}+c\braket{x}+\frac{1}{2}((\frac{(\frac{\partial
S}{\partial t})^{2}}{4S})+\frac{1}{4S}+(\frac{\partial\braket{x}(t)}{\partial
t})^{2}).$ (I17)

###  I2  Analytical Solution via EPES approach

EPES expression (Eq. 22) holds if/when the Hamiltonian expectation value is a
constant of motion in the quantum dynamics given by a suitable single-particle
Hamiltonian. Then, EPES approach is easier to work with by deriving the
second-order time-derivatives of the cumulants rather than the moments
directly. Note that the first cumulant is equivalent to the first-order
moment. We first relate the second cumulant (which is the variance) and their
time-derivatives to the moments by expressing them as

$\displaystyle S$ $\displaystyle\equiv\braket{x^{2}}-\braket{x}^{2},$ (I18)
$\displaystyle\dot{S}$
$\displaystyle=\dot{\braket{x^{2}}}-2\braket{x}\dot{\braket{x}},$
$\displaystyle\ddot{S}$
$\displaystyle=\ddot{\braket{x^{2}}}-2\dot{\braket{x}}^{2}-2\braket{x}\ddot{\braket{x}}.$

Using Eq. I18 and Eq. I17, the EPES equation of motion (Eq. 22) for the first-
order moment/cumulant can be written as

$\displaystyle\ddot{\braket{x}}$
$\displaystyle=-\frac{\partial\epsilon(x_{i},\dot{x}_{i})}{\partial\braket{x}}\left(\frac{\partial\epsilon(x_{i},\dot{x}_{i})}{\partial\dot{\braket{x}}}\right)^{-1}\dot{\braket{x}}$
(I19) $\displaystyle=-2\braket{x}-c.$

Rewriting Eq. 22 in terms of the cumulant via Eq. I18, the second-order time
derivative of the second cumulant can be rearrange as

$\ddot{S}=-\frac{\partial\epsilon(x_{i},\dot{x_{i}})}{\partial
S}\left(\frac{\partial\epsilon(x_{i},\dot{x_{i}})}{\partial\dot{S}}\right)^{-1}\dot{S}.$
(I20)

Then, substituting Eq. I17 in this expression yields

$\displaystyle\ddot{S}$
$\displaystyle=-(1-\frac{1}{2}(\frac{1+\dot{S}^{2}}{4S^{2}}))\frac{4S}{\dot{S}}\dot{S}$
(I21) $\displaystyle=-4S+\frac{1}{2}(\frac{1+\dot{S}^{2}}{S}).$

Converting this expression given in cumulants back to the one in moments by
using Eq. I18, we have

$\ddot{\braket{x^{2}}}=-4\braket{x^{2}}-2c\braket{x}+2(\frac{1+\dot{S}^{2}}{4S}+\dot{\braket{x}}^{2}),\\\
$ (I22)

which is the same as the one derived using the MPT without relying on the EPES
expression (see Section I1).

## References

  * Li  _et al._ [2020] X. Li, N. Govind, C. Isborn, A. E. I. DePrince, and K. Lopata, “Real-time time-dependent electronic structure theory,” Chemical Reviews 120, 9951–9993 (2020), pMID: 32813506, https://doi.org/10.1021/acs.chemrev.0c00223 .
  * Nascimento and DePrince [2016] D. R. Nascimento and A. E. I. DePrince, “Linear absorption spectra from explicitly time-dependent equation-of-motion coupled-cluster theory,” Journal of Chemical Theory and Computation 12, 5834–5840 (2016), pMID: 27779862, https://doi.org/10.1021/acs.jctc.6b00796 .
  * Koulias  _et al._ [2019] L. N. Koulias, D. B. Williams-Young, D. R. Nascimento, A. E. I. DePrince, and X. Li, “Relativistic real-time time-dependent equation-of-motion coupled-cluster,” Journal of Chemical Theory and Computation 15, 6617–6624 (2019), pMID: 31618584, https://doi.org/10.1021/acs.jctc.9b00729 .
  * Yabana and Bertsch [1996] K. Yabana and G. F. Bertsch, “Time-dependent local-density approximation in real time,” Phys. Rev. B 54, 4484–4487 (1996).
  * Yabana and Bertsch [1999] K. Yabana and G. F. Bertsch, “Time-dependent local-density approximation in real time: Application to conjugated molecules,” International Journal of Quantum Chemistry 75, 55–66 (1999).
  * Theilhaber [1992] J. Theilhaber, “Ab initio simulations of sodium using time-dependent density-functional theory,” Phys. Rev. B 46, 12990–13003 (1992).
  * Runge and Gross [1984] E. Runge and E. K. U. Gross, “Density-functional theory for time-dependent systems,” Phys. Rev. Lett. 52, 997–1000 (1984).
  * Shepard  _et al._ [2021] C. Shepard, R. Zhou, D. C. Yost, Y. Yao, and Y. Kanai, “Simulating electronic excitation and dynamics with real-time propagation approach to tddft within plane-wave pseudopotential formulation,” The Journal of Chemical Physics 155, 100901 (2021), https://doi.org/10.1063/5.0057587 .
  * Kononov  _et al._ [2022] A. Kononov, C.-W. Lee, T. P. dos Santos, B. Robinson, Y. Yao, Y. Yao, X. Andrade, A. D. Baczewski, E. Constantinescu, A. A. Correa, Y. Kanai, N. Modine, and A. Schleife, “Electron dynamics in extended systems within real-time time-dependent density-functional theory,” MRS Communications 12, 1002–1014 (2022).
  * Hekele  _et al._ [2021] J. Hekele, Y. Yao, Y. Kanai, V. Blum, and P. Kratzer, “All-electron real-time and imaginary-time time-dependent density functional theory within a numeric atom-centered basis function framework,” The Journal of Chemical Physics 155 (2021), 10.1063/5.0066753, 154801, https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0066753/14728668/154801_1_online.pdf .
  * Pemmaraju  _et al._ [2018] C. Pemmaraju, F. Vila, J. Kas, S. Sato, J. Rehr, K. Yabana, and D. Prendergast, “Velocity-gauge real-time tddft within a numerical atomic orbital basis set,” Computer Physics Communications 226, 30–38 (2018).
  * Sun  _et al._ [2021] J. Sun, C.-W. Lee, A. Kononov, A. Schleife, and C. A. Ullrich, “Real-time exciton dynamics with time-dependent density-functional theory,” Phys. Rev. Lett. 127, 077401 (2021).
  * Shepard, Yost, and Kanai [2023] C. Shepard, D. C. Yost, and Y. Kanai, “Electronic excitation response of dna to high-energy proton radiation in water,” Phys. Rev. Lett. 130, 118401 (2023).
  * Yao, Yost, and Kanai [2019] Y. Yao, D. C. Yost, and Y. Kanai, “$k$-shell core-electron excitations in electronic stopping of protons in water from first principles,” Phys. Rev. Lett. 123, 066401 (2019).
  * Yost, Yao, and Kanai [2020] D. C. Yost, Y. Yao, and Y. Kanai, “First-principles modeling of electronic stopping in complex matter under ion irradiation,” The Journal of Physical Chemistry Letters 11, 229–237 (2020), pMID: 31829604, https://doi.org/10.1021/acs.jpclett.9b02975 .
  * Petrone  _et al._ [2014] A. Petrone, D. B. Lingerfelt, N. Rega, and X. Li, “From charge-transfer to a charge-separated state: a perspective from the real-time tddft excitonic dynamics,” Phys. Chem. Chem. Phys. 16, 24457–24465 (2014).
  * Kümmel [2017] S. Kümmel, “Charge-transfer excitations: A challenge for time-dependent density functional theory that has been met,” Advanced Energy Materials 7, 1700440 (2017), https://onlinelibrary.wiley.com/doi/pdf/10.1002/aenm.201700440 .
  * Ratchford [2019] D. C. Ratchford, “Plasmon-induced charge transfer: Challenges and outlook,” ACS Nano 13, 13610–13614 (2019), pMID: 31809010, https://doi.org/10.1021/acsnano.9b08829 .
  * You  _et al._ [2021] P. You, C. Lian, D. Chen, J. Xu, C. Zhang, S. Meng, and E. Wang, “Nonadiabatic dynamics of photocatalytic water splitting on a polymeric semiconductor,” Nano Letters 21, 6449–6455 (2021), pMID: 34279962, https://doi.org/10.1021/acs.nanolett.1c01187 .
  * Miyamoto  _et al._ [2019] Y. Miyamoto, H. Zhang, X. Cheng, and A. Rubio, “Ab initio simulation of laser-induced water decomposition close to carbon nanotubes,” Phys. Rev. B 99, 165424 (2019).
  * Miyamoto  _et al._ [2017] Y. Miyamoto, H. Zhang, X. Cheng, and A. Rubio, “Modeling of laser-pulse induced water decomposition on two-dimensional materials by simulations based on time-dependent density functional theory,” Phys. Rev. B 96, 115451 (2017).
  * Makkonen  _et al._ [2021] E. Makkonen, T. P. Rossi, A. H. Larsen, O. Lopez-Acevedo, P. Rinke, M. Kuisma, and X. Chen, “Real-time time-dependent density functional theory implementation of electronic circular dichroism applied to nanoscale metal–organic clusters,” The Journal of Chemical Physics 154 (2021), 10.1063/5.0038904, 114102, https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0038904/9695019/114102_1_online.pdf .
  * Konecny  _et al._ [2019] L. Konecny, M. Repisky, K. Ruud, and S. Komorovsky, “Relativistic four-component linear damped response TDDFT for electronic absorption and circular dichroism calculations,” The Journal of Chemical Physics 151 (2019), 10.1063/1.5128564, 194112, https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5128564/15566486/194112_1_online.pdf .
  * Goings and Li [2016] J. J. Goings and X. Li, “An atomic orbital based real-time time-dependent density functional theory for computing electronic circular dichroism band spectra,” The Journal of Chemical Physics 144, 234102 (2016), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.4953668/15513453/234102_1_online.pdf .
  * Sun  _et al._ [2019] S. Sun, R. A. Beck, D. Williams-Young, and X. Li, “Simulating magnetic circular dichroism spectra with real-time time-dependent density functional theory in gauge including atomic orbitals,” Journal of Chemical Theory and Computation 15, 6824–6831 (2019), pMID: 31577442, https://doi.org/10.1021/acs.jctc.9b00632 .
  * Krumland  _et al._ [2020] J. Krumland, A. M. Valencia, S. Pittalis, C. A. Rozzi, and C. Cocchi, “Understanding real-time time-dependent density-functional theory simulations of ultrafast laser-induced dynamics in organic molecules,” The Journal of Chemical Physics 153 (2020), 10.1063/5.0008194, 054106, https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0008194/15577013/054106_1_online.pdf .
  * Rossi  _et al._ [2017] T. P. Rossi, M. Kuisma, M. J. Puska, R. M. Nieminen, and P. Erhart, “Kohn–sham decomposition in real-time time-dependent density-functional theory: An efficient tool for analyzing plasmonic excitations,” Journal of Chemical Theory and Computation 13, 4779–4790 (2017), pMID: 28862851, https://doi.org/10.1021/acs.jctc.7b00589 .
  * Jia and Lin [2019] W. Jia and L. Lin, “Fast real-time time-dependent hybrid functional calculations with the parallel transport gauge and the adaptively compressed exchange formulation,” Computer Physics Communications 240, 21–29 (2019).
  * Yost, Yao, and Kanai [2019] D. C. Yost, Y. Yao, and Y. Kanai, “Propagation of maximally localized Wannier functions in real-time TDDFT,” The Journal of Chemical Physics 150 (2019), 194113.
  * Gygi, Fattebert, and Schwegler [2003] F. Gygi, J.-L. Fattebert, and E. Schwegler, “Computation of maximally localized wannier functions using a simultaneous diagonalization algorithm,” Computer Physics Communications 155, 1–6 (2003).
  * Marzari  _et al._ [2012] N. Marzari, A. A. Mostofi, J. R. Yates, I. Souza, and D. Vanderbilt, “Maximally localized wannier functions: Theory and applications,” Rev. Mod. Phys. 84, 1419–1475 (2012).
  * Zhou, Yost, and Kanai [2021] R. Zhou, D. C. Yost, and Y. Kanai, “First-principles demonstration of nonadiabatic thouless pumping of electrons in a molecular system,” The Journal of Physical Chemistry Letters 12, 4496–4503 (2021), pMID: 33956458, https://doi.org/10.1021/acs.jpclett.1c01037 .
  * Zhou and Kanai [2023] R. Zhou and Y. Kanai, “Molecular control of floquet topological phase in non-adiabatic thouless pumping,” The Journal of Physical Chemistry Letters 14, 8205–8212 (2023), pMID: 37672485, https://doi.org/10.1021/acs.jpclett.3c01746 .
  * Carrasquilla [2020] J. Carrasquilla, “Machine learning for quantum matter,” Advances in Physics: X 5, 1797528 (2020), https://doi.org/10.1080/23746149.2020.1797528 .
  * Lin  _et al._ [2022] K. Lin, J. Peng, C. Xu, F. L. Gu, and Z. Lan, “Automatic evolution of machine-learning-based quantum dynamics with uncertainty analysis,” Journal of Chemical Theory and Computation 18, 5837–5855 (2022), pMID: 36184823, https://doi.org/10.1021/acs.jctc.2c00702 .
  * Luo  _et al._ [2021] D. Luo, Z. Chen, J. Carrasquilla, and B. K. Clark, “Autoregressive neural network for simulating open quantum systems via a probabilistic formulation,” (2021), arXiv:2009.05580 [cond-mat.str-el] .
  * Mohseni  _et al._ [2022] N. Mohseni, T. Fösel, L. Guo, C. Navarrete-Benlloch, and F. Marquardt, “Deep Learning of Quantum Many-Body Dynamics via Random Driving,” Quantum 6, 714 (2022).
  * Rodrí guez  _et al._ [2022] L. E. H. Rodrí guez, A. Ullah, K. J. R. Espinosa, P. O. Dral, and A. A. Kananenka, “A comparative study of different machine learning methods for dissipative quantum dynamics,” Machine Learning: Science and Technology 3, 045016 (2022).
  * Häse  _et al._ [2016] F. Häse, S. Valleau, E. Pyzer-Knapp, and A. Aspuru-Guzik, “Machine learning exciton dynamics,” Chemical Science 7, 5139–5147 (2016).
  * Ju  _et al._ [2021] C.-W. Ju, H. Bai, B. Li, and R. Liu, “Machine learning enables highly accurate predictions of photophysical properties of organic fluorescent materials: Emission wavelengths and quantum yields,” Journal of Chemical Information and Modeling 61, 1053–1065 (2021), pMID: 33620207, https://doi.org/10.1021/acs.jcim.0c01203 .
  * Chu, Saidi, and Prezhdo [2020] W. Chu, W. A. Saidi, and O. V. Prezhdo, “Long-lived hot electron in a metallic particle for plasmonics and catalysis: Ab initio nonadiabatic molecular dynamics with machine learning,” ACS Nano 14, 10608–10615 (2020), pMID: 32806073, https://doi.org/10.1021/acsnano.0c04736 .
  * Bogojeski  _et al._ [2020] M. Bogojeski, L. Vogt-Maranto, M. E. Tuckerman, K.-R. Müller, and K. Burke, “Quantum chemical accuracy from density functional approximations via machine learning,” Nature Communications 11 (2020), 10.1038/s41467-020-19093-1.
  * Westermayr  _et al._ [2021] J. Westermayr, M. Gastegger, K. T. Schütt, and R. J. Maurer, “Perspective on integrating machine learning into computational chemistry and materials science,” The Journal of Chemical Physics 154 (2021), 10.1063/5.0047760, 230903, https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0047760/15977579/230903_1_online.pdf .
  * Westermayr and Marquetand [2020] J. Westermayr and P. Marquetand, “Machine learning and excited-state molecular dynamics,” Machine Learning: Science and Technology 1, 043001 (2020).
  * Doerr  _et al._ [2021] S. Doerr, M. Majewski, A. Pérez, A. Krämer, C. Clementi, F. Noe, T. Giorgino, and G. De Fabritiis, “Torchmd: A deep learning framework for molecular simulations,” Journal of Chemical Theory and Computation 17, 2355–2363 (2021), pMID: 33729795, https://doi.org/10.1021/acs.jctc.0c01343 .
  * Secor, Soudackov, and Hammes-Schiffer [2021a] M. Secor, A. V. Soudackov, and S. Hammes-Schiffer, “Artificial neural networks as mappings between proton potentials, wave functions, densities, and energy levels,” The Journal of Physical Chemistry Letters 12, 2206–2212 (2021a).
  * Secor, Soudackov, and Hammes-Schiffer [2021b] M. Secor, A. V. Soudackov, and S. Hammes-Schiffer, “Artificial neural networks as propagators in quantum dynamics,” The Journal of Physical Chemistry Letters 12, 10654–10662 (2021b).
  * Yao  _et al._ [2022] Y. Yao, C. Cao, S. Haas, M. Agarwal, D. Khanna, and M. Abram, “Emulating quantum dynamics with neural networks via knowledge distillation,” (2022), arXiv:2203.10200 [quant-ph] .
  * Secor, Soudackov, and Hammes-Schiffer [2021c] M. Secor, A. V. Soudackov, and S. Hammes-Schiffer, “Artificial neural networks as propagators in quantum dynamics,” The Journal of Physical Chemistry Letters 12, 10654–10662 (2021c).
  * Prezhdo and Pereverzev [2000] O. V. Prezhdo and Y. V. Pereverzev, “Quantized Hamilton dynamics,” The Journal of Chemical Physics 113, 6557–6565 (2000), https://pubs.aip.org/aip/jcp/article-pdf/113/16/6557/10827502/6557_1_online.pdf .
  * Shigeta, Miyachi, and Hirao [2006] Y. Shigeta, H. Miyachi, and K. Hirao, “Quantal cumulant dynamics: General theory,” The Journal of Chemical Physics 125, 244102 (2006), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.2404677/15392483/244102_1_online.pdf .
  * Shen and Wang [2018] Y. Shen and L. Wang, “Semiclassical Moyal dynamics,” The Journal of Chemical Physics 149, 244113 (2018), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.5067005/15552634/244113_1_online.pdf .
  * Hall [2013] P. Hall,  _The bootstrap and Edgeworth expansion_ (Springer Science and Business Media, 2013).
  * Smith [1995] P. J. Smith, “A recursive formulation of the old problem of obtaining moments from cumulants and vice versa,” The American Statistician 49, 217–218 (1995).
  * Blinnikov and Moessner [1998] S. Blinnikov and R. Moessner, “Expansions for nearly gaussian distributions,” Astronomy and Astrophysics Supplement Series 130, 193–205 (1998).
  * Wallace [1958] D. L. Wallace, “Asymptotic approximations to distributions,” The Annals of Mathematical Statistics 29, 635–654 (1958).
  * Patarroyo [2019] K. Y. Patarroyo, “A digression on hermite polynomials,”  (2019).
  * Coffey, Wyatt, and Schieve [2011] T. M. Coffey, R. E. Wyatt, and W. C. Schieve, “Reconstruction of the time-dependent wave function exclusively from position data,” Phys. Rev. Lett. 107, 230403 (2011).
  * Ando [2018] K. Ando, “Potential energy surfaces for electron dynamics modeled by floating and breathing gaussian wave packets with valence-bond spin-coupling: An analysis of high-harmonic generation spectrum,” The Journal of Chemical Physics 148, 094305 (2018).
  * Shigeta  _et al._ [2015] Y. Shigeta, R. Harada, M. Kayanuma, and M. Shoji, “Quantal cumulant dynamics for real-time simulations of quantum many-body systems,” International Journal of Quantum Chemistry 115, 300–308 (2015), https://onlinelibrary.wiley.com/doi/pdf/10.1002/qua.24820 .
  * Andrews [2016] M. Andrews, “The evolution of oscillator wave functions,” American Journal of Physics 84, 270–278 (2016).
  * Newmark [1959] N. M. Newmark, “A method of computation for structural dynamics,” Journal of the Engineering Mechanics Division 85, 67–94 (1959), https://ascelibrary.org/doi/pdf/10.1061/JMCEA3.0000098 .
  * Crank and Nicolson [1947] J. Crank and P. Nicolson, “A practical method for numerical evaluation of solutions of partial differential equations of the heat-conduction type,” Mathematical Proceedings of the Cambridge Philosophical Society 43, 50–67 (1947).
  * Yabana  _et al._ [2006] K. Yabana, T. Nakatsukasa, J.-I. Iwata, and G. F. Bertsch, “Real-time, real-space implementation of the linear response time-dependent density-functional theory,” physica status solidi (b) 243, 1121–1138 (2006), https://onlinelibrary.wiley.com/doi/pdf/10.1002/pssb.200642005 .
  * Zangwill and Soven [1980] A. Zangwill and P. Soven, “Density-functional approach to local-field effects in finite systems: Photoabsorption in the rare gases,” Phys. Rev. A 21, 1561–1572 (1980).
  * Chen  _et al._ [2020] M. S. Chen, T. J. Zuehlsdorff, T. Morawietz, C. M. Isborn, and T. E. Markland, “Exploiting machine learning to efficiently predict multidimensional optical spectra in complex environments,” The Journal of Physical Chemistry Letters 11, 7559–7568 (2020), pMID: 32808797, https://doi.org/10.1021/acs.jpclett.0c02168 .
  * Unke  _et al._ [2021] O. T. Unke, S. Chmiela, H. E. Sauceda, M. Gastegger, I. Poltavsky, K. T. Schütt, A. Tkatchenko, and K.-R. Müller, “Machine learning force fields,” Chemical Reviews 121, 10142–10186 (2021), pMID: 33705118, https://doi.org/10.1021/acs.chemrev.0c01111 .
  * Botu  _et al._ [2017] V. Botu, R. Batra, J. Chapman, and R. Ramprasad, “Machine learning force fields: Construction, validation, and outlook,” The Journal of Physical Chemistry C 121, 511–522 (2017), https://doi.org/10.1021/acs.jpcc.6b10908 .
  * Zhang  _et al._ [2018] L. Zhang, J. Han, H. Wang, R. Car, and W. E, “Deep potential molecular dynamics: A scalable model with the accuracy of quantum mechanics,” Phys. Rev. Lett. 120, 143001 (2018).
  * Yao and Kanai [2020] Y. Yao and Y. Kanai, “Temperature dependence of nuclear quantum effects on liquid water via artificial neural network model based on SCAN meta-GGA functional,” The Journal of Chemical Physics 153, 044114 (2020), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0012815/14721884/044114_1_online.pdf .
  * Yao and Kanai [2021] Y. Yao and Y. Kanai, “Nuclear quantum effect and its temperature dependence in liquid water from random phase approximation via artificial neural network,” The Journal of Physical Chemistry Letters 12, 6354–6362 (2021), pMID: 34231366, https://doi.org/10.1021/acs.jpclett.1c01566 .
  * Han  _et al._ [2018] J. Han, L. Zhang, R. Car, and W. E, “Deep potential: A general representation of a many-body potential energy surface,” Communications in Computational Physics 23 (2018), 10.4208/cicp.oa-2017-0213.
  * Bonati, Zhang, and Parrinello [2019] L. Bonati, Y.-Y. Zhang, and M. Parrinello, “Neural networks-based variationally enhanced sampling,” Proceedings of the National Academy of Sciences 116, 17641–17647 (2019), https://www.pnas.org/doi/pdf/10.1073/pnas.1907975116 .
  * Behler and Parrinello [2007] J. Behler and M. Parrinello, “Generalized neural-network representation of high-dimensional potential-energy surfaces,” Phys. Rev. Lett. 98, 146401 (2007).
  * Chollet  _et al._ [2015] F. Chollet _et al._ , “Keras,” https://keras.io (2015).