  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.
