\section{Introduction}

The pseudogap phase of cuprate superconductors, in particular its link with the high temperature superconductivity, remains an enduring mystery of condensed matter physics.
While no clear phase transition at its characteristic temperature $T^{\star}$ is observed in transport and thermodynamic properties, the low temperature crossing of its critical doping $p^{\star}$ in the absence of superconductivity yields clear signatures of a transition~\cite{proust2019}.
The normal-state Hall coefficient of YBa$_2$Cu$_3$O$_y$ (YBCO)~\cite{badoux2016}, La$_{1.6-x}$Nd$_{0.4}$Sr$_{x}$CuO$_4$ (Nd-LSCO)~\cite{collignon2017} and Bi$_2$Sr$_{2-x}$La$_{x}$CuO$_{\rm 6+\delta}$ (Bi2201)~\cite{lizaire2020transport} shows that the carrier density $n$ changes abruptly from $n \simeq 1 + p$ to $n \simeq p$ when crossing $p^{\star}$ from above.
This is also observed in the electrical resistivity of La$_{2-x}$Sr$_x$CuO$_4$ (LSCO)~\cite{laliberte2016} and Nd-LSCO~\cite{collignon2017}, and the thermal conductivity (at $T\rightarrow 0$) of Nd-LSCO~\cite{michon2018}.
Clear evidence for a transformation of the Fermi surface of Nd-LSCO across $p^{\star}$ was recently obtained from angle-dependent magneto-resistance (ADMR)~\cite{fang2020}. The low temperature electronic specific heat \cite{michon2019,girod2021normal} shows a logarithmic divergence at $p^{\star}$, evidence for a quantum phase transition.

Recently, the thermopower was used as a complementary probe of the carrier density at low temperatures across $p^{\star}$ in Nd-LSCO~\cite{collignon2021thermopower}.
Indeed, in the $T = 0$ limit and for a single parabolic band, the Seebeck coefficient depends on only two parameters: the coefficient of electronic specific heat $\gamma = C_{\rm el}/T$ and the carrier density $n$ ($e$ is the electron charge):
\begin{equation}
\frac{S}{T} = \frac{\gamma}{ne}.
\label{eq:behnia_eq}
\end{equation}
According to this expression, the ratio $S/T$ can therefore be seen as the specific heat per carrier.
Simplified as it is, this physical picture was shown by Behnia \textit{et al.}~\cite{behnia2004}
to account empirically for the observed low-temperature value of the Seebeck coefficient for a great variety of materials that includes common metals and several strongly correlated materials.

In Nd-LSCO, Eq.~(\ref{eq:behnia_eq}) articulates that the sudden increase in the in-plane Seebeck
coefficient $S_{\rm a}/T$ as $p$ crosses below $p^{\star}$~\cite{collignon2021thermopower} is consistent with the drop in $n$ at $p^{\star}$ inferred from Hall effect, electrical resistivity, and thermal conductivity measurements. Although the relative change in $S_{\rm a}/T$ across $p^{\star}$ seems well captured by Eq. (\ref{eq:behnia_eq}), it is unable to explain the positive sign of the in-plane Seebeck coefficient, which is at odds with the electron-like band of Nd-LSCO above $p^{\star}$~\cite{verret2017}. To date, the sign of the Seebeck coefficient in overdoped cuprates remains a mystery.

At a more fundamental level, the Seebeck coefficient is actually a ratio of transport coefficients which
involves states immediately below and above the Fermi level, in contrast with electrical transport
that is only sensitive to the properties of the Fermi surface.
As a result, the Seebeck coefficient is controlled by the particle-hole asymmetry between
occupied and unoccupied states around the Fermi level \cite{behnia_2015, kondo2005}.
This asymmetry can originate both from the dispersion of electronic excitations (band structure)
and from the energy dependence of the scattering rate.
Accounting for these effects clearly goes beyond the simplified expression
Eq.~(\ref{eq:behnia_eq}).

