# Time-Dependent Self Consistent Harmonic Approximation:
Anharmonic nuclear quantum dynamics and time correlation functions

Lorenzo Monacelli  Francesco Mauri  Physics Department, University of Rome
“Sapienza”

###### Abstract

Most material properties of great physical interest are directly related to
nuclear dynamics, e.g. the ionic thermal conductivity, Raman/IR vibrational
spectra, inelastic X-ray, and Neutron scattering. A theory able to compute
from first principles these properties, fully accounting for the quantum
nature of the nuclei and the anharmonicity in the nuclear energy landscape
that can be implemented in systems with hundreds of atoms is missing. Here, we
derive an approximated theory for the quantum time evolution of lattice
vibrations at finite temperature. This theory introduces the time dynamics in
the Self-Consistent Harmonic Approximation (SCHA) and shares with the static
case the same computational cost. It is nonempirical, as pure states evolve
according to the Dirac least action principle and the dynamics of the thermal
ensemble conserves both energy and entropy. The static SCHA is recovered as a
stationary solution of the dynamical equations. We apply perturbation theory
around the static SCHA solution and derive an algorithm to compute efficiently
quantum dynamical correlation functions. Thanks to this new algorithm, we have
access to the response function of any general external time-dependent
perturbation, enabling the simulation of phonon spectra without following any
perturbative expansion of the nuclear potential or empirical methods. We
benchmark the method on the IR and Raman spectroscopy of high-pressure
hydrogen phase III, with a simulation cell of 96 atoms. Our work also explores
the nonlinear regime of the dynamical nuclear motion, providing a paradigm to
simulate the interaction with intense or multiple probes, as in pump-probe
spectroscopy, or chemical reactions involving light atoms, as the proton
transfer in biomolecules.

††preprint: APS/123-QED

##  I  Introduction

The computational power available for research dramatically increased in the
last decades paving the way to the birth of a new field of science: material
design. We can predict the physical properties of a material _in silico_ ,
anticipating the experimental data and enabling for an automatic search of
target compounds before their synthesis. An example of an impressive result is
the prediction of record-breaking high-temperature superconductivity in
hydrates, like hydrogen sulfide[1] and LaH10[2] that came together with
experimental result[3, 4, 5], thanks to the success of _ab initio_ crystal
structure prediction[6, 7, 8].

Material science can also assist in interpreting experimental data, as it is
possible, in principle, to anticipate the outcome of almost any experimental
technique by calculating dynamical correlation functions. For example, IR and
optical spectroscopy are obtained by simulating the dipole-dipole time
correlation functions, Raman signal with the polarizability-polarizability
dynamical correlation function, while Neutron and X-Ray scattering can be
simulated calculating the dynamical structure factor.

However, still, first-principle calculations are far from being perfect, and
they often miss the required precision. One of the most difficult aspect to
describe are ionic vibrations. Even if the energy of ionic vibration is orders
of magnitude smaller than the typical energy involved in chemical bonds, the
nuclear excited states are in the same range of energy as thermal excitation
available at room temperature. Therefore, ionic motion is responsible for
almost all properties of materials that depend on temperature as thermal
expansion and thermal conductivity.

While the development of time-dependent Hartree-Fock (TD-HF) and time-
dependent density functional theory[9, 10] (TD-DFT) paved the way to the
simulation of these quantum correlation functions fully _ab initio_ for
electrons, still a full quantum theory of nuclear motion that can be routinely
applied in realistic systems (with hundreds of atoms in the simulation cell)
is missing.

The usual approach to lattice dynamics is through perturbation theory, where
the ionic motion is expanded as a small perturbation around the static
position (minimum of the total electronic energy of the system with frozen
ions). Needless to say, this approach fails in many cases, as when ions are
located in a saddle point of the energy landscape and their position is
stabilized by entropy, or if a very strong anharmonicity is present.

Ionic time correlation functions can be computed also in the presence of
strong anharmonicity with _ab initio_ molecular dynamics (AIMD)[11, 12], but
it neglects nuclear quantum dispersion. This is critical for molecular
crystals or systems with light atoms, where the Debye temperature comparable
or above room temperature. In these cases, AIMD provides untrustable results.

The quantum equivalent of AIMD is path-integral molecular dynamics (PIMD), but
this is a static theory, in which time is replaced by an imaginary variable.
Thus, PIMD describes only equilibrium quantities. It is possible to obtain
dynamical correlation function by performing an analytical continuation to
real time[13], but this is computationally heavy and achievable only in
prototypical cases where the correlation function has few poles and if a
force-field is available to clean stochastic noise[14]. Empirical methods have
been developed[15, 16, 17, 18], but, still, their application is
computationally expensive as it requires a PIMD calculation.

Relevant cases where all these techniques fail are materials with light atoms
or close to structural instability, where ionic fluctuations are sizable. For
example, the simulation of high-pressure hydrates requires a full quantum and
anharmonic treatment of the nuclear motion. These materials are attracting a
lot of attention due to the discovery of room-temperature
superconductivity[19]. On the other side, systems close to a structural
instability play a fundamental role in technological and industrial
applications; among them, we have multiferroics, charge density waves in 2D
layers, and ferroelectrics, where the huge anharmonic phonon scattering is
exploited also to increase their thermoelectric efficiency. Moreover, quantum
nuclear dynamics beyond the linear regime is fundamental for studying chemical
reactions involving light atoms. A significant example is the proton transfer
in biomolecules.

In this work, we formulate a new dynamical theory of quantum nuclear motion at
finite temperature and derive an algorithm to calculate time-dependent
correlation functions fully _ab initio_. This method has a computational cost
of the same order of standard AIMD but fully accounts for the effect of
quantum fluctuations non empirically. Our theory is a time-dependent (TD)
extension of the Stochastic Self-Consistent Harmonic Approximation (SCHA),
that already proved to be very efficient in describing equilibrium properties
of many materials, as high-pressure hydrates[20, 21, 22, 23], hydrogen[24, 25,
26], charge density waves[27, 28, 29] and thermoelectric materials[30, 31,
32].

In Sec. II we introduce the theoretical framework of this paper. We revise the
general dynamical response to an external time-dependent perturbation and the
relationship between the nuclear time-correlation functions to some common
experimental techniques. In Sec. III we revise the static SCHA, redefining it
in a contest in which it is easier to link with the full time-dependent
theory. The Time-Dependent SCHA (TD-SCHA) is introduced in Sec. IV, where we
derive the master equation for the density matrix by imposing that a Gaussian
wave-packet must minimize the Dirac action. We also discuss generic features
of the time-evolution as energy and entropy conservation. We also discuss
steady-state solutions and the thermodynamic equilibrium (where the SCHA
solution is recovered). In Sec. V, we linearize the TD-SCHA equations of
motion around the equilibrium solution. We derive the explicit formula to
compute the quantum time-correlation functions and derive a new Lanczos based
algorithm to compute them. We show in Sec. V.4 that we recover the dynamical
_ansatz_ formulated by Bianco et al[33] for one-phonon Green functions. In
Sec. VII.2, we benchmark the method in the complex case of high-pressure
hydrogen phase III, where quantum dispersion and anharmonicity strongly affect
the vibrational spectra. We prove, comparing the simulation of IR and Raman
spectra to experiments, how the TD-SCHA can easily handle such a complex
system with 96 atoms in the simulation cell computing forces between ions
fully _ab initio_ within Density Functional Theory (DFT).

###  I.1 Notation

First of all, let us introduce the notation used in the paper.

We use bold fonts to indicate vectors or matrices; products between them are
the standard rows-by-columns product. Each index indicates both atomic and
Cartesian components, so, if not differently specified, it ranges from 1 to
$3N$ ($N$ is the number of atoms). Therefore, vectors and matrices have
dimension $3N$ and $3N\times 3N$, respectively. For quantum operators and
wave-functions, we use the Dirac notation: we distinguish quantum operators
with a hat $\hat{\cdot}$ from real numbers. So, for example, $R_{i}$ is the
position of the $i$-th atom/Cartesian component, $\hat{R}_{i}$ is quantum
position operator and $\ket{R_{i}}$ is the eigenvector of the position
operator with $R_{i}$ as eigenvalue:

$\hat{R}_{i}\ket{R_{i}}=R_{i}\ket{R_{i}}$ (1)

The bracket $\left\langle\cdot\right\rangle_{\rho}$ indicates an average on
the distribution specified in the lower index. The average is a standard
average on the nuclear probability distribution $\rho(\bm{R})$ if we do not
indicate the operator, otherwise it is a quantum average computed on the
density matrix:

$\left\langle O\right\rangle_{\rho}=\int d\bm{R}O(\bm{R})\rho(\bm{R})$ (2)
$\left\langle\hat{O}\right\rangle_{\hat{\rho}}=\textup{Tr}\left[\hat{\rho}\hat{O}\right]$
(3)

Moreover, if the density matrix, the wave-function and any parameter that
represent the status of the system does not have an explicit time dependency,
it represents the equilibrium solution, unless differently stated.

##  II The quantum nuclear evolution

Almost all experimental techniques used to characterize materials involve an
exchange in energy between the probe and the sample below hundreds of
$\text{\,}\mathrm{eV}$. At this energy scale, the only components of the
material that interact with the probe are valence electrons and ions, which
include the atomic nucleus plus the core electrons. These are the degrees of
freedom we are interested in studying. In the Born-Oppenheimer (BO)
approximation, electrons are fast moving with respect to ions. For this
reason, the electrons feel the ionic lattice frozen, and relax to their ground
state for each ionic configuration. Thus, the full electron-ion wave-function
$\ket{\Psi}$ can be factorized in a ionic and electronic part:

$\braket{\bm{R},{\bm{r}}}{\Psi}=\braket{\bm{R}}{\psi_{\text{ion}}}\braket{{\bm{r}}}{\psi_{\text{el}}[\bm{R}]}$
(4)

where ${\bm{r}}$ are the electron positions, $\ket{\psi_{\text{ion}}}$ is the
wave-function of ions and $\ket{\psi_{\text{el}}[\bm{R}]}$ is the electronic
ground state with the ions fixed in the $\bm{R}$ position.

In this work, we focus in physical properties that depends only on ions. For
this reason, we will drop the ${}_{\text{ion}}$ index from the wave-function,
and $\ket{\psi}$ always refers to the nuclear wave-function
$\ket{\psi_{\text{ion}}}$. Within BO approximation, the ionic wave-function
obeys to the Scroedinger equation with a BO Hamiltonian ${\hat{H}}$:

${\hat{H}}=\sum_{i=1}^{3N}\frac{\hat{p_{i}}^{2}}{2m_{i}}+V(\hat{\bm{R}})$ (5)

where $\hat{p}_{i}$ is the momentum operator of nucleus $i$, $m_{i}$ the mass
of the $i$-th atom. The BO potential $V(\bm{R})$ is the ionic energy
landscape, obtained as the ground state energy of the electronic problem with
ions fixed in positions $\bm{R}$. The BO energy landscape $V(\bm{R})$ is
usually obtained _ab initio_ by solving the electronic problem with fixed
nuclei within DFT or HF approximations.

This work aims to describe how (nuclear) physical properties change in time
when the system in equilibrium is perturbed by an external probe. This is the
typical setup for any experiment: when $t<t_{0}$, the system is in
equilibrium. At $t=t_{0}$, it starts interacting with an external time-
dependent perturbation $V^{(1)}(\bm{R},t)$. We are interested in the
expectation value of a generic nuclear observable $\hat{\mathcal{A}}$ at time
$t>t_{0}$.

At $T=$0\text{\,}\mathrm{K}$$, the nuclear wave-function in equilibrium
satisfies the static Schroedinger equation:

${\hat{H}}\ket{\psi}=E_{\text{GS}}\ket{\psi}$ (6)

The presence of a time-dependent perturbation for $t\geq t_{0}$ introduces a
time dependency in the nuclear wave-function $\ket{\psi(t)}$. The time-
dependent external perturbation is:

$V^{(1)}(\hat{\bm{R}},t)={\mathcal{B}}(\hat{\bm{R}})\mathcal{V}(t)$ (7)
$\mathcal{V}(t)=0\qquad\text{for}\qquad t<t_{0}$ (8)

where $\hat{\mathcal{B}}={\mathcal{B}}(\hat{\bm{R}})$ is the observable that
couples the external perturbation (the probe) with the ionic positions and
$\mathcal{V}(t)$ is the time envelope of our perturbation. Thus, the quantum
expectation value of the $\hat{\mathcal{A}}$ observable at time $t$ is:

${\mathcal{A}}(t)=\braket{\psi(t)}{\hat{\mathcal{A}}}{\psi(t)}$ (9)

The Schroedinger equation governs the time-evolution of $\ket{\psi(t)}$:

$i\hbar\frac{d}{dt}\ket{\psi(t)}={{\hat{H}}_{\text{td}}(t)}\ket{\psi(t)}.$
(10)

where the ${{\hat{H}}_{\text{td}}(t)}$ time-dependent Hamiltonian is:

${{\hat{H}}_{\text{td}}(t)}={\hat{H}}+V^{(1)}(\hat{\bm{R}},t),$ (11)

and Eq. (10) is solved with the initial condition that the wave-function at
$t=t_{0}$ is in the ground state:

$\ket{\psi(t_{0})}=\ket{\psi}.$ (12)

At finite temperature ($T>$0\text{\,}\mathrm{K}$$), we just need to replace
the wave-function with a density matrix $\hat{\rho}$ that can describe also
mixture of states. At equilibrium (when $t\leq t_{0}$), the density matrix is
${\hat{\rho}}$:

${\hat{\rho}}=\sum_{i}p_{i}\ket{\psi_{i}}\bra{\psi_{i}}.$ (13)

where the $\ket{\psi_{i}}$ are the eigenstates of the ${\hat{H}}$ Hamiltonian
with $E_{i}$ energy, and $p_{i}$ are the Boltzmann occupations ($k_{b}$ is the
Boltzmann constant and $T$ the temperature):

${\hat{H}}\ket{\psi_{i}}=E_{i}\ket{\psi_{i}}$ (14) $p_{i}=\frac{e^{-\beta
E_{i}}}{Z}\qquad Z=\sum_{i}e^{-\beta E_{i}}\qquad\beta=\frac{1}{k_{b}T}$ (15)

After $t_{0}$, the system interacts with the external perturbation, and the
density matrix depends on time $\hat{\rho}(t)$. This time evolution is simply
given by the evolution of each state in the mixture according to the
Schroedinger equation:

$\hat{\rho}(t)=\sum_{i}p_{i}\ket{\psi_{i}(t)}\bra{\psi_{i}(t)}.$ (16)

Where $\ket{\psi_{i}(t)}$ are the same eigenstates of the BO Hamiltonian
${\hat{H}}$ (Eq. 14) at $t=t_{0}$ evolving with the time dependent Hamiltonian
${{\hat{H}}_{\text{td}}(t)}$. In this work, we focus on isolated systems,
therefore $p_{i}$ does not depend on time, and the density matrix satisfy the
Liuoville-Von Neumann equation:

$i\hbar\frac{d}{dt}\hat{\rho}(t)={{\hat{H}}_{\text{td}}(t)}\hat{\rho}(t)-\hat{\rho}(t){{\hat{H}}_{\text{td}}(t)}$
(17)

In this case, the average of the observable is:

${\mathcal{A}}(t)=\left\langle\hat{\mathcal{A}}\right\rangle_{\hat{\rho}(t)}=\textup{Tr}\left[\hat{\rho}(t)\hat{\mathcal{A}}\right]$
(18)

A case of particular interest is the linear regime (when the interaction
between the probe and the system is small compared to the typical interaction
between the ions). In this case, the time-dependence of an observable is given
by the response function:

${\mathcal{A}}(t)={\mathcal{A}}(t_{0})+\int_{-\infty}^{\infty}\chi_{{\mathcal{A}}{\mathcal{B}}}(t-t^{\prime})\mathcal{V}(t^{\prime})\,dt^{\prime}$
(19)

The response function $\chi_{{\mathcal{A}}{\mathcal{B}}}(t)$ is directly
related to time-correlation function through the Kubo equation:

$\chi_{{\mathcal{A}}{\mathcal{B}}}(t)=-\frac{i}{\hbar}\left\langle[e^{\frac{i}{\hbar}\hat{H}t}\hat{\mathcal{A}}e^{-\frac{i}{\hbar}\hat{H}t},\hat{\mathcal{B}}]\right\rangle_{{\hat{\rho}}^{(0)}}\vartheta(t)$
(20)

Here, $e^{-\frac{i}{\hbar}\hat{H}t}$ is the time evolution operator in absence
of the external perturbation and $\vartheta(t)$ is the Heaviside function. In
particular, the convolution of Eq. (19) between
$\chi_{{\mathcal{A}}{\mathcal{B}}}$ and $\mathcal{V}$ becomes a simple product
in the frequency domain:

${\mathcal{A}}(\omega)-{\mathcal{A}}(t_{0})=\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)\mathcal{V}(\omega).$
(21)

The $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$ operators depend on the
particular experiment we want to simulate. For example, IR spectroscopy
studies the absorption of IR light across the sample. Absorbed light is
related to the imaginary part of the refractive index $n(\omega)$, that is
obtained from the dielectric tensor $\epsilon(\omega)$:

$n(\omega)=\sqrt{\epsilon(\omega)}$ (22)

The dielectric tensor of a material is directly related to the polarization-
polarization susceptibility $\chi_{MM}(\omega)$, for which we have an
electronic and ionic contribution. The electronic contribution is almost
always constant in the IR light for an insulator, and the features of IR
spectroscopies (peaks) are given by the ionic contribution to the dipole-
dipole response function:

$n(\omega)=\sqrt{1+4\pi\chi^{(el)}_{MM}(\omega)+4\pi\chi^{(ion)}_{MM}(\omega)}.$
(23)

Therefore, to simulate IR spectra, we need to compute the dynamical nuclear
correlation function $\chi_{MM}^{(ion)}(\omega)$, where
$\hat{\mathcal{A}}=\hat{\mathcal{B}}=M(\hat{\bm{R}})$, and $M(\hat{\bm{R}})$
is the dipole moment along the probe polarization of the system when nuclei
are located in $\bm{R}$. Analogously, the Raman signal is obtained measuring
the energy exchanged between the incoming and outcoming radiation. This energy
exchange occurs thanks to the polarizability tensor $\bm{\alpha}(\bm{R})$ of
the system caused by the displacements in the ionic position $\bm{R}$, induced
by the probe. In this case,
$\hat{\mathcal{A}}=\hat{\mathcal{B}}=\alpha_{xy}(\hat{\bm{R}})$, where the
$xy$ directions are the incoming and outcoming polarization of light[34]. In
the same way, the dynamical structure factor and the thermal conductivity are
dynamical nuclear correlation functions[35].

Indeed, we can also simulate experiments beyond linear regime, in which the
system interacts with multiple probes, as in impulsive vibrational
spectroscopy[36]. Here, the system is perturbed with two pulses: a first
optical pulse at $t=t_{0}$ brings the system out-of-equilibrium, then, the
probe pulse interacts with the sample at time $t_{1}>t_{0}$ measuring its
polarizability.

The numerical simulation of such an experiment requires the solution of Eq.
(17) and Eq. (18). This is computationally unfeasible for a system with more
than few atoms, as the memory required to store the density matrix goes as
$M^{(3N)^{2}}$, where $M$ is the dimension of the basis of the wave-function
and $N$ the number of atoms in the simulation cell. Also, the easier linear
regime requires the calculation of the time-correlation function in Eq. (20),
that involves the full diagonalization of the static interacting Hamiltonian
${\hat{H}}$ to calculate the free time-evolution
$e^{-\frac{i}{\hbar}\hat{H}t}$.

