1

## Atomistic simulations of out-of-equilibrium quantum nuclear dynamics

Francesco Libbi, 1, ∗ Anders Johansson, 1 Lorenzo Monacelli, 2, ∗ and Boris Kozinsky 1, 3, ∗

John A. Paulson School of Engineering and Applied Sciences, Harvard University, Cambridge, MA 02138, USA 2 Department of Physics, Sapienza University of Rome, Piazzale Aldo Moro 5, 00185 Rome, Italy 3 Robert Bosch LLC Research and Technology Center, Watertown, MA 02472, USA

The rapid advancements in ultrafast laser technology have paved the way for pumping and probing the out-of-equilibrium dynamics of nuclei in crystals. However, interpreting these experiments is extremely challenging due to the complex nonlinear responses in systems where lattice excitations interact, particularly in crystals composed of light atoms or at low temperatures where the quantum nature of ions becomes significant. In this work, we address the nonequilibrium quantum ionic dynamics from first principles. Our approach is general and can be applied to simulate any crystal, in combination with a first-principles treatment of electrons or external machine-learning potentials. It is implemented by leveraging the nonequilibrium time-dependent self-consistent harmonic approximation (TD-SCHA), with a stable, energy-conserving, correlated stochastic integration scheme that achieves an accuracy of O ( dt 3 ). We benchmark the method with both a simple one-dimensional model to test its accuracy and a realistic 40-atom cell of SrTiO 3 under THz laser pump, paving the way for simulations of ultrafast THz-Xray pump-probe spectroscopy like those performed in synchrotron facilities.

## I. INTRODUCTION

still producing predictions in very good agreement with experiments 2,3,7,25,26 .

Nuclear quantum effects often play a crucial role in determining properties of materials 1 , affecting their thermodynamic stability 2-4 , electronic structure 5,6 , and transport phenomena . 7 Understanding and accurately accounting for quantum contributions in tunneling and vibrational statistics requires simulations that go beyond classical nuclei approximation. Path-integral molecular dynamics (PIMD) is the most common approach to simulate nuclear quantum effects in complex anharmonic crystals 8-11 , especially with the ring-polymer formulation 12,13 , which is an exact framework for sampling the equilibrium nuclear density matrix. However, PIMD is rigorously formulated with an assumption of thermodynamic equilibrium as it evolves the trajectories in imaginary time. Simulating out-of-equilibrium quantum nuclear dynamics requires a theory for the realtime evolution, such as the path-integral quantum Monte Carlo (PIMC), 14-17 , which is extremely challenging due to the appearance of the so-called sign problem, in which different trajectories contribute with different signs, introducing noise in determining dynamical observables 18 . Therefore, most applications of real-time PIMC remain limited to systems composed of only a few degrees of freedom 19 .

Several techniques have been proposed to overcome the limits of real-time PIMC to study nuclear quantum effects in systems of realistic interest. Among these, the time-dependent self-consistent harmonic approximation (TD-SCHA 20,21 ) holds significant promise. The theory extends the stochastic self-consistent harmonic approximation (SSCHA) 22-24 , a well-established technique to simulate equilibrium thermodynamics of solids accounting for quantum nuclear fluctuations. The success of SSCHA lies in the adoption of approximations that are particularly effective for crystals, achieving a computational cost that is orders of magnitude lower than PIMD, while

TD-SCHA has already been employed in the linear response regime, where it enabled the prediction of Raman and IR spectra of metallic hydrogen with unprecedented accuracy 27 . However, up to now, no application of TD-SCHA beyond equilibrium has been attempted. The difficulty in obtaining accurate and stable dynamic solutions of TD-SCHA equations is challenging due to the need for the evaluation of ensemble averages of the nuclear potential energy landscape.

In this work we address this challenge by developing an algorithm to solve the dynamical TD-SCHA equations and simulate the out-of-equilibrium dynamics of nuclei in complex realistic systems. In Sec. II, we revise the TD-SCHA equations of motion for the nuclear density matrix. Sec. III introduces three different numerical algorithms to integrate the TD-SCHA equations, and their stability is discussed in Sec. IV. Sec. V presents the correlated sampling technique to perform the evaluation of the stochastic quantum averages. This formulation ensures both the efficient evaluation of ensemble averages and the numerical stability of the equations. Crucially, we show that such a correlated approach conserves energy in one dimensional problems independently of the number of stochastic configurations adopted. We benchmark the method in Sec. VI, where we compare the different numerical schemes on a one dimensional model system. Finally, in Sec. VII, we provide an example of the application of TD-SCHA to realistic systems by studying the quantum dynamics in SrTiO 3 (STO) when driven out of equilibrium by a strong laser pulse of THz frequency.

## II. TIME-DEPENDENT SELF-CONSISTENT HARMONIC APPROXIMATION

The TD-SCHA formulation leverages the Wigner formalism 28,29 . The Wigner transform of the nuclear quantum density matrix ˆ( ) ρ t describing the quantum state is defined as

$$\rho ( { \mathbf R }, { \mathbf P }, t ) = \int \frac { e ^ { - \frac { i } { h } { \mathbf P } \cdot { \mathbf R } ^ { \prime } } } { ( 2 \pi \hbar { ) } ^ { 3 N } } \langle { \mathbf R } + \frac { { \mathbf R } ^ { \prime } } { 2 } \Big | \hat { \rho } ( t ) \Big | { \mathbf R } - \frac { { \mathbf R } ^ { \prime } } { 2 } \rangle d { \mathbf R } ^ { \prime } \, \ ( 1 )$$

and it maps the quantum operator ˆ( ) ρ t into a function of positions and momenta ρ ( R P , , t ) that is analogous to the classical nuclear density. All atomic quantities are rescaled by mass to simplify the notation: R i = ˜ R i √ m i and P i = ˜ P / i √ m i , where the index i goes from 1 to 3 N , containing both the Cartesian and atomic index, and the tilde indicates the standard (not mass-rescaled) quantities. Dynamical averages of any quantum observable are obtained by tracing R and P on the density:

$$\langle O \rangle \left ( t \right ) = \int d { \mathbf R } d { \mathbf P } \rho ( { \mathbf R }, { \mathbf P }, t ) O ( { \mathbf R }, { \mathbf P } ). \quad ( 2 ) \quad \text{ for $\alpha$}$$

The TD-SCHA method is based on expressing the Wigner density matrix as a general Gaussian form in terms of ionic positions R and momenta P :

$$\rho ( { \mathbf R }, { \mathbf P }, t ) = \frac { 1 } { \mathcal { N } } e ^ { - \frac { 1 } { 2 } ( \delta { \mathbf R } \cdot \alpha \cdot \delta { \mathbf R } + \delta { \mathbf P } \cdot \beta \cdot \delta { \mathbf P } + \delta { \mathbf R } \cdot \gamma \cdot \delta { \mathbf P } ) } \. \quad ( 3 ) \quad \ \frac { \langle } { \mathcal { T } \cdot \cdot \cdot \cdot }$$

Here δ R ( ) t = R -R ( ) t and δ P ( ) t = P -P ( ) t and N is the normalization factor, where R ( ) t = ⟨ R ⟩ ( ) t represent average positions and P ( ) t = ⟨ P ⟩ ( ) t the average momenta. The α ( ), t β ( ) t and γ ( ) t matrices are related to, respectively, position-position, momentummomentum, and position-momentum covariances by the following relations:

$$\mathbf A ^ { - 1 } = \langle \delta R _ { i } \delta R _ { j } \rangle ^ { - 1 } = \alpha - \gamma \cdot \beta ^ { - 1 } \cdot \gamma ^ { T } \, \quad ( 4 ) \quad \text{is pr}$$

$$\mathbf B ^ { - 1 } = \langle \delta P _ { i } \delta P _ { j } \rangle ^ { - 1 } = - \gamma ^ { T } + \beta \cdot \gamma ^ { - 1 } \cdot \alpha \, \quad ( 5 ) \quad \text{and}$$

$$\mathbf r ^ { - 1 } = \langle \delta R _ { i } \delta P _ { j } \rangle ^ { - 1 } = \beta - \gamma ^ { T } \cdot \alpha ^ { - 1 } \cdot \gamma ^ { T } \. \quad ( 6 ) \quad \text{when}$$

The evolution of the density Wigner-space density matrix ρ ( R P , , t ) is determined by the propagation in time of the parameters R ( ), t P ( ), t A ( ), t B ( ) and t Γ ( ). t Analogously to the time-dependent Hartree-Fock or timedependent density functional theory for electrons, the time evolution is obtained by imposing the least action principle 21 , leading to the time-dependent self-consistent Liouville-von Neumann equation for the density matrix:

$$i \hbar { \frac { } { } } \frac { \partial \hat { \rho } } { \partial t } = [ \mathcal { H } [ \hat { \rho } ], \hat { \rho } ] \,, \quad \quad ( 7 ) \quad \begin{matrix} \text{di} \\ \text{en} \\ \text{ior} \end{matrix}$$

