\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}