In this work, we derive an approximated theory for the time evolution of the
nuclear density matrix. Our theory allows simulating the time evolution of
systems with hundreds of atoms with the BO energy landscape $V(\hat{\bm{R}})$
calculated _ab initio_ from the solution of the electronic problem. For the
case of linear perturbations, we also derive a new algorithm to calculate the
nuclear time-correlation functions $\chi_{{\mathcal{A}}{\mathcal{B}}}(t)$. We
present an _ab initio_ application of a system with 96 atoms in the unit cell.

In the next section, we revise the static SCHA theory, a theory for the
equilibrium density matrix that is the starting point of our derivation.

##  III Static Self-Consistent Harmonic Approximation

The SCHA is a mean-field theory developed to deal with interacting phonons and
to compute thermodynamic properties of solids.

Here, the interacting many-body ionic density matrix is replaced with a non-
interacting one. Starting from now through the rest of the paper,
${{\hat{\rho}}}$ is the SCHA equilibrium density matrix. As done in Hartree-
Fock for electrons, in the SCHA, the ${{\hat{\rho}}}$ density matrix minimizes
the free energy functional of a trial density matrix $\hat{\tilde{\rho}}$:

$F[{{\hat{\rho}}}]=\min_{\hat{\tilde{\rho}}}F[\hat{\tilde{\rho}}].$ (24)
$F[\hat{\tilde{\rho}}]=\left\langle\hat{H}\right\rangle_{\tilde{\hat{\rho}}}-TS[\hat{\tilde{\rho}}],$
(25)

where $S[\hat{\tilde{\rho}}]$ is the entropy functional:

$S[\hat{\tilde{\rho}}]=-k_{b}\textup{Tr}\left[\hat{\tilde{\rho}}\ln\hat{\tilde{\rho}}\right].$
(26)

The $\hat{\tilde{\rho}}$ is restricted only among Gaussian matrices:

$\displaystyle\braket{\bm{R}}{\hat{\tilde{\rho}}}{\bm{R}^{\prime}}={\mathcal{N}}\exp$
$\displaystyle\bigg{[}-\sum_{ab}\frac{\Theta_{ab}}{4}(R_{a}-\mathcal{R}_{a})(R_{b}-\mathcal{R}_{b})+$
$\displaystyle-\sum_{ab}\frac{\Theta_{ab}}{4}(R_{a}^{\prime}-\mathcal{R}_{a})(R_{b}^{\prime}-\mathcal{R}_{b})+$
$\displaystyle+\sum_{ab}A_{ab}(R_{a}-\mathcal{R}_{a})(R_{b}^{\prime}-\mathcal{R}_{b}^{\prime})\bigg{]}$
(27)

The parameters that uniquely determines $\hat{\tilde{\rho}}$ are the vector
$\bm{\mathcal{R}}$, that represent the average ionic positions, and the real
Hermitian matrices $\bm{\Theta}$ and $\bm{A}$, that represent the quantum and
thermal fluctuations around the average positions, respectively. The
${\mathcal{N}}$ factor is the normalization of the density matrix. The
$\bm{\Theta}$ and $\bm{A}$ matrices are not independent: they commute and are
constrained so that $\hat{\tilde{\rho}}$ must be normalized. This condition is
obtained if the real symmetric matrix $\bm{\Upsilon}$ has only positive
eigenvalues:

$\bm{\Upsilon}=\bm{\Theta}-2\bm{A}$ (28)

Eq. (27) is a bit more general than the one introduced in ref.[20]: it is the
density matrix of a quantum harmonic oscillator, where each normal mode is
thermally populated by a different temperature. However, as we prove later in
Sec. IV.5 and Appendix D, in the minimum of the free energy, all the modes are
always populated by the same temperature, and the result coincides with
ref.[20].

The SCHA is solved by substituting the expression of the trial density matrix
$\hat{\tilde{\rho}}$ (Eq. 27) into the free energy functional (Eq. 25), and
minimizing with respect to $\bm{\mathcal{R}}$, $\bm{\Theta}$ and $\bm{A}$.

Interestingly, the SCHA equilibrium density matrix ${{\hat{\rho}}}$ satisfies
the self-consistent equation:

${{\hat{\rho}}}=\frac{\exp\left(-\beta{\hat{\mathcal{H}}}[{\rho}]\right)}{Z[{\rho}]},$
(29)

where ${\hat{\mathcal{H}}}[{\rho}]$ is a harmonic Hamiltonian that depends
self-consistently on the nuclear equilibrium density ${\rho}(\bm{R})$:

${\rho}(\bm{R})=\braket{\bm{R}}{{{\hat{\rho}}}}{\bm{R}}$ (30)

${\hat{\mathcal{H}}}[{\rho}]=\sum_{i=1}^{3N}\frac{{\hat{p}_{i}}^{2}}{2m_{i}}+\sum_{ij}\left\langle\frac{\partial^{2}V}{\partial
R_{i}\partial
R_{j}}\right\rangle_{{\rho}}(\hat{R}_{i}-\mathcal{R}_{i}[{\rho}])(\hat{R}_{j}-\mathcal{R}_{j}[{\rho}]).$
(31)

The functional $\mathcal{R}_{i}[{\rho}]$ is the average position of the $i$-th
atom, and coincides with the $\mathcal{R}_{i}$ solution of the SCHA:

$\bm{\mathcal{R}}[{\rho}]=\int d\bm{R}\;{\rho}(\bm{R})\bm{R},$ (32)

and the averages are computed with the $\rho(\bm{R})$ probability density

$\left\langle\frac{\partial^{2}V}{\partial R_{a}\partial
R_{b}}\right\rangle_{{\rho}}=\int
d\bm{R}{\rho}(\bm{R})\frac{\partial^{2}V(\bm{R})}{\partial R_{a}\partial
R_{b}}.$ (33)

When $\hat{\tilde{\rho}}$ is the equilibrium solution ${{\hat{\rho}}}$, we
have:

$A_{ab}=\sqrt{m_{a}m_{b}}\sum_{\mu}\frac{2\omega_{\mu}n_{\mu}(n_{\mu}+1)}{\hbar(2n_{\mu}+1)}e_{\mu}^{a}e_{\mu}^{b}$
(34)
$\Upsilon_{ab}=\sqrt{m_{a}m_{b}}\sum_{\mu}\frac{2\omega_{\mu}}{\hbar(2n_{\mu}+1)}e_{\mu}^{a}e_{\mu}^{b}$
(35)

Where $\bm{e_{\mu}}$ and $\omega_{\mu}$ are the normal modes and frequencies
of the self-consistent harmonic Hamiltonian ${\hat{\mathcal{H}}}[\rho]$, and
$n_{\mu}$ is the Bose-Einstein occupation number:

$n_{\mu}=\frac{1}{e^{\beta\hbar\omega_{\mu}}-1}$ (36)

The SCHA in its stochastic implementation[20, 33, 37] proved to be a practical
tool for many equilibrium calculations thanks to its cheap numerical cost.
Both the kinetic energy and the entropy of Eq. (25) are analytical when a
Gaussian density matrix is employed. The only remaining term is the average of
the potential energy on ionic configurations. This can be efficiently
calculated with a Monte Carlo integral, sampling stochastically the ionic
configurations, distributed according to $\rho(\bm{R})$.

A particular case is the $T=$0\text{\,}\mathrm{K}$$ limit, when the
equilibrium SCHA density matrix is a pure state. Here, the SCHA equilibrium
pure state $\ket{\psi}$ is the ground state of the self-consistent Hamiltonian
${\hat{\mathcal{H}}}[\rho]$:

${\hat{\mathcal{H}}}[\rho]\ket{\psi}=E_{GS}\ket{\psi}$ (37)

and the density $\rho(\bm{R})$ is:

$\rho(\bm{R})=\braket{\bm{R}}{\psi}\braket{\psi}{\bm{R}}=\left|\braket{\bm{R}}{\psi}\right|^{2}$
(38)

The pure state is a Gaussian wave-packet:

$\braket{\bm{R}}{\psi}={\mathcal{N}}^{\frac{1}{2}}\exp\left[-\frac{1}{4}\sum_{ab}\Upsilon_{ab}(R_{a}-\mathcal{R}_{a})(R_{b}-\mathcal{R}_{b})\right]$
(39)

##  IV Time-Dependent Self-Consistent Harmonic Approximation

The SCHA introduced in Sec. III is a static theory: it cannot describe
dynamical properties, like phonons observed experimentally.

In this section, we derive a new theory, the Time-Dependent Self-Consistent
Harmonic Approximation (TD-SCHA), to describe correctly, without any empirical
approximation, the response of the system to any (small or not) external time-
dependent probe that interacts with ions.

We start from a pure state, i.e. the equilibrium solution at
$T=$0\text{\,}\mathrm{K}$$. At time $t=t_{0}$, we switch on a perturbation
$V^{(1)}(\bm{R},t)$ and the overall Hamiltonian becomes
${{\hat{H}}_{\text{td}}(t)}(t)$ defined in Eq. (11).

We constrain the wave-packet to the most general time-dependent Gaussian:

$\displaystyle\braket{\bm{R}}{\psi(t)}={\mathcal{N}}^{\frac{1}{2}}(t)\exp\bigg{(}i\sum_{a}Q_{a}(t)[R_{a}-\mathcal{R}_{a}(t)]+$
$\displaystyle-\sum_{ab}\left[\frac{\Theta_{ab}(t)}{4}-iC_{ab}(t)\right][R_{a}-\mathcal{R}_{a}(t)][R_{b}-\mathcal{R}_{b}(t)]\bigg{)}.$
(40)

In Eq. (40) we have two new parameters with respect to the static solution of
Eq. (38): $\bm{Q}$ and $\bm{C}$. They add a complex phase to our wave-packet,
and represent the momentum of the $\bm{R}$ and $\bm{\Theta}$ variables. The
time-dependency of the $\bm{Q}$, $\bm{\mathcal{R}}$, $\bm{\Theta}$, and
$\bm{C}$ parameters is found by minimizing the Dirac action along the time-
evolution path:

$A=\frac{1}{t-t_{0}}\int_{t_{0}}^{t}\braket{\psi(t^{\prime})}{\hat{H}_{\text{td}}(t^{\prime})-i\hbar\frac{d}{dt^{\prime}}}{\psi(t^{\prime})}\,dt^{\prime}$
(41)

The full equations of motion are derived in Appendix C.

It can be proved (as we show in Sec. IV.1) that the same equations of motion
are obtained if Eq. (40) is evolved by a self-consistent Schroedinger
equation:

$i\hbar\frac{d}{dt}\ket{\psi(t)}={\hat{\mathcal{H}}}[\rho(t)]\ket{\psi(t)}$
(42)

The self-consistent Hamiltonian that defines the time-evolution in the TD-SCHA
is:

$\displaystyle{\hat{\mathcal{H}}}[$
$\displaystyle\rho(t)]=\sum_{i=1}^{N}\frac{{\hat{p}_{i}}^{2}}{2m_{i}}+\sum_{i}\left\langle\frac{\partial{V^{(\text{tot})}}}{\partial
R_{i}}\right\rangle_{\rho(t)}(\hat{R}_{i}-\mathcal{R}_{i}[\rho(t)])+$
$\displaystyle+\sum_{ij}\left\langle\frac{\partial^{2}{V^{(\text{tot})}}}{\partial
R_{i}\partial
R_{j}}\right\rangle_{\rho(t)}(\hat{R}_{i}-\mathcal{R}_{i}[\rho(t)])(\hat{R}_{j}-\mathcal{R}_{j}[\rho(t)]),$
(43)

where

${V^{(\text{tot})}}(\bm{R},t)=V(\bm{R})+V^{(1)}(\bm{R},t)$ (44)

The ${\hat{\mathcal{H}}}[\rho(t)]$ Hamiltonian depends on the $\rho(\bm{R},t)$
probability distribution of finding the ions in the $\bm{R}$ configurations at
time $t$:

$\rho(\bm{R},t)=\left|\braket{\bm{R}}{\psi(t)}\right|^{2}$ (45)

The time-dependent self-consistent Hamiltonian (Eq. 43) has one extra linear
term in $\hat{\bm{R}}$ compared to the static one (Eq. 31). This extra linear
term is zero when the self-consistency of the equilibrium SCHA is achieved, as
the average of the derivative of the BO potential (forces) on the equilibrium
SCHA distribution is a necessary condition for the SCHA self-consistency[20].

The nuclear self-consistent Schroedinger equation (Eq. 42) has the same shape
of other mean-field theories for electrons, as TD-HF or TD-DFT. It is worth to
notice that Eq. (42) minimizes the action only if the wave-function is a
Gaussian wave-packet. Notably, as we show in the next section, if $\ket{\psi}$
is Gaussian, Eq. (42) is a closed equation: a Gaussian wave-packet evolving in
a general self-consistent harmonic Hamiltonian keeps its Gaussian form.

###  IV.1 Dynamics of a mixture of states

Eq. (42) describes the dynamics just of pure states. The equilibrium solution
of the SCHA is a pure state only if $T=$0\text{\,}\mathrm{K}$$.

We can extend the nuclear time-dependent evolution of a mixture of states by
replacing the time-dependent Scroedinger equation (Eq. 42) with the Liouville-
von Neumann equation (Eq. 17), as usually done in TD-DFT[38] and TD-HF:

$i\hbar\frac{d}{dt}\hat{\rho}(t)={\hat{\mathcal{H}}}[\rho(t)]\hat{\rho}(t)-\hat{\rho}(t){\hat{\mathcal{H}}}[\rho(t)]$
(46) $\rho(\bm{R},t)=\braket{\bm{R}}{\hat{\rho}(t)}{\bm{R}}$ (47)

Thanks to Eq. (46), we can describe the dynamics also of mixtures of states,
starting from the equilibrium SCHA solution at any temperature. Moreover, as
Eq. (42), this is a closed equation for a Gaussian wave-packet (as we show in
this section). We also prove in Sec. IV.2 and IV.3 that Eq. (46) correctly
conserves both the energy and the entropy, as expected from the correct
evolution of an isolated quantum system.

As we did for the pure state, we can get the explicit equation of motion for
the Gaussian mixed states wave-packet.

$\displaystyle\bra{\bm{R}^{\prime}}\hat{\rho}$
$\displaystyle(t)\ket{\bm{R}}={\mathcal{N}}(t)\exp\bigg{(}i\sum_{a}Q_{a}(t)(R_{a}^{\prime}-R_{a})+$
$\displaystyle-\sum_{ab}\left[\frac{\Theta_{ab}(t)}{4}-iC_{ab}(t)\right][R_{a}-\mathcal{R}_{a}(t)][R_{b}-\mathcal{R}_{b}(t)]+$
$\displaystyle-\sum_{ab}\left[\frac{\Theta_{ab}(t)}{4}+iC_{ab}(t)\right][R^{\prime}_{a}-\mathcal{R}_{a}(t)][R^{\prime}_{b}-\mathcal{R}_{b}(t)]$
$\displaystyle+\sum_{ab}A_{ab}(t)[R_{a}-\mathcal{R}_{a}(t)][R_{b}^{\prime}-\mathcal{R}_{b}(t)]\bigg{)}.$
(48)

With respect to the pure state introduced in Eq. (40), the time-dependent
density matrix (Eq. 48) has one more parameter: the $\bm{A}(t)$ complex
Hermitian matrix. $\mathcal{R}_{i}$ is the average position of the $i$-th atom
and coincides with the definition of Eq. (32). The variable $Q_{i}$ is a
linear phase modulation; multiplied by $\hbar$, it represents the momentum of
the $i$-th atom. In a flat potential (where ${\hat{\mathcal{H}}}$ only
contains the kinetic operator) the $i$-th atom average position drifts with
constant velocity $v_{i}=\frac{\hbar Q_{i}(t)}{m_{i}}$. The $\bm{Q}$ and
$\bm{\mathcal{R}}$ variables are similar to those of a classical molecular
dynamics. The $\bm{\Theta}$, $\bm{A}$ matrices describe the quantum and
thermal fluctuations. In particular, by looking at the diagonal elements of
the density operator (the density distribution, Eq. 45), we get the covariance
matrix of quantum-thermal fluctuations:

$\frac{\rho(\bm{R},t)}{{\mathcal{N}}(t)}=\exp\left\\{-\sum_{ab}\frac{1}{2}\Upsilon_{ab}(t)[R_{a}-\mathcal{R}_{a}(t)][R_{b}-\mathcal{R}_{b}(t)]\right\\}$
(49) $\Upsilon_{ab}(t)=\Theta_{ab}(t)-2\mathrm{Re}A_{ab}(t),$ (50)

where $\mathrm{Re}{}$ and $\mathrm{Im}{}$ identify the real and imaginary
part. $\bm{\Upsilon}$ is the inverse the covariance matrix of the Gaussian
distribution:

$\left(\bm{\Upsilon}^{-1}\right)_{ab}(t)=\left\langle[R_{a}-\mathcal{R}_{a}(t)][R_{b}-\mathcal{R}_{b}(t)]\right\rangle_{\rho(t)}.$
(51)

We give intuitive picture on the physical meaning of the parameters. The
$\bm{\Theta}$ matrix encodes pure quantum fluctuations and $\mathrm{Re}\bm{A}$
the thermal ones: if $\bm{A}=0$, Eq. (48) is a pure state and we recover Eq.
(40). $\bm{C}$ is a quadratic phase and represents the _chirp_ along the
quantum fluctuations. Its role is very similar to the chirp in signal
propagation, and represent a gradient in the speed of particles in different
positions in the wave-packet, as discussed in ref.[36]. On the other side, the
$\mathrm{Im}\bm{A}$ plays the role of the momentum for the thermal
fluctuations. It is nonzero only when there are more than 1 degree of freedom
($\mathrm{Im}\bm{A}$ is anti symmetric) and if the system is not in a pure
quantum state.

The ${\mathcal{N}}(t)$ is the density matrix normalization:

${\mathcal{N}}(t)=\sqrt{\frac{\det{\bm{\Upsilon}(t)}}{(2\pi)^{3N}}}$ (52)

Eq. (48) can be substituted in Eq. (46) to get the dynamical equations for the
parameters. For a convenient choice of the notation, it is better to express
the parameters rescaled by the masses. We indicate with a $\tilde{\cdot}$ the
rescaled matrices and vectors as:

$\tilde{C}_{ab}(t)=\frac{C_{ab}(t)}{\sqrt{m_{a}m_{b}}},\qquad\tilde{A}(t)=\frac{A_{ab}(t)}{\sqrt{m_{a}m_{b}}},$
(53a)
$\tilde{\Upsilon}_{ab}(t)=\frac{\Upsilon_{ab}(t)}{\sqrt{m_{a}m_{b}}},\qquad\tilde{\Theta}(t)=\frac{\Theta_{ab}(t)}{\sqrt{m_{a}m_{b}}},$
(53b)
$\tilde{\mathcal{R}}_{a}(t)=\sqrt{m}_{a}\mathcal{R}_{a}(t)\qquad\tilde{Q}_{a}(t)=\frac{Q_{a}(t)}{\sqrt{m}_{a}}$
(53c)

The final equation of motion are:

$\frac{d\tilde{\bm{\mathcal{R}}}}{dt}=\hbar\tilde{\bm{Q}}\qquad\frac{dQ_{a}}{dt}=\frac{\left\langle
f^{(\text{tot})}_{a}\right\rangle_{\rho(t)}}{\hbar}$ (54a)
$\frac{d\tilde{\bm{\Upsilon}}}{dt}=\hbar\left[\tilde{\bm{\Upsilon}}(2\tilde{\bm{C}}+\mathrm{Im}\tilde{\bm{A}})+(2\tilde{\bm{C}}-\mathrm{Im}\tilde{\bm{A}})\tilde{\bm{\Upsilon}}\right]$
(54b)
$\frac{d\mathrm{Re}\tilde{\bm{A}}}{dt}=\frac{\hbar}{2}\left(4\tilde{\bm{C}}\mathrm{Re}\tilde{\bm{A}}+4\mathrm{Re}\tilde{\bm{A}}\tilde{\bm{C}}-\tilde{\bm{\Theta}}\mathrm{Im}\tilde{\bm{A}}+\mathrm{Im}\tilde{\bm{A}}\tilde{\bm{\Theta}}\right)$
(54c)
$\frac{d\mathrm{Im}\tilde{\bm{A}}}{dt}=\frac{\hbar}{2}\left(4\tilde{\bm{C}}\mathrm{Im}\tilde{\bm{A}}+4\mathrm{Im}\tilde{\bm{A}}\tilde{\bm{C}}+\tilde{\bm{\Theta}}\mathrm{Re}\tilde{\bm{A}}-\mathrm{Re}\tilde{\bm{A}}\tilde{\bm{\Theta}}\right)$
(54d)
$\frac{d\tilde{\bm{C}}}{dt}=\frac{1}{2\hbar}\left\langle\frac{\partial^{2}{V^{(\text{tot})}}}{\partial\tilde{\bm{R}}\partial\tilde{\bm{R}}}\right\rangle_{\rho(t)}+\frac{\hbar}{2}\left[4\tilde{\bm{C}}^{2}-\frac{1}{4}\tilde{\bm{\Theta}}^{2}+\mathrm{Re}(\tilde{\bm{A}}\tilde{\bm{A}}^{\dagger})\right].$
(54e)

Here, we dropped the explicit time-dependency of these variables for brevity:
they represent the $\hat{\rho}(t)$ density matrix and not the equilibrium one.

The products of matrices is the standard rows-by-columns. The symbol
${\cdot}^{\dagger}$ after a matrix indicates the Hermitian conjugate, and
$\left\langle f^{(\text{tot})}_{a}\right\rangle_{\rho(t)}$ is the average of
BO forces on the $a$-th atom:

$f^{(\text{tot})}_{a}(\bm{R},t)=-\frac{\partial{V^{(\text{tot})}}(\bm{R},t)}{\partial
R_{a}},$ (55)

while the derivative with respect the $\tilde{\bm{R}}$ variable indicates the
mass rescale:

$\frac{\partial^{2}{V^{(\text{tot})}}}{\partial\tilde{R}_{a}\partial\tilde{R}_{b}}=\frac{1}{\sqrt{m_{a}m_{b}}}\frac{\partial^{2}{V^{(\text{tot})}}}{\partial
R_{a}\partial R_{b}}$ (56)

Eq. (54a) are the semi-classical equation of motion: they resemble the Newton
dynamics, but the force is averaged on the ionic probability distribution. The
dynamics preserves $\bm{\Upsilon}$, $\mathrm{Re}\bm{A}$ and $\bm{C}$ symmetric
and $\mathrm{Im}\bm{A}$ antisymmetric. The details of the derivation of Eq.
(54e) is reported in appendix A. These equations reduce to the evolution of
the pure state that minimizes the Dirac action, if we set $\bm{A}=0$ (see
Appendix C).

Notably, substituting the Gaussian wave-packet Eq. (48) into TD-SCHA equation
of motion Eq. (46), we get an extra condition on the time-dependency of the
${\mathcal{N}}(t)$ parameter.

$\frac{d{\mathcal{N}}(t)}{dt}=2\hbar{\mathcal{N}}(t)\textup{Tr}\left[\tilde{\bm{C}}(t)\right]$
(57)

This condition is automatically satisfied by Eq. (54e) if we substitute Eq.
(52) into Eq. (57). Since Eq. (57) does not depend explicitly on
$\hat{\bm{R}}$ operators and it is automatically satisfied by the equation of
motions, the Gaussian wave-packet evolution is closed with the TD-SCHA
equation (Eq. 46). In fact, the application of the Liouville operator on a
Gaussian density matrix:

${\hat{\mathcal{H}}}[\rho(t)]\hat{\rho}(t)-\hat{\rho}(t){\hat{\mathcal{H}}}[\rho(t)]$
(58)

with ${\hat{\mathcal{H}}}[\rho(t)]$ Harmonic, generates a polynomial of the
same order than the time derivative for the density matrix
$i\hbar\frac{d}{dt}\hat{\rho}(t)$. This means that, if the density matrix is
Gaussian at $t=t_{0}$ (equilibrium), it remains Gaussian for the whole time-
evolution, as Eq. (58) does not provide any term shifting
$i\hbar\frac{d}{dt}\hat{\rho}(t)$ from a Gaussian.

Very interestingly, the full dynamics is just obtained by standard rows-by-
columns product of small matrices (they are $3N\times 3N$, with $N$ the number
of atoms in the simulation cell). The only two quantities depending from the
physical system (the real BO Hamiltonian $\hat{H}$) are

$\left\langle\bm{f}^{(\text{tot})}\right\rangle_{\rho(t)}\qquad\left\langle\frac{\partial^{2}{V^{(\text{tot})}}}{\partial\tilde{\bm{R}}\partial\tilde{\bm{R}}}\right\rangle_{\rho(t)}.$
(59)

The calculation of these averages is also needed for a static SCHA calculation
and can be computed stochastically as described in [20, 33].

In particular, an efficient method to compute the average of the second
derivatives of the BO potential is obtained exploiting the methodology
introduced in ref.[33], that takes advantage of integration by parts:

$\left\langle\frac{\partial^{2}{V^{(\text{tot})}}}{\partial R_{a}\partial
R_{b}}\right\rangle_{\rho(t)}=\sum_{p}\Upsilon_{ap}\left\langle(R_{p}-\mathcal{R}_{p})f_{a}\right\rangle_{\rho(t)}$
(60)

In this way, only _ab initio_ forces are required. Therefore, the
implementation of the TD-SCHA equations has the same overall computational
cost as a static calculation.

Derivation of Gaussian wave-packet dynamics is discussed in literature
following different approaches[39, 40, 41, 42, 43, 44]. For example, refs.[42,
43, 44] project the real dynamics into the manifold of Gaussian states, while
refs.[39, 40, 41] exploit the Wentzel-Kramers-Brillouin (WKB) method to derive
semiclassical equations expanding the Schroedinger equation linearly around
$\hbar=0$.

###  IV.2 Energy conservation

Here we show that the TD-SCHA equation of motion (Eq. 46) satisfies energy
conservation when we switch off the time-dependent perturbation, as we expect
from a closed system. This is not trivial, as, for example, does not happen in
other methods for finite temperature dynamics of ions, as the finite
temperature multi-configurational time-dependent Hartree method[45].

The total energy is computed as the average of the BO Hamiltonian on the time-
dependent density matrix:

$E(t)=\left\langle\hat{H}\right\rangle_{\hat{\rho}(t)}=\operatorname{Tr}\left[\hat{\rho}(t)\hat{H}\right].$
(61)

By dividing into kinetic and potential contribution, the total energy is:

$E(t)=\left\langle\hat{K}\right\rangle_{\hat{\rho}(t)}+\left\langle\hat{V}\right\rangle_{\hat{\rho}(t)}$
(62)

The kinetic energy can be evaluated analytically directly from the density
matrix of Eq. (48):

$\displaystyle\left\langle\hat{K}\right\rangle_{\hat{\rho}(t)}=\frac{\hbar^{2}}{2}\bigg{[}$
$\displaystyle\operatorname{Tr}\bigg{(}4\tilde{\bm{C}}(\tilde{\bm{\Upsilon}})^{-1}\tilde{\bm{C}}-\mathrm{Im}\tilde{\bm{A}}(\tilde{\bm{\Upsilon}})^{-1}\mathrm{Im}\tilde{\bm{A}}\bigg{)}+$
$\displaystyle-2\operatorname{Tr}\bigg{(}\tilde{\bm{C}}(\tilde{\bm{\Upsilon}})^{-1}\mathrm{Im}\tilde{\bm{A}}\bigg{)}+$
$\displaystyle+\frac{1}{4}\operatorname{Tr}\tilde{\bm{\Upsilon}}+\operatorname{Tr}\tilde{\bm{A}}\bigg{]}+\sum_{i=1}^{n}\frac{\hbar^{2}Q_{i}^{2}}{2m_{i}}$
(63)

Also here, we dropped the exploit time dependency of the parameters that
represent $\hat{\rho}(t)$.

The last term is the kinetic energy of classical particles, the rest is the
contribution of the quantum and thermal fluctuations to the kinetic energy.

To prove that $\frac{d}{dt}E(t)=0$, we compute the derivative of the average
kinetic energy:

$\displaystyle\frac{d\left\langle\hat{K}\right\rangle_{\hat{\rho}(t)}}{dt}$
$\displaystyle=\frac{\hbar^{2}}{2}\bigg{[}\frac{1}{4}\operatorname{Tr}\frac{d\tilde{\bm{\Upsilon}}}{dt}+\sum_{i=1}^{3n}\frac{Q_{i}}{m_{i}}\frac{dQ_{i}}{dt}+$
$\displaystyle+2\operatorname{Tr}\left(4\tilde{\bm{C}}(\tilde{\bm{\Upsilon}})^{-1}\frac{d{\tilde{\bm{C}}}}{dt}-\tilde{\bm{C}}(\tilde{\Upsilon})^{-1}\frac{d\mathrm{Im}\tilde{\bm{A}}}{dt}\right)+$
$\displaystyle+\operatorname{Tr}\frac{d\mathrm{Re}\tilde{\bm{A}}}{dt}-2\operatorname{Tr}\left(\frac{d\tilde{\bm{C}}}{dt}(\tilde{\bm{\Upsilon}})^{-1}\mathrm{Im}\tilde{\bm{A}}\right)+$
$\displaystyle-2\operatorname{Tr}\left(\mathrm{Im}\tilde{\bm{A}}(\tilde{\bm{\Upsilon}})^{-1}\dot{\mathrm{Im}\tilde{\bm{A}}}+\tilde{\bm{C}}\frac{d(\tilde{\bm{\Upsilon}})^{-1}}{dt}\mathrm{Im}\tilde{\bm{A}}\right)+$
$\displaystyle+\operatorname{Tr}\left(4\tilde{\bm{C}}\frac{d(\tilde{\bm{\Upsilon}})^{-1}}{dt}\tilde{\bm{C}}-\mathrm{Im}\tilde{\bm{A}}\frac{d(\tilde{\bm{\Upsilon}})^{-1}}{dt}\mathrm{Im}\tilde{\bm{A}}\right)\bigg{]}.$
(64)

By substituting the equation of motion (Eq. 54e) we get the final derivative
of the kinetic energy:

$\displaystyle\frac{d\left\langle\hat{K}\right\rangle_{\hat{\rho}(t)}}{dt}$
$\displaystyle=\sum_{i}\frac{\hbar Q_{i}\left\langle
f_{i}\right\rangle_{\rho}}{m_{i}}+$
$\displaystyle+2\hbar\operatorname{Tr}\left[\left\langle\frac{\partial^{2}V}{\partial\tilde{\bm{R}}\partial\tilde{\bm{R}}}\right\rangle_{\rho(t)}(\tilde{\bm{\Upsilon}})^{-1}\tilde{\bm{C}}\right]+$
$\displaystyle-\hbar\operatorname{Tr}\left[\left\langle\frac{\partial^{2}V}{\partial\tilde{\bm{R}}\partial\tilde{\bm{R}}}\right\rangle_{\rho(t)}(\tilde{\bm{\Upsilon}})^{-1}\mathrm{Im}\tilde{\bm{A}}\right]$
(65)

In the same way, we can compute the derivative of the average potential:

$\displaystyle\frac{d\left\langle\hat{V}\right\rangle_{\hat{\rho}(t)}}{dt}=-\sum_{ab}\left\langle(R_{a}-\mathcal{R}_{a})(R_{b}-\mathcal{R}_{b})V\right\rangle_{\rho(t)}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}+$
$\displaystyle+\sum_{ab}\left\langle(R_{a}-\mathcal{R}_{a})V\right\rangle_{\rho(t)}\Upsilon_{ab}\frac{d\mathcal{R}_{b}}{dt}+\left\langle
V\right\rangle_{\rho(t)}\frac{1}{{\mathcal{N}}}\frac{d{\mathcal{N}}}{dt}$ (66)

By substituting the equation of motion and integrating by parts the averages,
it is possible to show that (see appendix B):

$\frac{d}{dt}\left\langle
V\right\rangle_{\hat{\rho}(t)}=-\frac{d}{dt}\left\langle
K\right\rangle_{\hat{\rho}(t)}$ (67)

###  IV.3 Entropy conservation

The equations of motion (54e) conserve also the entropy. This is a general
feature of any Hamiltonian dynamics, and the TD-SCHA makes no exception, even
if the Hamiltonian depends self-consistently from the density matrix. The
entropy defined on the many-body density matrix is:

$S[\hat{\rho}(t)]=-k_{b}\operatorname{Tr}\left[\hat{\rho}(t)\log\hat{\rho}(t)\right]$
(68)

for simplicity, we drop the explicit $t$ dependence of $\hat{\rho}(t)$.

$B(\hat{\rho})=\hat{\rho}\log\hat{\rho}$ (69)
$\frac{dS}{dt}=-k_{b}\operatorname{Tr}\left[\frac{d\hat{\rho}}{dt}\frac{dB(\hat{\rho})}{d\hat{\rho}}\right]$
(70)
$\frac{dS}{dt}=\frac{ik_{b}}{\hbar}\left[\operatorname{Tr}\left({\hat{\mathcal{H}}}[\rho]\hat{\rho}\frac{dB(\hat{\rho})}{d\hat{\rho}}\right)-\operatorname{Tr}\left(\hat{\rho}{\hat{\mathcal{H}}}[\rho]\frac{dB(\hat{\rho})}{d\hat{\rho}}\right)\right]$
(71)

By exploiting the cyclic permutation of the trace and the commutation between
$dB(\hat{\rho})/d\hat{\rho}$ and $\hat{\rho}$, Eq. (71) is zero.

The entropy conservation derives from the unitary time evolution: we are
evolving a closed quantum system and there is no de-coherence in the dynamics.
In other words, the dynamics is reversible, as if we change the initial sign
of $\mathrm{Im}\bm{A}$, $\bm{C}$, and $\bm{Q}$, the evolution proceeds
backward in time.

###  IV.4 Steady-state solutions

Thanks to Eq. (46), the steady state solution of the dynamical equations
occurs when the density matrix $\hat{\rho}$ commutes with the self-consistent
Hamiltonian ${\hat{\mathcal{H}}}[\rho]$. This means that there is a basis that
simultaneously diagonalizes both the Hamiltonian ${\hat{\mathcal{H}}}[\rho]$
and the density matrix. This condition can also be inferred from the equation
of motions Eq. (54e), imposing that the time-derivatives are zero. In this
case we have:

$C_{ab}=0\qquad\mathrm{Im}A_{ab}=0\qquad Q_{a}=0$ (72)

This is quite intuitive, as discussed in Sec. IV.1, these variables are
related to the instantaneous average momentum and chirp.

Substituting Eq. (72) inside Eq. (54e) we get the following conditions:

$\left\langle f_{a}\right\rangle_{\rho}=0$ (73a)
$\frac{1}{2\hbar}\left\langle\frac{\partial^{2}V}{\partial\tilde{R}_{a}\partial\tilde{R}_{b}}\right\rangle_{\rho}=\frac{\hbar}{2}\left(\frac{1}{4}\tilde{\Theta}_{ab}^{2}-\mathrm{Re}{\tilde{A}}^{2}_{ab}\right)$
(73b)

Eq. (73a) is a necessary condition for the equilibrium SCHA: the average of
the BO forces is the gradient of the SCHA free energy versus the centroids
positions $\bm{\mathcal{R}}$[20]. Moreover, from Eq. (54d) we have that
$\tilde{\bm{\Theta}}$ and $\mathrm{Re}\tilde{\bm{A}}$ commute in the static
solution, from which also $\tilde{\bm{\Upsilon}}$ and
$\mathrm{Re}\tilde{\bm{A}}$ commute (that is a direct consequence of the
commutation between $\hat{\rho}$ and ${\hat{\mathcal{H}}}[\rho]$). Since it is
better to express our quantity as a function of the total dispersion
$\tilde{\bm{\Upsilon}}$ (Eq. 50) we have:

$\frac{\tilde{\bm{\Upsilon}}^{2}}{4}+\tilde{\bm{\Upsilon}}\mathrm{Re}\tilde{\bm{A}}=\frac{1}{\hbar^{2}}\left\langle\frac{\partial^{2}V}{\partial\tilde{\bm{R}}\partial\tilde{\bm{R}}}\right\rangle_{\rho}.$
(74)

If we express Eq. (74) in the basis that diagonalizes both
$\tilde{\bm{\Upsilon}}$ and $\mathrm{Re}\tilde{\bm{A}}$, the same basis must
diagonalize also
$\braket{\frac{\partial^{2}V}{\partial\tilde{\bm{R}}\partial\tilde{\bm{R}}}}$.
Let us define $\bm{e_{\mu}}$ the eigenvector of
$\braket{\frac{\partial^{2}V}{\partial\tilde{\bm{R}}\partial\tilde{\bm{R}}}}$,
$\tilde{\bm{\Upsilon}}$ and $\mathrm{Re}\tilde{\bm{A}}$ and
$\omega_{\mu}^{2}$, $\tilde{\Upsilon}_{\mu}$ and $\mathrm{Re}\tilde{A}_{\mu}$
the corresponding eigenvalues:

$\sum_{b}\left\langle\frac{\partial^{2}V}{\partial\tilde{R}_{a}\partial\tilde{R}_{b}}\right\rangle_{\rho}e_{\mu}^{b}=\omega_{\mu}^{2}e_{\mu}^{a},$
(75)

Eq. (74) becomes:

$\frac{\tilde{\Upsilon}^{2}_{\mu}}{4}+\tilde{\Upsilon}_{\mu}\mathrm{Re}\tilde{A}_{\mu}=\frac{\omega_{\mu}^{2}}{\hbar^{2}}.$
(76)

Without loss of generality, if we change variable introducing a new parameter
$n_{\mu}$ and define:

$\tilde{\Upsilon}_{\mu}=\frac{2\omega_{\mu}}{\hbar(2n_{\mu}+1)},$ (77a) from
Eq. (76), we get:
$\mathrm{Re}\tilde{A}_{\mu}=\frac{2\omega_{\mu}n_{\mu}(n_{\mu}+1)}{\hbar(2n_{\mu}+1)}.$
(77b)

In these equations, $n_{\mu}$ is a free parameter for each frequency: the TD-
SCHA equations are stationary for any choice of $n_{\mu}$. In particular, we
can write $n_{\mu}$ as the Bose-Eninstein occupation number with a temperature
that depends on $\mu$:

$n_{\mu}=\frac{1}{e^{\beta_{\mu}\hbar\omega_{\mu}}-1}$ (78)