where H [ ˆ] ρ is a self-consistent harmonic Hamiltonian whose parameters depend on the anharmonic potential and the density matrix ˆ and the square brackets indiρ cate the quantum commutator (more details in Ref. 20 ). Notably, in TD-SCHA, H [ ˆ] ρ is local in time, so the time evolution depends only on the current quantum state. Expressing Eq. 7 in the Wigner formalism and substituting the Gaussian form for the density matrix leads to the set of differential equations

$$\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \ delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \ dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots
\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \ delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \delta \ delta
\dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots
\begin{cases} \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } } = \mathcal { P } } \\ \dot { \mathcal { R } \end{cases} \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \dots \end{cases}$$

where the dot over a tensor ˙ ◦ indicates the timederivative, the product between tensors is the standard rows-by-columns contraction among all indices, and the dagger symbol indicates the matrix transposition operation O † ij = O ji . Here, the atomic potential energy landscape (PES), V ( R , t ), enters the quantum averages of forces ⟨ f ⟩ and the average curvature tensor ⟨ ∂ 2 V ⟩ , defined as

$$\langle f _ { a } \rangle ( t ) = - \int d { \mathbf R } d { \mathbf P } \frac { \partial V } { \partial R _ { a } } ( t ) \rho ( { \mathbf R }, { \mathbf P }, t ), \quad ( 9 )$$

$$\langle \partial _ { a b } ^ { 2 } V \rangle ( t ) = \int d { \mathbf R } d { \mathbf P } \frac { \partial ^ { 2 } V } { \partial R _ { a } \partial R _ { b } } ( t ) \rho ( { \mathbf R }, { \mathbf P }, t ). \quad ( 1 0 )$$

The solution of Eq. (8) provides the quantum state ρ ( R P , , t ), enabling the direct computation of the time envelope of any quantum observable.

The stationary solution of Eq. (8) coincides with the equilibrium fixed-volume state that minimizes the Helmholtz free energy and can be obtained with the standard SSCHA algorithm 21,22 .

When simulating a pump-probe experiment, the system is prepared at equilibrium and perturbed with a radiation pulse modeled as a time-dependent external potential V ext ( R , t ). The overall potential that enters in Eq. (9) and (10) is

$$V ( { \mathbf R }, t ) = V _ { \text{BO} } ( { \mathbf R } ) + V _ { \text{ext} } ( { \mathbf R }, t ), \quad \ \ ( 1 1 )$$

where V BO ( R ) is the instantaneous interaction potential of nuclei within the Born-Oppenheimer approximation that depends only on the nuclear positions.

## III. NUMERICAL INTEGRATION OF TD-SCHA EQUATIONS

Numerical solutions of the TD-SCHA equations of motion (8) have so far been limited to simple onedimensional models 21 and linear-response calculations, enabled by an efficient Lanczos algorithm 20 . The major challenge to applying TD-SCHA in the out-ofequilibrium regime is associated with the cost of calculating the averages ⟨ f ⟩ and particularly ⟨ ∂ 2 V ⟩ , Eq. (10), as

sampling the second derivatives of the potential is computationally expensive. This section introduces a finitedifference scheme to integrate the TD-SCHA equations with an error scaling as O ( dt 3 ) that requires the compu- tation of Eq. (10) only once per time step. Expanding the time evolution of TD-SCHA quantities in the Taylor series to second order in the time step increment dt , we derive the following expressions:

$$\begin{cases} \mathcal { R } _ { t + d t } = \mathcal { R } _ { t } + \mathcal { P } _ { t } d t + \frac { 1 } { 2 } \left \langle \mathfrak { E } _ { t } \right \rangle d t ^ { 2 } + \mathcal { O } ( d t ^ { 3 } ) \\ \mathcal { P } _ { t + d t } = \mathcal { P } _ { t } + \left \langle \mathfrak { E } _ { t } \right \rangle d t + \frac { 1 } { 2 } \mathcal { P } _ { t ^ { \prime \prime } } d t ^ { 2 } + \mathcal { O } ( d t ^ { 3 } ) \\ A _ { t + d t } = A _ { t } + \left ( \Gamma + \Gamma ^ { 1 } \right ) _ { t } d t + \frac { 1 } { 2 } \left ( B - A \left ( \partial ^ { 2 } V \right ) \right ) _ { t } d t ^ { 2 } + \frac { 1 } { 2 } \left ( B - \left ( \partial ^ { 2 } V \right ) A \right ) _ { t } d t ^ { 2 } + \mathcal { O } ( d t ^ { 3 } ) \\ B _ { t + d t } = B _ { t } - \left ( \left \langle \partial ^ { 2 } V \right ) \Gamma + \Gamma ^ { 1 } \left \langle \partial ^ { 2 } V \right ) _ { t } d t + \frac { 1 } { 2 } \mathcal { B } _ { t ^ { \prime \prime } } d t ^ { 2 } + \mathcal { O } ( d t ^ { 3 } ) \\ \Gamma _ { t + d t } = \Gamma _ { t } + \left ( B - A \left \langle \partial ^ { 2 } V \right ) _ { t } d t + \frac { 1 } { 2 } \Gamma _ { t ^ { \prime \prime } } d t ^ { 2 } + \mathcal { O } ( d t ^ { 3 } ) \end{cases}$$

Notably, while the equations for the evolution of R and A are explicit up to O ( dt 3 ), we need the values of the second derivatives of B t , Γ t and P t .

To preserve the O ( dt 3 ) error of the time propagation we use the central difference formula to approximate those second derivatives (see Appendix B):

Euler method. Since Γ is related to the derivatives of A and B (see Eqs. 8), the semi-implicit Euler (SIE) scheme consists in updating Γ first, then computing A t + dt and B t + dt by using Γ t + dt instead of Γ t :

$$F _ { t + d t } = F _ { t } + \frac { 1 } { 2 } ( F _ { t } ^ { \prime } + F _ { t + d t } ^ { \prime } ) d t + \mathcal { O } ( d t ^ { 3 } ) \, \quad ( 1 3 ) \quad \begin{cases} \mathbf B _ { t } \end{cases}$$

where F represents a generic variable. This expression only requires the knowledge of its first derivatives at times t and t + dt , but not of the second derivative. We can rely on two observations: (i) The calculation of ⟨ f ⟩ and ⟨ ∂ 2 V ⟩ depends only on R and A as the potential is a function of only the positions (see Eqs. 9-10 and Appendix A); (ii) both R and A can be integrated explicitly with accuracy O ( dt 3 ). The GV algorithm we devise comprises the following steps: (I) calculate R t + dt and A t + dt with accuracy O ( dt 3 ) using the first and the third of Eqs. 12 respectively; (II) use R t + dt and A t + dt to determine ⟨ f ⟩ t + dt and ⟨ ∂ 2 V ⟩ t + dt ; (III) determine P , B and Γ using Eq. 13

$$\begin{cases} \mathcal { P } _ { t + d t } = \mathcal { P } _ { t } + \langle \mathfrak { f } \rangle _ { \bar { t } } d t + \mathcal { O } ( t ^ { 3 } ) & \mathbf R \\ \mathbf B _ { t + d t } = \mathbf B _ { t } - \left ( \langle \partial ^ { 2 } V \rangle \mathbf \Gamma + \mathbf \Gamma ^ { \dagger } \langle \partial ^ { 2 } V \rangle \right ) _ { \bar { t } } d t + \mathcal { O } ( d t ^ { 3 } ) & \mathbf H _ { t } \\ \mathbf \Gamma _ { t + d t } = \mathbf \Gamma _ { t } + \left ( \mathbf B - \mathbf A \langle \partial ^ { 2 } V \rangle \right ) _ { \bar { t } } d t + \mathcal { O } ( d t ^ { 3 } ) \. & \text{Her} \\ \text{schem} \\ \text{Here we use the shorthand $F_{\bar{t}$ to indicate $\frac{1}{2}(F_{t}+F_{t+dt}$)}. } & \text{and $\mathbf \Gamma$} \end{cases}$$

Here we use the shorthand F ¯ t to indicate 1 2 ( F t + F t + dt ). Even though Eq. 14 is implicit in the variables B and Γ , it requires the calculation of ⟨ f ⟩ and ⟨ ∂ 2 V ⟩ just once per time step (due to observation (i)), and it evolves the parameters with accuracy O ( dt 3 ). Furthermore, the equations for integrating R and P coincide with the familiar velocity Verlet scheme 30 .

We now compare the above GV integration scheme with alternatives, namely explicit Euler (EE) and semiimplicit Euler (SIE) algorithms, showing that the GV is more accurate than both EE and SIE and stable for larger values of the time-step dt . Alternative integration strategies for Eqs. (12) are founded on the semi-implicit