In the present article, we present measurements of the out-of-plane Seebeck coefficient $S_{\rm c}$. We report a sudden qualitative change of $S_{\rm c}$ in Nd-LSCO across $p^{\star}$. We show that in order to successfully describe the behaviour of the Seebeck coefficient along both directions, it is crucial to  take into account the  particle-hole asymmetry (`skewness') in the energy dependence of the scattering rate.

Nd-LSCO is a single-layer, tetragonal cuprate superconductor with low critical temperature $T_c$ and field $H_{c2}$, making it an ideal candidate to study the field-induced normal state down to low temperatures.
Its phase diagram is shown in Fig.~\ref{fig:diagram}, where the pseudogap temperature $T^{\star}$ extracted from resistivity measurements~\cite{collignon2021thermopower} is displayed and seen to be in agreement with angle-resolved photoemission spectroscopy (ARPES) measurements~\cite{matt2015}.

Above $p^{\star}$, at $p = 0.24$, our measurements reveal that $S_{\rm a}$ and $S_{\rm c}$ are positive and essentially equal at low temperature, both ending with a $\log(1/T)$ dependence below about 10~K. Although it is consistent with Eq. (\ref{eq:behnia_eq}), the isotropy of the Seebeck coefficient at $p=0.24$ as $T\rightarrow 0$ cannot be understood in terms of the band structure since it predicts the wrong magnitude of $S_{\rm c}$, and is of the wrong sign for $S_{\rm a}$. We show that taking into account the angle dependence and temperature dependence of the scattering rate deduced from ADMR~\cite{grissonnanche2020b} and adding a linear energy dependence that is particle-hole asymmetric reconciles the calculated and measured Seebeck coefficients along both directions.

Upon crossing into the pseudogap phase at $p^{\star}$, we find that $S_{\rm c}$ becomes negative at low temperature, in contrast to $S_{\rm a}$ which remains positive. This contrasts with Eq. (\ref{eq:behnia_eq}) that predicts that the drop in $n$ at $p^{\star}$ should also be observed in the interlayer $c$ axis component $S_{\rm c}$.
The negative $S_{\rm c}$ reflects a profound transformation of the electronic structure across $p^{\star}$.

\begin{figure}[t!]
\includegraphics[width=0.45\textwidth]{Fig1.pdf}
\caption{Temperature-doping phase diagram of the cuprate superconductor Nd-LSCO. Red points represent the temperature below which the pseudogap phase appears, $T$$^{\star}$, obtained from resistivity \cite{daou2009,collignon2017} (full points show the dopings studied here) and ARPES~\cite{matt2015} (full square) measurements. The red shaded area represents the pseudogap (PG) phase, which disappears at the critical doping $p^{\star} = 0.23 \pm 0.01$ (red diamond). The full line is a guide to the eye. The superconducting phase (in zero field) is bounded by $T_{\rm c}$ (dashed line).
}
\label{fig:diagram}
\end{figure}

\section{Methods}

Single crystals of Nd-LSCO were grown by a traveling solvent/floating zone technique in an image furnace, with a Nd content of 0.4 and doping $p$ = 0.20, 0.21, 0.22, 0.23 and 0.24 (for more details see Collignon et al. \cite{collignon2017}).
These were subsequently cut into bar-shaped samples with typical dimensions of 1 mm x 0.5 mm x 0.2 mm.
For each doping, separate samples with their length respectively along the $a$ axis and the $c$ axis were cut.

We measured the Seebeck coefficient using an AC-technique (Alternative Current) that was originally developed for thin films~\cite{wang2019} and that we have adapted for measuring bulk materials.
An AC thermal excitation was generated through the sample by sending an electric current at a frequency $\omega \sim 1$~Hz to a 5k$\Omega$ strain gauge used as a heater. This generates a longitudinal thermal gradient $\Delta T_{\rm AC}$ along the length of the sample, measured at a frequency $2\omega$ using two absolute type E thermocouples.
In response to that thermal gradient, a Seebeck voltage $\Delta V_{\rm AC}$ is measured at a frequency $2\omega$ with phosphor-bronze wires using the same contacts as for $\Delta T_{\rm AC}$, which eliminates uncertainties associated with the geometric factor.
The Seebeck coefficient is then given by $S = - \Delta V_{\rm AC} / \Delta T_{\rm AC}$.
Our method differs from previous AC Seebeck measurements by accounting for both the modulus and the phase of the thermal and thermoelectric signals.
This technique can be in principle extended to all orders in frequencies.
The advantages of using an AC method over DC is that the DC technique requires to measure the Seebeck coefficient at constant temperature steps and wait the necessary time to reach equilibrium, which can take minutes and even hours closer to room temperature, it is also limited in resolution to few milliKelvins fluctuations for the thermal gradient and few nanoVolts for the Seebeck voltage.
The AC Seebeck technique on the other hand does not require for the sample to reach equilibrium, which makes the measurements quasi-instantaneous, and comes down to a resolution in the $\sim 10$~$\mu$K range because of the use of lock-in amplifiers.
However, it takes a real time measurement of the sample temperature, which comes with its own set of challenges.

The thermocouple and Seebeck voltages were amplified using EM Electronics A10 preamplifiers and picked-up using SR830 lock-in amplifiers at the thermal excitation frequency $2\omega$.
Our method allows us to measure $S(T)$ continuously from 2 to 300K within a few hours and with a much better signal to noise ratio compared to a standard steady-state DC technique. It was carefully benchmarked against a steady-state DC method on several samples.\\

\begin{figure}[t!]
\includegraphics[width=0.45\textwidth]{Fig2.pdf}
\caption{In-plane ($S_{\rm a}$, blue) and out-of-plane ($S_{\rm c}$, red) Seebeck coefficient of Nd-LSCO, plotted as $S/T$ vs T at $H=16$~T for (\textbf{a}) $p = 0.24>p^{\star}$ and (\textbf{b}) $p~=~0.22~<~p^{\star}$.}
\label{fig:SoT_0p24_0p22}
\end{figure}

\section{Results}

\subsection{Outside the pseudogap phase $p > p^{\star}$}

In Fig.~\ref{fig:SoT_0p24_0p22}(a) we show $S_{\rm a}/T$ and $S_{\rm c}/T$ as a function of temperature in Nd-LSCO just outside the pseudogap phase, at $p = 0.24$, in the normal state induced by an applied field of $H = 16$~T ($H_{\rm c2} = 9$~T in Nd-LSCO and Eu-LSCO at $p = 0.24$ \cite{michon2019}).
We observe that $S_{\rm a}$ is negative at room temperature, similar to what is observed in overdoped Bi2201~\cite{kondo2005}, for example.
Upon cooling, $S_{\rm a}/T$ increases and changes sign around 250 K, reaching a positive value of about 0.45 $\mu$V/K$^{2}$ at the lowest measured temperature $T \simeq 2$~K.
Conversely, $S_{\rm c}/T$ starts from a positive value at 300~K, it increases until it reaches a plateau around 60 K, but at lower temperature starts to increase again.
Below about 10~K a remarkable isotropy appears: $S_{\rm a}/T$ and $S_{\rm c}/T$ are essentially identical despite the quasi two-dimensional character of the
electronic structure and associated anisotropy in the electrical resistivity (in Nd-LSCO at $p = 0.24$ the $c$ axis resistivity $\rho_c$ is approximately 250 times larger than $a$ axis resistivity $\rho_a$~\cite{daou2009}).

\begin{figure}[t!]
\includegraphics[width=0.45\textwidth]{Fig3.pdf}
\caption{In-plane ($S_{\rm a}$, blue) and out-of-plane ($S_{\rm c}$, red) Seebeck coefficients of Nd-LSCO at $p = 0.24> p^{\star}$ plotted as $S/T$ vs $\log(T)$ at $H=16$~T. Note the logarithmic dependence below $\sim 10$~K, where $S$ is essentially isotropic ($S_{\rm a} \simeq S_{\rm c}$).}
\label{fig:logT_0p24}
\end{figure}

As seen in Fig. \ref{fig:logT_0p24}, zooming in on the low temperature data reveals that $S_{\rm a}/T$ and $S_{\rm c}/T$ both display a $\log(1/T)$ dependence below $\sim 10$~K, as previously reported for $S_{\rm a}/T$ in Nd-LSCO~\cite{daou2009S}, Eu-LSCO~\cite{laliberte2011} at $p~=~0.24$ and in Bi2201 near $p^{\star}$~\cite{lizaire2020transport}.
This behavior is also observed in the electronic specific heat of Nd-LSCO and Eu-LSCO at $p = 0.24$ whereby $C_{\rm el} / T \sim \log(1/T)$ below $\sim 10$~K is interpreted as a signature of quantum criticality in the vicinity of the pseudogap critical point~\cite{michon2019,girod2021normal}.
Another typical signature of quantum criticality is a $T$-linear dependence of the electrical resistivity as $T \to 0$ observed in Nd-LSCO at $p = 0.24$ both in the plane and along the $c$ axis~\cite{daou2009}.

\begin{figure}[t]
\includegraphics[width=0.42\textwidth]{Fig4.pdf}
\caption{(\textbf{a}) $S_{\rm a}/T$ and (\textbf{b}) $S_{\rm c}/T$ as a function of temperature, at $H=16$~T, in Nd-LSCO at dopings as indicated, on both sides of the pseudogap critical point $p^\star = 0.23$. Note the qualitative change in $S_{\rm c}$ upon crossing below $p^\star$, whereby $S_{\rm c}$ suddenly becomes negative at low $T$.}
\label{fig:SoT_doping}
\end{figure}

\subsection{Inside the pseudogap phase $p < p^{\star}$}

Now going inside the pseudogap phase, we see in Fig.~\ref{fig:SoT_doping}(a) that $S_{\rm a}/T$ undergoes a large enhancement at low temperature for $p < 0.23$.
This was previously reported and discussed in detail in ref.~\cite{collignon2021thermopower} and here we show the continuous $T$-dependent curves of $S_{\rm a}/T$ in $H = 16$~T obtained using our AC technique (our AC data are in excellent agreement with the previous DC data \cite{collignon2021thermopower}).
As seen in the Hall effect and electrical resistivity of Nd-LSCO~\cite{collignon2017}, the enhancement of $S_{\rm a}/T$ occurs below $T^{\star}$ and reaches a maximum as $T \rightarrow 0$ that corresponds roughly to a 5-fold increase between $p = 0.24$ and $p = 0.20$.
In agreement with Eq. (\ref{eq:behnia_eq}), this roughly matches the change in carrier density from $n \simeq 1 + p = 1.24$ at $p = 0.24$ to $n \simeq p = 0.20$ at $p = 0.20$, suggesting that it is the principal cause for the enhancement of these three transport coefficients ($R_{\rm H}$, $\rho$, $S$) inside the pseudogap phase.

Turning to $S_{\rm c}/T$ below $p^{\star}$, at $p = 0.22$, we see in Fig.~\ref{fig:SoT_0p24_0p22} that it tracks $S_{\rm c}/T$ at $p = 0.24$ down to about 100~K, but then drops upon further cooling to reach negative values below $T \sim 20$~K.
This behavior is seen for all the measured dopings at and below $p^{\star}$ (Fig.~\ref{fig:SoT_doping}(b)), showing that it is a property of the pseudogap phase.
(We stress that the field dependence of $S_{\rm c}$ is very weak, as seen in Fig. S1(b), where the curves taken at $H = 0$ and 16~T essentially overlap, establishing that the sole effect of the field is to suppress superconductivity and reveal $S_{\rm c}/T$ at low temperatures, not to induce the negative $S_{\rm c}/T$.)

As we discuss below, the highly contrasting behavior between $S_{\rm a}/T$ and $S_{\rm c}/T$ below $p^{\star}$ contradicts Eq.~(\ref{eq:behnia_eq}) that predicts an isotropy between the two current directions, and is a puzzle for standard theories of the Seebeck effect, likely carrying significant hints for the underlying nature of the pseudogap phase.
Whatever mechanism causes the drop in carrier density from $n = 1 + p$ to $n = p$ at $p^{\star}$ and the associated Fermi surface transformation observed by ADMR measurements~\cite{fang2020}, it seems to also affect the asymmetry of the band dispersion around the Fermi level in an anisotropic fashion, as attested by the negative $S_{\rm c}/T$.
$S_{\rm c}$ therefore adds to the list of transport properties undergoing a dramatic change upon crossing $p^{\star}$, such as the negative thermal Hall effect~\cite{grissonnanche2019}, attributed to chiral phonons~\cite{grissonnanche2020}, that suddenly appears in the pseudogap phase of Nd-LSCO.

\section{Discussion}

\subsection{$p > p^{\star}$: particle-hole asymmetric energy-dependent scattering rate}

Previous studies of the Seebeck coefficient of cuprates have largely focused on the behavior above $T_c$, revealing that the in-plane $S(T)$ decreases with increasing doping~\cite{obertelli1992} and is negative at high temperature and high doping, as seen for instance in Bi$_{\rm 2}$Sr$_{\rm 2}$CaCu$_{\rm 2}$O$_{\rm 8+x}$ (Bi2212)~\cite{munakata1992}, LSCO~\cite{elizarova2000}, Bi2201~\cite{kondo2005}, and HgBa$_{\rm 2}$CuO$_{\rm 4+\delta}$~\cite{yamamoto2000}.
Nd-LSCO follows this trend as a function of doping at $T$~=~300~K, as showed in Fig.~S2, and as a function of temperature in Fig.~\ref{fig:SoT_0p24_0p22}(a).
In Nd-LSCO above~$p^{\star}$, at $p = 0.24$, $S_{\rm a}(T)$ is slightly negative at 300~K but becomes positive below about $T \sim 250$~K and keeps increasing upon further cooling (Fig.~\ref{fig:SoT_0p24_0p22}(a)).
This behavior contrasts with overdoped Bi2201 where $S_{\rm a}(T)$ is negative
at all temperatures at dopings close to $p^{\star}$~\cite{lizaire2020transport}.

\begin{figure}[t]
\includegraphics[width=0.4\textwidth]{Fig5.pdf}
\caption{Calculations of the Seebeck coefficient ($S/T$, full line, left axis) and conductivity ($\sigma$, dashed line, right axis) as a function of doping from the bare band dispersion of Nd-LSCO, as measured by ARPES \cite{matt2015,horio2018}, with an energy-independent scattering rate $1/\tau$. (\textbf{a}) In-plane coefficients at $T=6$~K; (\textbf{b}) Out-of-plane coefficients at $T=6$~K. We can see from the figures that the sign of $S/T$ is set by the slope of the conductivity $\sigma$ as a function of doping. This reflects the Mott formula in Eq.~(\ref{eq:Mott_formula}). The minimum in $S_{\rm a}/T$ and maximum in $S_{\rm c}/T$ coincide approximately with the van Hove point ($p_{\rm vHs}$, gray band) where the Fermi surface goes from hole-like (below $p_{\rm vHs}$) to electron-like (above $p_{\rm vHs}$).}
\label{fig:calc_vs_p}
\end{figure}

\begin{figure*}[t]
\includegraphics[width=0.65\textwidth]{Fig6.pdf}
\caption{Left panels: Seebeck coefficient of Nd-LSCO $p~=~0.24$ for both $a$-axis (blue) and $c$-axis (red). Right panels: Corresponding calculated Seebeck coefficients for a momentum-dependent but energy-independent scattering rate $1/\tau_{\rm 0}(\vecck) = A + B |\cos (2\phi)|^\nu$ (see Appendix~\ref{app:seebeck_calculations} for more details) inferred from the elastic scattering rate extracted from ADMR measurements \cite{grissonnanche2020b} (dashed line) and a momentum- and energy-dependent skewed scattering rate of the form $1/\tau_{\rm 0}(\vecck) + 1/\tau_{\rm SMFL}(\epsilon, T)$ (full line).
}
\label{fig:calc_anisotropic_vs_T}
\end{figure*}

In order to understand the observed behavior of $S_{\rm a}/T$ and $S_{\rm c}/T$, we turn to
Boltzmann transport theory in a relaxation time approximation, using the experimentally determined (quasiparticle) band dispersion $E(\vecck)$
of Nd-LSCO at $p = 0.24$ obtained from ADMR~\cite{fang2020,grissonnanche2020b} and ARPES measurements~\cite{matt2015,horio2018} (see Appendix~\ref{app:seebeck_calculations} for details).
It is certainly not immediately obvious that the Boltzmann formalism can be successfully used to describe the non-Fermi liquid regime considered in the present work.
The justification of this procedure and a discussion of its limitations, including a comparison to a full Kubo calculation,
are provided in Appendix~\ref{app:kubo_vs_boltz}.
We start by considering a constant scattering rate $1/\tau$, independent of energy and momentum.
In this case, the Seebeck coefficient does not depend on $\tau$ because it is a ratio of two transport coefficients, the Peltier coefficient $\alpha$ and the conductivity $\sigma$, and it is given by the Mott formula at low temperature:
\begin{equation}
S_{\rm Mott} = - \frac{\pi^2 k_B^2 T}{3e}\frac{\partial \ln \sigma_{\rm ii}(\epsilon) }{\partial \epsilon}\Big\vert_{\epsilon=0} = - \frac{\pi^2 k_B^2 T}{3e}\frac{\sigma_{\rm ii}^{\prime}(\epsilon)}{\sigma_{\rm ii}(\epsilon)}\Big\vert_{\epsilon=0},
\label{eq:Mott_formula}
\end{equation}
In this expression, $e$ is the electron charge, $k_B$ the Boltzmann constant, $i = x, z$ and
the Fermi level has been set at $\epsilon=0$.
$\sigma(\epsilon)$ is the energy-dependent transport function which only depends on the band structure of the material for a constant $\tau$. Hence, the sign and magnitude of $S$ are entirely determined by the particle-hole asymmetry of the band structure in this case.
In Fig. \ref{fig:calc_vs_p} we display the calculated $S$ for both directions
as a function of doping.
These calculations are performed by allowing for the momentum dependence of the scattering rate along the Fermi surface, inferred from the ADMR measurements \cite{grissonnanche2020b}.
We see that the sign of $S$ is simply set by the slope of $\sigma(\epsilon=0)$ vs $p$ (the inverse of the resistivity $\rho$). In Fig. \ref{fig:calc_vs_p}(a), the in-plane electrical conductivity $\sigma_{\rm a})$ decreases monotonically with doping. This results in a negative $S_{\rm a}$ predicted for all overdoped cuprates. The sign of $S_{\rm a}$ does not change in the considered doping range, whereas the sign of $S_{\rm c}$ flips (becomes negative) at high $p$. One can see that crossing the van Hove singularity around $p = p_{\rm vHs}$ has a large effect on the out-of-plane coefficient, but a modest one for the in-plane coefficient. This is because above $p_{\rm vHs}$ the Fermi surface loses states in the antinodal regions, where the $c$ axis dispersion is largest.

We turn to the temperature dependence of the Seebeck coefficients $S_{\rm a}/T$ and $S_{\rm c}/T$ in Nd-LSCO at $p = 0.24$.
Fig.~\ref{fig:calc_anisotropic_vs_T} shows a direct comparison between our experimental data (left column) and the calculations (right column).
We see from Fig.~\ref{fig:calc_anisotropic_vs_T}(b) and (d), that the calculated $S(T)$ assuming an energy-independent (but momentum-dependent) scattering rate (dashed lines) disagrees with the data in magnitude for both directions.
Furthermore, it yields the wrong sign for the in-plane component, $S_{\rm a}$,
as also found in Ref.~\cite{verret2017}.
To summarize, Fig.~\ref{fig:calc_anisotropic_vs_T} shows that neither the sign of $S_{\rm a}$ nor the $T$ dependence and magnitude of $S_{\rm a}$ and $S_{\rm c}$ can be explained when assuming that particle-like and hole-like excitations have equal scattering rates, even when taking into account the momentum dependence of this scattering rate along the Fermi surface.

We now consider the effect of a scattering rate which depends on energy, allowing for this energy dependence to reflect an asymmetry between particle and holes.
The energy dependence of the inelastic scattering rate has long been recognized to be unconventional in cuprates.
It is often described by the marginal Fermi liquid (MFL) form~\cite{varma1989,varma2020}:
$1/\tau_{\rm MFL} = \sqrt{(a \epsilon)^2 + (\alpha \frac{k_{\rm B} T}{\hbar})^2}$.
As written, this expression is symmetric (even) under $\epsilon\rightarrow - \epsilon$, corresponding to equal scattering rates for particles and holes.
Here, we introduce a modification of the MFL scattering rate with a particle-hole asymmetric energy dependence (`skewed' marginal Fermi liquid - SMFL), namely:
\begin{equation}
1/\tau_{\rm SMFL}(\epsilon, T)=\sqrt{(a_{\pm}\epsilon)^2 + (\alpha \frac{k_{\rm B} T}{\hbar})^2}
\label{eq:SMFL}
\end{equation}
in which $a_{+}$ applies to $\epsilon>0$ and $a_{-}$ to $\epsilon<0$. The particle-hole asymmetry is encoded in the difference between the coefficients $a_{+}$ and $a_{-}$.
We justify the use of the marginal Fermi liquid ansatz for Nd-LSCO at $p=0.24$ by the observation of perfectly $T$-linear resistivity (along both the $a$ axis and the $c$ axis) below $T \simeq 50$~K~\cite{daou2009,legros2019}.
The dimensionless coefficient $\alpha = 1.2 \pm 0.4$ extracted from ADMR~\cite{grissonnanche2020b} obeys the so-called Planckian limit~\cite{bruin2013,legros2019}.
In addition to this SMFL inelastic scattering rate, we also take into account, as above, the momentum-dependent
elastic scattering rate $1/\tau_{\rm 0}(\vecck)$ obtained from the ADMR experiment~\cite{grissonnanche2020b}
for Nd-LSCO at $p=0.24$. Hence, we use in our calculations the total scattering rate:
\begin{equation}
1/\tau(\epsilon, \vecck, T) = 1/\tau_{\rm 0}(\vecck) + 1/\tau_{\rm SMFL}(\epsilon, T).
\label{eq:total_scattering}
\end{equation}
Note that the ADMR study yields the remarkable finding that the inelastic scattering rate is isotropic, \textit{i.e.} $k$-independent. Correspondingly, we assume the energy dependence to also be isotropic, \textit{i.e.} we take $a_+$ and $a_-$ to be $k$-independent.
In Fig.~\ref{fig:calc_anisotropic_vs_T}(b) and (d), we see that the asymmetric scattering rate (plotted in Fig.~S3 at $T=25$~K) changes the calculated $S(T)$ dramatically for both directions and matches well the experimental data.
Crucially, the calculated $S_{\rm a}(T)$ now has the correct sign. Furthermore, $S_{\rm c}(T)$ is reproduced well, including the nontrivial plateau behavior between 15~K and 60~K.
The same model and the exact same parameters lead to calculated in-plane and out-of-plane resistivities that also reproduces the data well in agreement with Ref.~\cite{grissonnanche2020b}.
It is important to note that the particle-hole asymmetry of the scattering rate has a radical effect on the Seebeck coefficient, but no effect on the resistivity.

In Fig. \ref{fig:calc_isotropic_vs_T} we reveal the importance of incorporating the anisotropic elastic scattering rate $1/\tau_{\rm 0}(\vecck)$, inferred from the ADMR \cite{grissonnanche2020b}, in the calculations of the Seebeck coefficient.
We see that using an isotropic scattering rate does not affect $S_{\rm a}$ much, it greatly enhances the amplitude of $S_{\rm c}$ and does not exhibit the experimentally observed plateau which was reproduced with an anisotropic scattering rate.
It is not surprising that $S_{\rm c}$ is more affected by a large anisotropic scattering rate, as most of the $c$ axis dispersion is located closer to the antinodal region of the Fermi surface, the same region whose contribution gets dramatically reduced by the anisotropy of the scattering rate.
It therefore appears now that the plateau in $S_{\rm c}/T$ below 60~K (Fig.~\ref{fig:SoT_0p24_0p22})results from a competition between two effects that grow upon cooling: first the skewness increases upon cooling (see Fig.~S3) which increases $S_{\rm c}/T$, second the anisotropy of the scattering rate along the Fermi surface that decreases $S_{\rm c}/T$.

We now provide insight into the temperature-dependence and sign of the Seebeck coefficient
found from this calculation.
The SMFL expression (3) is a particular example of a general class of ‘skewed’ non-Fermi
liquids recently discussed by two of us in Ref.~\cite{georges2021}.
In that work, it is shown that an inelastic scattering rate obeying an $\epsilon/T$ scaling with an intrinsic particle-hole asymmetry of the scaling function remarkably leads to a modification of the low-temperature value of the Seebeck coefficient, even in the presence of elastic scattering.
This is in strong contrast to Fermi liquids, in which the dominant term in the inelastic scattering rate is particle-hole symmetric and particle-hole asymmetric subdominant corrections affect the Seebeck coefficient only at high temperature but not in the low-$T$ regime dominated by elastic scattering~\cite{georges2021,haule_2009}.
In the case of the SMFL, following a simplification that is discussed in the Appendix~\ref{app:kubo_vs_boltz} and detailed in Ref.~\cite{georges2021}, in the $T\rightarrow 0$ limit we obtain:
\begin{equation}
\frac{S}{T}\bigg|^{\mathrm{SMFL}}_{T=0} \simeq
\frac{1}{Z(T)}\left(\frac{S}{T}\right)_{{\rm Mott}} \, +  c_{\rm a}\, \frac{k_B^2}{e \hbar} \alpha\,\tau_{\rm 0}
\label{eq:Seebeck_slope_lowT_Planck}
\end{equation}
In this expression, the Mott value corresponds to expression (\ref{eq:Mott_formula}) which ignores the energy dependence of the scattering rate and whose sign is entirely determined by band structure.
$Z(T)$ is a mass enhancement renormalization to which we return below.
Here, we want to put the emphasis on the second term in (\ref{eq:Seebeck_slope_lowT_Planck}), which
involves a coefficient $c_{\rm a}$ determined  by the particle-hole asymmetry (`skew') of the scattering rate ($c_{\rm a}= \int_0^\infty dx x/[4 \cosh(x/2)^2] [\sqrt{1 + (a_+ \hbar x /\alpha)^2}-\sqrt{1 + (a_- \hbar x /\alpha)^2}]  \approx 1.1 \hbar  (a_+ - a_-) / \alpha \approx 0.4$ for the values of parameters reported in Appendix~\ref{app:seebeck_calculations}), the dimensionless coupling constant $\alpha$ measuring the strength of the inelastic scattering, and the elastic scattering time $\tau_{\rm 0}$.
This unconventional contribution to the Seebeck coefficient implies that the particle-hole asymmetry of the inelastic scattering rate can affect the magnitude and the sign of the Seebeck coefficient even in the low temperature limit where the magnitude of the inelastic scattering is small in comparison with the elastic scattering. This effect applies when the scattering rate is non-Fermi liquid and explains the observed sign change in the calculation above.
We also note that, in contrast to the band (Mott) term, this contribution does not depend on any details of the band structure and is, in particular, isotropic. This may provide a hint into the observed isotropy of our experimental data at low temperature.

\begin{figure}[t!]
\includegraphics[width=0.5\textwidth]{Fig7.pdf}
\caption{Left panels: Seebeck coefficient of Nd-LSCO $p~=~0.24$ for both $a$-axis (blue) and $c$-axis (red). Right panels: Corresponding calculated Seebeck coefficients for a momentum-independent but energy-dependent skewed scattering rate $1/\tau_{\rm 0} + 1/\tau_{\rm SMFL}(\epsilon, T)$, with $B=0$ in Eq.~S5 for the elastic scattering rate (isotropic, dashed line) and a momentum and energy-dependent skewed scattering rate of the form $1/\tau_{\rm 0}(\vecck) + 1/\tau_{\rm SMFL}(\epsilon, T)$ (anisotropic, full line), same as in Fig. \ref{fig:calc_anisotropic_vs_T}.
}
\label{fig:calc_isotropic_vs_T}
\end{figure}

We now discuss the renormalization $Z(T)$. By Kramers-Kronig, the ansatz for the scattering rate (imaginary part of the self energy) with a linear energy/temperature dependence  implies a logarithmic term in the low-energy slope of the real part of self energy, that is, a mass renormalization
$1/Z(T)=m^*/m \propto (1 + \log(\Lambda/T))$, where $\Lambda$ is a cutoff.
Indeed, a logarithmic enhancement of the specific heat coefficient has been observed experimentally for Nd-LSCO at $p=0.24$~\cite{michon2019}.
We see that this renormalization in turn implies a logarithmic behavior of $S/T$ at low temperatures, as seen in the first term of Eq.~\ref{eq:Seebeck_slope_lowT_Planck}.
It is important to notice however that this renormalization multiplies the band (Mott) part of the Seebeck coefficient.
As a consequence, when the logarithmic mass term is included in our modeling its effect is to increase $S_{\rm c}$ (improving the agreement with the experiment), but also to
diminish $S_{\rm a}$ (worsening the agreement somewhat for this direction). A further understanding of the mass renormalization effect is left for future work.
Note that logarithmic terms in the $T$-dependence of the Seebeck coefficient have been
discussed in \cite{paul_kotliar_2001} in relation to the proximity of a quantum critical point.

Recently, Jin \textit{et al.}~\cite{jin2021} reported measurements of in-plane Seebeck coefficient for LSCO $p=0.33$. At this doping, $S_{\rm a}$ displays also a positive sign down to the lowest temperature, while it shows a $T^2$ resistivity characteristic of a Fermi liquid. Because particle-hole asymmetry in a Fermi liquid cannot cause a change of sign in the low-temperature regime \cite{georges2021}, this could point to the influence of non-Fermi liquid corrections at the antinodes even in the overdoped regime~\cite{chang2013}.

We conclude this section by contrasting our analysis with previous calculations of the Seebeck coefficient of other families of cuprate compounds for $p>p^\star$.
In overdoped Bi2201 and Bi2212, in which the in-plane Seebeck coefficient is mostly negative, calculations by Kondo et al.~\cite{kondo2005} successfully reproduced the $T$-dependence of $S$ by using the electronic structure measured by ARPES and assuming that the scattering time depends on momentum in such a way as to maintain a constant scattering length $l= v \tau$, that is  $\tau(\vecck,\epsilon) \propto 1 / v(\vecck,\epsilon)$ (the same assumption was made in Ref.~\cite{storey2013}).
In this expression, $\vecck$ denotes a momentum on the Fermi surface and $\epsilon$ the energy of an excitation when moving away from the Fermi surface. The asymmetry between $\epsilon>0$ (particles) and $\epsilon<0$ (holes) in this model originates from the nearby van Hove singularity for momenta near the antinodes.
Although this model was able to capture the temperature dependence and sign of the data, its validity is uncertain.
Indeed, in this model the scattering rate is the largest at the nodes, where the velocity is maximum. This is at odds with most of the literature on cuprates \cite{abdeljawad2006, chang2013, grissonnanche2020b,yoshida_2007}.

Finally, it is interesting to note that, although quasi-1D organic superconductors can have very different physics from cuprates at the pseudogap critical point, calculations on a microscopic model have shown that incorporating the energy dependence of the scattering rate in those materials causes the Seebeck coefficient to change sign when close to the magnetic quantum critical point \cite{shahbazi2016}.

\subsection{$p < p^{\star}$: Fermi surface transformation}
In Nd-LSCO, reducing the doping below $p^{\star}$ induces two main anomalies in the Seebeck coefficient at low temperatures: a large enhancement of $S_{\rm a}$ and a negative $S_{\rm c}$.
Both are unambiguously connected with the onset of the pseudogap phase: they occur upon crossing $p^{\star}$, and at low temperatures below the pseudogap temperature $T^{\star}$.
The sudden sign change in the out-of-plane Seebeck coefficient necessarily reflects a change in particle-hole asymmetry upon entering the pseudogap phase, either through a transformation of the electronic structure, and therefore the Fermi surface, or the scattering rate.
As discussed previously~\cite{collignon2021thermopower}, the enhancement in $S_{\rm a}$ below $p^{\star}$ reflects the change in $n$ from $n \simeq 1 + p$ above $p^{\star}$ to $n \simeq p$ below $p^\star$ associated with the Fermi surface transformation at $p^{\star}$~\cite{fang2020}.
The origin of the negative $S_{\rm c}$, however, is not clear and is the focus of the remaining discussion.

In Nd-LSCO, the Fermi surface undergoes two changes at $p^{\star}$.
First, the large Fermi surface goes from hole-like to electron-like when the van Hove singularity crosses the Fermi level, which occurs at a doping $p_{\rm vHs}$ that coincides with $p^{\star}$~\cite{matt2015}.
As we show in Fig. \ref{fig:calc_vs_p}, the bare band calculations predict no sign change in $S_{\rm a}$ across $p_{\rm vHs}$, and a change to a positive sign in $S_{\rm c}$, in contrast to our results.
In Fig.~SS2, we plot $S_{\rm a}/T$ and $S_{\rm c}/T$ at 300~K as a function of doping and see that both quantities decrease smoothly from $p = 0.20$ to $p = 0.24$, with no anomaly.
Similarly, we observe that $S_{\rm a}/T$ and $S_{\rm c}/T$ show the same smooth increase upon cooling from 300~K down to about 100~K (Figs.~\ref{fig:SoT_doping}), irrespective of doping.
At those high temperatures above $T^{\star}$, only the vHs has an impact on the band structure, and since there is no clear difference in the behaviors of $S_{\rm a}/T$ and $S_{\rm c}/T$ we conclude that the vHs alone is not the main cause for the negative $S_{\rm c}/T$.
Second, there is a transformation of the Fermi surface associated with the drop in carrier density as deduced from Hall effect, resistivity~\cite{collignon2017}, thermal conductivity~\cite{michon2018}, and in-plane Seebeck~\cite{collignon2021thermopower} measurements.
This Fermi surface transformation was recently detected directly by ADMR measurements~\cite{fang2020}, which found that the Fermi surface below $p^{\star}$ is likely made of four nodal hole pockets, consistent with a transformation by a $Q = (\pi,\pi)$ wavector.
ARPES measurements see Fermi arcs that are consistent with the side of these pockets residing in the 1st Brillouin zone \cite{matt2015}. A model is needed to explain how such a Fermi surface transformation can account for a negative $S_{\rm c}$.

The scattering rate, on the other hand, is not known to undergo a significant change at $p^{\star}$.
Resistivity measurements on either side of $p^{\star}$~\cite{collignon2017} find that the normal-state magnetoresistance (MR) is comparable at $p = 0.22$ and $p = 0.24$.
In the weak-field limit, the magnetoresistance varies as MR~$\propto(\omega_c \tau)^2$, where $\omega_c$ is the cyclotron frequency and $\tau$ is the scattering time, suggesting that $\tau$ does not change significantly across $p^{\star}$.
Similarly, in ADMR measurements on Nd-LSCO~\cite{fang2020} the magnitude of the oscillations seen on each side of $p^{\star}$, at $p = 0.21$ and $p = 0.24$, is similar, implying that the amplitude of $\tau$ is roughly the same.
So there is no indication, a priori, that the sign change in $S_{\rm c}$ is correlated with a change in $\tau$ at the Fermi level.
Finally, the fact that the change of $S_{\rm a}$ through $p^{\star}$ roughly matches the ratio of the carrier densities further indicates that any change in $\tau$ is small.

A change in Fermi surface, and specifically in the asymmetry of the band near the Fermi level is therefore the likely cause of the negative $S_{\rm c}$.
It seems clear that without the pseudogap, which appears at temperatures below $\sim 100$~K over this doping range, $S_{\rm a}/T$ and $S_{\rm c}/T$ would keep evolving in tandem down to low temperatures, likely to match in the $T \rightarrow 0$ limit as they do at $p = 0.24$. (This was actually shown to be the case by application of pressure to move $p^\star$ down below $p = 0.22$: the low-$T$ rise in $S_{\rm a}/T$ at $p = 0.22$ seen at ambiant pressure is gone under pressure~\cite{gourgout2021pressure}.)
Consequently, we attribute the negative $S_{\rm c}$ to a transformation of the band structure and therefore of the Fermi surface caused by the pseudogap phase.
Based on ARPES~\cite{matt2015} and recent ADMR measurements~\cite{fang2020}, the Fermi surface below $p^{\star}$ is truncated at the antinodes. Now this antinodal region is where $c$ axis dispersion is largest, thereby dominating, lending the bulk of the $c$ axis transport. From there, it is not hard to imagine that out-of-plane Seebeck would be most sensitive to a transformation of the Fermi surface in the ($\pi$, 0) direction caused by the pseudogap phase. Further work  is needed to model the sign change of $S_{\rm c}$ below $p^*$.

\section{Summary}

In summary, we have measured the Seebeck coefficient of the cuprate superconductor Nd-LSCO at dopings close to the critical doping $p^{\star}$\xspace~$= 0.23$ where the pseudogap phase ends at $T=0$.
In particular, we examine the $c$ axis component, $S_{\rm c}$, for a heat current normal to the CuO$_2$ planes, in the presence of a magnetic field sufficient to suppress superconductivity and thus track the normal-state behaviour down to low temperature.

At $p >$~$p^{\star}$\xspace, for $p = 0.24$, we find that $S_{\rm c}(T)$ and the in-plane component $S_{\rm a}(T)$ are both positive below 250~K, and become roughly equal below 10~K. This is in contrast to what calculations based on the well-characterized band structure predict, namely that $S_{\rm a}$ should be negative at all temperatures. We show that a good quantitative description of (and the correct sign for) both $S_{\rm c}(T)$ and $S_{\rm a}(T)$ is obtained if we add a linear asymmetric energy dependence to the scattering rate $1/\tau$ previously extracted in Nd-LSCO as a function of angle and temperature via ADMR measurements~\cite{grissonnanche2020b}.
This suggests that $1/\tau$ is Planckian not only in its $T$ dependence ($1/\tau \simeq k_{\rm B}T / \hbar$) but also in its energy dependence ($1/\tau \propto \epsilon$),
with an intrinsic particle-hole asymmetry of the $\epsilon/T$ scaling function -- as expected for a skewed Planckian metal~\cite{georges2021}.

At $p <$~$p^{\star}$\xspace, we find that $S_{\rm c}(T)$ undergoes a dramatic change at low temperature, suddenly becoming negative as soon as $p$ falls below $p^{\star}$\xspace. This is a striking new experimental signature of the pseudogap phase, which we attribute to a transformation of the Fermi surface. Further work is needed to identify what modification of the Fermi surface is responsible for this sign change.

\section*{Acknowledgements}

We thank S.~Fortier for his assistance with the experiments.
L.T. acknowledges support from the Canadian Institute for Advanced Research
(CIFAR) as a CIFAR Fellow
and funding from
the Institut Quantique,
the Natural Sciences and Engineering Research Council of Canada (PIN:123817),
the Fonds de Recherche du Qu\'ebec -- Nature et Technologies (FRQNT),
the Canada Foundation for Innovation (CFI),
and a Canada Research Chair.
J.M. is supported by the Slovenian Research Agency (ARRS) under Program No. P1-0044 and Project No. J1-2458, N1-0088, and J1-2455-1.
J.S.Z. was supported by NSF MRSEC under Cooperative Agreement No. DMR-1720595.
This research was undertaken thanks in part to funding from the Canada First Research Excellence Fund
and the Gordon and Betty Moore Foundation's EPiQS Initiative (Grant GBMF5306 to L.T.).
The Flatiron Institute is a division of the Simons Foundation.


\section{Introduction}

Whether classical and quantum mechanical, the transport of statistical density is our primary means of making predictions of macroscopic, nonequilibrium behavior from microscopic dynamics~\cite{zwanzig2001nonequilibrium}.
Classically, Jacobi's form of Liouville's equation of motion for the phase space density of mechanical systems is the foundation of classical statistical mechanics~\cite{tolman1979principles}.
Its usefulness derives largely from its many forms and approximations, including the Boltzmann equation, the Vlasov approximation, the Bogoliubov-Born-Green-Kirkwood-Yvon hierarchy, that underlie applications across physics and chemistry~\cite{zwanzig2001nonequilibrium}.
In quantum mechanics, the Liouville-von Neumann equation describes the evolution of the density operator~\cite{Neu1927}; it is the fundamental equation of quantum statistical mechanics and a main ingredient in quantum computing, tomography, and decoherence~\cite{Blu2012}.
To translate between the classical and quantum mechanical Liouville equations, one can use Dirac's rule~\cite{dirac1981principles} of replacing Poisson brackets by commutators.
Here, we establish a density matrix formalism for classical systems, with not only a density matrix but also a classical commutator, supplanting Dirac's heuristic with a more direct correspondence between these physical theories.

Operator-theoretic methods, such as Frobenius-Perron and its dual Koopman formalism~\cite{Koopman315}, give a formal analogy to quantum mechanics by lifting the description of classical systems to infinite dimensions~\cite{neumann_zur_1932,*neumann_zusatze_1932}.
They preserve global nonlinear features and guarantee exact linearization of the dynamics, providing useful connections between classical dynamical systems and statistical physics~\cite{dorfman1999introduction,gaspard2005chaos}.
However, they can be difficult to apply to systems under active external control and to find observables representing the nonlinear system in the lifted linear space~\cite{lusch_deep_2018}.
Symmetries can make the calculation of the Koopman operator approximation and its spectral properties more efficient~\cite{salova_koopman_2019}, but, in practice~\cite{rowley2009,budisic_applied_2012}, the number of variables must be truncated to finite-dimensions (e.g., through extended~\cite{williams_datadriven_2015,korda_convergence_2018} or kernel~\cite{williamsrowly2015} dynamic mode decomposition~\cite{schmid_2010}).
While reminiscent of quantum mechanics, these are classical theories that add weight to the question of whether other formulations of quantum mechanics might have classical counterparts that are fundamental to statistical physics~\cite{BudRoh2017}.

While Liouville's equation is the foundation of nonequilibrium statistical physics, many theories avoid, approximate, or subject Liouville's equation to model specific solutions~\cite{zwanzig2001nonequilibrium}.
Here, we construct a classical density matrix formulation of dynamical systems on the local stability of nonlinear dynamics~\cite{Eckmann1985Ergodic} -- Lyapunov exponents and vectors~\cite{PikovskyP16}.
The infinitesimal perturbations, Lyapunov vectors, we use to define the density matrix have been used to analyze rare trajectories~\cite{TailleurK07}, jamming~\cite{banigan_chaotic_2013}, nonequilibrium self-assembly~\cite{GreenCGS13}, equilibrium and nonequilibrium fluids~\cite{evans1990statistical, Bosetti_2014,das_self-averaging_2017}, and critical phenomena~\cite{dascritical2019}.
From these finite-dimensional vectors, we derive a classical analogue of the von Neumann equation for the density matrix dynamics.
We show this classical density matrix is similar to the dual metric tensor and that its determinant evolves according to a generalized Liouville equation and satisfies a generalized Liouville theorem.
And, by imposing a norm-preserving dynamics with Lyapunov exponents not only normalizes the density matrix, it reinstates the form of the usual Liouville equation for generic, non-Hamiltonian dynamical systems.
All together, these results establish a computationally-tractable basis for nonequilibrium statistical mechanics grounded in dynamical systems theory.

The structure and this formalism suggests the possibility of other classical counterparts to elements of quantum theory, including uncertainty relations and speed limits; In Ref.~\cite{das2021speed}, we use the classical density matrix theory here to derive a family of time-information uncertainty relations that set speed limits on the evolution of (observables of) arbitrary dynamical systems.

\begin{figure*}[!htb]
\vspace{-2cm}
\includegraphics[width=0.49\linewidth, height=0.3\textheight]{Lorenz_with_LVs_on_panel1.pdf}\hspace{-0.5cm}
\includegraphics[width=0.49\linewidth, height=0.3\textheight]{Lorenz_with_LVs_on_panel2.pdf}
\vspace*{-2cm}
\caption{\label{fig:Lorenz_attractors}Snapshots of the strange attractor and an initially random, unit perturbation transported along a chaotic solution of the Lorenz system. Parameters are those originally used by Lorenz~\cite{Lorenz63}: $\sigma = 10$, $\beta = 8/3$, $\rho = 28$.}

\end{figure*}

To start, we define an unnormalized density matrix from the linearization of the classical dynamics, Sec.~\ref{sec:dynamics_unnorm}.
The properties of this density matrix lead to a generalization of Liouville's theorem and equation, Sec.~\ref{sec:GLE/GLT}.
For Hamiltonian dynamics, Sec.~\ref{sec:PB}, we show the reduction to the usual Liouville theorem and equation and establish a connection to Poisson brackets.
The dynamics of a normalized density matrix, Sec.~\ref{sec:dynamics_norm}, transform the generalized Liouville's theorem and equation to the usual form, Sec.~\ref{sec:GLE_norm}.
The basis chosen determines the representation of the density matrix, Sec.~\ref{sec:basis}, which is a common consideration in dynamical systems theory.
This formalism is illustrated analytically for the linear and damped harmonic oscillator in Sec.~\ref{sec:case_studies}.

\section{Dynamics of the classical unnormalized density matrix}\label{sec:dynamics_unnorm}
\noindent Consider a generic dynamical system with state-space variables $\{x^i\}$.
At any moment in time, these variables together mark a point ${\bm{x}}(t) := [x^1(t),x^2(t),\ldots,x^n(t)]^\top$ in an $n$-dimensional state space $\mathcal{M}$ that evolves according to: $\dot{{\bm{x}}} = \boldsymbol{F}[{\bm{x}}(t)]$.
Perturbations to the system will also evolve under the flow of the dynamics.
Because of their analytical and computational tractability, infinitesimal perturbations $\ket{\delta{\bm{x}}(t)}:= [\delta x^1(t), \delta x^2(t), \ldots, \delta x^n(t)]^\top\in T\mathcal{M}$ and their linearized dynamics are a well established means of analyzing the stability of nonlinear dynamical systems~\cite{PikovskyP16}.
These perturbations represent an uncertainty about the initial condition and stretch, contract, and rotate over time,
\begin{equation}
  \ket{\delta\dot{{\bm{x}}}(t)} = {\bm{A}}[{\bm{x}} (t)]\ket{\delta {\bm{x}}(t)},
  \label{equ:EOM_per_ket2}
\end{equation}
with the phase point. Their time evolution is governed by the local stability matrix ${\bm{A}}:={\bm{A}}[{\bm{x}}(t)] = \grad\boldsymbol{F}$ with elements
$({\bm{A}})^i_{j}=\partial \dot{x}^i(t)/\partial x^j(t)$.
Figure~\ref{fig:Lorenz_attractors} shows a unit perturbation vector as it is transported across the Lorenz attractor.
Surrounding the phase point ${\bm{x}}(t)$ is an infinitesimal $k$-dimensional phase space volume that transforms its shape over time as it evolves.
The volume we take to be spanned by a finite set $\{\ket{\delta{\bm{\psi}}_i}\}$ of $k\leq n$ basis vectors $\ket{\delta{\bm{\psi}}_i} \in T\mathcal{M}$ ($i=1,2,\ldots, k$) that also obey the linearized dynamics (e.g., Lyapunov vectors~\cite{f_ginelli_characterizing_2007,Wolfe2007}).
If the dynamics are Hamiltonian, then according to Liouville's theorem, the volume spanned by these tangent states is conserved at all times.

Comparing to quantum mechanics, the evolution of classical, finite-dimensional, real tangent space vectors according to Eq.~\ref{equ:EOM_per_ket2} is analogous to the evolution of infinite-dimensional, complex Hilbert space vectors according to Schr\"odinger's equation~\cite{dirac1981principles}.
An alternative representation of quantum states, however, is the density operator formalism, which provides a basis for quantum technology, dissipative decoherence, and statistical mechanics~\cite{Fan1957,Haa1961,Blu2012}.
Because of the power of this formulation of quantum mechanics, we break from traditional classical dynamical systems, constructing this formulation by defining the classical state with a density operator:
\begin{equation}
  {\bm{\xi}}(t) := \sum_{i=1}^k\dyad{\delta {\bm{\psi}}_i(t)}.
\end{equation}
This \textit{unnormalized} matrix form of this operator represents an alternative state of a classical dynamical system.
Expressed here in the $\{\delta{\bm{\psi}}_i\}$ basis, it is the outer product of perturbation vectors or what Gibbs called the dyadic~\cite{gibbs1901vector} product.

Partitioning the stability matrix ${\bm{A}} = {\bm{A}}_+ + {\bm{A}}_-$ into its symmetric and anti-symmetric parts, ${\bm{A}}_\pm = \frac{1}{2}({\bm{A}} \pm {\bm{A}}^\top)$, the time evolution of ${\bm{\xi}}$ is:
\begin{equation}\label{eq:EOM_xi}
  \frac{d{\bm{\xi}}}{dt} = \{{\bm{A}}_+,{\bm{\xi}}\}+[{\bm{A}}_-,{\bm{\xi}}].
\end{equation}
Its solution,
\begin{align}\label{eq:op_M}
  {\bm{\xi}}(t) = {\bm{M}}(t,t_0){\bm{\xi}}(t_0){\bm{M}}^\top(t,t_0),
\end{align}
is in terms of the non-normal propagator, $({\bm{M}})^i_{j}=\partial x^i(t)/\partial x^j(t_0)$ (Appendix~\ref{sec:EOM_LV}) when ${\bm{\xi}}$ is built from the tangent vectors $\{\ket{\delta {\bm{x}}^i}\}$ evolved by ${\bm{M}}$.

The generally non-symmetric stability matrix ${\bm{A}}$ plays the role of the quantum mechanical Hamiltonian in the \textit{classical} commutator $[\boldsymbol{X},\boldsymbol{Y}]=\boldsymbol{Y}\boldsymbol{X}-\boldsymbol{X}\boldsymbol{Y}$ and anti-commutator $\{\boldsymbol{X},\boldsymbol{Y}\}=\boldsymbol{Y}\boldsymbol{X}+\boldsymbol{X}\boldsymbol{Y}$.
Two facts are striking about this solution: it is entirely computable from standard methods in dynamical systems theory~\cite{PikovskyP16} and it is a purely classical analogue of the Liouville-von Neumann equation in quantum dynamics.
Unlike its quantum counterpart, the determinant of ${\bm{\xi}}$ is directly related to Liouville's equation and theorem with this equation of motion, as we show next.

\subsection{Generalized Liouville's theorem/equation}\label{sec:GLE/GLT}

\noindent Liouville's theorem and equation are the foundation for nonequilibrium statistical mechanics~\cite{zwanzig2001nonequilibrium,dorfman1999introduction,gaspard2005chaos} and the point at which statistical mechanics departs from classical dynamics.
Here, we derive generalizations of both from the determinant of the unnormalized density matrix. To make this connection, consider the arbitrary set of basis vectors $\{\delta {\bm{\psi}}_i\}$ that span the entire $n$-dimensional phase space volume, $d\mathcal{V}$, with the perturbation state ${\bm{\xi}}$.
The (square of the) phase space volume is determined by the determinant $|{\bm{\xi}}|$,
which has the equation of motion (Appendix~\ref{sec:GLE}):
\begin{equation}\label{eq:gen_Liouville}
\frac{1}{2}\frac{d}{dt}\ln|{\bm{\xi}}(t)| = \grad\cdot \dot{\bm{x}} = \Tr{\bm{A}}_+ = \Lambda.
\end{equation}
Both this equation of motion and its solution,
\begin{equation}\label{eq:liouville_unnorm_den_mat}
  |{\bm{\xi}}(t)| = |{\bm{\xi}}(t_0)|\,e^{2\int_{t_0}^t\Lambda(t') \,dt'},
\end{equation}
depend on the divergence of the phase space velocity $\dot{{\bm{x}}}$ or the phase space volume contraction/expansion rate $\Lambda = \Tr{\bm{A}}_+ = \Tr{\bm{A}}$.
This rate $\Lambda$ is (twice) the sum of the Lyapunov exponents~\cite{dorfman1999introduction}, which is related to physical quantities.
For example, the phase space contraction rate is related to the entropy production rate in fluid transport~\cite{dorfman1999introduction}.

From these results for the density matrix determinant, we can generalize Liouville's theorem and equation.
The geometric interpretation of Liouville's theorem is that the velocity field $\dot{{\bm{x}}}$ has zero divergence: $\grad\cdot\dot{{\bm{x}}} = \Tr({\bm{\mathcal{H}}}) = 0$,
and, consequently, the ``phase fluid'' flow is incompressible and phase space volumes are conserved $d\mathcal{V}(t)=d\mathcal{V}(t_0)$.
Here, ${\bm{\mathcal{H}}}$ is the stability matrix for Hamiltonian dynamics with Hamiltonian $H$. The phase space volume element $d\mathcal{V}(t)$ spanned by the basis vectors has a coordinate transformation: $d\mathcal{V}(t)=|{\bm{M}}(t,t_0)|\,d\mathcal{V}(t_0)$. Combining this fact with the determinant of Eq.~\ref{eq:op_M}, we obtain a generalization of Liouville's theorem for general dynamical systems:
\begin{align}\label{eq:invariant_measure}
  |{\bm{\xi}}(t)|^{-\frac{1}{2}} d\mathcal{V}(t) &= |{\bm{\xi}}(t_0)|^{-\frac{1}{2}} d\mathcal{V}(t_0).
\end{align}
Any dynamics conserves the measure, $|{\bm{\xi}}|^{-\frac{1}{2}} d\mathcal{V}$ or $e^{-\int_{t_0}^t\Lambda(t')\,dt'}d\mathcal{V}$.
For dissipative systems with $\Lambda < 0$, volumes contract at a rate $\Tr{\bm{A}}$.
In the Lorenz model, for example, ${\bm{A}}$ is constant, so $|{\bm{\xi}}|$ decays linearly on a semi-log scale with a slope proportional to $2\Tr{\bm{A}}$ as shown in Fig.~\ref{fig:Lorenz_trace_det}.
But, when the dynamics are Hamiltonian, we recover the conventional form of Liouville's theorem~\cite{tolman1979principles} because $\Lambda = 0$ and $|{\bm{\xi}}|$ is a constant of motion.

Equation~\ref{eq:invariant_measure} also appears in Riemannian geometry as the transformation of a metric determinant on a Riemannian manifold of an arbitrary curvature endowed with a covariant metric tensor $\boldsymbol{g}_{\bm{\xi}}$.
Metrics have been considered previously in nonlinear dynamics~\cite{greeneIntroductionMetricTensor1989,thiffeaultCovariantTimeDerivatives2001}.
But, as we shown in Appendix~\ref{sec:GLE}, the (covariant) metric tensor is similar to the inverse of the unnormalized density matrix, ${\bm{\xi}}^{-1}$.
With this identification, a number of results follow.
Most immediate is that Eq.~\ref{eq:invariant_measure} becomes the transformation of the metric determinant: $\sqrt{g_{\bm{\xi}}(t)}d\mathcal{V}(t)=\sqrt{g_{\bm{\xi}}(t_0)}d\mathcal{V}(t_0)$ with $g_{\bm{\xi}}=|{\bm{\xi}}^{-1}|$.
However, what traditionally follows from Liouville's theorem is Liouville's equation, a formally exact equation of motion for the probability density in phase space.
This equation derives from another statement of Liouville's theorem: the density of representative points in the phase space is conserved along the trajectories of Hamiltonian systems, $d_t\rho({\bm{x}})=0$~\cite{tolman1979principles}.
The density matrix also yields a generalization of this form of Liouville's equation (Appendix~\ref{sec:GLE}).

With the similarity of the unnormalized density matrix and the metric tensor, the flow compressibility accounts for the metric's compatibility with the dynamics~\cite{Ezra2004},
$-d_t\ln\sqrt{g_{\bm{\xi}}} = (n/2)d_t\ln\Tr\boldsymbol{g}^{-1}_{\bm{\xi}} = \grad\cdot\dot{{\bm{x}}}$.
By identifying the density matrix as similar to the dual metric, $\boldsymbol{g}_{\bm{\xi}}^{-1}$, we can also find this compatibility condition as the equation of motion for the metric determinant, Eq.~\ref{eq:gen_Liouville}, and $\Tr{\bm{A}}_+$.
Therefore, $|{\bm{\xi}}|^{-\frac{1}{2}}$ obeys the generalized Liouville's equation:
\begin{equation}\label{eq:GLE_xi}
  \frac{\partial}{\partial t}|{\bm{\xi}}|^{-\frac{1}{2}} + \grad\cdot(|{\bm{\xi}}|^{-\frac{1}{2}}\dot{{\bm{x}}}) = 0,
\end{equation}
which is related to other generalizations~\cite{Tuckerman_1999,TucLiuGioMar2001,Ramshaw_2002,Ezra2004}.
It applies to non-Hamiltonian dynamics but gives the usual Liouville equation if the dynamics are Hamiltonian.
For Hamiltonian dynamics, the metric determinant $|{\bm{\xi}}|^{-1}$ is time-independent, the divergence of the flow vanishes, $d\mathcal{V}(t)=d\mathcal{V}(t_0)$ from Eq.~\ref{eq:invariant_measure}, and the generalized Liouville equation,
\begin{equation}\label{eq:GLE_xi_Hamiltonian}
\frac{\partial}{\partial t}|{\bm{\xi}}|^{-\frac{1}{2}} + \dot{{\bm{x}}}\cdot\grad|{\bm{\xi}}|^{-\frac{1}{2}} = 0,
\end{equation}
becomes the usual Liouville equation.

\begin{figure}[t!]
	\includegraphics[width=1\columnwidth]{Trace_determinant_Lorenz.pdf}
	\caption{\label{fig:Lorenz_trace_det} The trace of the unnormalized density matrix, ${\bm{\xi}}(t)$, as a function of time for 100 random perturbations drawn from a uniform distribution. The inset shows the evolution of pure states $|{\bm{\xi}}(t)|$ for a set of 100 local perturbations. Data for the vector evolving on the Lorenz attractor in Fig.~\ref{fig:Lorenz_attractors} are shown in blue. The trace and determinant of the normalized density matrix ${\bm{\varrho}}$ (inset) are time-invariant (dashed).}
	\vspace{-.2in}
\end{figure}

\subsection{Poisson brackets}\label{sec:PB}
\noindent So far, we have analyzed the properties of the determinant of the unnormalized density matrix. However, its trace also appears in the compatibility condition. Analyzing the trace, we find connections to classical dynamics, well-known quantities in dynamical systems, and similarities with quantum mechanics.

Because the dynamics of ${\bm{\xi}} = \sum_{i=1}^{k}\dyad{\delta{\bm{\psi}}_i}{\delta{\bm{\psi}}_i}$ are not trace preserving, the rate of change of
$\Tr{\bm{\xi}}$ governs the rate of change of the perturbation,
\begin{equation}
  \frac{1}{2}\frac{d}{dt}\Tr{\bm{\xi}} = \frac{1}{2}\Tr\{{\bm{A}}_+,{\bm{\xi}}\}
= \Tr({\bm{\xi}} {\bm{A}}_+) = \langle {\bm{A}}_+\rangle_{{\bm{\xi}}}.
\end{equation}
Averages over the unnormalized density matrix are defined akin to quantum mechanical averages; here, the trace evolves at a rate determined by the instantaneous Lyapunov exponents~\cite{PikovskyP16}, measures of local (in)stability, $\langle {\bm{A}}_+\rangle_{{\bm{\xi}}_i}$. Defining $\langle {\bm{A}}_+\rangle = \langle {\bm{A}}_+\rangle_{{\bm{\xi}}}/\Tr{\bm{\xi}} $, the solution to this equation of motion is:
\begin{equation}\label{eq:trace_xi_sol}
  \Tr{\bm{\xi}}(t) = \Tr{\bm{\xi}}(t_0)e^{2\int^t_{t_0}\langle {\bm{A}}_+(t')\rangle\,dt'}.
\end{equation}
As numerical verification of this result, and all our others, we simulated the dynamics of Hamiltonian and dissipative dynamical systems. As a prototypical dissipative system, we chose the Lorenz model. Figure~\ref{fig:Lorenz_trace_det} shows the time evolution of the trace for a chaotic orbit for 100 random perturbation states drawn from a uniform distribution.
The rapid increase (on the semi-log scale) in $\Tr{\bm{\xi}}$ is indicative of the chaotic nature of the chosen orbit.

For Hamiltonian systems, the equation of motion for the trace can be expressed as a Poisson bracket.
In classical statistical mechanics, a dynamical variable, $f(\boldsymbol{q},\boldsymbol{p})$, can be expressed as $\dot{f}=\{f, H\}_P$, in terms of the Poisson bracket, $\{.\}_P$.
Combined with our result above, this fact gives a correspondence,
\begin{equation}\label{eq:poisson_connection}
\Tr\dot{{\bm{\xi}}} = 2\langle{\bm{\mathcal{H}}}_+\rangle_{\bm{\xi}} = \{\Tr{\bm{\xi}}, H\}_P,
\end{equation}
between the Poisson bracket and the tangent-space average of the symmetric part of the stability matrix ${\bm{\mathcal{H}}}_+$.
For example, for an arbitrary perturbation $(u,v)^\top$ in the phase space of the linear harmonic oscillator, $\Tr\,{\bm{\xi}} =u^2+v^2$ and $\Tr\dot{{\bm{\xi}}}=\{\Tr{\bm{\xi}}, H\}_P = uv(1-\omega^2)$ where $\omega$ is the oscillation frequency.

\begin{figure*}[t]
	\centering
		\resizebox{1.8\columnwidth}{!}
		{
	\begin{minipage}{0.3\textwidth}
		\hspace{-2cm}
		\includegraphics[]{HH_trace_equ_pt_lines.pdf}
	\end{minipage}
	\begin{minipage}{.3\textwidth}
		\vspace*{0.7cm}\hspace*{2cm}\small(b)

		\hspace{-1.5cm}
		\vspace{0.5cm}
		\resizebox{0.92\columnwidth}{!}
		{
			\begin{tikzpicture}
			\definecolor{S_DBlue}{RGB}{3,78,123}
			\definecolor{blue2}{RGB}{32,99,155}
			\definecolor{S_DGrey}{RGB}{166,166,166}
			\definecolor{S_Gold}{RGB}{255,192,0}
			\begin{axis}
			[
			view={60}{50},
			axis line style={draw=none},
			tick style={draw=none},
			tick label style={color=white}
			]
			\addplot3[
			surf,
			colormap/blackwhite,
			shader=interp,
			domain=-1.2:1.2,
			domain y=-1.13:1.13,
			]
			{exp(-.15*x^2-.15*y^2)};
			\addplot3[fill=blue2, opacity=0.5]
			coordinates{
				(-1,-1,1)
				(-1,1,1)
				(1,1,1)
				(1,-1,1)
			};
			\end{axis}
			\draw[black,very thick,-latex] (3.6,3.52) -- ++(94:2.5) node[below left] {\large $\grad H$};
			\draw[S_Gold,very thick,-latex] (3.6,3.52) to[bend left] (6.05,2.77);
			\draw[black,very thick,-latex] (3.6,3.52) -- ++(14:2.1);

			\node at (5.1,4.19) {\large\color{black}$\dot{{\bm{x}}}$};
			\node at (3.52,3.18) {\large$P$};
			\node at (4.8,2.2) {\large$\mathcal{M}$};
			\node at (1.25,4.0) {\large$T\mathcal{M}$};
			\node at (4.8,3.2) {\large\color{S_Gold}$\Gamma(t)$};
			\end{tikzpicture}
		}
	\end{minipage}%
	\begin{minipage}{0.35\textwidth}
		\hspace{-0.8cm}\includegraphics[]{HH_ILEs_example.pdf}
	\end{minipage}
}
\caption{\label{fig:manifold}(a) The trace of the unnormalized density matrices for a regular (black) and a chaotic orbit (blue) for the H\'enon-Heiles system.
The trace of the normalized density matrices is time invariant (gray).
The inset shows the equipotential surface (and lines) for the potential $U = \frac{1}{2}(x^2+y^2+2x^2y-\frac{2}{3}y^3)$.
Line colors in (a) correspond to those of the equipotential lines.
(b) Schematic illustration of the neighborhood of a point $P$ along a trajectory $\Gamma(t)$ through the phase space $\mathcal{M}$ with the associated tangent space $\mathcal{T}_\mathcal{M}$ and conjugate vectors, $\grad H$ and $\dot{\boldsymbol{x}}$.
(c) Pair of instantaneous Lyapunov exponents (ILE) for the regular with $E=0.0833$ (black) and a chaotic orbit with $E=0.1667$ (blue).}
\end{figure*}

\section{Dynamics of the classical normalized density matrix}\label{sec:dynamics_norm}
\noindent The connections between the density matrix and its properties so far have immediate implications for classical statistical physics. The key results differ in some respects from their analogues in quantum mechanics. For example, the dynamics of ${\bm{\xi}}$ are not norm-preserving, while the norm-preserving dynamics of Hilbert state vectors are a postulate of many formulations of quantum mechanics.

To derive a norm-preserving dynamics for the density matrix of general, classical dynamical systems, we consider a unit perturbation $\ket{\delta {\bm{u}}} = \ket{\delta {\bm{x}}}/\|\delta {\bm{x}}\|$, where $\|.\|$ is the $\ell_2$-norm.
Its equation of motion,
\begin{eqnarray}\label{equ:EOM_unit_per}
  \frac{d}{dt}\ket{\delta {\bm{u}}}
&=& ({\bm{A}}+{\bm{A}}_-)\,\ket{\delta {\bm{u}}}- r\,\ket{\delta {\bm{u}}},
\end{eqnarray}
contains a source/sink term with the instantaneous rate:
$r:=r(t) = \bra{\delta {\bm{u}}}{\bm{A}}_+\ket{\delta {\bm{u}}}= d_t\ln\|\delta {\bm{x}}(t)\|$.
This rate is the instantaneous Lyapunov exponent for a linearized dynamics, related to the finite-time Lyapunov exponent, $\lambda(t) := \lambda(t,t_0) = |t-t_0|^{-1}\int_{t_0}^tr(t)\,dt$,
and the Lyapunov exponent, $\lambda_i=\lim_{t\to\infty}\lambda_i(t)$. The maximum instantaneous Lyapunov exponent is also referred to as \textit{reactivity} -- the maximum amplification rate over all perturbations, immediately  following the perturbation~\cite{NeuCas1997}.

As before, we represent the state of the dynamical system $\dot{{\bm{x}}}=\boldsymbol{F}$ as a density matrix. But now, we express it in terms of a unit tangent-space basis $\{\ket{\delta{\bm{\phi}}_i\}}$.
Normalizing each $\ket{\delta{\bm{\psi}}_i} =c_i\ket{\delta{\bm{\phi}}_i}$, we can define the pure states ${\bm{\varrho}}_i(t) := \dyad{\delta {\bm{\phi}}_i(t)}{\delta {\bm{\phi}}_i(t)}$ with the expected properties: ${\bm{\varrho}}_i^2 = {\bm{\varrho}}_i$, $\Tr{\bm{\varrho}}_i = 1$, $\Tr{\bm{\varrho}}_i^2=1$, symmetric, ${\bm{\varrho}}_i \succeq 1$, i.e., ${\bm{\varrho}}_i$ is positive semi-definite.
Proving these properties requires the dynamics of ${\bm{\varrho}}_i$ be norm-preserving.
The perturbation state is the \textit{normalized} density matrix:
\begin{equation}\label{equ:norm_unnorm_den_mat_rel}
  {\bm{\varrho}}(t) = \frac{{\bm{\xi}}(t)}{\Tr{\bm{\xi}}(t)} = K^{-1}\sum_{i=1}^kc_i^2\dyad{\delta{\bm{\phi}}_i}{\delta{\bm{\phi}}_i},
\end{equation}
which is directly related to the unnormalized ${\bm{\xi}}$ state with $K=\sum_i c_i^2$.
If the state is maximally mixed, with $c_i=1$ $\forall i$ then $\Tr{\bm{\xi}}=k$.
The tangent-space average $\expval{{\bm{\varrho}}}{\delta{\bm{\phi}}}$ over the unit basis is analogous to the quantum-mechanical probability of finding the system at $\ket{\delta{\bm{\phi}}}$ given that its state is ${\bm{\varrho}}$; for example, if  all tangent space directions contribute equally $\expval{{\bm{\varrho}}}{\delta{\bm{\phi}}}=k^{-1}$.

The normalized state evolves in time according to:
\begin{equation}\label{equ:EOM_rho}
  \frac{d{\bm{\varrho}}}{dt} = \{{\bm{A}}_+,{\bm{\varrho}}\} + [{\bm{A}}_-,{\bm{\varrho}}] - 2r{\bm{\varrho}},
\end{equation}
another classical analogue of the von Neumann equation in quantum mechanics.
The equation indicates that the average of ${\bm{A}}_+$ at each instant of time, i.e., the instantaneous Lyapunov exponent is crucial to norm preservation.
It offsets the stretching and contraction of $\ket{\delta {\bm{u}}}$ due to ${\bm{A}}$ along a given trajectory.
Solving this equation of motion, we find the density matrix ${\bm{\varrho}}$ at time $t_0$ and at a later time $t$ are similar:
\begin{align}\label{eq:rho_sim_trans}
  {\bm{\varrho}}(t)   = \tilde{\bm{M}}(t,t_0){\bm{\varrho}}(t_0)\tilde {\bm{M}}^{-1}(t,t_0),
\end{align}
when ${\bm{\varrho}}$ is composed of vectors $\{\ket{\delta {\bm{u}}^i}\}$ evolving in time under $\tilde {\bm{M}}$.
Regardless of the dynamics, the norm-preserving evolution operator $\tilde{{\bm{M}}} := {\bm{M}}\boldsymbol{\Gamma}$ is orthogonal (Appendices~\ref{sec:EOM_LV} and~\ref{sec:proof}) and, so, $|\tilde{\bm{M}}|=1$; the matrix $\boldsymbol{\Gamma}$ has the inverse expansion factors $\|\delta {\bm{x}}^i(t_0)\|/\|\delta {\bm{x}}^i(t)\|=e^{-\lambda_i(t) (t-t_0)}$ on the diagonal.
Expressing the unnormalized density matrix ${\bm{\xi}}'=\sum_{i=1}^k\dyad{\delta{\bm{x}}_i(t)}$ (in terms of the set of vectors that evolve under ${\bm{M}}$), all $\ket{\delta{\bm{x}}_i(t)}$ collapse on the most rapidly expanding direction.
This collapse motivates well-known methods for computing tangent basis sets~\cite{Benettin1976Kolmogorov,greeneCalculationLyapunovSpectra1987,Wolfe2007,f_ginelli_characterizing_2007} that can also be used to construct density matrices. However, the formalism is not limited to these bases, and, as we will show, can be used to analyze dynamics with other (observable) representations.

\subsection{Generalized Liouville's theorem/equation}\label{sec:GLE_norm}

\noindent Imposing a norm-preserving dynamics not only normalizes the density matrix, it also reinstates the original form of Liouville's equation for non-Hamiltonian dynamics.
To see this, we again take the basis to span the phase space $k\to n$.
As must be the case for a norm-preserving dynamics, the magnitude of perturbations do not evolve with time and $\Tr{\bm{\varrho}}(t) = \Tr{\bm{\varrho}}(t_0) = 1$.
The preservation of the trace $\operatorname{Tr}{\bm{\varrho}}$ follows from Eq.~\ref{equ:EOM_rho}, which shows $\Tr\dot{{\bm{\varrho}}}=0$ because $\Tr \{{\bm{A}}_+,{\bm{\varrho}}\} = 2r$.
However, both the trace and the determinant
\begin{equation}\label{eq:Liouville_rho}
  |{\bm{\varrho}}(t)| = |{\bm{\varrho}}(t_0)|,
\end{equation}
are similarity-invariant constants of motion, Figs.~\ref{fig:Lorenz_trace_det} and~\ref{fig:manifold}.
Defining the metric tensor $\boldsymbol{g}_{\bm{\varrho}}$ such that its determinant $g_{\bm{\varrho}} = |{\bm{\varrho}}^{-1}|$,
the conservation of the determinant is another form of the generalized Liouville's equation,
\begin{align}

  \frac{\partial}{\partial t}|{\bm{\varrho}}|^{-\frac{1}{2}} + \dot{{\bm{x}}}\cdot \grad |{\bm{\varrho}}|^{-\frac{1}{2}} &= 0.
\end{align}
in terms of ${\bm{\varrho}}$. The compatibility condition of the metric is
$\partial_t \sqrt{g_{\bm{\varrho}}} +  \dot{{\bm{x}}}\cdot \grad \sqrt{g_{\bm{\varrho}}} = 0$.

Imposing a norm-preserving dynamics and normalizing the density matrix, however, does not change the form of the generalized Liouville theorem.
With these norm-preserving dynamics, the generalized Liouville theorem in Eq.~\ref{eq:invariant_measure} becomes: $|\boldsymbol{\Gamma}|d\mathcal{V}(t)=|{\bm{M}}||\boldsymbol{\Gamma}| d\mathcal{V}(t_0) =|\tilde{\bm{M}}| d\mathcal{V}(t_0) = d\mathcal{V}(t_0)$ and we can recognize $|\boldsymbol{\Gamma}|^2=|{\bm{\xi}}(t_0)|/|{\bm{\xi}}(t)|$.
Geometrically, the generalized Liouville theorem here is: scaled phase space volumes are conserved under the norm-preserving evolution of a perturbation state with a basis that spans the $n$-dimensional phase space.
This generalization of Liouville's theorem is not limited to Hamiltonian dynamics, however.
In phase space of non-Hamiltonian systems, any part of the initial volume lost (or gained) in course of the dynamics is continually and entirely compensated for by the stretching/contracting of the volume. As a result, the scaled volume $|\boldsymbol{\Gamma}|d\mathcal{V}$ is conserved for the $n$-dimensional phase space of any dynamical system.
For both the harmonic and the damped harmonic oscillator, if the initial volume is set by $|{\bm{\xi}}(t_0)| = 1$ and $\Tr{\bm{\xi}}(t_0)=2$, then $|{\bm{\varrho}}(t)|=1/4\,\forall t$.
Figure~\ref{fig:manifold}(a) shows $\operatorname{Tr}{\bm{\xi}}$ for the classical H\'enon Heiles system on a regular and chaotic orbit.
\begin{table*}[ht]
	\caption{Equation of motion and instantaneous Lyapunov exponent for the pure state ${\bm{\varrho}}_i = \dyad{\delta{\bm{\phi}}_i}{\delta{\bm{\phi}}_i}$ in the eigenbases of ${\bm{A}}_+$, ${\bm{A}}_-$, and ${\bm{A}}$, for Hamiltonian and dissipative dynamics.}
	\label{table:summary}

	\begin{tabular}{cccc}
		\hline
		\hline
		& ${\bm{A}}_+$ & ${\bm{A}}_-$ & ${\bm{A}}$ \\
		\hline
		\multirow{2}{6.5em}{\textit{Hamiltonian}} & $\dot{{\bm{\varrho}}_i}= [{\bm{\mathcal{H}}}_-,{\bm{\varrho}}_i]$ & $\dot{{\bm{\varrho}}_i}=\{{\bm{\mathcal{H}}}_+,{\bm{\varrho}}_i\}$& $\dot{{\bm{\varrho}}_i}=0$\\
		& $r^{{\bm{\mathcal{H}}}_+}_i = \operatorname{Eig}{\bm{\mathcal{H}}}_{+,i}$ & $r_i^{{\bm{\mathcal{H}}}_-} = 0$ & $r_i^{{\bm{\mathcal{H}}}}=\Re(\operatorname{Eig}{\bm{\mathcal{H}}}_i)$ \\
		\hline
		\multirow{2}{6em}{\textit{Dissipative}} & $\dot{{\bm{\varrho}}_i}= [{\bm{A}}_-,{\bm{\varrho}}_i]$ & $\dot{{\bm{\varrho}}_i}=\{{\bm{A}}_+,{\bm{\varrho}}_i\}-2r^{{\bm{A}}_-}_i{\bm{\varrho}}_i$& $\dot{{\bm{\varrho}}_i}=0$ \\
		& $r_i^{{\bm{A}}_+} = \operatorname{Eig}{\bm{A}}_{+,i}$ & $r_i^{{\bm{A}}_-} = \Tr({\bm{A}}_+{\bm{\varrho}}_i)$ & $r^{{\bm{A}}}_i=\Re(\operatorname{Eig}{\bm{A}}_i)$\\
		\hline
		\hline
	\end{tabular}
\end{table*}

\subsection{Basis representation} \label{sec:basis}
\noindent In quantum statistical mechanics, the choice of basis states provides an explicit matrix representation of the quantum state.
Here, there are also different basis sets that can be used for perturbation states: ${\bm{\varrho}}$ and ${\bm{\xi}}$ (Appendix~\ref{sec:basis_rep}).
In dynamical systems theory, it is common to analyze Lyapunov vectors, such as Gram-Schmidt vectors~\cite{Benettin1976Kolmogorov,shimada_numerical_1979}.
More recently, however, there has been an interest in covariant Lyapunov vectors~\cite{f_ginelli_characterizing_2007}.
Based on early work~\cite{DellagoPH96}, Lyapunov vectors with small, but finite, exponents are hydrodynamic modes that characterize macroscopic transport~\cite{McNamaraM01,YangR05}.

The choice of basis has consequences for the form of the equation of motion for the density matrix (Appendix~\ref{sec:basis_rep} and Table~\ref{table:summary}).
A natural basis that spans the phase space is the eigenvectors of ${\bm{A}}_+$ because the instantaneous Lyapunov exponents are the associated eigenvalues.
In this basis, the equation of motion for each ${\bm{\varrho}}_i=\dyad{\delta{\bm{\phi}}_i}{\delta{\bm{\phi}}_i}$ reduces to the commutator $\dot{\bm{\varrho}}_i = [{\bm{A}}_-,{\bm{\varrho}}_i]$.
However, choosing the eigenvectors of ${\bm{A}}_-$ as the basis commutator vanishes and $\dot{\bm{\varrho}}_i = \{{\bm{A}}_+,{\bm{\varrho}}_i\} - 2r_i^{{\bm{A}}_-}{\bm{\varrho}}_i$.
The eigenbasis of ${\bm{A}}$ appears to be unique in that ${\bm{\varrho}}$ is time independent.
Hamiltonian dynamics with $H(q^i, p_i)\in\mathbb{R}$ lead to interesting consequences for these basis-dependent results (Appendix~\ref{sec:DMD_Hamiltonian}).
For example, in the ${\bm{\mathcal{H}}}_-$ basis: $\dot{\bm{\varrho}}_i = \{{\bm{\mathcal{H}}}_+,{\bm{\varrho}}_i\}$.
But, in the ${\bm{\mathcal{H}}}$ basis, $\Tr{\bm{\xi}}_i$ is a constant of motion, $\{\Tr{\bm{\xi}}_i,H\}_P=0$.

Hamiltonian systems also have special tangent space directions associated with conserved quantities. On a constant energy manifold, for example, there are two conjugate tangent directions: the phase velocity $\dot {\bm{x}}$ and the gradient of the Hamiltonian $\gradH$, Fig.~\ref{fig:manifold}(b).
They are related $\dot {\bm{x}}=\bm{\Omega} \gradH$ through the Poisson matrix $\bm{\Omega}$, orthogonal to each other, $\gradH\cdot \bm{\Omega} \grad H = 0$, and have equal magnitude $\|\dot {\bm{x}}\|=\|\grad H\|$ through Hamilton's equations.
The vector $\gradH$ is also orthogonal to the constant energy manifold and used to define the invariant measure~\cite{dorfman1999introduction}.
In general, the vector sets defining the density matrices here need not span the whole phase space.
So, forming density matrices from each conjugate vector, we find the pure states, ${\bm{\varrho}}_{\dot {\bm{x}}}$ and ${\bm{\varrho}}_{\bm{\Omega}\grad H}$, both have unit trace and are related by ${\bm{\varrho}}_{\dot {\bm{x}}} = \bm{\Omega}{\bm{\varrho}}_{\grad H}\bm{\Omega}^\top$ (Appendix~\ref{sec:DMD_Hamiltonian}). The conjugate pure states in the tangent space for a 2D-Hamiltonian system are:
\begin{align}
\|\dot {\bm{x}}\|^2{\bm{\varrho}}_{\dot {\bm{x}}} &=
\begin{pmatrix}
\dot q^2 & \dot p\dot q \\[2pt]
\dot p\dot q & \dot p^2\\
\end{pmatrix}, \nonumber\\
\|\gradH\|^2{\bm{\varrho}}_{\gradH}
&=
\begin{pmatrix}
\dot p^2 & -\dot p\dot q\\
-\dot p\dot q & \dot q^2
\end{pmatrix}.\label{eq:tangent_pure_states}
\end{align}
These states are formed from the outer product of the unit tangent vectors, $\ket{\delta{\bm{\phi}}_{\dot{{\bm{x}}}}} = \|\dot{{\bm{x}}}\|^{-1}(q,p)^\top$ and $\ket{\delta{\bm{\phi}}_{\grad{H}}} = \|\gradH\|^{-1}(-\dot{p},\dot{q})^\top$,
where $\|\dot {\bm{x}}\|^2 = \|\gradH\|^2 = \dot{p}^2 + \dot{q}^2$.

These particular density matrices have an interesting (lower dimensional) parallel with Liouville's theorem. Liouville's theorem can be thought of as an equivalence of the divergence of the phase flow, the trace of the Jacobian, and the intrinsic rate:
$d_t\ln \delta\mathcal{V}=\grad\cdot\dot{{\bm{x}}} = n\Tr[{\bm{\varrho}}{\bm{\mathcal{H}}}(\boldsymbol{x})]$,
for ${\bm{\varrho}} = n^{-1}\sum_{i=1}^{n}\dyad{\delta{\bm{\phi}}_i}{\delta{\bm{\phi}}_i}$. There is a similar equivalence for the instantaneous Lyapunov exponents of $\dot{\bm{x}}$ and $\grad H$.
For a $2D$-Hamiltonian system, these take the form:
\begin{align}\nonumber
\frac{d\ln\|\dot {\bm{x}}\|}{d t} \nonumber &= \frac{\dot q\cdot\dot p}{\|\dot {\bm{x}}\|^2}\grad\cdot
\begin{pmatrix} \dot p \\ \dot q \end{pmatrix} = \Tr ({\bm{\varrho}}_{\dot x}{\bm{\mathcal{H}}}_+), \\
\frac{d \ln\|\grad H\|}{d t}  &= \frac{-\dot q\cdot\dot p}{\|\grad H\|^2}\grad\cdot
\begin{pmatrix} \dot p \\ \dot q \end{pmatrix} = \Tr({\bm{\varrho}}_{\grad H}{\bm{\mathcal{H}}}_+),
\end{align}
denoting the phase point ${\bm{x}}=(q,p)$.
Unlike, the intrinsic rate of the volume element $\delta\mathcal{V}=\delta q\delta p$ that appears in Liouville's theorem, these instantaneous Lyapunov exponents are not zero -- they are related to the divergence in a common direction $(\dot p, \dot q)^\top$ that is a reflection about $\dot p= \dot q$.
Instead, they are conjugate, so they sum to zero and span an area that is conserved.

\section{Case studies}\label{sec:case_studies}

To illustrate this density matrix formalism and show that the main ingredients of the theory are computable, we apply it analytically to the linear and damped harmonic oscillator and numerically to the H\'enon-Heiles and Lorenz models (Appendix~\ref{sec:case_examples}).
Below, we explicitly derive the instantaneous Lyapunov exponents of the damped and simple harmonic oscillator in special tangent space directions.
Analytical expressions for these exponents in the eigenvector bases of the stability matrix ${\bm{A}}$ are summarized in Table~\ref{table:HOsummary}.

\begin{table*}[ht]
	\caption{Instantaneous Lyapunov exponents for the linear and damped harmonic oscillator in the eigenvector basis of the stability matrix ${\bm{A}}$ and its (anti)symmetric parts. The harmonic oscillator with frequency $\omega$ has the Hamiltonian: $H = p^2 + \omega^2q$; in this case, $\{\Tr {\bm{\varrho}},H\} = 0$. The equations of motion for the damped harmonic oscillator with the dissipation parameter $\gamma$ are $\dot{q} = p$ and $\dot{p} = -\omega^2q - \gamma p$.}
	\label{table:HOsummary}

	\begin{tabular}{cccc}
		\hline
		\hline
		Harmonic oscillator & ${\bm{A}}_+$ & ${\bm{A}}_-$ & ${\bm{A}}$\\
		\hline
		\vspace{2pt}
		\textit{Linear} & $r_i^{{\bm{\mathcal{H}}}_+} = \pm (1-\omega^2)/2$&
		\vspace{2pt}
		$r_i^{{\bm{\mathcal{H}}}_-} = 0$&
		$r_i^{{\bm{\mathcal{H}}}} = 0$\\[2pt]
		\textit{Damped} & $r_i^{{\bm{A}}_+} = \frac{1}{2}(-\gamma \pm \sqrt{\gamma^2 + (1-\omega^2)^2})$ &
		$r_i^{{\bm{A}}_-} = -\frac{1}{2}\gamma$ &
		$r_i^{{\bm{A}}} = \frac{1}{2}(-\gamma \pm \Re (\sqrt{\gamma^2 - 4\omega^2}))$\\[2pt]
		\hline
		\hline
	\end{tabular}
\end{table*}
\medskip

\noindent\textit{a. Linear harmonic oscillator.--} The simple one-dimensional harmonic oscillator for unit mass and frequency $\omega$ has the Hamiltonian $H = \frac{1}{2}(p^2 + \omega^2 q^2)$. From the usual equations of motion, $\dot{q} = p$, and $\dot{p} = -\omega^2\,q$, the stability matrix ${\bm{\mathcal{H}}}$ is time independent,
\begin{align}
{\bm{\mathcal{H}}} &= \left(
\begin{array}{cc}
0 & 1\\
-\omega^2 & 0\\
\end{array}\right),
\end{align}
with symmetric and anti-symmetric parts:
\begin{align}
{\bm{\mathcal{H}}}_+&= \frac{1}{2}\left(
\begin{array}{cc}
0 & 1-\omega^2\\
1-\omega^2 & 0\\
\end{array}\right), \nonumber\\
{\bm{\mathcal{H}}}_-&= \frac{1}{2}\left(
\begin{array}{cc}
0 & 1+\omega^2 \\
-1-\omega^2 & 0\\
\end{array}\right).
\end{align}
In the basis of ${\bm{\mathcal{H}}}$ and ${\bm{\mathcal{H}}}_-$, the instantaneous Lyapunov exponents vanish, Appendix~\ref{sec:DMD_Hamiltonian}.
Only the instantaneous Lyapunov exponents in the ${\bm{\mathcal{H}}}_+$ eigenvector basis is nonzero, Table~\ref{table:HOsummary}. Because the dynamics are Hamiltonian, there is a conjugate pair of exponents with a magnitude that depends on the oscillator frequency $\omega$.
Appendix Figs.~\ref{fig:SI_HO}(a) and~\ref{fig:SI_HO}(b) show the time evolution of the instantaneous Lyapunov exponents and $\Tr{\bm{\xi}}$ for an arbitrary pure state (formed from a random unit perturbation vector) on a sampled trajectory.

The instantaneous Lyapunov exponents in the conjugate tangent space directions $\dot{\bm{x}}$ and $\grad H$ for the linear harmonic oscillator are also a conjugate pair. From Eq.~\ref{eq:tangent_pure_states}, the pure states are:
\begin{align}
{\bm{\varrho}}_1 &=  \|\dot {\bm{x}}\|^{-2}
\begin{pmatrix}
p^2 &	-\omega^2pq\\
-\omega^2pq &	\omega^4q^2
\end{pmatrix}\nonumber \\
{\bm{\varrho}}_2 &=  \|\gradH\|^{-2}\begin{pmatrix}
\omega^4q^2 &	\omega^2pq\\
\omega^2pq & p^2\\
\end{pmatrix},
\end{align}
where $\|\dot {\bm{x}}\|^{2} = \|\gradH\|^2 = p^2 + \omega^4q^2$. The instantaneous Lyapunov exponents for these states,
\begin{align}\label{eq:ILE_HO}
r_{1,2} & = \Tr({\bm{\mathcal{H}}}_{+}{\bm{\varrho}}_{1,2})= \pm (1-\omega^2)\omega^2\frac{pq}{p^2+\omega^4q^2},
\end{align}
are functions of time, as shown in Appendix Fig.~\ref{fig:SI_HO}(c) using a frequency of $\omega = 0.5$.
In the eigenbasis of ${\bm{\mathcal{H}}}_+$, the instantaneous Lyapunov exponents are (from Appendix Eq.~\ref{eq:ILE_symA}): $r^{{\bm{\mathcal{H}}}_+}_{1,2} = \pm \frac{1}{2}(1-\omega^2)$.
Therefore, the instantaneous Lyapunov exponents in these special directions are the eigenvalues of ${\bm{\mathcal{H}}}_+$ scaled by a function of the state.

These results can be generalized to higher dimensions, as we show for the H\'enon-Heiles system in Appendix~\ref{sec:case_examples}. Figure~\ref{fig:manifold}(c) shows the conjugate pair of instantaneous Lyapunov exponents for the H\'enon Heiles system.

\smallskip

\noindent\textit{b. Damped harmonic oscillator.--} For the damped harmonic oscillator and a velocity dependent dissipation parameterized by $\gamma$, the equations of motion are:
\begin{align}\label{eq:DHO}
\dot{q} &= p  \nonumber\\
\dot{p} &= -\omega^2 q - \gamma p.
\end{align}
Conservative dynamics are recovered for $\gamma = 0$.
The stability matrix of this system,
\begin{align}
{\bm{A}} =\left(
\begin{array}{cc}
0 & 1\\
-\omega^2 & -\gamma\\
\end{array}\right),
\end{align}
has symmetric and anti-symmetric parts:
\begin{align}\nonumber
{\bm{A}}_+ &=\frac{1}{2}\left(
\begin{array}{cc}
0 & 1-\omega^2\\
1-\omega^2 & -2\gamma\\
\end{array}\right), \\
{\bm{A}}_- &=\frac{1}{2}\left(
\begin{array}{cc}
0 & 1+\omega^2\\
-1-\omega^2 & 0\\
\end{array}\right).
\end{align}
Appendix Figs.~\ref{fig:SI_DHO}(a) and~\ref{fig:SI_DHO}(b) show the time evolution of $\Tr{\bm{\xi}}$ and instantaneous Lyapunov exponent for a pure state given by a random unit perturbation vector.

The instantaneous Lyapunov exponents in the conjugate tangent space directions  for the damped harmonic oscillator are found, as before, from the pure states in Eq.~\ref{eq:tangent_pure_states}:
\begin{align}
{\bm{\varrho}}_1 &=  \|\dot {\bm{x}}\|^{-2}
\begin{pmatrix}
p^2 &	-(\omega^2q +\gamma p)p\\
-(\omega^2q +\gamma p)p &	(\omega^2q + \gamma p)^2
\end{pmatrix}\nonumber \\
{\bm{\varrho}}_2 &=  \|\gradH\|^{-2}\begin{pmatrix}
(\omega^2q + \gamma p)^2 &	(\omega^2q +\gamma p)p\\
(\omega^2q +\gamma p)p & p^2\\
\end{pmatrix},
\end{align}
where $\|\dot {\bm{x}}\|^{2} = \|\gradH\|^2 = p^2 + (\omega^2q + \gamma p)^2$. The instantaneous Lyapunov exponent for these states,
\begin{align}
r_{1} & = \Tr({\bm{A}}_{+}{\bm{\varrho}}_{1})\nonumber\\
&= \frac{-(\omega^2q + \gamma p)}{p^2+(\omega^2q + \gamma p)^2}\Big[(1-\omega^2)p+ \gamma(\omega^2q+\gamma p)\Big], \nonumber\\[5pt]
r_{2} & = \Tr({\bm{A}}_{+}{\bm{\varrho}}_{2})\nonumber\\
& = \frac{\omega^2p}{p^2+(\omega^2q + \gamma p)^2}\Big[(1-\omega^2)q - \gamma p\Big].
\end{align}
are shown in Appendix Fig.~\ref{fig:SI_DHO}(c) using a frequency of $\omega = 0.5$ and $\gamma=0.05$.
Note that Eq.~\ref{eq:ILE_HO} is recovered when $\gamma = 0$.

We can also compute these instantaneous Lyapunov exponents in the other eigenbases. For instance, in the ${\bm{A}}_+$-basis, its eigenvalues are:
\begin{align}
r^{{\bm{A}}_+}_{1,2} &= \frac{1}{2}(-\gamma \pm \sqrt{\gamma^2 + (1-\omega^2)^2}).
\end{align}
In the ${\bm{A}}$ eigenbasis, these are $r_i^{{\bm{A}}} = \Re(\lambda^{\bm{A}}_{i})$ (using Appendix  Eq.~\ref{eq:ILE_A}),
\begin{align}
r_i^{{\bm{A}}} = \Re \lambda^{\bm{A}}_i = \Re \frac{1}{2}(-\gamma \pm \sqrt{\gamma^2 - 4\omega^2}),
\end{align}
where $\lambda^{\bm{A}}_i$ indicates the eigenvalue of ${\bm{A}}$.
The real part of $\lambda^{\bm{A}}_i$ gives the instantaneous Lyapunov exponent in the ${\bm{A}}$ basis. These exponents are unequal when $\gamma^2>4\omega^2$.
In the ${\bm{A}}_-$ basis, we can compute the instantaneous Lyapunov exponents as an expectation value over the normalized density matrix, $r_i^{{\bm{A}}_-} = \Tr ({\bm{A}}_-{\bm{\varrho}}_i)$, where each ${\bm{\varrho}}_i$ is the eigenvectors of ${\bm{A}}_-$ (see Appendix Eq.~\ref{eq:basis_states_A_anti}).
Doing so, we find the instantaneous Lyapunov exponents are directly proportional to the damping coefficient $r_i^{{\bm{A}}_-} = -\gamma/2$.

These results for a dissipative system can be also generalized to higher dimensions, as we show for the Lorenz model in Appendix~\ref{sec:case_examples}.
Appendix Fig.~\ref{fig:SI_Lorenz} shows the instantaneous Lyapunov exponents for the Lorenz model.

\section{Conclusions}\label{sec:conclusions}

\noindent Liouville's equation and theorem are the foundation of statistical mechanics established by Gibbs, Maxwell, and Boltzmann.
Boltzmann, for example, approximated Liouville's equation to derive his $H$-theorem for irreversible processes, making an assumption of ``molecular chaos''.
Here, we have established a density matrix formulation of dynamical systems that explicitly and quantitatively accounts for measures of local instability and chaos -- Lyapunov exponents.
Through this connection, we could derive generalizations Liouville's theorem/equation for any differentiable dynamical system.
And, when the dynamics are Hamiltonian, these generalizations reduce to the traditional forms of the Liouville theorem and Liouville's equation.
We have shown they derive from the properties of classical density matrices, which themselves evolve under an equation of motion akin to the von Neumann equation at the foundation of quantum statistical mechanics.
From these results, the generalized Liouville equation becomes numerically computable and, thus, a new basis for analyzing classical speed limits on observables~\cite{das2021speed}, the spread of perturbations, and the transport of statistical density in the living, synthetic, and engineered dynamical systems across physics.

\begin{acknowledgments}
This material is based upon work supported by the National Science Foundation under Grant No. 2124510 and 1856250.
This publication was also made possible, in part, through the support of a grant from the John Templeton Foundation.
\end{acknowledgments}


\section*{\emph{Essay written for the Gravity Research Foundation 2021 Awards for Essays on Gravitation.}}
\thispagestyle{plain}
\pagestyle{plain}
\newpage

In General Relativity (GR), there has been a focus on scalar curvature invariants (hereafter called curvature invariants) which are scalars constructed from the Riemann, Weyl and Ricci tensors and their covariant derivatives \cite{Stephani}. One can find extensive studies on the mathematical structure of these scalars and the way they can be constructed \cite{Stephani,Ref1,Ref2,Cherubini:2003nj,Coley:2009eb}, as well as case studies for different spacetimes \cite{Mattingly1,Mattingly2,Mattingly3,Mattingly4,Overduin,Ghodrati:2020mtx,Chrusciel:2019xuf,Henry:1999rm,Aksteiner:2018vze,Abdelqader:2014vaa,Page:2015aia,McNutt:2017gjg,Coley:2017vxb,Tavlayan,McNutt:2017paq}. As there are an infinite number of invariants that can be generated, one may ask: what is the minimum number of algebraically independent invariants required in order to generate every other invariant for any given metric? Such a finite set is called a complete set.

In the case of algebraic curvature invariants, which are built from only the Ricci tensor and Weyl tensor, it was conjectured that at most 14 independent scalars are needed to cover all possible sub-spaces \cite{P1,P2,P3,P4,P5,P6}. These independent scalars were categorized into four different types: the Ricci scalar, two complex invariants $I$ and $J$ formed from the Weyl spinor, and three real invariants $I_6$, $I_7$ and $I_8$ formed from the trace-free Ricci spinor and three complex invariants $K$, $L$ and $M$. Carminati-McLenaghan (CM) investigated the aforementioned proposal and showed that the 14 previously suggested invariants fail to describe solutions admitting a perfect fluid \cite{Ref1}. Thus, CM came up with 16 necessary independent invariants to cover all sub-spaces including the degenerate cases such as the Einstein-Maxwell system and the perfect fluid. Later, Zakhary-McIntosh (ZM) \cite{Ref2} constructed  the first complete set of algebraic invariants for all possible type of metrics. The tensorial and spinorial forms of ZM invariants are given in Refs. \cite{Overduin,Eichhorn2} and \cite{Ref2}, respectively. Less is known about the higher order curvature invariants, which involve the covariant derivatives of the Ricci tensor and Weyl tensor up to a particular order. Most efforts here have focused on ``brute-force'' calculations to determine a complete set of invariants up to a particular order \cite{martin2008invar}.

Curvature invariants can give insight into the physical properties of solutions in GR and we will provide here a brief, yet incomplete, list of their application in gravity theories. The CM invariants have been used to study the behaviour of the Riemann tensor for Lorentzian traversable wormholes and warp drives \cite{Mattingly1,Mattingly2,Mattingly3,Mattingly4}, and the ZM invariants have been used to study charged and rotating black holes \cite{Overduin}. Curvature invariants have also been shown to detect the horizon for all stationary black holes and a large class of dynamical black holes \cite{Abdelqader:2014vaa,Page:2015aia,McNutt:2017gjg,Coley:2017vxb,Tavlayan,McNutt:2017paq} and distinguish this hypersurface from a wormhole throat \cite{McNutt}.

In this essay, we will discuss the role of curvature invariants in three important aspects of gravitational theories. We will discuss how curvature invariants can be used to study the nature of singularities \cite{Mattingly1,Mattingly2,Mattingly3,Mattingly4,Overduin,Ghodrati:2020mtx,Chrusciel:2019xuf,Henry:1999rm,Aksteiner:2018vze}, their relationship to gravitational entropy of a spacetime through the Weyl curvature conjecture \cite{Gron1,Gron2,Gron3,Romero,Li,shakerin}, and their use in quantum gravity in dynamical singularity-resolution \cite{Eichhorn,Calmet:2020vuh,Falls} and in loop quantum cosmology \cite{Singh}.

It has been proven that a large class of spacetime metrics necessarily must admit singularities \cite{Hawking1,Hawking2,Senovilla:2006db}. In GR, when a metric has a singularity, it can be very difficult to determine whether this is a physical singularity or merely an artifact of the coordinate system one chooses. The latter are known as coordinate singularities and they can be removed using a coordinate transformation. In the case of the Schwarzschild metric, after rewriting the Schwarzschild solution in terms of the Kruskal-Szekeres coordinates \cite{Kruskal,Szekeres}, it becomes clear that $r = 2M$ (for $c = G =1$) is a coordinate singularity, not a physical one. In most recent GR textbooks, the Kretschmann scalar, $R^{\alpha\beta\mu\nu}R_{\alpha\beta\mu\nu} = \frac{48M^2}{r^6}$, is introduced to deal with the problem. Since this is not infinite at the event horizon of $r = 2M$, we can easily conclude that $r = 2M$ is just a coordinate
singularity and it can be removed via a proper coordinate transformation. Similarly, the singularity at $r = 0$ is present in the Kretschmann scalar and must be a real singularity. But the question is still open, what about the cases where finding a proper transformation is challenging to rewrite the metric in a conformally equivalent way to the initial metric and studying the resulting Penrose diagrams? Is there
any way one could study the nature of the underlying metric singularities? A solution is to work with quantities that are independent of the choice of coordinates and here the curvature invariants are ideally suited to this task.

In two long technical essays, Penrose \cite{Penrose1,Penrose2} raised some concerns about a conceptual and fundamental conflict between the second law of thermodynamics and cosmology. In an attempt to summarize Penrose's points, P.C.W. Davies \cite{Davies} sent a letter to Nature to address the issue. According to Davies, in the big-bang model, material content of the universe began in a thermodynamic equilibrium. The universe today is highly ordered and far from equilibrium. In other words, a gravity-dominated evolution may violate the second law of thermodynamics as long as only thermodynamic entropy is taken into account. Penrose suggested that this problem could be solved by assigning entropy to the gravitational field itself. Penrose also conjectured that  the Weyl curvature tensor could be used as a measure of the gravitational entropy. One can express this conjecture in terms of the ratio of Weyl and Ricci curvature invariants \cite{Wainwright}, $P^2=\frac{C_{\alpha\beta\mu\nu}C^{\alpha\beta\mu\nu}}{R_{\mu\nu}R^{\mu\nu}}$, where $P^2$ represents the gravitational entropy and vanishes at the initial singularity of the universe. This is called Penrose’s Weyl Curvature Conjecture, according to which
the Ricci part of the curvature dominates over the Weyl part at the initial singularity of the universe.

GR admits solutions that include essential, physical singularities which are regions where some scalar curvature diverges at a point in question. Physical singularities indicate a breakdown of
GR, but it is expected that a quantum theory of gravity will somehow cure this. Investigators have been examining proposed corrections to GR that bound certain curvature invariants to values less than appropriate powers of the Planck length within the context of cosmologies and low-dimensional black holes. Loop Quantum Gravity (LQG) uses curvature invariants to probe
physical singularities in homogeneous spacetimes, black hole spacetimes, and inhomogeneous Gowdy spacetimes \cite{Singh}. There is confirmation from phenomenological, analytical,
and numerical investigations that singularities are resolved for various homogeneous spacetimes in Loop Quantum Cosmology (LQC) wherein cosmological spacetimes are quantized using the methods of LQG \cite{Singh}. There’s also evidence for singularity resolution
for the various isotropic and anisotropic Bianchi spacetimes using the LQG methods. Singularities are resolved via non-perturbative gravitational effects because of LQG’s discrete
quantum geometry which bounds the curvature invariants by the Planck curvature scale. This leads to the replacement of the big bang by a big bounce close to the Planck scale which is a
robust result in LQC \cite{Singh}.

Recent developments in computer algebra systems \cite{MacCallum:2018csx} have provided researchers with opportunities to perform lengthy or complex calculations more easily than ever before. The calculation of the  ZM or CM curvature invariants for a given metric are examples of such lengthy and time consuming calculation, and this is now achievable on most desk-top computers. As we have emphasized throughout this essay, calculation of these invariants is beneficial for many areas of theoretical physics such as GR, black hole thermodynamics and quantum gravity. It is worth mentioning that the calculation and analysis of curvature invariants can still present research opportunities for talented graduate students who are interested in the interface between theoretical physics and computer programming.

\section*{References}