The stationary density matrix identified by Eq. (77b) and (72) is the product
of equilibrium densities matrices of harmonic oscillators of frequencies
$\omega_{\mu}$ and temperatures $\beta_{\mu}$.

If we fix all $\beta_{\mu}=\beta$, we recover the standard finite temperature
SCHA solution (Eq. 34, 35, and 36), with $\beta$ the Boltzman factor[20].

$\beta=\frac{1}{k_{b}T}$ (79)

However, in principle, also solutions with different temperatures in different
modes are a stationary point of the TD-SCHA equations. This is true also for
the real dynamics: a perfectly harmonic oscillator is stationary even if each
normal mode is populated with a different temperature (since there is no
coupling and heat transfer between modes). In general, any mixture of
eigenstates of the time-independent Hamiltonian is a stationary density
matrix, even if the population is not the equilibrium one.

###  IV.5 Equilibrium and thermodynamics

In the previous section, we showed that infinite steady states exist. The
steady-state solutions are equal to the product of equilibrium noninteracting
quantum Harmonic oscillators, where each normal mode has a thermal occupation
number $n_{\mu}$ with a temperature that depends on the mode. This is not,
indeed, the equilibrium solution of the SCHA, which requires all normal modes
populated by the same temperature. However, not all steady-states are
equilibrium solutions: if we prepare the system in a mixture of state with
eigenstates of the Hamiltonian, the exact time-dependent density matrix is
stationary even if the occupation probabilities are not the Boltzmann factors.

The system reaches equilibrium if we introduce an interaction with an external
bath (or a dephasing mechanism). In an isolated system, the equilibrium
solution is the one that maximizes the entropy among all possible steady
states at fixed energy. Maximizing the entropy fixing the energy is equivalent
in minimizing the Helmholtz free energy. This is the starting point for the
static SCHA. In Appendix D we prove that this condition correspond imposing an
uniform $\beta$ on each mode.

Therefore, we recover the SCHA as the static solution of the TD-SCHA that
maximizes the entropy. It is worth noting that the TD-SCHA will not
spontaneously evolve into the SCHA solution, unless the equations are modified
to account for the coupling with a reservoir that provides a mechanism for
quantum decoherence[46], allowing the entropy of the subsystem to increase.
The extension of TD-SCHA to describe the dynamics of an open quantum system is
beyond the scope of the current work.

##  V Linear response theory

Almost all experimental data are collected by probing the response of the
system to a time-dependent external perturbation. This perturbation could be
either electromagnetic radiation (static electric fields, IR, optical light,
X-ray) or particles like electrons and neutrons. If the perturbation does not
provide enough energy to heat the system, we are in the linear regime[47, 48].
Typical experiments that involve interactions with ionic degrees of freedom
are Raman and IR spectroscopy, neutron, and X-Ray scattering. In this section,
we present the dynamical linear response of the TD-SCHA equations on top of
the static SCHA solution. This enables the computation of the response
function for any experiment probing the nuclear motion fully _ab initio_ and
considering both quantum/thermal fluctuations and anharmonicity beyond
perturbation theory.

Since the probe heating of the sample (Joule effect) enter at higher orders in
the perturbation, the linear response does not depend on the coupling with the
thermal bath[47, 48]. Therefore, even if the TD-SCHA introduced in this work
describes closed quantum systems, the results we derive in this section are
general, and apply also to systems coupled with a bath.

We start from the thermodynamic equilibrium (the SCHA solution) and then we
add “small” time-dependent external potential $V^{(1)}(\bm{R},t)$ that acts on
the nuclei for $t\geq t_{0}$. We study how the parameters of our density
matrix changes under the effect of the external perturbation:

$\mathcal{R}_{i}(t)=\mathcal{R}_{i}+\mathcal{R}_{i}^{(1)}(t)\qquad
Q_{i}(t)=Q_{i}^{(1)}(t)$ (80a)
$\Upsilon_{ab}(t)=\Upsilon_{ab}+\Upsilon_{ab}^{(1)}(t)\qquad
C_{ab}(t)=C_{ab}^{(1)}(t)$ (80b) $A_{ab}(t)=A_{ab}+A_{ab}^{(1)}(t)$ (80c)
${V^{(\text{tot})}}(\bm{R},t)=V(\bm{R})+V^{(1)}(\bm{R},t)$ (80d)

where the (1) indicates the small perturbation around the equilibrium SCHA
solution. We restrict to cases in which $\hat{\mathcal{A}}$ and
$\hat{\mathcal{B}}$ depend only on nuclear positions. For this reason, the
only dynamics we are interested on is the diagonal elements of the density
matrix, i.e the ionic probability distribution $\rho(\bm{R},t)$. This depends
only on $\bm{\Upsilon}(t)$ and $\bm{\mathcal{R}}(t)$. Substituting Eq. (80d)
into the equations of motion (Eq. 54e) and keeping only linear terms in the
perturbation, we get a system of linear differential equations. Since we are
interested in $\bm{\Upsilon}(t)$ and $\bm{\mathcal{R}}(t)$, the minimal linear
system that describe their dynamics is the following:

$\frac{d^{2}{\tilde{\bm{\Upsilon}}}^{(1)}}{dt^{2}}={\bar{\bm{X}}}{\tilde{\bm{\Upsilon}}}^{(1)}+{\bar{\bm{Y}}}{\mathrm{Re}\tilde{\bm{A}}}^{(1)}+{\bar{\bm{Z}}}{\tilde{\bm{\mathcal{R}}}}^{(1)}+\bm{f}^{(1)}_{\Upsilon},$
(81a)
$\frac{d^{2}{\mathrm{Re}\tilde{\bm{A}}}^{(1)}}{dt^{2}}={\bar{\bm{X}}}^{\prime}{\tilde{\bm{\Upsilon}}}^{(1)}+{\bar{\bm{Y}}}^{\prime}{\mathrm{Re}\tilde{\bm{A}}}^{(1)}+{\bar{\bm{Z}}}^{\prime}{\tilde{\bm{\mathcal{R}}}}^{(1)}+\bm{f}^{(1)}_{\mathrm{Re}A},$
(81b)
$\frac{d^{2}{\tilde{\bm{\mathcal{R}}}}^{(1)}}{dt^{2}}={\bar{\bm{X}}}^{\prime\prime}{\tilde{\bm{\Upsilon}}}^{(1)}+{\bar{\bm{Z}}}^{\prime\prime}{\tilde{\bm{\mathcal{R}}}}^{(1)}+\bm{f}_{\mathcal{R}}^{(1)}.$
(81c)

Even if $\mathrm{Re}\tilde{\bm{A}}(t)$ does not affect the average of the
observable directly, we need to keep it, as $\bm{\Upsilon}(t)$ depends
explicitly on its dynamics. In Eq. (81c), the bar $\bar{\cdot}$ over a symbol
indicates a tensor. $\bar{\bm{X}}$, $\bar{\bm{Y}}$, $\bar{\bm{X}^{\prime}}$,
and $\bar{\bm{Y}^{\prime}}$ are 4-rank tensors, $\bar{\bm{Z}}$,
$\bar{\bm{Z}^{\prime}}$, and $\bar{\bm{X}^{\prime\prime}}$ are 3-rank tensors,
while $\bar{\bm{Z}^{\prime\prime}}$ is a 2-rank tensor. We recall that
$\tilde{\bm{\Upsilon}}^{(1)}$ and $\mathrm{Re}\tilde{\bm{A}}^{(1)}$ are 2-rank
tensor, as well as $\bm{f}_{\Upsilon}^{(1)}$ and
$\bm{f}_{\mathrm{Re}A}^{(1)}$, while $\bm{\mathcal{R}}^{(1)}$ and
$\bm{f}_{\mathcal{R}}^{(1)}$ are vectors (1-rank).

The product between tensors is the defined by the operator on the left. 4-rank
tensors $\bar{\bm{X}}$, $\bar{\bm{Y}}$, $\bar{\bm{X}^{\prime}}$, and
$\bar{\bm{Y}^{\prime}}$ are contracted on the last two indices:

$\left(\bar{\bm{X}}{\tilde{\bm{\Upsilon}}}^{(1)}\right)_{ab}=\sum_{cd}{\bar{X}}_{abcd}{\tilde{\Upsilon}}^{(1)}_{cd},$
(82a) also 3-rank tensors $\bar{\bm{Z}}$, $\bar{\bm{Z}^{\prime}}$ are
contracted on the last two indices:
$\left(\bar{\bm{Z}}{\tilde{\bm{\Upsilon}}}^{(1)}\right)_{a}=\sum_{bc}\bar{Z}_{abc}{\tilde{\Upsilon}}^{(1)}_{bc}.$
(82b) The 3-rank tensor $\bar{\bm{X}^{\prime\prime}}$ is contracted only on
the last index:
$\left(\bar{\bm{X}^{\prime\prime}}{\tilde{\bm{\mathcal{R}}}}^{(1)}\right)_{ab}=\sum_{bc}\bar{X}^{\prime\prime}_{abc}{\tilde{\mathcal{R}}}^{(1)}_{c}.$
(82c)

The 2-rank tensor $\bar{\bm{Z}^{\prime\prime}}$ is contracted on the last
index, as the standard matrix rows-by-columns product.

The details on the derivation of Eq. (81c) are reported in Appendix E, and the
explicit expression of the tensors in Appendix F.

It is convenient to express Eq. (81c) in matrical form, as:

$\frac{d^{2}}{dt^{2}}\begin{pmatrix}{\tilde{\bm{\Upsilon}}}^{(1)}\\\
\mathrm{Re}{\tilde{\bm{A}}}^{(1)}\\\
{\tilde{\bm{\mathcal{R}}}}^{(1)}\end{pmatrix}=\begin{pmatrix}\bar{\bm{X}}&\bar{\bm{Y}}&\bar{\bm{Z}}\\\
\bar{\bm{X}^{\prime}}&\bar{\bm{Y}^{\prime}}&\bar{\bm{Z}^{\prime}}\\\
\bar{\bm{X}^{\prime\prime}}&0&\bar{\bm{Z}^{\prime\prime}}\end{pmatrix}\begin{pmatrix}{\tilde{\bm{\Upsilon}}}^{(1)}\\\
\mathrm{Re}{\tilde{\bm{A}}}^{(1)}\\\
{\tilde{\bm{\mathcal{R}}}}^{(1)}\end{pmatrix}+\begin{pmatrix}\bm{f}_{\Upsilon}^{(1)}\\\
\bm{f}_{\mathrm{Re}A}^{(1)}\\\ \bm{f}_{\mathcal{R}}^{(1)}\end{pmatrix}$ (83)

The tensors introduced in Eq. (81c) account for the free time-evolution of the
system with the full anharmonic interaction: they are defined by the static
unperturbed Hamiltonian ${\hat{\mathcal{H}}}[\rho]$. In particular, they
depend on phonon scattering vertexes: the 4-phonon scattering tensor
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ and 3-phonon scattering
tensor $\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$, and the free
evolution $\overset{\scriptscriptstyle{\text{(2)}}}{\bm{D}}$

$\overset{\scriptscriptstyle{\text{(4)}}}{D}_{abcd}=\left\langle\frac{\partial^{4}V}{\partial\tilde{R}_{a}\partial\tilde{R}_{b}\partial\tilde{R}_{c}\partial\tilde{R}_{d}}\right\rangle_{{\rho}},$
(84)
$\overset{\scriptscriptstyle{\text{(3)}}}{D}_{abc}=\left\langle\frac{\partial^{3}V}{\partial\tilde{R}_{a}\partial\tilde{R}_{b}\partial\tilde{R}_{c}}\right\rangle_{{\rho}},$
(85)
$\overset{\scriptscriptstyle{\text{(2)}}}{D}_{ab}=\left\langle\frac{\partial^{2}V}{\partial\tilde{R}_{a}\partial\tilde{R}_{b}}\right\rangle_{{\rho}}$
(86)

where the averages are computed over the unperturbed ionic probability
distribution $\rho(\bm{R})$. Due to the Gaussian constrain on the density
matrix, the TD-SCHA does not account directly for higher-order phonon
scattering processes. However,
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$,
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$, and
$\overset{\scriptscriptstyle{\text{(2)}}}{\bm{D}}$ depend “self-consistently”
on higher-order anharmonicities, as they are averaged on the equilibrium
distribution. In fact, they are temperature dependent.

In particular, $\bar{\bm{X}}$ and $\bar{\bm{X}^{\prime}}$ depend on
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$, while $\bar{\bm{Z}}$,
$\bar{\bm{Z}^{\prime}}$ and $\bar{\bm{X}^{\prime\prime}}$ depend on
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$. The $\bar{\bm{X}}$,
$\bar{\bm{X}^{\prime}}$, $\bar{\bm{Y}}$, $\bar{\bm{Y}^{\prime}}$ and
$\bar{\bm{Z}^{\prime\prime}}$ contain terms of the free evolution, that are
non zero even if the system is a perfect harmonic crystal. Table 1 summarizes
these dependencies.

|  $\overset{\scriptscriptstyle{\text{(2)}}}{\bm{D}}$ |  $\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$ |  $\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$
---|---|---|---
$\bar{\bm{X}}$ |  $\bm{\circ}$ |  |  $\bm{\circ}$
$\bar{\bm{Y}}$ |  $\bm{\circ}$ |  |
$\bar{\bm{Z}}$ |  |  $\bm{\circ}$ |
$\bar{\bm{X}^{\prime}}$ |  $\bm{\circ}$ |  |  $\bm{\circ}$
$\bar{\bm{Y}^{\prime}}$ |  $\bm{\circ}$ |  |
$\bar{\bm{Z}^{\prime}}$ |  |  $\bm{\circ}$ |
$\bar{\bm{X}^{\prime\prime}}$ |  |  $\bm{\circ}$ |
$\bar{\bm{Z}^{\prime\prime}}$ |  $\bm{\circ}$ |  |
Table 1: Dependency of the coefficients of the linear response system (Eq.
81c) on the free evolution
($\overset{\scriptscriptstyle{\text{(2)}}}{\bm{D}}$), on anharmonic coupling
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$ (Eq. 85) and
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ (Eq. 84). A $\bm{\circ}$ in
the grid indicates that the tensor on the first column depends on the
corresponding term in the first row. This table helps to see which terms can
be set to zero if we neglect
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$ or
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$, and to visualize how
anharmonicity couples different degrees of freedom

The $\bm{f}_{\Upsilon}^{(1)}$, $\bm{f}_{\mathrm{Re}A}^{(1)}$, and
$\bm{f}_{\mathcal{R}}^{(1)}$ represent how the bare perturbation
$V^{(1)}(\bm{R},t)$ enters in the equation of the motion. In the polarization
basis (i.e. the basis that diagonalizes ${\tilde{\bm{\Upsilon}}}^{(0)}$ and
$\mathrm{Re}{\tilde{\bm{A}}}^{(0)}$ introduced in Sec. IV.4), they are:

${f_{\Upsilon}^{(1)}}_{\mu\nu}=\frac{1}{\hbar}\left(\frac{2\omega_{\mu}}{2n_{\mu}+1}+\frac{2\omega_{\nu}}{2n_{\nu}+1}\right)\left\langle\frac{\partial^{2}V^{(1)}}{\partial\tilde{R}_{\mu}\partial\tilde{R}_{\nu}}\right\rangle_{{\rho}},$
(87a) $\displaystyle{f_{\mathrm{Re}A}^{(1)}}_{\mu\nu}=\frac{1}{\hbar}$
$\displaystyle\bigg{(}\frac{2(n_{\mu}+1)n_{\mu}\omega_{\mu}}{2n_{\mu}+1}+$
$\displaystyle+\frac{2(n_{\nu}+1)n_{\nu}\omega_{\nu}}{2n_{\nu}+1}\bigg{)}\left\langle\frac{\partial^{2}V^{(1)}}{\partial\tilde{R}_{\mu}\partial\tilde{R}_{\nu}}\right\rangle_{{\rho}}$
(87b) ${f^{(1)}_{\mathcal{R}}}_{\mu}=-\left\langle\frac{\partial
V^{(1)}}{\partial\tilde{R}_{\mu}}\right\rangle_{{\rho}}$ (87c)

If we rewrite the perturbation in frequency domain, we can define the Green
function:

$\bm{G}(\omega)=-\begin{pmatrix}\bar{\bm{X}}+\omega^{2}&\bar{\bm{Y}}&\bar{\bm{Z}}\\\
\bar{\bm{X}^{\prime}}&\bar{\bm{Y}^{\prime}}+\omega^{2}&\bar{\bm{Z}^{\prime}}\\\
\bar{\bm{X}^{\prime\prime}}&0&\bar{\bm{Z}^{\prime\prime}}+\omega^{2}\end{pmatrix}^{-1}$
(88)

Given any possible time-dependent perturbation one nuclei
$V^{(1)}(\bm{R},\omega)$ we can compute the perturbation vector
$\bm{f}^{(1)}_{\Upsilon|\mathrm{Re}A|\mathcal{R}}$ and obtain the response as:

$\begin{pmatrix}{\tilde{\bm{\Upsilon}}}^{(1)}(\omega)\\\
\mathrm{Re}{\tilde{\bm{A}}}^{(1)}(\omega)\\\
{\tilde{\bm{\mathcal{R}}}}^{(1)}(\omega)\end{pmatrix}=\bm{G}(\omega)\begin{pmatrix}\bm{f}^{(1)}_{\Upsilon}(\omega)\\\
\bm{f}^{(1)}_{\mathrm{Re}A}(\omega)\\\
\bm{f}^{(1)}_{\mathcal{R}}(\omega)\end{pmatrix}.$ (89)

###  V.1 The general response function

In this section, we derive the general expression of the response function
$\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)$ within the TD-SCHA theory for any
couple of ionic time-independent observables $\hat{\mathcal{A}}$ and
$\hat{\mathcal{B}}$.

The $\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)$ describes how an external
perturbation, interacting with the ions through $\hat{\mathcal{B}}$, affects
the average of $\hat{\mathcal{A}}$, as described in Sec. II.

$V^{(1)}(\hat{\bm{R}},t)={\mathcal{B}}(\hat{\bm{R}})\mathcal{V}(t)$ (90)
$\left\langle{\mathcal{A}}(\bm{R})\right\rangle_{\rho^{(1)}(\omega)}=\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)\mathcal{V}(\omega)$
(91)

Since we are in linear regime, we expand the average of $\hat{\mathcal{A}}$ at
first order around the equilibrium solution:

$\displaystyle\left\langle{\mathcal{A}}(\bm{R})\right\rangle_{\rho^{(1)}(\omega)}$
$\displaystyle=\sum_{\mu\nu}\frac{\partial\left\langle{\mathcal{A}}(\bm{R})\right\rangle_{\rho(\omega)}}{\partial{\tilde{\Upsilon}}_{\mu\nu}}{\tilde{\Upsilon}}^{(1)}_{\mu\nu}+$
$\displaystyle+\sum_{\mu\nu}\frac{\partial\left\langle{\mathcal{A}}(\bm{R})\right\rangle_{\rho(\omega)}}{\partial\mathrm{Re}{\tilde{A}}_{\mu\nu}}\mathrm{Re}{\tilde{A}}^{(1)}_{\mu\nu}+$
$\displaystyle+\sum_{\mu}\frac{\partial\left\langle{\mathcal{A}}(\bm{R})\right\rangle_{\rho(\omega)}}{\partial{\tilde{\mathcal{R}}}_{\mu}}{\tilde{\mathcal{R}}}^{(1)}_{\mu}$
(92)