$$\underset { ( 1 3 ) } { \text{$\mathcal{Sigma}$} } \quad \mathcal { T } _ { t + a t } \, \mathcal { Y } _ { t } \, \text{$\mathcal{Sigma}$} \, \text{$\mathcal{T}_{t+a t }$} \, \text{$\mathcal{V}_{t}$} \, \text{$\mathcal{V}_{t}$} } \, \text{$\mathcal{V}_{t}$} } \,. \underset { ( 1 3 ) } { \text{$\mathcal{A}$} } \quad \begin{cases} \mathbf A _ { t + d t } = \mathbf A _ { t } + \left ( \mathbf \Gamma + \mathbf \Gamma ^ { \dagger } \right ) _ { t + d t } d t + \mathcal { O } ( d t ^ { 2 } ) \\ \mathbf B _ { t + d t } = \mathbf B _ { t } - \left ( \langle \partial ^ { 2 } V \rangle _ { t } \mathbf \Gamma _ { t + d t } + \mathbf \Gamma _ { t + d t } ^ { \dagger } \langle \partial ^ { 2 } V \rangle _ { t } \right ) d t + \mathcal { O } ( d t ^ { 2 } ) \\ \text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text	$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text}$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\font$$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\tt$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text
$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text{$\text                                                                                                                                                                                                         $ \text{                                                                                                                                                                                                        \text{                                                                                                                                                                                                     $ \text$                                                                                                                                                                                                        }                                                                                                                                                                                                       
                                                                                                                                                                                                        |                                                                                                                                                                                                        -                                                                                                                                                                                                        +$$

Instead, R and P are evolved according to Verlet. These integration schemes can be compared to the simplest approach, the explicit Euler (EE) scheme, where all the parameters are evolved simultaneously:

$$\intertext { a n e } \colon \quad & \text{p. a. c. v o v a m u a c o v a m u a c o v a n y } \text{d. } \text{p. } \text{allicitly } \begin{cases} \mathbf A _ { t + d t } = \mathbf A _ { t } + \left ( \mathbf \Gamma + \mathbf \Gamma ^ { \dagger } \right ) _ { t } d t + \mathcal { O } ( d t ^ { 2 } ) \\ \mathbf B _ { t + d t } = \mathbf B _ { t } - \left ( \langle \partial ^ { 2 } V \rangle _ { t } \mathbf \Gamma _ { t } + \mathbf \Gamma ^ { \dagger } _ { t } \langle \partial ^ { 2 } V \rangle _ { t } \right ) d t + \mathcal { O } ( d t ^ { 2 } ) \end{cases}. \\ \intertext { f. e q s. } \text{rmine } \text{using}$$

## IV. STABILITY OF THE INTEGRATION SCHEMES

Here, we investigate the stability of the integration schemes introduced in the previous section. Particularly, we are interested in the dynamics of the variables A B , , and Γ , which are not present in classical Newtonian nuclear dynamics equations.

Let us consider a 1D wave packet evolving in a Harmonic potential. Thanks to the constant curvature κ of the PES, A B , and Γ do not depend on the average coordinates (centroids) R and P since

$$\langle \partial ^ { 2 } V \rangle = \kappa. \quad \text{(1 7)}$$

The TD-SCHA equations for these variables reduce to

$$\begin{pmatrix} \dot { A } \\ \dot { B } \\ \dot { \Gamma } \end{pmatrix} = \begin{pmatrix} 0 & 0 & 2 \\ 0 & 0 & - 2 \kappa \\ - \kappa & 1 & 0 \end{pmatrix} \begin{pmatrix} A \\ B \\ \Gamma \end{pmatrix} \. \quad \ \ ( 1 8 )$$

.

By rescaling the parameters as

$$\begin{cases} A ^ { \prime } = \sqrt { \frac { \kappa } { 2 } } A \\ B ^ { \prime } = \frac { 1 } { \sqrt { 2 \kappa } } B & \quad ( 1 9 ) \quad \text{de} \\ t ^ { \prime } = \sqrt { 2 \kappa } t & \quad \text{ler} \\ \quad \text{less} \end{cases}$$

we get a generalized equation that does not depend on the PES curvature κ

$$\begin{pmatrix} \dot { A } ^ { \prime } \\ \dot { B } ^ { \prime } \\ \dot { \Gamma } ^ { \prime } \end{pmatrix} = \begin{pmatrix} 0 & 0 & 1 \\ 0 & 0 & - 1 \\ - 1 & 1 & 0 \end{pmatrix} \begin{pmatrix} A ^ { \prime } \\ B ^ { \prime } \\ \Gamma ^ { \prime } \end{pmatrix} \. \quad \ \ ( 2 0 ) \quad \ k n o \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

The above equation can be written in a compact notation as

$$\dot { x } = M \cdot x \quad \quad \ \ ( 2 1 ) \quad p ($$

where the propagation matrix M is skew-symmetric. This symmetry imposes that the norm of x is conserved

$$A ^ { \prime 2 } + B ^ { \prime 2 } + \Gamma ^ { \prime 2 } = \text{const }. \quad \quad ( 2 2 )$$

In the following, we omit the prime symbol to maintain a cleaner notation and introduce the integer step n as n = t/dt . The stability of the methods is investigated by calculating the step transformation matrix S ( dt ), which connects the degrees of freedom at time step n +1 with those at step n :

$$\mathbf x _ { n + 1 } = \mathbf S ( d t ) \mathbf x _ { n } \. \quad \quad \quad ( 2 3 )$$

Iterating Eq. (23), we obtain

$$\mathbf x _ { n + 1 } = \mathbf S ( d t ) ^ { n } \mathbf x _ { 0 } \. \quad \ \ ( 2 4 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

The stability condition is achieved if the propagator S ( dt ) n remains finite for arbitrary large powers n . This is equivalent to requiring that all its eigenvalues λ are such that | λ | ≤ 1. We calculate the step transformation matrices for EE, SIE, and GV schemes. The details of the derivation are reported in Appendix B. For the EE method we find that

$$\lambda _ { \max } = 1 + 2 d t ^ { 2 } > 1 \quad \forall d t \, \quad \quad ( 2 5 ) \quad,$$

meaning that the EE method is unconditionally unstable . The stability condition obtained for the SIE method is instead

$$d t _ { S I E } \leq \frac { 1 } { \omega } \, \quad \ \ \ ( 2 6 ) \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \
 \quad d t _ { S I E } \leq \frac { 1 } { \omega }$$

where ω = √ κ is the frequency of the harmonic oscillator, while for the GV method we obtain

$$\ m e n u { u } \, w e \, \chi \varpropto \Gamma \Gamma \\ d t _ { G V } & \leq \frac { \sqrt { 2 } } { \omega } \. & \quad ( 2 7 ) & \quad \langle \\ & \quad \rho.$$

Thus, both the SIE and GV are stable for sufficiently small dt , with the GV method having a larger stability range. The integration step of the GV algorithm must thus be shorter than approximately 1/5th of the shortest period of vibrational motion of an atomic system.

## V. STOCHASTIC FORMULATION

The ensemble averages of the potential energy and its derivatives are multidimensional integrals that are challenging to calculate. One strategy for addressing this issue involves expanding the potential energy in a Taylor series centered at a high-symmetry point of the structure , which allows analytically computing the thermodynamic averages. This approach relies on the analytic knowledge of Gaussian integrals, and it is at the basis of the self-consistent phonon (SCP) approach 31 . The alternative consists in evaluating the integrals through a stochastic Monte Carlo algorithm, as exploited by the SSCHA approach 22 . Here, we introduce a stochastic formulation for the TD-SCHA. The ensemble average of the potential energy on the nuclear density in Eq. (3) can be calculated as:

$$\langle V \rangle _ { \mathcal { D } } = \frac { 1 } { N _ { c } } \sum _ { i = 1 } ^ { N _ { c } } V ( \mathcal { R } ( t ) + \mathbf J ( t ) \cdot \mathbf Y _ { i } ( t ) ) \, \quad ( 2 8 )$$

where y i ( ) are i.i.d. t normal random variables , N c is the number of stochastic configurations, and the subscript D stands for the discrete evaluation of the ensemble average (a more detailed introduction is in Appendix C). J ( ) is t the principal square root (one of many possibilities) of the position-position covariance A ( ): t

$$J _ { a b } = \sum _ { \mu } \sqrt { \lambda _ { \mu } } e _ { \mu a } e _ { \mu b }$$

where λ µ and e µ are respectively eigenvalues and eigenvectors of A . The TD-SCHA equations only require the averages of the first and second derivatives of the potential. In our formulation, both averages require only the calculation of forces, which can be obtained either from first principles or from surrogate machine-learning force field models. The average of the first derivative of the potential simplifies to

$$\left \langle \frac { \partial V } { \partial R _ { a } } \right \rangle _ { \mathcal { D } } = - \frac { 1 } { N _ { c } } \sum _ { i = 1 } ^ { N _ { c } } f _ { a } ( \mathcal { R } + \mathbf J \cdot \mathbf y _ { i } ) \, \quad ( 3 0 )$$

where the Cartesian force component is

$$f _ { a } = - \frac { \partial V } { \partial R _ { a } } \. \quad \quad \quad ( 3 1 )$$

(Here, we have omitted the time dependence for clarity in notation). The calculation of the ensemble average of the second derivatives leverages integration by parts to solely utilize the forces 23

$$\left \langle \frac { \partial ^ { 2 } V } { \partial R _ { a } \partial R _ { b } } \right \rangle _ { \mathcal { D } } = - \sum _ { c d } A _ { a c } ^ { - 1 } \sum _ { i = 1 } ^ { N _ { c } } J _ { c d } \, y _ { d i } f _ { b } ( \mathcal { R } + \mathbf J \cdot y _ { i } ) \,. \ ( 3 2 )$$

Eq. 32 is symmetric in the Cartesian indexes a and b only in the limit N c →∞ . For a finite number of configurations, it must be symmetrized:

$$\stackrel { \dots } { \ r t e s t } \quad \left \langle \frac { \partial ^ { 2 } V } { \partial R _ { a } \partial R _ { b } } \right \rangle _ { \mathcal { D } } ^ { s y m } = \frac { 1 } { 2 } \left \langle \frac { \partial ^ { 2 } V } { \partial R _ { a } \partial R _ { b } } \right \rangle _ { \mathcal { D } } + \frac { 1 } { 2 } \left \langle \frac { \partial ^ { 2 } V } { \partial R _ { a } \partial R _ { b } } \right \rangle _ { \mathcal { D } } ^ { T } \ ( 3 3 )$$

The stochastic evaluation of these integrals is characterized by Gaussian noise, which decreases as 1 /N c . If the random displacements y i are sampled in uncorrelated way at each time step, this implies the presence of stochastic noise as input in the TD-SCHA differential equations, which can significantly affect their accuracy and stability. As demonstrated in Ref. 21 , the TD-SCHA equations conserve the total energy in the absence of external potentials acting on the system:

$$\frac { d } { d t } \sum _ { a } \left [ \frac { B _ { a a } + \mathcal { P } _ { a } ^ { 2 } } { 2 } \right ] + \frac { d } { d t } \left \langle V \right \rangle = 0 \quad \ \ ( 3 4 ) \quad \ \ r a t e \, \\ \quad \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \$$

However, the total energy is conserved only in the limit for N c → ∞ . We can show that these issues can be fixed by using the same random configurations ¯ y i in different time step evaluations, which we refer to as correlated sampling. This choice introduces a systematic bias but eliminates the stochastic noise of the ensemble averages of forces and curvatures, making the time evolution smooth. Moreover, we can demonstrate that for onedimensional problems, energy conservation holds true for any finite number of configurations when a constant ¯ y i is employed as dt → 0 (see Appendix D). For higherdimensional problems, energy conservation remains dependent on the number of configurations due to arbitrariness in the definition of J . Nevertheless, the correlated sampling approach drastically improves the energy conservation for a given number of configurations N , as c shown numerically in the following sections. The possibility of choosing a gauge for J that allows for energy conservation independent of the number of configurations in higher-dimensional problems is discussed in the Appendix D and will be the subject of future research.

FIG. 1: The solid curve represents the potential energy surface of the model, while the dashed Gaussian corresponds to the equilibrium nuclear distribution.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00902v2-figure-ab3b9a05-4d75-4ea4-9ef5-7b402a2caf30.jpg)

## VI. TESTS

We test the integration schemes on a one-dimensional model with potential energy

$$V ( u ) = \frac { 1 } { 2 } \left ( - a u ^ { 2 } - b u ^ { 3 } + c u ^ { 4 } \right ) \,, \quad \ \ ( 3 5 )$$

where a = 1 00eV ˚ , . / A 2 b = 1 00eV ˚ , . / A 3 and c = 1 00 eV . / ˚ . A 4 The ionic mass is 1 amu. As illustrated in Fig. 1, this potential exhibits two local minima separated by a barrier of approximately 0 58 eV. . The initial conditions for the parameters R , A B Γ , , correspond to thermodynamic equilibrium, determined by solving the SSCHA equations 22 at 100 K. In such equilibrium, the nuclear density is centered at the lowest minimum, with a spread of about 0.5 ˚ due to the ion's light mass. A The initial momentum of the oscillator is set to P √ m = 0 075 . √ eV. We solve the TD-SCHA equations using the stochastic formulation, comparing uncorrelated and correlated random displacements with N c = 100 random configurations and a time-step of 1 fs for 400 fs with the GV scheme. The results, reported in FIG. 2, are compared with the solution of the TD-SCHA equations in the N c →∞ limit (numerically sampled via trapezoidal integration on a dense grid). The uncorrelated sampling algorithm quickly deviates from the N c → ∞ solution. The correlated sampling algorithm, instead, remains stable throughout the dynamics due to the suppression of the stochastic noise across different time steps. FIG. 3 report the energy conservation over time for the uncorrelated GV, correlated GV, and correlated SIE methods. The uncorrelated sampling approach fails to conserve energy, although reducing the simulation time step by half partially mitigates this issue. The SIE algorithm with correlated sampling suppresses the energy oscillations, but it suffers from a uniform energy drift that decreases with a reduction in the simulation time step. In contrast, the correlated GV method demonstrates flawless energy conservation for both time steps, exhibiting no energy fluctuations or drift. Indeed, the quality of energy conservation increases with the number of configurations for the uncorrelated sampling algorithm, while it is unaffected when employing the correlated sampling, as shown in FIG. 4a (see Appendix D for the formal proof).

## VII. DYNAMICS IN SRTIO 3

In this section, we showcase a realistic application of TD-SCHA by investigating the out-of-equilibrium quantum dynamics in SrTiO 3 (STO) that follows a resonanceexciting short pulse of infrared light. STO is a prototypical quantum paraelectric 32-34 , where nuclear quantum fluctuations suppress the ferroelectric order at low temperatures. STO has been extensively studied due to the emerging phenomena occurring when driven out of equilibrium by strong electric field pulses 35-40 . Notably,

FIG. 2: The panels in the upper row compare the N c →∞ solution to the TD-SCHA equations (blue line) with the finite sampling solution using uncorrelated displacement (red line). The lower row of panels compares the exact solution to the TD-SCHA equation with the numerical solution using correlated sampling (green line). The quantities represented are (moving from left to right) the position R , position-position correlation A , average force ⟨ f ⟩ , and average curvature ⟨ ∂ V 2 ⟩ . For both the calculations with correlated and uncorrelated displacements, we employed N c = 100 and a time-step of 1 fs.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00902v2-figure-882f041f-41dc-4202-ba11-9f197595161b.jpg)

the anharmonic coupling between them.

FIG. 3: Energy conservation for different integration schemes and time-step.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00902v2-figure-ad32069f-8543-44b1-9d52-0fe5b85c8fde.jpg)

irradiating STO with a THz-frequency pulse at low temperatures induces a long-lasting second harmonic generation signal 35,36 , suggesting the occurrence of a lightinduced ferroelectric phase transition; however, this interpretation is still debated 41 . Furthermore, recent studies have demonstrated the possibility of transferring energy from lower frequency phonons, pumped by the optical excitation, to higher frequency phonons in an outof-equilibrium process called upconversion 37 , enabled by

In our simulation, a 40-atom supercell of STO originally equilibrated at 100 K through a static SSCHA calculation is excited by an infrared pulse with an amplitude of 833 kV/cm, which is resonant with the soft phonon mode (SPM, represented in Fig. 5a of STO. We account for the light-matter interaction in the dipole approximation through the Born effective charges. Details on the coupling with the electric field and the atomic energy landscape calculation are discussed in appendix E. We integrate the TD-SCHA equations using the GV scheme in the correlated formulation, adopting a time step of 1 fs and sampling the potential energy landscape with N c =4000. The system, originally in equilibrium at 100 K, interacts at t = 0fs with an external pulse of oscillating electric field, triggering a non-equilibrium evolution of the density matrix. Fig. 5b shows the motion of the SPM phonon coordinate

$$Q _ { \mu } = \sum _ { a x } e _ { \mu a x } ( \mathcal { R } _ { a x } - \mathcal { R } _ { a x } ^ { e q } )$$

as a function of the time delay after the pulse. Here R eq ax represents the equilibrium centroid position of the atom a in the direction x , and e µax is the equilibirum eigenvector of the soft phonon mode µ . The irradiation of STO with resonant pulses drives large oscillations of the SPM, which slowly decay due to the interaction with other phonon modes. The blue area in the figure corresponds to the quantum uncertainty in the position of

FIG. 4: (a) Energy conservation of the GV method as a function of the number of configurations N c in the uncorrelated formulation. (b) Energy conservation in the correlated formulation. Note that the y-scale of the two panels differ by more than two orders of magnitudes.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00902v2-figure-2356ae05-e75c-43b8-9ccf-8b66cb2087f4.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00902v2-figure-bf9b9c1a-749b-4858-84aa-089b2da8701a.jpg)

the SPM. It is computed as ± √ A µµ , which is equal to √ ⟨ Q 2 µ ⟩ . The large extent of this uncertainty relative to the motion of Q µ highlights the fundamental importance of quantum effects in the dynamics of STO. An extensive discussion of the relevance of the simulation for the physics of STO goes beyond the scope of this work, and is subject of a separate publication 42 .

## VIII. CONCLUSIONS

In this work, we introduced the first approach to simulate nonequilibrium quantum nuclear dynamics using the TD-SCHA. We derived an integration scheme, the Generalized Verlet, which allows for the evolution of the

FIG. 5: (a) The STO unit cell. The arrows indicate the displacement pattern of the SPM excited resonantly by the THz impulsive pump. (b) Dynamics of the SPM as a function of the time delay after the pulse. The blue area represents the quantum uncertainty, which is comparable to the amplitude of the oscillations, witnessing the importance of nuclear quantum effects in the dynamics.

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00902v2-figure-159fdf90-5a73-42ce-ad20-e1902b6724c4.jpg)

![image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00902v2-figure-ca9f77b9-c4b5-4316-a61e-ec4cf8eaece7.jpg)

equations with an error of O ( dt 3 ), demonstrating that its conditional stability is consistent and of the same order as the Nyquist sampling rate. Additionally, we introduced a stochastic formulation of the TD-SCHA, which enables efficient calculation of ensemble averages while ensuring the stability of the evolution. Finally, we showcased the method's potential, proving it is well-suited for simulating quantum nonequilibrium processes in pumpprobe setups on the scale of hundreds of atoms.

## ACKNOWLEDGEMENTS

This research was partially funded by the Swiss National Science Foundation (SNSF, mobility fellowship P500PT 217861), the Department of Navy award N00014-20-1-2418 issued by the Office of Naval Research and Robert Bosch LLC. L. M. thanks the European

Union under the program Horizon 2020 for the award and funding of the MSCA individual fellowship (grant number 101018714). Computational resources were provided by the FAS Division of Science Research Computing Group at Harvard University.

## CONTRIBUTIONS

F.L. conducted the theoretical derivations, implemented the methods, and carried out the atomistic sim-

## Appendix A: Notes on the Wigner formulation of TD-SCHA

The nuclear density introduced in Eq. 3 corresponds to the most general Gaussian form in the positions R and momenta P variables. It is parametrized by the vectors R and P and the matrices A B , , and Γ . The physical meaning of these parameters is immediately clear after noting that