where the derivatives in Eq. (92) are evaluated with the equilibirum SCHA
density matrix. In particular, since the $\hat{\mathcal{A}}$ observable
depends only on the ionic positions $\hat{\bm{R}}$, only the
${\tilde{\bm{\Upsilon}}}^{(1)}$ and ${\tilde{\bm{\mathcal{R}}}}^{(1)}$ affects
its average.

By exploiting the integration by parts as illustrated in ref.[33], we get (in
the polarization basis of Eq. (76)):

$\frac{\partial\left\langle{\mathcal{A}}(\bm{R})\right\rangle_{\rho(\omega)}}{\partial{\tilde{\Upsilon}}_{\mu\nu}}=-\hbar^{2}\frac{(2n_{\mu}+1)(2n_{\nu}+1)}{8\omega_{\mu}\omega_{\nu}}\left\langle\frac{\partial^{2}{\mathcal{A}}}{\partial\tilde{R}_{\mu}\partial\tilde{R}_{\nu}}\right\rangle_{{\rho}}$
(93a)
$\frac{\partial\left\langle{\mathcal{A}}(\bm{R})\right\rangle_{\rho(\omega)}}{\partial\mathrm{Re}{\tilde{A}}_{\mu\nu}}=0$
(93b)
$\frac{\partial\left\langle{\mathcal{A}}(\bm{R})\right\rangle_{\rho(\omega)}}{\partial{\tilde{\mathcal{R}}}_{\mu}}=\left\langle\frac{\partial{\mathcal{A}}}{\partial\tilde{R}_{\mu}}\right\rangle_{{\rho}}$
(93c)

In the same way, we employ Eq. (87c) to express how the perturbation
$\hat{\mathcal{B}}$ affects the density matrix dynamics:

${f_{\Upsilon}^{(1)}}_{\mu\nu}(\omega)=\frac{1}{\hbar}\left(\frac{2\omega_{\mu}}{2n_{\mu}+1}+\frac{2\omega_{\nu}}{2n_{\nu}+1}\right)\left\langle\frac{\partial^{2}{\mathcal{B}}}{\partial\tilde{R}_{\mu}\partial\tilde{R}_{\nu}}\right\rangle_{{\rho}}\mathcal{V}(\omega),$
(94a) $\displaystyle{f_{\mathrm{Re}A}^{(1)}}_{\mu\nu}(\omega)=\frac{1}{\hbar}$
$\displaystyle\bigg{(}\frac{2(n_{\mu}+1)n_{\mu}\omega_{\mu}}{2n_{\mu}+1}+$
$\displaystyle+\frac{2(n_{\nu}+1)n_{\nu}\omega_{\nu}}{2n_{\nu}+1}\bigg{)}\left\langle\frac{\partial^{2}{\mathcal{B}}}{\partial\tilde{R}_{\mu}\partial\tilde{R}_{\nu}}\right\rangle_{{\rho}}\mathcal{V}(\omega)$
(94b)
${f^{(1)}_{\mathcal{R}}}_{\mu}(\omega)=-\left\langle\frac{\partial{\mathcal{B}}}{\partial\tilde{R}_{\mu}}\right\rangle_{{\rho}}\mathcal{V}(\omega)$
(94c)

Since $\bm{f}^{(1)}$ is proportional to $\mathcal{V}$, it is convenient to
define $\textbf{f}^{(1)}$ as:

$\bm{f}_{x}^{(1)}(\omega)=\textbf{f}_{x}^{(1)}\mathcal{V}(\omega)\qquad
x=\tilde{\Upsilon},\mathrm{Re}\tilde{A},\tilde{\mathcal{R}}$ (95)

We can express the response function with a standard linear-algebra matrix-
vector multiplications. We introduce the response vector $\bm{p}$ and the
perturbation vector $\bm{q}$ as:

$\bm{p}=\begin{pmatrix}\frac{\partial\left\langle\hat{\mathcal{A}}\right\rangle_{\hat{\rho}(\omega)}}{\partial{\tilde{\bm{\Upsilon}}}}&\frac{\partial\left\langle\hat{\mathcal{A}}\right\rangle_{\hat{\rho}(\omega)}}{\partial\mathrm{Re}{\tilde{\bm{A}}}}&\frac{\partial\left\langle\hat{\mathcal{A}}\right\rangle_{\hat{\rho}(\omega)}}{\partial{\tilde{\bm{\mathcal{R}}}}}\end{pmatrix},$
(96) $\bm{q}=\begin{pmatrix}{\textbf{f}_{\Upsilon}^{(1)}}\\\
{\textbf{f}_{\mathrm{Re}A}^{(1)}}\\\
{\textbf{f}^{(1)}_{\mathcal{R}}}\end{pmatrix}.$ (97)

In this way, the average of $\hat{\mathcal{A}}$ (Eq. 92) is a simple scalar
product:

$\left\langle\hat{\mathcal{A}}\right\rangle_{\rho^{(1)}(\omega)}=\bm{p}\cdot\begin{pmatrix}{\tilde{\bm{\Upsilon}}}^{(1)}(\omega)\\\
\mathrm{Re}{\tilde{\bm{A}}}^{(1)}(\omega)\\\
\tilde{\bm{\mathcal{R}}}^{(1)}(\omega)\end{pmatrix}$ (98)

The last vector in Eq. (98) is the result of the linear response system (Eq.
89), and it is the product between the Green function $\bm{G}(\omega)$ to the
$\bm{q}$ vector:

$\begin{pmatrix}{\tilde{\bm{\Upsilon}}}^{(1)}(\omega)\\\
\mathrm{Re}{\tilde{\bm{A}}}^{(1)}(\omega)\\\
\tilde{\bm{\mathcal{R}}}^{(1)}(\omega)\end{pmatrix}=\bm{G}(\omega)\bm{q}\mathcal{V}(\omega).$
(99)

Combining Eq. (98) with Eq. (99), we get the expression of the response
function in the TD-SCHA linear response formalism:

$\left\langle\hat{\mathcal{A}}\right\rangle_{\rho^{(1)}(\omega)}=\bm{p}\bm{G}(\omega)\bm{q}\mathcal{V}(\omega)$
(100)

Comparing Eq. (100) with the definition of the response function
$\chi_{{\mathcal{A}}{\mathcal{B}}}$ (Eq. 91) we get:

$\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)=\bm{p}\bm{G}(\omega)\bm{q},$ (101)

where $\bm{G}(\omega)$ is defined in Eq. (88), $\bm{p}$ in Eq. (96) and
$\bm{q}$ in Eq. (97).

###  V.2 Interacting one-phonon Green function

As a simple illustration, here, we compute the expression for the interacting
one-phonon Green function. This is the response function when
$\hat{\mathcal{A}}=\sqrt{m}_{a}\hat{R}_{a}$ and
$\hat{\mathcal{B}}=\sqrt{m}_{b}\hat{R}_{b}$.

$\left\langle\sqrt{m}_{a}\hat{R}_{a}\right\rangle_{\rho^{(1)}(\omega)}=\tilde{\mathcal{R}}^{(1)}_{a}(\omega)$
(102)

Thus, the response $\bm{p}$ vector (Eq. 96) is:

$\bm{p}=\begin{pmatrix}0&0&\bm{\delta}_{a}\end{pmatrix}$ (103)

where we indicate with $\bm{\delta}_{a}$ the vector with 1 in the $a$-th
atom/Cartesian index and zeros elsewhere.

On the other side, the perturbation vector $\bm{q}$ can be obtained
substituting $\hat{\mathcal{B}}$ into Eq. (94c):

$\textbf{f}_{\Upsilon}^{(1)}=0\qquad\textbf{f}_{\mathrm{Re}A}^{(1)}=0$ (104)
${\text{f}_{\mathcal{R}}^{(1)}}_{c}(\omega)=-\delta_{cb}.$ (105)

Then, from Eq. (97), we get:

$\bm{q}=-\begin{pmatrix}0\\\ 0\\\ \bm{\delta}_{b}\end{pmatrix}$ (106)

From which we get the interacting one-phonon Green function, as the last block
of the Green function:

$\chi_{\sqrt{m}_{a}R_{a},\sqrt{m}_{b}R_{b}}(\omega)=\mathcal{G}_{ab}(\omega)=-\begin{pmatrix}0&0&\bm{\delta}_{a}\end{pmatrix}{\bm{G}}(\omega)\begin{pmatrix}0\\\
0\\\ \bm{\delta}_{b}\end{pmatrix}$ (107)

The one-phonon Green function describes lattice excitations. Its poles are the
energy of the physical phonons in the system, and its imaginary part is
proportional to the intensity of the fully anharmonic phonon branches probed
in X-ray or neutron scattering experiments.

###  V.3 Physical phonons when cubic anharmonicity is negligible

The full expressions for the Green function is quite complex. However, it
becomes easy if we consider the special case for:

$\overset{\scriptscriptstyle{\text{(3)}}}{D}_{abc}=\left\langle\frac{\partial^{3}V^{(0)}}{\partial\tilde{R}_{a}\partial\tilde{R}_{b}\partial\tilde{R}_{c}}\right\rangle_{{\rho}}=0.$
(108)

This condition is satisfied in some points of the Brillouin zone in highly
symmetric crystals[49], or perfectly harmonic systems.

If this is the case, the $\bm{\mathcal{R}}^{(1)}$ block of the total Green
function $\bm{G}(\omega)$ is decoupled from the others as
$\bar{\bm{X}^{\prime\prime}}$, $\bar{\bm{Z}}$ and $\bar{\bm{Z}^{\prime}}$ are
all equal to zero (see Table 1).

Thanks to this block diagonal form of $\bm{G}(\omega)$, the inversion needed
to compute the one-phonon Green function requires only inverting the
$\bm{\mathcal{R}}^{(1)}$ block. In particular:

${\bm{\mathcal{G}}}(\omega)=(\omega^{2}+\bar{\bm{Z}^{\prime\prime}})^{-1}$
(109)

If we write this in the basis of polarization modes substituting the
expression of $\bar{\bm{Z}^{\prime\prime}}$ (see Appendix F), we get:

${\mathcal{G}}_{\mu\nu}(\omega)=\frac{\delta_{\mu\nu}}{\omega^{2}-\omega_{\mu}^{2}}.$
(110)

This expression resembles the Green function of harmonic phonons. However,
this is not true, as the $\omega_{\mu}$ frequencies do not coincide with the
harmonic result. In particular, $\omega_{\mu}$ frequencies are the SCHA
frequencies, obtained as the average of the second derivative of the BO
potential in the equilibrium distribution (Eq. 75). Thanks to this average,
the $\omega_{\mu}$ frequencies are affected by the
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ and, thus, they are
temperature dependent. Therefore, the static SCHA frequencies $\omega_{\mu}$
coincides with the poles of the dynamical Green function only if
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}=0$. However, this is not
true in the most general case, and to find the poles of the Green function, we
must invert Eq. (88).

Notably, this is a general misunderstanding in empirical mean-field methods to
deal with strong ionic anharmonicity, as TDEP[50], ALAMODE[51] and the static
SCHA[20]. The self-consistent frequencies $\omega_{\mu}$ extracted from these
methods are not the phonon frequencies probed by dynamical experiments. The
reason is that the $\omega_{\mu}$s extracted from this methods do not coincide
with the poles of the interacting one-phonon Green function, unless the cubic
anharmonicity is negligiable.

###  V.4 TD-SCHA self-energy

In Sec. V.3, we derived the expression of the one-phonon Green function if
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}=0$. Following the standard
many-body approach to Green functions, we define the TD-SCHA self-energy
$\bm{\Sigma}(\omega)$ as:

$\bm{{\mathcal{G}}}^{-1}(\omega)=\left(\bar{\bm{Z}^{\prime\prime}}+\omega^{2}\right)-\bm{\Sigma}(\omega),$
(111)

where $({\bm{Z}^{\prime\prime}}-\omega^{2})$ is the inverse of the one-phonon
Green function with $\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}=0$ (Eq.
109).

Due to the importance in calculating the interacting one-phonon Green
function, the authors of ref.[33] formulated a dynamical _ansatz_ for
$\bm{\Sigma}(\omega)$, proposing an analytical continuation at finite
frequencies of the static self-energy calculated from equilibrium SCHA.
Ref.[33] proved the correctness of this _ansatz_ in the limit of small
anharmonicity (perturbation theory) and for the static case $\omega=0$.

In this section, we derive the true expression of the self-energy
$\bm{\Sigma}(\omega)$ within the TD-SCHA, which includes the strong anharmonic
regime and high-frequency (as far as BO approximation holds).

The complete one-phonon green function is given by Eq. (107).

${\mathcal{G}}_{\mu\nu}(\omega)=\begin{pmatrix}0&0&\bm{\delta_{\mu}}\end{pmatrix}{\begin{pmatrix}\bar{\bm{X}}+\omega^{2}&\bar{\bm{Y}}&\bar{\bm{Z}}\\\
\bar{\bm{X}^{\prime}}&\bar{\bm{Y}^{\prime}}+\omega^{2}&\bar{\bm{Z}^{\prime}}\\\
\bar{\bm{X}^{\prime\prime}}&0&\bar{\bm{Z}^{\prime\prime}}+\omega^{2}\end{pmatrix}^{-1}}\begin{pmatrix}0\\\
0\\\ \bm{\delta_{\nu}}\end{pmatrix}.$ (112)

When $\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}\neq 0$,
$\bar{\bm{X}^{\prime\prime}}$ is non zero and the one-phonon Green function
interacts with the rest of the big $\bm{G}(\omega)$ matrix through
$\bar{\bm{X}^{\prime\prime}}$.

Thus, the self-energy is

$\bm{\Sigma}(\omega)={\bar{\bm{X}}}^{\prime\prime}\frac{\partial{\tilde{\bm{\Upsilon}}}^{(1)}}{\partial{\tilde{\bm{\mathcal{R}}}}^{(1)}}(\omega)$
(113)

and $\bar{\bm{X}^{\prime\prime}}$ is contracted with the indices of
${\tilde{\bm{\Upsilon}}}^{(1)}$. In particular, this term indicates how a
perturbation in the quantum fluctuations ${\tilde{\bm{\Upsilon}}}^{(1)}$
affects the average positions ${\tilde{\bm{\mathcal{R}}}}^{(1)}$. If
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}\neq 0$, the dynamics in
${\tilde{\bm{\mathcal{R}}}}^{(1)}$ affects ${\tilde{\bm{\Upsilon}}}^{(1)}$ and
$\mathrm{Re}{\tilde{\bm{A}}}^{(1)}$ through $\bar{\bm{Z}}$ and
$\bar{\bm{Z}^{\prime}}$ (Eq. 81a and 81b). They evolve freely (harmonic
propagation) and interacting anharmonically through
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$. Finally,
${\tilde{\bm{\Upsilon}}}^{(1)}$ affects back $\tilde{\bm{\mathcal{R}}}^{(1)}$
through $\bar{\bm{X}^{\prime\prime}}$ (Eq. (81c)).

To calculate $\bm{\Sigma}(\omega)$, we need to get an explicit expression from
the dependency of ${\tilde{\bm{\Upsilon}}}^{(1)}$ and
${\tilde{\bm{\mathcal{R}}}}^{(1)}$. It is convenient to remove the
$\mathrm{Re}{\tilde{\bm{A}}}^{(1)}(\omega)$ firstly (Eq. 81b):

$\bar{\bm{X}^{\prime}}{\tilde{\bm{\Upsilon}}}^{(1)}+(\bar{\bm{Y}^{\prime}}+\omega^{2})\mathrm{Re}{\tilde{\bm{A}}}^{(1)}+\bar{\bm{Z}^{\prime}}{\tilde{\bm{\mathcal{R}}}}^{(1)}=0$
(114)

$\mathrm{Re}{\tilde{\bm{A}}}^{(1)}=-(\bar{\bm{Y}^{\prime}}+\omega^{2})^{-1}\left(\bar{\bm{X}^{\prime}}{\tilde{\bm{\Upsilon}}}^{(1)}+\bar{\bm{Z}^{\prime}}{\tilde{\bm{\mathcal{R}}}}^{(1)}\right).$
(115)

We substitute it in the ${\tilde{\bm{\Upsilon}}}^{(1)}(\omega)$ equation (Eq.
81a):

$\displaystyle\left[(\bar{\bm{X}}+\omega^{2})-\bar{\bm{Y}}(\bar{\bm{Y}^{\prime}}+\omega^{2})^{-1}\bar{\bm{X}^{\prime}}\right]{\tilde{\bm{\Upsilon}}}^{(1)}+$
$\displaystyle+\left[\bar{\bm{Z}}-\bar{\bm{Y}}(\bar{\bm{Y}^{\prime}}+\omega^{2})^{-1}\bar{\bm{Z}^{\prime}}\right]{\tilde{\bm{\mathcal{R}}}}^{(1)}=0$
(116)

From which we get the relationship between
${\tilde{\bm{\Upsilon}}}^{(1)}(\omega)$ and
${\tilde{\bm{\mathcal{R}}}}^{(1)}(\omega)$. Here, products between tensors
follows the same convention as Eq. (82c), where the number of indices to be
contracted is determined by the tensor on the left. The inversion of a 4-rank
tensor is equivalent of inverting a matrix where we group the frist two and
last two indices:

$\bar{X}_{abcd}=\bar{X}_{(ab)(cd)}$ (117)

Now we have all the ingridients to write the interacting one-phonon self-
energy:

$\displaystyle\bm{\Sigma}(\omega)=-\bar{\bm{X}^{\prime\prime}}$
$\displaystyle\left[(\bar{\bm{X}}+\omega^{2})-\bar{\bm{Y}}(\bar{\bm{Y}^{\prime}}+\omega^{2})^{-1}\bar{\bm{X}^{\prime}}\right]^{-1}\cdot$
$\displaystyle\cdot\left[\bar{\bm{Z}}-\bar{\bm{Y}}(\bar{\bm{Y}^{\prime}}+\omega^{2})^{-1}\bar{\bm{Z}^{\prime}}\right]$
(118)

Interestingly, $\bar{\bm{X}^{\prime\prime}}$, $\bar{\bm{Z}}$ and
$\bar{\bm{Z}^{\prime}}$ depend on
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$. In Appendix G we prove
that this self-energy is equal to the one proposed as dynamical ansatz in ref.
[33]:

$\bm{\Sigma}(\omega)=\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}\bm{\Lambda}(\omega)\left[\mathbbm{1}-\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}\bm{\Lambda}(\omega)\right]^{-1}\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$
(119)

with $\mathbbm{1}$ the identity 4-rank tensor and

$\bm{\Lambda}(\omega)=-\frac{1}{2}\bm{\chi}(\omega)$ (120)

With $\bm{\chi}(\omega)$ the non interacting two-phonon Green function. We
calculate it in Eq. (135).

Contemporary to this work, another independent study[52] proved the dynamical
_ansatz_ of the SCHA, underlying the relevance of the topic.

###  V.5 IR and Raman response

In this section, we derive the explicit expression of the IR and Raman
response. To compute the general time-correlation function of the two
observables $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$ in the TD-SCHA
formalism, we follow the procedure introduced in Sec. V.1: we get the $\bm{p}$
and $\bm{q}$ vectors of Eq. (96) and (97) from the observables
$\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$, and compute the response
$\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)$ with Eq. (101).

Let us analyze the IR response first. This is related to the ionic dipole-
dipole correlation function:

$\hat{\mathcal{A}}=M_{\alpha}(\hat{\bm{R}})\qquad\hat{\mathcal{B}}=M_{\alpha^{\prime}}(\hat{\bm{R}})$
(121)

where $M_{\alpha}(\hat{\bm{R}})$ is the $\alpha$ Cartesian component of the
net dipole moment when the ions are displaced in the $\bm{R}$ position.

To compute the IR response we replace ${\mathcal{A}}=M_{\alpha}$ and
${\mathcal{B}}=M_{\alpha^{\prime}}$ inside the expressions of vectors $\bm{p}$
(Eq. 96) and $\bm{q}$ (Eq. 97).

So the quantities we need to compute to obtain $\bm{p}$ and $\bm{q}$ are the
averages of the dipole derivatives:

$\left\langle\frac{\partial M_{\alpha}(\bm{R})}{\partial
R_{a}}\right\rangle_{{\rho}}=\left\langle\mathcal{Z}_{\alpha
a}(\bm{R})\right\rangle_{{\rho}}$ (122)

where $\mathcal{Z}(\bm{R})$ is the effective charge of the system when ions
are displaced along $\bm{R}$, and, by exploiting integration by parts as
introduced in ref.[33],

$\left\langle\frac{\partial^{2}M_{\alpha}(\bm{R})}{\partial R_{a}\partial
R_{b}}\right\rangle_{{\rho}}=\sum_{c}\Upsilon_{ac}\left\langle(R_{c}-\mathcal{R}_{c})\mathcal{Z}_{\alpha
b}(\bm{R})\right\rangle_{{\rho}}.$ (123)

Both Eq. (122) and Eq. (123) can be computed by averaging the effective
charges in a random ensemble extracted according to the equilibrium
distribution $\rho(\bm{R})$. Since most standard _ab initio_ codes calculate
effective charges, the full IR response is accessible fully _ab initio_.

In general, the IR response involves the full interacting Green function. A
particular case is when effective charges do not depend on the ionic position.
In this case, the calculation of IR response is much easier, as Eq. (122)
becomes:

$\left\langle\mathcal{Z}_{\alpha
a}(\bm{R})\right\rangle_{{\rho}}=\mathcal{Z}_{\alpha a}(\bm{\mathcal{R}}),$
(124)

while Eq. (123):

$\sum_{c}\Upsilon_{ac}\left\langle(R_{c}-\mathcal{R}_{c})\mathcal{Z}_{\alpha
b}(\bm{R})\right\rangle_{{\rho}}=0$ (125)

In this simple case, only the $\bm{\mathcal{R}}$ block of $\bm{p}$ and
$\bm{q}$ is different from zero, thus we can link the IR response to the
interacting one-phonon Green function:

$\chi_{M_{\alpha}M_{\alpha^{\prime}}}(\omega)=\sum_{ab}\frac{\mathcal{Z}_{\alpha
a}\mathcal{Z}_{\alpha^{\prime}b}}{\sqrt{m_{a}m_{b}}}{\mathcal{G}}_{ab}(\omega)$
(126)

Notably, this equation, that is the standard approximation for IR signal, is
only valid for effective charges that do not depend on the ionic displacement.
TD-SCHA allows to compute the IR response even in the general scenario.

The same procedure holds also for the Raman response. Here, the
${\mathcal{A}}$ and ${\mathcal{B}}$ observables are the polarizability
$\alpha_{\alpha\beta}(\hat{\bm{R}})$ of the sample along $\alpha$ and $\beta$
direction (the incoming and outcoming polarization of the light) when ions are
displaced in the $\bm{R}$ position.

${\mathcal{A}}=\alpha_{\alpha\beta}(\hat{\bm{R}})\qquad{\mathcal{B}}=\alpha_{\alpha^{\prime}\beta^{\prime}}(\hat{\bm{R}})$
(127)

The procedure of deriving $\bm{p}$ and $\bm{q}$ is the same as for the IR. If
we define the Raman tensor $\bm{{\Xi}}$ as:

${\Xi}_{\alpha\beta
a}(\bm{R})=\frac{\partial\alpha_{\alpha\beta}(\bm{R})}{\partial R_{a}}$ (128)

Then we have:

$\left\langle\frac{\partial\alpha_{\alpha\beta}(\bm{R})}{\partial
R_{a}}\right\rangle_{{\rho}}=\left\langle{\Xi}_{\alpha\beta
a}(\bm{R})\right\rangle_{{\rho}}$ (129)

$\left\langle\frac{\partial^{2}\alpha_{\alpha\beta}(\bm{R})}{\partial
R_{a}\partial
R_{b}}\right\rangle_{{\rho}}=\sum_{c}\Upsilon_{ac}\left\langle(R_{c}-\mathcal{R}_{c}){\Xi}_{\alpha\beta
b}(\bm{R})\right\rangle_{{\rho}}.$ (130)

Also in this case, the full response function can be obtained by simply
calculating the raman tensor $\bm{{\Xi}}$ in a randomly distributed ensemble
of ionic configurations according to $\rho(\bm{R})$. The Raman tensor can be
calculated _ab initio_ efficiently with the method introduced by ref.[34]. If
the Raman tensor depends on nuclear positions, Eq. (130) is different from
zero, and we have a contribution to the response function from the complete
interacting Green function. However, if we neglect the dependence of
$\bm{{\Xi}}$ from the nuclear position, we can express the Raman response only
from the one-phonon Green function, as we did for the IR:

$\chi_{\alpha_{\alpha\beta}\alpha_{\alpha^{\prime}\beta^{\prime}}}(\omega)=\sum_{ab}\frac{{\Xi}_{\alpha\beta
a}{\Xi}_{\alpha^{\prime}\beta^{\prime}b}}{\sqrt{m_{a}m_{b}}}{\mathcal{G}}_{ab}(\omega).$
(131)

In next section we give a schematic overview of the processes neglected by
considering the response only due to the interacting one-phonon Green
function.

###  V.6 Mixed one-two and two-phonons Green function

The one-phonon Green function does not provide the response to any general
external perturbation $V^{(1)}(\bm{R},t)$, but only to those that depend
linearly on the ionic displacements $\bm{R}$, as we show for the specific
cases of IR and Raman (Sec. V.5). If the perturbation or the observable we
probe are nonlinear in $\bm{R}$, then the vectors $\bm{p}$ and $\bm{q}$ that
determine the response function $\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)$
(Eq. 101) have a non zero contribution also in the
${\tilde{\bm{\Upsilon}}}^{(1)}$ and $\mathrm{Re}{\tilde{\bm{A}}}^{(1)}$ sector
(Eq. 96 and 97).

This contribution is important if ionic fluctuations are sizable (e.g. in
presence of light atoms or close to a second-order phase-transition), or if
the linear term in $\bm{R}-\bm{\mathcal{R}}$ of ${\mathcal{A}}(\hat{\bm{R}})$
and ${\mathcal{B}}(\hat{\bm{R}})$ is zero by symmetry. In these conditions,
quadratic terms in the coupling between the probe and ionic displacement may
become important. It is fundamental to explain IR/Raman spectra of ice VII and
X[53] and liquid water[54]. In these works, the authors employed molecular
dynamics to go beyond the one-phonon Green function (neglecting quantum
fluctuations). Another class of materials with Raman spectra arising from a
nonlinear coupling between the probe and the ionic displacement are high-
symmetry structures like diamond[55]. Here, the authors calculated the Raman
signal within the harmonic approximation. A full treatment with both
anharmonicity and quantum fluctuations with a nonlinear probe is missing in
the literature.

In Figure 1 we report a diagrammatic representation of the response function.
Let time flow from left to right; when the probe interacts with the sample, it
can either excite a single phonon, with the linear term of ${\mathcal{B}}$
(panels c,d) or two phonons with the quadratic dependency of ${\mathcal{B}}$
on ionic positions (panels a,b). Then, the excited phonons evolve interacting
through anharmonicity. At the end of the anharmonic evolution, we may either
have a single phonon (panels b,d) or two phonons (panels a,c). Then, the
excited phonons give a response in the ${\mathcal{A}}$ observable; the average
of ${\mathcal{A}}(\hat{\bm{R}})$ is affected by two phonons proportionally to
its quadratic dependence on ionic displacement (panels a,c), while a single
phonon affects the ${\mathcal{A}}$ average proportionally to linear term of
${\mathcal{A}}(\hat{\bm{R}})$ (panels b,d). The overall process is divided in
(interacting) two-phonon Green functions (a), mixed one-two phonon Green
functions (b,c) and one-phonon Green functions (d). Harmonic systems do not
have a contribution from mixed Green functions, as phonons cannot decay or
scatter during the free propagation (panel e,f).

Figure 1: Schematic representation of all processes captured by the TD-SCHA
response function. The circles represent the interaction between the probe and
the ionic displacements, that occurs trhough observable ${\mathcal{A}}$. We
have two possible interaction, the linear one give rise to a single phonon
(dark red), the quadratic one excites two-phonons (light cyan). Panels a,b,c,d
represent the TD-SCHA processes, where the excited phonons propagate with the
full anharmonic interaction. We have a two-phonon signal (panel a), mixed one-
two phonons (panels b,c) and one-phonon (panel d). Panels e,f show the
possible processes with harmonic propagation.

The expression of the IR (Eq. 126) and Raman (Eq. 131) in the approximation of
linear coupling with the probe correspond to take only the contribution of the
process in panel d, neglecting all the others.

Interestingly, also in the pure harmonic case there is an additional
contribution due to non-linear couplings between the probe and the system,
that is reported in panel e. For example, the harmonic two-phonon IR signal
is:

$\chi^{(2ph)}_{M_{\alpha}M_{\alpha^{\prime}}}(\omega)=\sum_{abcd\mu\nu}\frac{e_{\mu}^{a}e_{\nu}^{b}e_{\mu}^{c}e_{\nu}^{d}}{\sqrt{m_{a}m_{b}m_{c}m_{d}}}\frac{d\mathcal{Z}_{\alpha
a}}{dR_{b}}\chi_{\mu\nu\mu\nu}(\omega)\frac{d\mathcal{Z}_{\alpha^{\prime}c}}{dR_{d}}$
(132)

where $\chi_{\mu\nu\mu\nu}(\omega)$ is the harmonic two-phonon Green function.
We can compute the expression of the $\bm{\chi}(\omega)$ tensor as the time
correlation between the nuclear quadratic displacements:

$\hat{\mathcal{A}}=\sqrt{m_{a}m_{b}}(\hat{R}_{a}-\mathcal{R}_{a})(\hat{R}_{b}-\mathcal{R}_{b})$
(133)
$\hat{\mathcal{B}}=\sqrt{m_{a}m_{b}}(\hat{R}_{c}-\mathcal{R}_{c})(\hat{R}_{d}-\mathcal{R}_{d})$
(134)

Since, in second-quantization formalism, each position operator $\hat{\bm{R}}$
is proportional to the creation-annihilation of one phonon,
$\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$ contain creations and
annihilations of two phonons; this is the reason why it is called the two-
phonon Green function.

We compute the $\bm{p}$ and $\bm{q}$ vectors (Eq. 96 and 97) from
$\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$ and insert them into Eq. (101) to
get the interacting two-phonon Green functions.

The non-interacting (harmonic) two-phonon Green function is obtained when both
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$ and
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ equal 0. In appendix J we
report the details of the calculation. The final result is:

$\displaystyle\chi_{\mu\nu\mu\nu}(\omega)=\frac{1}{2\omega_{\mu}\omega_{\nu}}$
$\displaystyle\bigg{[}\frac{(\omega_{\mu}+\omega_{\nu})(n_{\mu}+n_{\nu}+1)}{(\omega_{\mu}+\omega_{\nu})^{2}-\omega^{2}}+$
$\displaystyle-\frac{(\omega_{\mu}-\omega_{\nu})(n_{\mu}-n_{\nu})}{(\omega_{\mu}-\omega_{\nu})^{2}-\omega^{2}}\bigg{]}.$
(135)

Indeed, Eq. (135) is the well known two-phonon propagator for harmonic
systems, and it coincides with the one obtained with the standard many-body
approach.

Notably, since usually the non-linear interaction between the probe and the
ionic position is small, the diagrams b,c in Figure 1 provide a lower order
signal than the harmonic (e) diagram.

Indeed, all the processes represented in panels a,b,c,d are automatically
included in the response function of the TD-SCHA introduced in Sec. V.1. In
principle, also higher-order processes exist, where the probe interacts with
more than two phonons. However, the Gaussian constrain on the TD-SCHA density
matrix does not allow these excitations, but they are accounted for with a
mean-field approach affecting the average of the derivatives of
${\mathcal{A}}$ and ${\mathcal{B}}$ observables on the ionic probability
distribution, providing a temperature dependency for the static SCHA
frequencies and vertices in the diagrams reported in Figure 1.

##  VI Lanczos algorithm to compute the response function

While a numerical implementation to calculate the dynamical one-phonon Green
function has been presented[33], its computational cost diverges quickly for
systems with more than 10 atoms, as the self-energy of Eq. (119) requires, for
any frequency to probe, the inversion of the matrix

$\left[\mathbbm{1}-\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}\bm{\Lambda}(\omega)\right]^{-1}.$
(136)

This matrix has a dimension of $(3N)^{2}\times(3N)^{2}$; its numerical
inversion is a heavy computational task. For example, just to store in memory
the matrix of Eq. (136) in a system with 100 atoms with 64-bit floating-point
precision, more than 60 Gb are needed. The inversion of Eq. (136) requires a
LU decomposition of a matrix of size $n=9N^{2}$. The LU decomposition scales a
$n^{3}$; the overall scaling for computing the inversion is $N^{6}$, which
quickly diverges for realistic systems with more than 10 atoms. For this
reason, the application of the full dynamical one-phonon Green function has
been performed in realistic systems always under the assumption that
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$.

Moreover, as discussed in Sec. V.4, the one-phonon Green function does not
provide the most general response to the experimental probe, but it is limited
to probes interacting linearly with atomic displacements.

In this section, we derive a very efficient algorithm that allows computing
the elements of the response function to any general external perturbation,
that is computationally achievable in systems of hundreds of atoms even with
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}\neq 0$.

The response function to any generic probe can be obtained from the total
Green function (Eq. 88).

By defining the matrix $\mathcal{L}$ as:

$\mathcal{L}=-\begin{pmatrix}\bar{\bm{X}}&\bar{\bm{Y}}&\bar{\bm{Z}}\\\
\bar{\bm{X}^{\prime}}&\bar{\bm{Y}^{\prime}}&\bar{\bm{Z}^{\prime}}\\\
\bar{\bm{X}^{\prime\prime}}&0&\bar{\bm{Z}^{\prime\prime}}\end{pmatrix},$ (137)

the full Green function matrix is obtained as:

$\bm{G}(\omega)=(\mathcal{L}-\omega^{2})^{-1}.$ (138)

However, we are usually interested only in the response function to two
specific observables $\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$, and not the
full Green function in Eq. (138). In this case, instead of fully inverting the
big $\mathcal{L}-\omega^{2}$ matrix, we need just to compute the matrix
element:

$\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)=\bm{p}(\mathcal{L}-\omega^{2})^{-1}\bm{q}$
(139)

where $\bm{p}$ and $\bm{q}$ are defined in Eq. (96) and Eq. (97) from
$\hat{\mathcal{A}}$ and $\hat{\mathcal{B}}$.

The algorithm we discuss in this section computes only the matrix element of
Eq. (139) for any value of $\omega$. This is an easy task if we use a basis
where $\mathcal{L}$ is tridiagonal (as we show later). Therefore, the purpose
of the algorithm is to find the basis in which $\mathcal{L}$ is tridiagonal.

The bi-conjugate Lanczos algorithm does exactly this job: it builds with an
iterative procedure the basis-changing $Q$ matrix (non-singular) so that:

$Q^{-1}\mathcal{L}Q={\mathcal{T}},$ (140)

where $\mathcal{T}$ is tridiagonal:

${\mathcal{T}}=\begin{pmatrix}\alpha_{1}&\gamma_{1}&&\cdots&0\\\
\beta_{1}&\alpha_{2}&\ddots&&\vdots\\\ &\ddots&\ddots&\ddots&\\\
\vdots&&\ddots&\ddots&\gamma_{n-1}\\\
0&\cdots&&\beta_{n-1}&\alpha_{n}\end{pmatrix}$ (141)

The basis (non orthonormal) in which $\mathcal{L}$ is tridiagonal is
represented by the columns of the $Q$ matrix:

$Q=\begin{pmatrix}{\bm{q}_{1}}&\cdots&{\bm{q}_{n}}\end{pmatrix}$ (142)

The rows of the $Q^{-1}$ matrix defines the conjugate vectors:

$Q^{-1}=\begin{pmatrix}{\bm{p}_{1}}\\\ \vdots\\\ {\bm{p}_{n}}\end{pmatrix}.$
(143)

The coefficients of the $\mathcal{T}$ matrix (the $\alpha_{k}$, $\beta_{k}$,
$\gamma_{k}$), the $\bm{q}$ and $\bm{p}$ vectors are found with the iterative
bi-conjugate Lanczos algorithm[56]:

$\alpha_{k}=\bm{p_{k}}\cdot{\mathcal{L}}\bm{q_{k}}$ (144a)
$\beta_{k}\bm{q_{k+1}}=\bm{r_{k}}=(\mathcal{L}-\alpha_{k})\bm{q_{k}}-\gamma_{k-1}\bm{q_{k-1}}$
(144b)
$\gamma_{k}\bm{p_{k+1}}=\bm{s}_{k}=({\mathcal{L}}^{\dagger}-\alpha_{k})\bm{p_{k}}-\beta_{k-1}\bm{p_{k}}$
(144c) $\beta_{k}=|\bm{r_{k}}|$ (144d)
$\gamma_{k}=\frac{\bm{s_{k}}\cdot\bm{r_{k}}}{\beta_{k}}$ (144e)

This recursion formally ends either when either $\bm{q}_{k}$ or $\bm{p}_{k}$
are linear combinations of the previous vectors or if
$\bm{p}_{k}\cdot\bm{q}_{k}=0$. Unless the system is perfectly harmonic, this
condition is usually never reached in practical runs, and the algorithm is
truncated after a maximum number of steps $N_{\text{max}}$.

To facilitate the calculation of the Green function, we initialize the
algorithm in the following way:

$\bm{q_{1}}=\frac{\bm{q}}{|\bm{q}|}$ (145)
$\bm{p_{1}}=\bm{p}\frac{|\bm{q}|}{\bm{p}\cdot\bm{q}}$ (146)

The response function can be rewritten as:

$\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)=({\bm{p}}\cdot{\bm{q}})\;\bm{p_{1}}\cdot(\mathcal{L}-\omega^{2})^{-1}\bm{q_{1}}$
(147)
$\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)=({\bm{p}}\cdot{\bm{q}})\;\bm{p_{1}}QQ^{-1}\cdot(\mathcal{L}-\omega^{2})^{-1}QQ^{-1}\bm{q_{1}}$
(148)

Thanks to the definitions in Eq. (142) and (143) we have

$Q^{-1}\bm{q_{1}}=\begin{pmatrix}1\\\ 0\\\ \vdots\\\ 0\end{pmatrix}$ (149)
$\bm{p_{1}}Q=\begin{pmatrix}1&0&\cdots&0\end{pmatrix}$ (150)

Thus, the response function is the first element of the inverse matrix in the
tridiagonal basis:

$\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)=({\bm{p}}\cdot{\bm{q}})\left[\left(\mathcal{T}-\omega^{2}\right)^{-1}\right]_{11}$
(151)