$$\mathcal { R } ( t ) = \left \langle \mathbf R \right \rangle ( t )$$

$$\mathcal { P } ( t ) = \left \langle \mathbf P \right \rangle ( t )$$

$$A _ { i j } ( t ) & = \langle \delta R _ { i } \delta R _ { j } \rangle$$

$$B _ { i j } ( t ) & = \langle \delta P _ { i } \delta P _ { j } \rangle$$

$$\Gamma _ { i j } ( t ) = \langle \delta R _ { i } \delta P _ { j } \rangle \.$$

These properties can be easily shown using the definition of ensemble averages given in Eq. 2. According to Eqs. A1 and A2, the parameters R and P correspond to the expected values for the positions and momenta, which are analogous to the classical positions and momenta. The variables A B , and Γ , instead, correspond to the positionposition, momentum-momentum and position-momentum correlation matrices, as suggested by Eqs. A3-A5. If an observable is a function of only the position O ( R ), the momenta degrees of freedom can be integrated out in Eq. 2. After some tedious algebra, it is possible to show that in such case

$$\langle O \rangle = \int O ( \mathcal { R } ( t ) + \delta \mathbf R ) \, \sqrt { \frac { 1 } { ( 2 \pi ) ^ { 3 N } \det A } } e ^ { - \frac { 1 } { 2 } \delta \mathbf R \cdot A ^ { - 1 } ( t ) \cdot \delta \mathbf R } \, d ( \delta \mathbf R ) \.$$

We can clearly see that ⟨ O ⟩ ( ) is a function of t R and A only.

## Appendix B: Numerical Stability of Integation Schemes

First of all, we prove an identity for the central finite difference approximation of a derivative. We first write the Taylor expansion of a function F up to the third order centered at t :

$$F _ { t + d t } = F _ { t } + F _ { t } ^ { \prime } d t + \frac { 1 } { 2 } F _ { t } ^ { \prime \prime } d t ^ { 2 } + \mathcal { O } ( d t ^ { 3 } )$$

we can write the same expression by centering it in t + dt ,

$$F _ { t } = F _ { t + d t } - F ^ { \prime } _ { t + d t } d t + \frac { 1 } { 2 } F ^ { \prime \prime } _ { t + d t } d t ^ { 2 } + \mathcal { O } ( d t ^ { 3 } )$$

Adding up these equations we get

$$F _ { t + d t } = F _ { t } + \frac { 1 } { 2 } ( F _ { t } ^ { \prime } + F _ { t + d t } ^ { \prime } ) d t + \frac { 1 } { 2 } ( F _ { t } ^ { \prime \prime } - F _ { t + d t } ^ { \prime \prime } ) d t ^ { 2 } + \mathcal { O } ( d t ^ { 3 } ) \.$$

ulations. A.J. trained the machine-learned potential. L. M. guided the numerical implementation and algorithminc development. L.M. and B.K. supervised the full project.

However, from the following equation we get

leading to the result

$$F _ { t + d t } = F _ { t } + \frac { 1 } { 2 } ( F _ { t } ^ { \prime } + F _ { t + d t } ^ { \prime } ) d t + \mathcal { O } ( d t ^ { 3 } )$$

which corresponds to Eq. 13 of the main text. This equation indicates that we can evolve a function F with accuracy up to third-order by utilizing the average of the derivative at time t and that at time t + dt . Eq. 13 can be employed to derive the Verlet algorithm:

$$\begin{cases} \mathcal { P } _ { t + d t } = \mathcal { P } _ { t } + \frac { 1 } { 2 } \Big ( \langle \mathfrak { f } \rangle _ { t } + \langle \mathfrak { f } \rangle _ { t + d t } \Big ) d t + \mathcal { O } ( d t ^ { 3 } ) \\ \mathcal { R } _ { t + d t } = \mathcal { R } _ { t } + \mathcal { P } _ { t } d t + \frac { 1 } { 2 } \langle \mathfrak { f } \rangle _ { t } d t ^ { 2 } + \mathcal { O } ( d t ^ { 3 } ) \end{cases}.$$

In the following, we study the stability of the explicit Euler, semi-implicit Euler and generalized Verlet methods.

## 1. Explicit Euler method

At each step of the explicit Euler method, the variables are updated according to the following rule:

$$\mathbf x _ { n + 1 } = \begin{pmatrix} 1 & 0 & d t \\ 0 & 1 & - d t \\ - d t & d t & 1 \end{pmatrix} \mathbf x _ { n } = \mathbf S ( d t ) \mathbf x _ { n }$$

leading to

$$x _ { n + 1 } & = \mathbf S ( d t ) ^ { n } \mathbf x _ { 0 }$$

The exponentiation of such a matrix requires the calculation of its eigenvalue, which are

$$\begin{cases} \lambda _ { 1 } = 1 \\ \lambda _ { 2 } = 1 - i \sqrt { 2 } d t \\ \lambda _ { 3 } = 1 + i \sqrt { 2 } d t \end{cases}$$

The stability condition for preventing divergence of the power sequence in Eq. B9 is

$$| 1 \pm i \sqrt { 2 } d t | \leq 1$$

which is never satisfied. The explicit Euler method is thus unconditionally unstable .

## 2. Semi-implicit Euler method

At each step, first Γ is updated:

and then A and B

$$F _ { t + d t } ^ { \prime \prime } = F _ { t } ^ { \prime \prime } + F _ { t } ^ { \prime \prime \prime } d t + \mathcal { O } ( d t ^ { 2 } )$$

$$( F ^ { \prime \prime } _ { t + d t } - F ^ { \prime \prime } _ { t } ) d t ^ { 2 } = F ^ { \prime \prime \prime } _ { t } d t ^ { 3 } + \mathcal { O } ( d t ^ { 4 } )$$

$$\Gamma _ { n + 1 } & = \Gamma _ { n } + ( B _ { n } - A _ { n } ) d t$$

$$\begin{cases} A _ { n + 1 } = A _ { n } + \Gamma _ { n + 1 } d t \\ B _ { n + 1 } = B _ { n } - \Gamma _ { n + 1 } d t \end{cases}$$

Substituting Eq. B12 into B13, we obtain

$$\begin{cases} A _ { n + 1 } = A _ { n } + \Gamma _ { n } d t + ( B _ { n } - A _ { n } ) d t ^ { 2 } \\ B _ { n + 1 } = B _ { n } - \Gamma _ { n } d t - ( B _ { n } - A _ { n } ) d t ^ { 2 } \end{cases}$$

or equivalently

$$\begin{pmatrix} A _ { n + 1 } \\ B _ { n + 1 } \\ \Gamma _ { n + 1 } \end{pmatrix} = \begin{pmatrix} 1 - d t ^ { 2 } & d t ^ { 2 } & d t \\ d t ^ { 2 } & 1 - d t ^ { 2 } & - d t \\ - d t & d t & 1 \end{pmatrix} \begin{pmatrix} A _ { n } \\ B _ { n } \\ \Gamma _ { n } \end{pmatrix} \.$$

The eigenvalues of the step matrix are

$$\begin{cases} \lambda _ { 1 } = 1 & \\ \lambda _ { 2 } = 1 - d t ^ { 2 } - d t \sqrt { d t ^ { 2 } - 2 } & \\ \lambda _ { 3 } = 1 - d t ^ { 2 } + d t \sqrt { d t ^ { 2 } - 2 } & \end{cases}$$

The eigenvalues λ 2 , λ 3 are real for dt ≥ √ 2, or complex conjugated otherwise. We easily note that

$$\lambda _ { 2 } < - 1 \quad \forall d t \geq \sqrt { 2 }$$

meaning that this method is unstable in such range. For dt &lt; √ 2, λ 2 3 , are complex, with modulus

$$| \lambda _ { 2, 3 } | ^ { 2 } = ( 1 - d t ^ { 2 } ) ^ { 2 } + 2 d t ^ { 2 } - d t ^ { 4 } = 1$$

therefore the method is stable . Remembering the definition 19, the stability condition for this method is

$$t \leq \frac { 1 } { \sqrt { \kappa } } = \frac { 1 } { \omega }$$

where ω is the frequency of the harmonic oscillator (we remind that κ = 1 m ∂R ∂ 2 V 2 due to the mass rescaling convention adopted).

## 3. Generalized Verlet

The same analylisis can be performed on the GV algorithm. In this case, the variables are updated according to

$$\begin{cases} A _ { n + 1 } = A _ { n } + \Gamma _ { n } d t + \frac { 1 } { 2 } ( B _ { n } - A _ { n } ) d t ^ { 2 } \\ B _ { n + 1 } = B _ { n } - \frac { 1 } { 2 } ( \Gamma _ { n } + \Gamma _ { n + 1 } ) d t \\ \Gamma _ { n + 1 } = \Gamma _ { n } + \frac { 1 } { 2 } ( B _ { n } - A _ { n } + B _ { n + 1 } - A _ { n + 1 } ) d t \end{cases}$$

After some tedious algebra, it is possible to write explicitly the transformation in matrix form as

$$\begin{pmatrix} A _ { n + 1 } \\ B _ { n + 1 } \\ \Gamma _ { n + 1 } \end{pmatrix} = \begin{pmatrix} 1 - \frac { d t ^ { 2 } } { 2 } & \frac { d t ^ { 2 } } { 2 } & d t \\ \frac { 2 d t ^ { 2 } + 4 d t ^ { 2 } } { 2 d t ^ { 2 } + 8 } & \frac { d t ^ { 3 } - 4 d t } { 2 d t ^ { 2 } + 4 } & \frac { d t ^ { 3 } - 4 d t } { d t ^ { 2 } + 4 } \\ \frac { d t ^ { 3 } - 4 d t } { d t ^ { 2 } + 4 } & \frac { - d t ^ { 3 } + 4 d t } { d t ^ { 2 } + 4 } & \frac { - 3 d t ^ { 2 } + 4 } { d t ^ { 2 } + 4 } \end{pmatrix} \begin{pmatrix} A _ { n } \\ B _ { n } \\ \Gamma _ { n } \end{pmatrix} \,,$$

with eigenvalues

$$\begin{cases} \lambda _ { 1 } = 1 & \\ \lambda _ { 2 } = \frac { 4 - 3 d t ^ { 2 } - 2 \sqrt { 2 } d t \sqrt { d t ^ { 2 } - 4 } } { d t ^ { 2 } + 4 } & \quad. & & ( B 2 2 ) \\ \lambda _ { 3 } = \frac { 4 - 3 d t ^ { 2 } + 2 \sqrt { 2 } d t \sqrt { d t ^ { 2 } - 4 } } { d t ^ { 2 } + 4 } & \end{cases}$$

The eigenvalues are complex for dt &lt; 2, with modulus

$$| \lambda _ { 2, 3 } | ^ { 2 } = 1$$

and the method is stable, whereas for dt &gt; 2

$$\lambda _ { 2 } < - 1$$

thus the method is unstable. To conclude, the stability condition for such an algorithm is

$$\Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \Phi _ { 0 } \cdot \\ d t & < \sqrt { \frac { 2 } { \kappa } } = \frac { \sqrt { 2 } } { \omega } \.$$

The stability range is thus larger than that of the SIE.

## Appendix C: Importance sampling

In the following section we define some properties which will be useful for the development of the theory reported in the main text. The ensemble average of the potential energy at time t is given by

$$\langle V \rangle = \int V ( \mathcal { R } ( t ) + \delta \mathbf R ) \, \sqrt { \frac { 1 } { ( 2 \pi ) ^ { 3 N } \det A } } e ^ { - \frac { 1 } { 2 } \delta \mathbf R \cdot A ^ { - 1 } ( t ) \cdot \delta \mathbf R } \, d ( \delta \mathbf R )$$

where the integration is carried over the variables δR a whose number is 3 N . We now perform a change of variables which is meant to turn the Gaussian in Eq. C1 into a normal distribution. We start from the modal decomposition of the matrix A

$$A _ { a b } = \sum _ { \mu } \lambda _ { \mu } e _ { \mu a } e _ { \mu b }$$

and express the displacement δR a in normal coordinates:

$$\delta R _ { a } = \sum _ { b } J _ { a b } y _ { b } \.$$

The determinant of the Jacobian matrix is

$$\det J = \prod _ { \mu } \sqrt { \lambda _ { \mu } } = \sqrt { \det A } \.$$

The change of variables of Eq. C3 turns the integral C1 into

$$\langle V \rangle = \int V ( \mathcal { R } _ { a } ( t ) + \sum _ { b } J _ { a b } ( t ) y _ { b } ) \ \prod _ { b } \frac { e ^ { - \frac { 1 } { 2 } y _ { b } ^ { 2 } } } { \sqrt { 2 \pi } } d y _ { b } \.$$

Eq. C5 provides the starting point for the stochastic calculation of the ensemble averages, expressed in Eq. 28 of the main text.

## Appendix D: Energy conservation

Here, we outline the main steps introduced in Ref. 21 to derive energy conservation for infinitely many configurations. We then show the conditions under which such conservation holds regardless of the number of configurations, which represents the result of the current work. The time derivative of the quantum kinetic energy reads

$$\frac { d } { d t } \sum _ { i } \langle \frac { P _ { i } ^ { 2 } } { 2 } \rangle = \frac { d } { d t } \frac { \sum _ { i } \mathcal { P } _ { i } ^ { 2 } + \text{Tr} \text{ B} } { 2 } = \sum _ { i } \left ( \mathcal { P } _ { i } \dot { P } _ { i } + \frac { \dot { B } _ { i i } } { 2 } \right )$$

Combining the first and second of Eqs. 8, we obtain

$$\sum _ { i } \mathcal { P } _ { i } \dot { \mathcal { P } } _ { i } = - \sum _ { i } \langle \frac { \partial V } { \partial R _ { i } } \rangle \dot { \mathcal { R } } _ { i }$$

while the third and fourth of Eqs. 8 give

$$\sum _ { i } \frac { \dot { B } _ { i i } } { 2 } = - \frac { 1 } { 2 } \sum _ { i j } \langle \frac { \partial ^ { 2 } V } { \partial R _ { i } \partial R _ { j } } \rangle \dot { A } _ { i j } \.$$

Importantly, these equations also hold when the discrete expressions for forces (Eq. 30) and curvatures (Eq. 33) are used to drive the dynamics

$$\sum _ { i } \mathcal { P } _ { i } \dot { \mathcal { P } } _ { i } = - \sum _ { i } \left \langle \frac { \partial V } { \partial R _ { i } } \right \rangle _ { \mathcal { D } } \dot { \mathcal { R } } _ { i }$$

and

$$\frac { 1 } { 2 } \langle \frac { \partial ^ { 2 } V } { \partial R _ { i } R _ { j } } \rangle = \frac { \partial \langle V \rangle } { \partial \mathcal { A } _ { i j } } \.$$

Inserting Eqs. D6 and D7 into Eq. D1 we obtain

$$\frac { 1 } { 2 } \frac { d } { d t } \left ( \sum _ { i } \mathcal { P } _ { i } ^ { 2 } + \text{Tr} \, \text{B} \right ) = - \frac { \partial \left \langle V \right \rangle } { \partial \mathcal { R } _ { i } } \dot { \mathcal { R } } _ { i } - \frac { \partial \left \langle V \right \rangle } { \partial \mathcal { A } _ { i j } } \dot { A } _ { i j } = - \frac { d } { d t } \left \langle V \right \rangle$$

This proves the conservation of energy in absence of external forces. However, for this energy conservation to hold for any configuration, we need to demonstrate that Eqs. D6 and D7 also hold for the discrete expression of forces and curvatures, namely:

$$\left \langle \frac { \partial V } { \partial R _ { i } } \right \rangle _ { \mathcal { D } } = \frac { \partial \left \langle V \right \rangle _ { \mathcal { D } } } { \partial \mathcal { R } _ { i } }$$

and

$$\frac { 1 } { 2 } \langle \frac { \partial ^ { 2 } V } { \partial R _ { i } R _ { j } } \rangle _ { \mathcal { D } } = \frac { \partial \langle V \rangle _ { \mathcal { D } } } { \partial \mathcal { A } _ { i j } } \.$$

Equation D9 can be easily derived by applying the chain rule to differentiate Equation 28:

$$\left \langle \frac { \partial V } { \partial R _ { a } } \right \rangle _ { \mathcal { D } } = - \frac { 1 } { N _ { c } } \sum _ { i = 1 } ^ { N _ { c } } f _ { a } ( \mathcal { R } + J \cdot y _ { i } ) = \frac { 1 } { N _ { c } } \sum _ { i = 1 } ^ { N _ { c } } \frac { \partial V } { \partial \mathcal { R } _ { a } } ( \mathcal { R } + J \cdot y _ { i } ) = \frac { \partial \left \langle V \right \rangle _ { \mathcal { D } } } { \partial \mathcal { R } _ { a } }$$

The derivation of Eq. D10 is more complex, and does not hold in the most general case. From Eqs. 32 and 33, we know that

$$\left \langle \frac { \partial ^ { 2 } V } { \partial R _ { a } \partial R _ { b } } \right \rangle _ { \mathcal { D } } ^ { s y m } = - \frac { 1 } { 2 } \sum _ { c d } A _ { a c } ^ { - 1 } \sum _ { i = 1 } ^ { N } J _ { c d } y _ { d i } f _ { b } ( \mathcal { R } + J \cdot y _ { i } ) \ - \frac { 1 } { 2 } \sum _ { c d } A _ { b c } ^ { - 1 } \sum _ { i = 1 } ^ { N } J _ { c d } y _ { d i } f _ { a } ( \mathcal { R } + J \cdot y _ { i } ) \.$$

At the same time, the derivative of ⟨ V ⟩ D with the respect to the element A ab is

∂

1

N

c

⟨

V

⟩

=

-

∑∑

∂J

∂A

ab

N

c

D

cd

i

=1

cd

∂A

ab

In order to prove Eq. D10, it would be enough to show that

$$\frac { \partial J _ { c d } } { \partial A _ { a b } } = \frac { 1 } { 4 } \delta _ { a c } J _ { b d } ^ { - 1 } + \frac { 1 } { 4 } \delta _ { b c } J _ { a d } ^ { - 1 } \.$$

In one dimensional problems, identity D14 reduces to

$$\frac { \partial J } { \partial A } = \frac { 1 } { 2 \sqrt { A } }$$

which is satisfied since J = √ A . Therefore, in one dimensional problems, the conservation of energy holds independently on the number of configurations, provided that the time step is reasonably small.

c

i

f

(

R

+

J

·

y

)

$$\sum _ { i } \frac { \dot { B } _ { i i } } { 2 } = - \frac { 1 } { 2 } \sum _ { i j } \left \langle \frac { \partial ^ { 2 } V } { \partial R _ { i } \partial R _ { j } } \right \rangle _ { \mathcal { D } } \dot { A } _ { i j } \,.$$

Following Ref. 21 , it is possible to show that

$$\left \langle \frac { \partial V } { \partial R _ { i } } \right \rangle = \frac { \partial \left \langle V \right \rangle } { \partial \mathcal { R } _ { i } }$$

y

d

(D13)

For higher dimensional problems, the validity of the identity D 14 depends on the definition of J . In fact, there are infinite many ways to obtain the square root of a symmetric matrix, defined as the matrix J such that

$$A = J J ^ { T } \.$$

This is because multiplying the matrix J by any orthogonal matrix O results in a matrix that still satisfies Eq. D16. The definition of J in Eq. 29 does not satisfy the identity D14. In fact, using the rules derived in Appendix 5 of Ref. 23 , we easily obtain

$$\sum _ { a b } \frac { \partial J _ { c d } } { \partial A _ { a b } } \dot { A } _ { a b } = \sum _ { a b } \sum _ { \mu \nu } \frac { e _ { a \mu } e _ { b \nu } e _ { c \mu } e _ { d \nu } } { \sqrt { \lambda _ { \mu } } + \sqrt { \lambda _ { \nu } } } \dot { A } _ { a b } \.$$

The problem of finding an expression for J that satisfies identity D14 is very important and will be the subject of future research efforts.

## Appendix E: Atomic energy landscape

The TD-SCHA simulations on STO utilize machine-learned potentials (MLP) to model the potential energy surface. We opted for FLARE 43 for its active learning capabilities, enabling efficient data generation and rapid inference times 44 . Through active learning, we explored various temperatures and volumes to develop a broadly applicable potential.

Simulations were conducted at 100, 300, and 500 K, each for 200 ps, at the DFT-relaxed lattice parameter, as well as at ± 2% strain, amounting to a total of 1.8 ns of dynamics. A timestep of 2 fs and a thermostat damping time of 200 fs were employed, utilizing the default Nos´-Hoover thermostat in LAMMPS e 45 .

DFT calculations were performed at the PBE level of theory using Quantum ESPRESSO 46 . We used a plane-wave cutoff of 80 Ry and a k-point grid of 6 × 6 × 8 for the 20-atom cell, adopting the pseudopotentials recommended by the SSSP efficiency library 47 .

The forces on the atom i by the electric field is obtained as

$$\mathfrak { f } _ { i } = \frac { 1 } { \varepsilon _ { e f f } } \, \mathbf Z _ { i } \cdot \mathcal { E } \.$$

Here E is the external electric field, ε eff is the dielectric constant and Z i are the Born effective charge tensors. The effective charges are computed through DFPT, using the same parameters as above. Their value for the different atomic species is reported in Tab. I. Here we employ the screening model proposed in Refs. 38,48

|     |   Z ∗ xx |   Z ∗ yy |   Z ∗ zz |
|-----|----------|----------|----------|
| Ti  |    7.338 |    7.338 |    7.338 |
| Sr  |    2.549 |    2.549 |    2.549 |
| O 1 |   -2.024 |   -5.845 |   -2.024 |
| O 2 |   -2.024 |   -2.024 |   -5.845 |
| O 3 |   -5.845 |   -2.024 |   -2.024 |

Table I: Born effective charges for the cubic STO unitcell, computed through DFPT.

$$\varepsilon _ { e f f } = \frac { 1 + \sqrt { \varepsilon _ { D F P T } } } { 2 }$$

with ε DFPT = 6 31. .

∗ Corresponding authors: libbi@g.harvard.edu, lorenzo.monacelli@uniroma1.it, bkoz@g.harvard.edu

1 T. E. Markland and M. Ceriotti, Nature Reviews Chemistry 2 , 0109 (2018).

2 I. Errea, M. Calandra, C. J. Pickard, J. R. Nelson, R. J. Needs, Y. Li, H. Liu, Y. Zhang, Y. Ma, and F. Mauri, Nature 532 , 81 (2016).

3 I. Errea, F. Belli, L. Monacelli, A. Sanna, T. Koretsune, T. Tadano, R. Bianco, M. Calandra, R. Arita, F. Mauri, and J. A. Flores-Livas, Nature 578 , 66 (2020).