We can compute the first element of the inverse of a tridiagonal matrix as a
continued fraction recursion:

$\chi_{{\mathcal{A}}{\mathcal{B}}}(\omega)=\frac{({\bm{p}}\cdot{\bm{q}})}{\alpha_{1}-\omega^{2}-\frac{\gamma_{1}\beta_{1}}{\alpha_{2}-\omega^{2}-\frac{\gamma_{2}\beta_{2}}{\ddots}}}$
(152)

Indeed, the application of this method is limited to response function where
the perturbation and the response are not orthogonal, i.e. the product
$\bm{p}\cdot\bm{q}\neq 0$. The condition $\bm{p}\cdot\bm{q}=0$ is met, for
example, by off-diagonal elements of the interacting one-phonon Green
function. Luckily, the spectral function depends only on the trace of the
interacting one-phonon Green function, so only diagonal elements are required
(for which $\bm{p}\cdot\bm{q}=1$).

When the recursion is truncated without reaching the stopping condition, we
can extend the $\alpha$,$\beta$ and $\gamma$ above $N_{\text{max}}$ by
assuming them as constant:

$k\geq N_{\text{max}}$
$\alpha_{k}=\alpha_{N_{\text{max}}}\qquad\beta_{k}=\beta_{N_{\text{max}}}\qquad\gamma_{k}=\gamma_{N_{\text{max}}}.$
(153)

Under this hypothesis, the last part of the continued fraction satisfy the
self-consistent equation

$g(\omega)=\frac{1}{\alpha_{N_{\text{max}}}-\omega^{2}-\gamma_{N_{\text{max}}}\beta_{N_{\text{max}}}g(\omega)},$
(154)

that is solved:

$g(\omega)=\frac{\alpha_{N_{\text{max}}}-\omega^{2}}{2\gamma_{N_{\text{max}}}\beta_{N_{\text{max}}}}-\frac{\sqrt{(\alpha_{N_{\text{max}}}-\omega^{2})^{2}-4\gamma_{N_{\text{max}}}\beta_{N_{\text{max}}}}}{2\gamma_{N_{\text{max}}}\beta_{N_{\text{max}}}}.$
(155)

This introduces an imaginary term to the Green function when the argument
inside the square root becomes negative (branch cut). This termination of the
recursion is similar to the proposed one for the turbo Lanczos algorithm in
TD-DFT[57]. With this method, we compute the response function at any
frequency, with a single tri-diagonalization procedure: the application of
$\mathcal{L}$ to $\bm{p}_{i}$ and $\bm{q}_{i}$ in Eq. (144b) and (144c) is the
expensive part of the algorithm. The $\alpha_{i},\beta_{i},\gamma_{i}$
coefficients depend on the ${\mathcal{A}}$ and ${\mathcal{B}}$ observables (on
the vectors $\bm{p}$ and $\bm{q}$ used to start the Lanczos algorithm) but not
on the frequency $\omega$. The application of $\mathcal{L}$ to a vector does
not require any matrix inversion; it is efficient and parallelizable. In
particular, $\mathcal{L}$ is a square matrix of dimension $3N(1+6N)\sim N^{2}$
($3N$ the $\bm{\mathcal{R}}$ block and $9N^{2}$ both the $\bm{\Upsilon}$ and
$\mathrm{Re}\bm{A}$ blocks). Therefore, the application of the $\mathcal{L}$
matrix to the $\bm{p}_{i}$ and $\bm{q}_{i}$ vectors scales as $N^{4}$, much
better than the inversion of Eq. (136) ($N^{6}$). To calculate the full
response function, we need to apply the $\mathcal{L}$ matrix for each
iteration. Thus the overall computational cost scales as $N^{4}\cdot
N_{\text{steps}}$, where $N_{\text{steps}}$ is the number of iterations of the
Lanczos algorithm.

The heaviest part of applying $\mathcal{L}$ to a vector is the calculation of
the $\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ tensor inside
$\bar{\bm{X}}$ and $\bar{\bm{X}^{\prime}}$. Since
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ is computed as a Monte
Carlo integral over randomly displaced ionic configurations, as discussed in
ref.[33], its calculation scales as $N^{4}\cdot N_{\text{config}}$, where
$N_{\text{config}}$ is the number of random ionic configurations in which used
to evaluate the average of $\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$
(see Appendix H). The number of configurations needed to converge
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ average depends on the
system and on the symmetries; usually $N_{\text{config}}$ is of the order of
thousands.

On the other side, as we show in Sec. VII.2, the Lanczos continued fraction
converges quite well already with $N_{\text{steps}}=40$ iterations. Therefore,
the calculation of the $\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$
matrix has a bigger pre-factor than the full Lanczos procedure, and most of
the computational cost goes in the calculation of
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$.

However, the memory required to store
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ diverges quickly. In big
systems, it may be convenient to recompute
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ on-the-fly for each
iteration of the Lanczos, to avoid storing it. In this case the overall
computational cost of the Lanczos algorithm scales as $N^{4}\cdot
N_{\text{steps}}N_{\text{config}}$. Since $N_{\text{steps}}$ is small, we can
compute the response function for systems of hundreds of atoms even on a
standard desktop PC, without requiring a prohibitive amount of RAM to store
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ with a moderate additional
computational cost.

##  VII Applications

We now illustrate two applications of the TD-SCHA. First, as a simple
benchmark, we show how the TD-SCHA performs in a toy model composed by a
single one-dimensional particle in an external anharmonic potential, comparing
the results with the exact (numerical) solution of the problem (Sec. VII.1).

Then, we apply the TD-SCHA in linear response regime to a real physical
system, the phase III of solid high-pressure hydrogen, and simulate the IR and
Raman spectroscopy. This system is challenging for any numerical technique
that aims to calculate quantum nuclear time-correlation functions _ab initio_
: it is strongly anharmonic, quantum effects are important and we have many
atoms in the simulation cell ($N=96$).

###  VII.1 1D particle in an external potential

To benchmark the TD-SCHA, we compare it with the numerical solution of the
Schroedinger equation. This is possible only for very small systems. In this
section, we simulate a one-dimensional particle in a strongly anharmonic
external potential. The external potential $V(R)$ plays the role of Born-
Oppenheimer energy landscape, its shape is reported in Figure 2 (panel a). In
Hartree atomic units it is:

$V(R)=3R^{4}+\frac{1}{2}R^{3}-3R^{2}$ (156)

This example is performed at $T=$0\text{\,}\mathrm{K}$$, where the density
matrix is a pure state. Our particle has the mass of an electron. We compare
the probability distribution of the ground state wave-funciton
$\ket{\psi_{GS}}$ (obtained diagonalizing the real Hamiltonian) with the SCHA
one in panel a of Figure 2. The SCHA solution is a Gaussian, while the exact
ground state is a complex function. In panel b of the same figure, we present
the linear response to a dynamical perturbation. We plot the one-phonon
spectral function multiplied by $\omega$, defined as:

$S(\omega)=-\frac{\omega}{\pi}\mathrm{Im}\mathcal{G}(\omega).$ (157)

In this way, $S(\omega)$ satisfy the sum rule that its integral is
proportional to the number of phonons[22] (the intensity of a peak does not
depend on the phonon frequency). We compare the results computed with the
harmonic approximation and the TD-SCHA to the exact (numerical) solution.

The exact solution is calculated with the Lehmann representation of the Green
function, where we performed the full diagonalization on the Hamiltonian
$\hat{H}$ to obtain the excited sates $\ket{\psi_{i}}$:

$\displaystyle{\mathcal{G}}(\omega)$
$\displaystyle=\sum_{i=1}^{\infty}\frac{\left|\braket{\psi_{i}}{\sqrt{m}(\hat{R}-\mathcal{R})}{\psi_{GS}}\right|^{2}}{\hbar\omega-
E_{i}+E_{GS}+i0^{+}}+$
$\displaystyle+\sum_{i=1}^{\infty}\frac{\left|\braket{\psi_{i}}{\sqrt{m}(\hat{R}-\mathcal{R})}{\psi_{GS}}\right|^{2}}{\hbar\omega+E_{i}-E_{GS}+i0^{+}}$
(158) $\mathcal{R}=\braket{\psi_{GS}}{\hat{R}}{\psi_{GS}}$ (159)

Since we have a 1D system, we can calculate the TD-SCHA response analytically.
In particular, we can write the $\mathcal{L}$ matrix explicitly (Eq. 137):

$\mathcal{L}=\begin{pmatrix}-\frac{\overset{\scriptscriptstyle{\text{(4)}}}{D}}{2\omega_{s}}-2\omega_{s}^{2}&4\omega_{s}^{2}&4\omega_{s}\overset{\scriptscriptstyle{\text{(3)}}}{D}\\\
0&0&0\\\
\frac{\overset{\scriptscriptstyle{\text{(3)}}}{D}}{8\omega_{s}^{2}}&0&-\omega_{s}^{2}\end{pmatrix}$
(160)

where $\omega_{s}$ is the frequency of the self-consistent harmonic
Hamiltonian of the equilibrium SCHA. Since we are at
$T=$0\text{\,}\mathrm{K}$$ the line of $\mathcal{L}$ that represent the
evolution of $\mathrm{Re}A^{(1)}$ is zero, as the system is in a pure state.
Therefore, $\mathrm{Re}A$ does not affect the dynamics and we can remove it:

$\mathcal{L}=\begin{pmatrix}-\frac{\overset{\scriptscriptstyle{\text{(4)}}}{D}}{2\omega_{s}}-2\omega_{s}^{2}&4\omega_{s}\overset{\scriptscriptstyle{\text{(3)}}}{D}\\\
\frac{\overset{\scriptscriptstyle{\text{(3)}}}{D}}{8\omega_{s}^{2}}&-\omega_{s}^{2}\end{pmatrix}$
(161)

The full anharmonic Green function is:

$\bm{G}(\omega)=-(\mathcal{L}-\omega^{2})^{-1}$ (162)

The one-phonon Green function ${\mathcal{G}}(\omega)$ is the element in the
$\bm{\mathcal{R}}$ block (the last) of the full Green function
$\bm{G}(\omega)$. Performing the inversion analytically we can compute the
self-energy:

${\mathcal{G}}(\omega)^{-1}=\omega^{2}-\omega_{s}^{2}-\Sigma(\omega)$ (163)
$\Sigma(\omega)=-\frac{\overset{\scriptscriptstyle{\text{(3)}}}{D}^{2}/(2\omega_{s})}{\omega^{2}-4\omega_{s}^{2}-\overset{\scriptscriptstyle{\text{(4)}}}{D}/(2\omega_{s})}$
(164)

As we proved for the general case in Sec. V.4, this self-energy is the same as
the one obtained by exploiting the dynamical ansatz formulated in ref.[33]:

$\Sigma(\omega)=\overset{\scriptscriptstyle{\text{(3)}}}{D}\Lambda(\omega)\left[1-\overset{\scriptscriptstyle{\text{(4)}}}{D}\Lambda(\omega)\right]^{-1}\overset{\scriptscriptstyle{\text{(3)}}}{D}$
(165)

where in 1D at $T=$0\text{\,}\mathrm{K}$$ we have:

$\Lambda(\omega)=-\frac{1}{4\omega_{s}^{2}}\frac{1}{4\omega_{s}^{2}-\omega^{2}}$
(166)

As shown in Figure 2, the strong anharmonicity downshifts the energy of the
peak by 50 % of the harmonic result. Here, the TD-SCHA displays a relative
error with the exact solution on the phonon energy of 5 %; an important
improvement from the 100 % of the harmonic approximation. The TD-SCHA still
slightly overestimates the vibrational energy. This is a general feature of
the method, as the TD-SCHA wave-function is more rigid than the real one, as
we constrained its Gaussian form.

Figure 2: One dimensional anharmonic model. Panel a: The Born-Oppenheimer
energy landscape $V(R)$ is plotted with the exact ground state $\psi(R)$
probability density (full diagonalization) and the SCHA equilibrium
distribution $\rho(R)$. Panel b: the dynamical spectral function of the model.
Comparison between the exact result (full diagonalization), the linear
response from TD-SCHA, and the Harmonic approximation. The finite line-width
arise from a smearing we introduced for presentation purposes.

Moreover, the TD-SCHA spectral function displays an overtone, that arises from
the anharmonic coupling between the oscillations of the average position and
the fluctuations. This overtone is produced by the pole of the self-energy in
Eq. (164). We plot, in Figure 3, the spectral function zoomed in the energy
region where the overtone appears. In this example, it is 100 times smaller
than the principal peak. This is a consequence that the overtone is off-
resonant: there is not a one-phonon excitation close to its frequency
triggering a Fermi resonance. The TD-SCHA correctly reproduces the overtone
intensity but with a relative error on its energy of a 60%.

Figure 3: Zoom in the spectral function at high frequencies. We plot the
comparison between the exact overtone and the TD-SCHA overtone.

Interestingly, the energy of the overtone is not twice the energy of the
fundamental phonon mode, as there is a very strong anharmonicity. The TD-SCHA
is unable to simulate third or higher overtones, as we have only two poles
from the Green function that are the zeros of the determinant of the inverse
of $\bm{G}(\omega)$.

Finally, we can simulate the wave-function dynamics adding a time-dependent
external potential in a nonlinear regime. For this purpose, we integrate the
TD-SCHA equation of motion comparing the result with the exact evolution. We
introduce, at $t_{0}=0$, a perturbation of the form:

$V^{(1)}(R,t)=E_{0}R\sin(\omega_{0}t)$ (167)

where $E_{0}=1\,\text{Ha/Bohr}$, and $\omega_{0}=1$ Ha. The results are shown
in Figure 4.

Figure 4: Time evolution of the wave-function in a time-dependent external
potential. The simulation starts from the equilibrium configuration of Figure
2. Panel a: time evolution of the modulus square of the exact wave-function.
Panel b: Time evolution of the TD-SCHA Gaussian distribution. Panel c: average
position of the particle as a function of time. Panel d: quantum dispersion
(mean square displacement).

The TD-SCHA time evolution is very close to the exact dynamics in the first
oscillations, where the wave-function is well localized. The two solutions
deviate when the exact wave-function becomes delocalized: The TD-SCHA
probability distribution is more localized during the whole dynamics, as shown
by the mean square displacement reported in panel d. However, the average
position of the particle matches very nicely the exact result (panel c). The
TD-SCHA can reproduce well the nuclear dynamics even in the nonlinear regime.

###  VII.2 High-pressure hydrogen

In this section, we employ the linear response theory of the TD-SCHA to
calculate the Raman and spectra of high-pressure hydrogen phase III. The
simulation of vibrational spectra in high-pressure phases of molecular
hydrogen is very important to dissect the crystal geometry, as it is the only
experimental signature directly related to the lattice. In fact, both X-ray
spectroscopy and neutron scattering are extremely challenging for the small
samples of hydrogen available. Thus, the structure identification is possible
only by comparing results from _ab initio_ simulations with experimental data.
The physical interpretation of the results and their comparison with
experimental data has been presented in ref.[26]. Here, instead, we focus on
the details of the calculation of the Raman response at
$260\text{\,}\mathrm{GPa}$ ($T=$0\text{\,}\mathrm{K}$$). Hydrogen phase III is
a monoclinic crystal of symmetry group _C2/c_ , with 24 atoms in the primitive
cell.

We computed only the one-phonon contribution to the full response function: we
approximate the Raman tensor as independent from the ionic displacement equal
to the value on the SCHA equilibrium positions, as described in Sec. V.5 and
V.6.

The Lanczos algorithm introduced in Sec. VI enables the calculation of the
dynamical response accounting also for the contribution from 4-phonon
scattering processes non perturbatively
($\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$). In Figure 5 we report
the result of the Raman spectra by progressively switching on anharmonicity.
In panel a, we report the harmonic Raman spectra, computed within perturbation
theory. As an illustration, we report in panel b the Raman spectra from the
SCHA equilibrium self-consistent harmonic Hamiltonian (Eq. 31). Here, phonons
have an infinite lifetime, as in the harmonic case. However, the peak
positions and intensities are already strongly affected by anharmonicity
through the equilibrium SCHA self-consistency. This response is what we would
obtain when cubic anharmonicity is negligible as discussed in Sec. V.3, and is
equivalent to the phonon spectra presented in the original derivation of the
static SCHA[20, 58, 59].

Indeed, phonons defined from the static SCHA self-consistent harmonic
Hamiltonian are not the true dynamical response. To get the real Raman
spectrum, one has to calculate the dynamical response function, as introduced
in Sec. V.1. In panel c,d, we report the anharmonic Raman spectrum calculated
within TD-SCHA and the Lanczos algorithm introduced in Sec. VI. Panel c
presents the results with
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$, which is equivalent to
add a “bubble diagram” to the interacting one-phonon self-energy applied to
the equilibrium SCHA Green function. This approximation has been employed for
the dynamical response within the SCHA in most studies[33, 22, 30, 31]. In
panel d, we report the full anharmonic Raman spectra, also accounting for the
full expression of the self-energy, including also the four-phonon scattering
process explicitly $\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$. Notice
that nonperturbative quartic anharmonicity is also accounted for within the
self-consistency of the SCHA Hamiltonian in panels b,c.

Figure 5: Raman spectra of high-pressure hydrogen phase III
($260\text{\,}\mathrm{GPa}$, $0\text{\,}\mathrm{K}$). a: Harmonic spectra. b:
Anharmonic spectra of static SCHA quasi-particles. It is obtained neglecting
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}$ and
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ in the TD-SCHA response
function (Sec. V.3) c: Anharmonic spectra within TD-SCHA, where we set
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$ in the calculation of the
response function Eq. (101). d: Full anharmonic spectra within TD-SCHA. Here,
the full response function was computed with
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}\neq 0$.

The downshift of the vibron peak (the strongest peak in the spectrum at high
frequencies) from the Harmonic (panel a) to the full anharmonic spectrum
(panel d) occurs already when considering equilibrium SCHA auxiliary phonons.
The vibron acquires a broad shape and a finite lifetime when we calculate the
response function with the Lanczos considering
$\overset{\scriptscriptstyle{\text{(3)}}}{\bm{D}}\neq 0$. The contribution of
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}$ here is not negligible, as
it further reduces the lifetime, deviating from the Lorentzian shape.

In Figure 6 we report how the Raman spectrum depends on the number of
iterations in the Lanczos algorithm, showing the convergence of the algorithm.

Figure 6: Convergence of the Raman signal with the number of steps in the
Lanczos chain. In panel a we report the Raman signal with
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$ and 41 Lanczos
iterations. In panel b, we report the same signal after 200 iterations.

From the comparison between Figure 6, we show how the spectrum becomes more
refined as we increase the number of Lanczos iterations. Most of the changes
in the spectra occur in the low-frequency region. On the other side, the
algorithm describes the vibron shape very well already with few iterations.
The good convergence with a limited number of iteration of the algorithm is
mainly due to the use of a terminator, as discussed in Sec. VI and Appendix H.
In Figure 7, we report the spectra obtained with 200 iterations with and
without the terminator, displaying how it improves the convergence of the
algorithm. Further details are reported in Appendix H.

Figure 7: Raman signal computed within the TD-SCHA Lanczos algorithm and
$\overset{\scriptscriptstyle{\text{(4)}}}{\bm{D}}=0$ after 200 iterations,
with and without the terminator of Eq. (155).

Since the energy of the vibron is more than $3000\text{\,}{\mathrm{cm}}^{-1}$,
the temperature needed to populate excited states is above
$4000\text{\,}\mathrm{K}$. For this reason, there is no hope to simulate this
lattice vibration with standard AIMD, that neglects quantum fluctuations. For
this reason, previous studies on high-pressure hydrogen that neglected quantum
effects reported a much more modest anharmonicity in the Raman signal[60, 61].

These results indicate that linear response of TD-SCHA can tackle complex open
physical problems from first principles, providing an unprecedented precision
on the evaluation of dynamical correlation functions that is sufficient to
compare with experimental results[26].

##  VIII Conclusions

In this work, we introduced a new time-dependent theory for lattice dynamics
based on the self-consistent harmonic approximation that can be applied to
real systems with a first-principle treatment of the electrons. We discussed
the linear response on the static SCHA equilibrium solution of the new
equations and derived an efficient algorithm that computes the response
function at any frequency with a single calculation. This algorithm, with the
simple knowledge of two observables of ionic positions, can compute
efficiently time-correlation function, with the inclusion of thermal and
quantum effects. Notably, this result is not empirical but derived from a
grounded basis: the least action principle. We benchmark the theory both in a
simple strongly anharmonic one-dimensional case and on phase III of solid
hydrogen, a real complex molecular crystal, with _ab initio_ treatment of
electrons.

The theory proved to be able to accurately reproduce complex phonon line-width
that arises from strong anharmonic coupling. The TD-SCHA theory paves the way
to predict from first-principles the outcome of most experiments that probe
the lattice dynamics, like IR and Raman spectroscopy, Neutron, and X-Ray
scattering. The theory can be employed in strongly anharmonic regimes, even
when quantum fluctuations are dominant and where perturbative approach fails,
as in cryogenic conditions, when light ions are present, or in systems close
to a second order structural phase transition, as multiferroics, charge
density waves, ferroelectrics and thermoelectric materials. Moreover, TD-SCHA
describes also out-of-equilibrium dynamics, enabling the study of chemical
reactions with light atoms, as the proton transfer in biomolecules, and the
simulation of pump-probe spectroscopies.

##  Appendix A Equations of motion

In this section, we report the details of the derivation of the dynamical
equation of motion in Eq. (54e).

We need to cast the density matrix definition Eq. (48) into the TD-SCHA
equation (Eq. 46).

First, we compute the left-hand side of the dynamical equation of motion:

$i\hbar\frac{d}{dt}\hat{\rho}$

where $\hat{\rho}$ is the general Gaussian density matrix defined in Eq. (48).
We write the density matrix in real space $\rho(\bm{R},\bm{R}^{\prime})$.

$\displaystyle\frac{1}{\rho}\frac{d\rho}{dt}$
$\displaystyle=\sum_{ab}(R_{a}-\mathcal{R}_{a})(R_{b}-\mathcal{R}_{b})\left[i\frac{d}{dt}C_{ab}-\frac{1}{4}\frac{d}{dt}\Theta_{ab}\right]+$
$\displaystyle+\sum_{ab}(R_{a}^{\prime}-\mathcal{R}_{a})(R_{b}^{\prime}-\mathcal{R}_{b})\left[-i\frac{d}{dt}C_{ab}-\frac{1}{4}\frac{d}{dt}\Theta_{ab}\right]+$
$\displaystyle+\sum_{ab}(R_{a}^{\prime}-\mathcal{R}_{a})(R_{b}-\mathcal{R}_{b})\frac{d}{dt}A_{ab}+$
$\displaystyle-\sum_{ab}(R-\mathcal{R}_{a})\left[2iC_{ab}-\frac{1}{2}\Theta_{ab}+A_{ab}\right]\frac{d}{dt}\mathcal{R}_{b}+$
$\displaystyle-\sum_{ab}(R^{\prime}-\mathcal{R}_{a})\left[-2iC_{ab}-\frac{1}{2}\Theta_{ab}+A_{ab}\right]\frac{d}{dt}\mathcal{R}_{b}+$
$\displaystyle-i\sum_{a}(R_{a}-\mathcal{R}_{a})\frac{dQ_{a}}{dt}+i\sum_{a}(R^{\prime}_{a}-\mathcal{R}_{a})\frac{dQ_{a}}{dt}+\frac{1}{N}\frac{dN}{dt}$
(168)

This is a polynomial in $(\bm{R}-\bm{\mathcal{R}})$ and
$(\bm{R}^{\prime}-\bm{\mathcal{R}})$. On the right-hand side of the TD-SCHA
equation (Eq. 46) we have the Liouville operator:

${\hat{\mathcal{H}}}[\rho]\hat{\rho}-\hat{\rho}{\hat{\mathcal{H}}}[\rho],$
(169)

where the Hamiltonian is:

$\displaystyle{\hat{\mathcal{H}}}[\rho]$
$\displaystyle=\frac{1}{2}\sum_{ab}(R_{a}-\mathcal{R}_{a})\left\langle\frac{d^{2}{V^{(\text{tot})}}}{dR_{a}dR_{b}}\right\rangle_{\rho(t)}(R_{b}-\mathcal{R}_{b})+$
$\displaystyle-\sum_{a}\left\langle
f_{a}\right\rangle_{\rho(t)}(R_{a}-\mathcal{R}_{a})+\sum_{a}\frac{p_{a}^{2}}{2m_{a}}.$
(170)

We compute first the commutator with the kinetic operator:

$p_{a}^{2}\rightarrow-\hbar^{2}\frac{\partial^{2}}{\partial R_{a}^{2}}$

We start computing the first derivative of the $R_{a}$ variable:

$\displaystyle\frac{1}{\rho}\frac{d\rho}{dR_{a}}$
$\displaystyle=\sum_{b}(R_{b}-\mathcal{R}_{b})\left[2iC_{ab}-\frac{\Theta_{ab}}{2}\right]+$
$\displaystyle+\sum_{b}(R_{b}^{\prime}-\mathcal{R}_{b})A_{ab}-iQ_{a}$ (171)

Then, we perform the second derivative:

$\displaystyle\frac{1}{\rho}\frac{d^{2}\rho}{dR_{a}^{2}}$
$\displaystyle=\sum_{bc}(R_{b}-\mathcal{R}_{b})(R_{c}-\mathcal{R}_{c})\left[-4C_{ab}C_{ac}+\frac{\Theta_{ab}\Theta_{ac}}{4}-iC_{ab}\Theta_{ac}-iC_{ac}\Theta_{ab}\right]+$
$\displaystyle+\sum_{bc}(R_{b}^{\prime}-\mathcal{R}_{b})(R_{c}^{\prime}-\mathcal{R}_{c})A_{ab}A_{ac}+\sum_{bc}(R_{b}-\mathcal{R}_{b})(R^{\prime}_{c}-\mathcal{R}_{c})\left(4iC_{ab}A_{ac}-\Theta_{ab}A_{ac}\right)+$
$\displaystyle+\sum_{b}(R_{b}-\mathcal{R}_{b})(4C_{ab}Q_{a}+iQ_{a}\Theta_{ab})-\sum_{b}(R^{\prime}_{b}-\mathcal{R}_{b})2iA_{ab}Q_{a}-Q_{a}^{2}+2iC_{aa}-\frac{1}{2}\Theta_{aa}$
(172)

If we change the derivative of $R$ with $R^{\prime}$ (that is obtained from
$\hat{\rho}{\hat{\mathcal{H}}}[\rho]$), we get the complex conjugate of Eq.
(172) where $R$ and $R^{\prime}$ variables exchanged.

Therefore, we collect all the terms in Eq. (168) and Eq. (169) quadratic in
$\bm{R}$, i.e. the coefficients of
$\sum_{ab}(R_{a}-\mathcal{R}_{a})(R_{b}-\mathcal{R}_{b})$).

$\displaystyle
i\hbar\left(i\frac{dC_{ab}}{dt}-\frac{1}{4}\frac{d\Theta_{ab}}{dt}\right)$
$\displaystyle=-\frac{1}{2}\left\langle\frac{d^{2}{V^{(\text{tot})}}}{dR_{a}dR_{b}}\right\rangle_{\rho}+$
$\displaystyle-\sum_{c}\frac{\hbar^{2}}{2m_{c}}\bigg{[}4C_{ac}C_{cb}-\frac{1}{4}\Theta_{ac}\Theta_{cb}+$
$\displaystyle+2iC_{ac}\Theta_{cb}+A_{ac}A_{bc}\bigg{]}$ (173)

We can split the imaginary and real part to get the first two equation of
motion for $\bm{C}$ and $\bm{\Theta}$:

$\displaystyle\frac{dC_{ab}}{dt}=\frac{1}{2\hbar}$
$\displaystyle\left\langle\frac{d^{2}V}{dR_{a}dR_{b}}\right\rangle_{\rho}+\sum_{c}\frac{\hbar}{2m_{c}}\bigg{[}4C_{ac}C_{cb}+$
$\displaystyle-\frac{1}{4}\Theta_{ac}\Theta_{cb}+\mathrm{Re}\left[A_{ac}A_{bc}\right]\bigg{]}$
(174)
$\displaystyle\frac{d\Theta_{ab}}{dt}=\sum_{c}\frac{2\hbar}{m_{c}}\left(C_{ac}\Theta_{cb}+C_{bc}\Theta_{ca}+\mathrm{Im}\left[A_{ac}A_{bc}\right]\right)$
(175)

We now consider the equality between Eq. (168) and (169) of the coefficient
that multiplies $(\bm{R}-\bm{\mathcal{R}})(\bm{R}^{\prime}-\bm{\mathcal{R}})$.

$\displaystyle i\hbar\frac{dA_{ab}}{dt}=\sum_{c}\frac{\hbar^{2}}{2m_{c}}$
$\displaystyle\big{[}4iC_{ac}A_{cb}+4iC_{bc}A_{ac}+$
$\displaystyle-\Theta_{ac}A_{cb}+\Theta_{bc}A_{ac}\big{]}$ (176)

We can split also this equation in real and imaginary part:

$\displaystyle\mathrm{Re}\frac{dA_{ab}}{dt}=\sum_{c}\frac{\hbar}{2m_{c}}$
$\displaystyle\big{(}4C_{ac}\mathrm{Re}A_{cb}+4C_{bc}\mathrm{Re}A_{ca}+$
$\displaystyle-\Theta_{ac}\mathrm{Im}A_{cb}+\Theta_{bc}\mathrm{Im}A_{ac}\big{)}$
(177)

$\displaystyle\mathrm{Im}\frac{dA_{ab}}{dt}=\sum_{c}\frac{\hbar}{2m_{c}}$
$\displaystyle\big{[}4C_{ac}\mathrm{Im}A_{cb}+4C_{bc}\mathrm{Im}A_{ac}+$
$\displaystyle+\Theta_{ac}\mathrm{Re}A_{cb}-\Theta_{bc}\mathrm{Re}A_{ac}\big{]}$
(178)

The last two equations for $\bm{Q}$ and $\bm{\mathcal{R}}$ are obtained from
the real and imaginary part of the terms linear in $\bm{R}$

$\displaystyle-i\hbar$
$\displaystyle\sum_{b}\left(2iC_{ab}-\frac{1}{2}\Theta_{ab}+A_{ab}\right)\frac{d\mathcal{R}_{b}}{dt}-\hbar\frac{dQ_{a}}{dt}=$
$\displaystyle=\sum_{c}\frac{\hbar^{2}}{2m_{c}}\left(4C_{ca}Q_{c}+iQ_{c}\Theta_{ca}-2iA_{ca}Q_{c}\right)-\braket{f_{a}^{(tot)}}$
(179)

Separating real and imaginary parts we get:

$\displaystyle\sum_{b}$
$\displaystyle\left(2C_{ab}+\mathrm{Im}A_{ab}\right)\frac{d\mathcal{R}_{b}}{dt}-\frac{dQ_{a}}{dt}=$
$\displaystyle=\sum_{b}\frac{\hbar}{m_{b}}\left(2C_{ab}+\mathrm{Im}A_{ab}\right)Q_{b}-\frac{\braket{f^{(tot)}_{a}}}{\hbar}$
(180) $\displaystyle\sum_{b}\hbar$
$\displaystyle\left(\frac{1}{2}\Theta_{ab}-\mathrm{Re}A_{ab}\right)\frac{d\mathcal{R}_{b}}{dt}=$
$\displaystyle=\sum_{b}\frac{\hbar^{2}}{2m_{b}}\left(\Theta_{ab}-2\mathrm{Re}A_{ab}\right)Q_{b}$
(181)

Simplifying we get the last equations:

$\frac{d\mathcal{R}_{a}}{dt}=\frac{\hbar Q_{a}}{m_{a}}$ (182)
$\frac{dQ_{a}}{dt}=\frac{\braket{f_{a}}}{\hbar}$ (183)

The last term of the equivalence between Eq. (168) and (169) is the term that
does depend in neither $(\bm{R}-\bm{\mathcal{R}})$ nor
$(\bm{R}^{\prime}-\bm{\mathcal{R}})$. This last equation expresses the time
evolution for the norm $N(t)$, and it is just the conservation of the density
matrix normalization.

The final equations of motion are obtained substituting the expression of
$\bm{\Upsilon}$ (Eq. 50) and dividing for the masses to obtain the mass-
rescaled matrices (Eq. 53c). The convenience of rescaling the masses is that
we can express all the summations as standard rows-by-columns products.

##  Appendix B Energy conservation

Here, we calculate the time derivative of the average potential (Eq. 66),
showing that it balances the time derivative of the average kinetic energy,
proving that the total energy is conserved by the equation of motions.

For this purpose, we use formalism introduced by Bianco et al[33] (appendix C,
Eq. C1): given an observable $O(\bm{R})$, the average of its derivative can be
written as:

$\left\langle\frac{dO}{dR_{a}}\right\rangle_{\rho}=\sum_{b}\Upsilon_{ab}\left\langle(R_{b}-\mathcal{R}_{b})O\right\rangle_{\rho}$
(184)

For this reason, we can write:

$\sum_{ab}\frac{d\mathcal{R}_{b}}{dt}\Upsilon_{ab}\left\langle(R_{a}-\mathcal{R}_{a})V\right\rangle_{\rho}=-\sum_{b}\frac{d\mathcal{R}_{b}}{dt}\left\langle
f_{a}(\bm{R})\right\rangle_{\rho}$ (185)

For simplicity, we define $\bm{u}=(\bm{R}-\bm{\mathcal{R}})$. The other term
of Eq. (66) is

$\displaystyle\sum_{ab}\left\langle
u_{a}u_{b}V\right\rangle_{\rho}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}$
$\displaystyle=\sum_{abcd}\left\langle
u_{c}u_{b}V\right\rangle_{\rho}\Upsilon_{cd}\Upsilon^{-1}_{da}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}$
$\displaystyle=\sum_{abd}\left\langle\frac{d(u_{b}V)}{dR_{d}}\right\rangle_{\rho}\Upsilon^{-1}_{da}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}$
(186)

The derivative gives:

$\displaystyle\sum_{ab}\left\langle
u_{a}u_{b}V\right\rangle_{\rho}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}$
$\displaystyle=\sum_{ab}\left\langle
V\right\rangle_{\rho}\Upsilon^{-1}_{ba}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}+$
$\displaystyle+\sum_{abd}\left\langle
u_{b}\frac{dV}{dR_{d}}\right\rangle_{\rho}\Upsilon^{-1}_{da}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}$
(187)

The first term is zero, because:

$\sum_{ab}\Upsilon^{-1}_{ba}\frac{d\Upsilon_{ab}}{dt}=\frac{1}{2}\frac{d}{dt}\textup{Tr}\left[\bm{\Upsilon}^{-1}\bm{\Upsilon}\right]=0$
(188)

While for the second term we can proceed again:

$\sum_{ab}\left\langle
u_{a}u_{b}V\right\rangle_{\rho}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}=\sum_{abcde}\left\langle
u_{e}\frac{dV}{dR_{d}}\right\rangle_{\rho}\Upsilon_{ec}\Upsilon^{-1}_{cb}\Upsilon^{-1}_{da}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}$
(189)

$\sum_{ab}\left\langle
u_{a}u_{b}V\right\rangle_{\rho}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}=\sum_{abcd}\left\langle\frac{d^{2}V}{dR_{c}dR_{d}}\right\rangle_{\rho}\Upsilon^{-1}_{cb}\Upsilon^{-1}_{da}\frac{1}{2}\frac{d\Upsilon_{ab}}{dt}$
(190)

Introducing the mass re-scaled for $\bm{\Upsilon}^{-1}$:

${\tilde{\Upsilon}}^{-1}_{ab}=\sqrt{m_{a}m_{b}}\Upsilon^{-1}_{ab},$ (191)

we get the final equivalence:

$\frac{1}{2}\textup{Tr}\left[\left\langle\bm{u}\otimes\bm{u}V\right\rangle_{\rho}\frac{d\bm{\Upsilon}}{dt}\right]=\frac{1}{2}\textup{Tr}\left[\left\langle\frac{d^{2}V}{d\tilde{\bm{R}}d\tilde{\bm{R}}}\right\rangle_{\rho}{\tilde{\bm{\Upsilon}}}^{-1}{\tilde{\bm{\Upsilon}}^{-1}}\frac{d\tilde{\bm{\Upsilon}}}{dt}\right].$
(192)

Substituting Eq. (185) and (192) into Eq. (66), and substituting the equation
of motion, it is easy to show that

$\frac{d}{dt}\left\langle\hat{K}\right\rangle_{\rho(t)}=-\frac{d}{dt}\left\langle\hat{V}\right\rangle_{\rho(t)}$

##  Appendix C Least action

Here, we prove that for a Gaussian pure state expressed in Eq. (40), the
action (Eq. 41) with the TD-SCHA dynamics is stationary.

First we break the action in three parts: kinetic, potential, and time:

$A_{1}=\frac{1}{t_{2}-t_{1}}\int_{t_{1}}^{t_{2}}dt\,\braket{\psi(t)}{\hat{K}}{\psi(t)}$
(193)
$A_{2}=\frac{1}{t_{2}-t_{1}}\int_{t_{1}}^{t_{2}}dt\,\braket{\psi(t)}{\hat{V^{(\text{tot})}}}{\psi(t)}$
(194)
$A_{3}=-\frac{i}{\hbar}\frac{1}{t_{2}-t_{1}}\int_{t_{1}}^{t_{2}}dt\,\braket{\psi(t)}{\frac{d}{dt}}{\psi(t)}$
(195)

by substituting Eq. (40) in Eq. (193), (194) and (195), we get:

$A_{1}=\frac{\hbar^{2}}{t_{2}-t_{1}}\int_{t_{1}}^{t_{2}}dt\left\\{\sum_{a}\frac{Q_{a}^{2}}{2m_{a}}+\textup{Tr}\left[\frac{\tilde{\bm{\Upsilon}}}{8}+2{\tilde{\bm{C}}}_{ab}{\tilde{\bm{\Upsilon}}}^{-1}\tilde{\bm{C}}\right]\right\\}$
(196)

From which the derivatives of the action are:

$\frac{\delta
A_{1}}{\delta\tilde{\bm{C}}}=2\hbar^{2}{\tilde{\bm{\Upsilon}}}^{-1}\tilde{\bm{C}}+2\hbar^{2}\tilde{\bm{C}}{\tilde{\bm{\Upsilon}}}^{-1}$
(197) $\frac{\delta A_{1}}{\delta
Q_{a}}=\frac{\hbar^{2}Q_{a}}{m_{a}}\qquad\frac{\delta
A_{1}}{\delta\bm{\mathcal{R}}}=0$ (198) $\frac{\delta