- 4 M. Takesada, M. Itoh, and T. Yagi, Phys. Rev. Lett. 96 , 227602 (2006).
- 5 M. Zacharias, C. E. Patrick, and F. Giustino, Phys. Rev. Lett. 115 , 177401 (2015).
- 6 G. Antonius, S. Ponc´, P. Boulanger, M. Cˆt´, e o e and X. Gonze, Phys. Rev. Lett. 112 , 215501 (2014).
- 7 U. Aseginolaza, R. Bianco, L. Monacelli, L. Paulatto, M. Calandra, F. Mauri, A. Bergara, and I. Errea, Phys. Rev. Lett. 122 , 075901 (2019).
- 8 D. Marx and M. Parrinello, The Journal of Chemical Physics 104 , 4077 (1996), https://pubs.aip.org/aip/jcp/articlepdf/104/11/4077/19179556/4077 1 online.pdf.
- 9 D. M. Ceperley, Rev. Mod. Phys. 67 , 279 (1995).
- 10 M. Parrinello and A. Rahman, The Journal of Chemical Physics 80 , 860 (1984), https://pubs.aip.org/aip/jcp/articlepdf/80/2/860/18946078/860 1 online.pdf.
- 11 D. Chandler and P. G. Wolynes, The Journal of Chemical Physics 74 , 4078 (1981), https://pubs.aip.org/aip/jcp/articlepdf/74/7/4078/18928433/4078 1 online.pdf.
- 12 I. R. Craig and D. E. Manolopoulos, The Journal of Chemical Physics 121 , 3368 (2004), https://pubs.aip.org/aip/jcp/articlepdf/121/8/3368/19312704/3368 1 online.pdf.
- 13 I. Miller, Thomas F. and D. E. Manolopoulos, The Journal of Chemical Physics 123 , 154504 (2005), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/1.2074967/13408998/154504 1 online.pdf.
- 14 B. Chen, M. A. Gomez, M. Sehl, J. D. Doll, and D. L. Freeman, The Journal of Chemical Physics 105 , 9686 (1996), https://pubs.aip.org/aip/jcp/article-pdf/105/21/9686/19172206/9686 1 online.pdf.
- 15 D. Kim, J. D. Doll, and J. E. Gubernatis, The Journal of Chemical Physics 106 , 1641 (1997), https://pubs.aip.org/aip/jcp/article-pdf/106/4/1641/19148728/1641 1 online.pdf.
- 16
- G. Krilov and B. J. Berne, The Journal of Chemical Physics 111 , 9147 (1999), https://pubs.aip.org/aip/jcp/articlepdf/111/20/9147/19249140/9147 1 online.pdf.
- 17 G. Baym and N. D. Mermin, Journal of Mathematical Physics 2 , 232 (1961), https://pubs.aip.org/aip/jmp/articlepdf/2/2/232/19141881/232 1 online.pdf.
- 18 A. Alexandru, G. Ba ˚ ar, P. F. Bedaque, A and N. C. Warrington, Reviews of Modern Physics 94 (2022), 10.1103/revmodphys.94.015006.
- 19 L. M¨ uhlbacher and E. Rabani, Phys. Rev. Lett. 100 , 176403 (2008).
- 20 L. Monacelli and F. Mauri, Phys. Rev. B 103 , 104305 (2021).
- 21 A. Siciliano, L. Monacelli, G. Caldarelli, and F. Mauri, Phys. Rev. B 107 , 174307 (2023).
- 22 L. Monacelli, R. Bianco, M. Cherubini, M. Calandra, I. Errea, and F. Mauri, Journal of Physics: Condensed Matter 33 , 363001 (2021).
- 23 R. Bianco, I. Errea, L. Paulatto, M. Calandra, and F. Mauri, Phys. Rev. B 96 , 014111 (2017).
- 24 L. Monacelli, I. Errea, M. Calandra, and F. Mauri, Phys. Rev. B 98 , 024106 (2018).
- 25 R. Bianco, L. Monacelli, M. Calandra, F. Mauri, and I. Errea, Phys. Rev. Lett. 125 , 106101 (2020).
- 26 U. Aseginolaza, J. Diego, T. Cea, R. Bianco, L. Monacelli, F. Libbi, M. Calandra, A. Bergara, F. Mauri, and I. Errea, Nature Physics (2024), 10.1038/s41567-024-02441-z.
- 27 L. Monacelli, I. Errea, M. Calandra, and F. Mauri, Nature Physics 17 , 63 (2021).
- 28 E. Wigner, Phys. Rev. 40 , 749 (1932).
- 29 K. ¨ A mre, ° E. ˜ zizmir, A M. Rosenbaum, and P. F. Zweifel, Journal of Mathematical Physics 8 , 1097 (1967), https://pubs.aip.org/aip/jmp/article-pdf/8/5/1097/19085677/1097 1 online.pdf.
- 30 L. Verlet, Phys. Rev. 159 , 98 (1967).
- 31 T. Tadano and S. Tsuneyuki, Phys. Rev. B 92 , 054301 (2015).
- 32 D. Shin, S. Latini, C. Sch¨fer, S. A. Sato, U. De Giovannini, H. H¨bener, a u and A. Rubio, Phys. Rev. B 104 , L060103 (2021).
- 33 C. Verdi, L. Ranalli, C. Franchini, and G. Kresse, Phys. Rev. Mater. 7 , L030801 (2023).
- 34 T. Esswein and N. A. Spaldin, Phys. Rev. Res. 4 , 033020 (2022).
- 35 T. F. Nova, A. S. Disa, M. Fechner, and A. Cavalleri, Science 364 , 1075 (2019), https://www.science.org/doi/pdf/10.1126/science.aaw4911.

36

X.

Li,

T.

Qiu,

J.

Zhang,

E.

Baldini,

J.

Lu,

A.

M.

Rappe, and

K.

A.

Nelson,

Science

364

,

1079

(2019), https://www.science.org/doi/pdf/10.1126/science.aaw4913.

- 37 M. Kozina, M. Fechner, P. Marsik, T. van Driel, J. M. Glownia, C. Bernhard, M. Radovic, D. Zhu, S. Bonetti, U. Staub, and M. C. Hoffmann, Nature Physics 15 , 387 (2019).
- 38 D. Shin, S. Latini, C. Sch¨fer, S. A. Sato, E. Baldini, U. De Giovannini, H. H¨bener, a u and A. Rubio, Phys. Rev. Lett. 129 , 167401 (2022).
- 39 I. Katayama, H. Aoki, J. Takeda, H. Shimosato, M. Ashida, R. Kinjo, I. Kawayama, M. Tonouchi, M. Nagai, and K. Tanaka, Phys. Rev. Lett. 108 , 097401 (2012).
- 40 G. Orenstein, V. Krapivin, Y. Huang, Z. Zhan, G. de la Pena Munoz, R. A. Duncan, Q. Nguyen, J. Stanton, S. Teitelbaum, H. Yavas, T. Sato, M. C. Hoffmann, P. Kramer, J. Zhang, A. Cavalleri, R. Comin, M. P. M. Dean, A. S. Disa, M. Forst, S. L. Johnson, M. Mitrano, A. M. Rappe, D. Reis, D. Zhu, K. A. Nelson, and M. Trigo, 'Observation of polarization density waves in srtio3,' (2024), arXiv:2403.17203 [cond-mat.mtrl-sci].
- 41 B. Cheng, P. L. Kramer, Z.-X. Shen, and M. C. Hoffmann, Phys. Rev. Lett. 130 , 126902 (2023).
- 42 F. Libbi, A. Johansson, B. Kozinsky, and L. Monacelli, 'Ultrafast quantum dynamics in srtio 3 under impulsive thz radiation,' (2024), arXiv:2408.12421 [cond-mat.mtrl-sci].

- 43 J. Vandermause, S. B. Torrisi, S. Batzner, Y. Xie, L. Sun, A. M. Kolpak, and B. Kozinsky, npj Computational Materials 6 , 20 (2020).
- 44 A. Johansson, Y. Xie, C. J. Owen, J. S. Lim, L. Sun, J. Vandermause, and B. Kozinsky, (2022), 10.48550/ARXIV.2204.12573, arXiv:2204.12573 [physics.comp-ph].
- 45 A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. In't Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, et al. , Computer Physics Communications 271 , 108171 (2022).
- 46 P. Giannozzi, S. Baroni, N. Bonini, M. Calandra, R. Car, C. Cavazzoni, D. Ceresoli, G. L. Chiarotti, M. Cococcioni, I. Dabo, A. D. Corso, S. de Gironcoli, S. Fabris, G. Fratesi, R. Gebauer, U. Gerstmann, C. Gougoussis, A. Kokalj, M. Lazzeri, L. Martin-Samos, N. Marzari, F. Mauri, R. Mazzarello, S. Paolini, A. Pasquarello, L. Paulatto, C. Sbraccia, S. Scandolo, G. Sclauzero, A. P. Seitsonen, A. Smogunov, P. Umari, and R. M. Wentzcovitch, Journal of Physics: Condensed Matter 21 , 395502 (2009).
- 47 G. Prandini, A. Marrazzo, I. E. Castelli, N. Mounet, and N. Marzari, npj Computational Materials 4 , 72 (2018), http://materialscloud.org/sssp.
- 48 K. Yabana, T. Sugiyama, Y. Shinohara, T. Otobe, and G. F. Bertsch, Phys. Rev. B 85 , 045134 (2012).
Francesco Libbi Anders Johansson Lorenzo Monacelli Boris Kozinsky
Atomistic simulations of out-of-equilibrium quantum nuclear dynamics
The rapid advancements in ultrafast laser technology have paved the way for
pumping and probing the out-of-equilibrium dynamics of nuclei in crystals.
However, interpreting these experiments is extremely challenging due to the
complex nonlinear responses in systems where lattice excitations interact,
particularly in crystals composed of light atoms or at low temperatures where
the quantum nature of ions becomes significant. In this work, we address the
nonequilibrium quantum ionic dynamics from first principles. Our approach is
general and can be applied to simulate any crystal, in combination with a
first-principles treatment of electrons or external machine-learning
potentials. It is implemented by leveraging the nonequilibrium time-dependent
self-consistent harmonic approximation (TD-SCHA), with a stable,
energy-conserving, correlated stochastic integration scheme that achieves an
accuracy of $\mathcal{O}(dt^3)$. We benchmark the method with both a simple
one-dimensional model to test its accuracy and a realistic 40-atom cell of
SrTiO3 under THz laser pump, paving the way for simulations of ultrafast
THz-Xray pump-probe spectroscopy like those performed in synchrotron
facilities.