# Maximum speed of dissipation

Swetamber Das  Jason R. Green<EMAIL_ADDRESS>Department of Chemistry,
University of Massachusetts Boston, Boston, Massachusetts 02125, USA
Department of Physics, University of Massachusetts Boston, Boston,
Massachusetts 02125, USA

(September 3, 2024)

###### Abstract

We derive statistical-mechanical speed limits on dissipation from the
classical, chaotic dynamics of many-particle systems. In one, the rate of
irreversible entropy production in the environment is the maximum speed of a
deterministic system out of equilibrium, $\bar{S}_{e}/k_{B}\geq 1/2\Delta t$,
and its inverse is the minimum time to execute the process, $\Delta t\geq
k_{B}/2\bar{S}_{e}$. Starting with deterministic fluctuation theorems, we show
there is a corresponding class of speed limits for physical observables
measuring dissipation rates. For example, in many-particle systems interacting
with a deterministic thermostat, there is a trade-off between the time to
evolve between states and the heat flux, $\bar{Q}\Delta t\geq k_{B}T/2$. These
bounds constrain the relationship between dissipation and time during
nonstationary process, including transient excursions from steady states.

Introduction.— When a piece of macroscopic matter is set in motion, it will
irreversibly lose energy due to frictional or viscous heat generation [1]. If
the system undergoes these dissipative transitions on finite timescales, it
will expel energy as heat, wasting free energy and producing entropy. But, how
quickly does a, potentially finite-size, classical system of particles incur
these thermodynamic costs as it evolves and interacts with its environment?
Are there limits on the speed at which these observables accumulate during
finite time processes? Answering these questions is important to both
synthetic and biological systems. For example, the microscopic machinery of
the cell includes nanoscale motors, such as kinesin, that perform essential
functions on finely-tuned timescales [2]. Given one second to power motion and
haul molecular cargo, kinesin will dissipate about $650k_{B}T$ of energy to
its environment [3], where $k_{B}$ is Boltzmann’s constant and $T$ is the
temperature. One might then hypothesize that the dissipation of energy and the
production of entropy could constrain the time it takes for such systems to
transition between states [4, 5].

Recently, there have been stochastic thermodynamic predictions of the speed of
irreversible processes driven by gradients in temperature, pressure, and
concentration. These thermodynamic speed limits [6, 7, 8, 9, 10, 11] are
giving new understanding of the trade-off between time and dissipation. In
some forms, these purely stochastic results are analogues of quantum speed
limits [12], which bound the speed at which physical systems (or their
observables) evolve between two distinguishable states. While thermodynamic
speed limits on fluxes of energy and entropy are an emerging feature of
stochastic thermodynamics, another important feature is fluctuation theorems.
However, fluctuation theorems were originally derived for mechanical systems
out of equilibrium [13] with the techniques of dynamical systems theory [14,
15]. Still open is the question of whether thermodynamic speed limits on
dissipation have counterparts that derive from the physical dynamics of
mechanical systems. That is, despite progress in classical speed limits and
the classical limit of quantum speed limits [16, 17, 18, 19], there are not
speed limits on dissipation that derive from atomistic dynamics and apply to
many observables.

In this Letter, we derive a classical speed limit with entropy production,
dissipative flux, energy dissipated as heat, and transport coefficients for
dissipative, continuous-time dynamical systems. We show the mean entropy
production puts a bound on the time for a classical many-body system to
transition between two nonequilibrium states as defined by the integrated
deterministic fluctuation theorem. This trade-off between the time to evolve
between two distinguishable states and the entropy production predicts how
much entropy or energy will be dissipated by a physical process over a finite
period of time or how much time will it take for a process to dissipate an
amount of entropy or energy. Like early fluctuation theorems [20, 21, 22],
this bound can be expressed in terms of Lyapunov exponents, the fundamental
quantities characterizing deterministic chaos. These quantities have been used
to analyze rare trajectories [23], jamming [24], nonequilibrium self-assembly
[25], equilibrium and nonequilibrium fluids [26, 27, 28], and critical
phenomena [29]. The finite system size, finite-time behaviors in these
examples are subject to speed limits on dissipation.

Figure 1: (a) A typical trajectory $\Omega$ evolving in the $n$-dimensional
phase space of a dissipative system from an initial time $t_{0}$ to a final
$t$. An infinitesimal volume surrounding the initial condition contracts at a
rate $\bar{\Lambda}$ over the time interval $\Delta t=t-t_{0}$. The
contraction rate caused by dissipation satisfies the inequality $\Delta
t\bar{\Lambda}\geq(1-\eta)/(1+\eta)$ (see text for details). (b) Time
evolution of the contraction rate $\Lambda$ at $\epsilon=0.5$ of a heat-
conducting oscillator, the 0532 model [30, 31]. A nonzero value of $\epsilon$
imposes a temperature gradient and the system relaxes to a nonequilibrium
steady state. The gray dashed lines mark time windows that follow the speed
limit $\Delta t\bar{\Lambda}\geq 1/2$ for $\eta=1/3$ (Eq. 2). The largest time
window $\Delta t$ in the time series is highlighted.

Dissipation from phase space volume contraction.— Equilibrium statistical
physics has successful methods for counting indistinguishable microscopic
states for a given set of physical conditions [1]. There, the thermodynamic
entropy has long been a measure of the number of states,
$S=k_{B}\ln\mathcal{V}$ using phase space volume as a surrogate for the raw
count. However, away from equilibrium, dissipative forces cause the mechanical
contraction of phase space [14, 15]. To construct measures of
distinguishability between physical states along deterministic trajectories
away from equilibrium, we can take a complementary approach using the
relationship between the dissipation of energy and the contraction rate of the
phase space volume $\mathcal{V}$ [14, 15]. With this approach, we can derive
speed limits on dissipation in nonequilibrium processes by measuring the
number of distinct physical states visited by a classical trajectory during a
time $\Delta t=t-t_{0}$.

The rate of phase space contraction is a theoretical link between dynamics and
distinguishable states that are necessary for speed limits on physical
observables.Consider a $d$-dimensional system of $N$ classical particles at a
point, $\bm{x}\in\mathcal{M}$, in phase space $\mathcal{M}$. As the system
evolves under chosen physical conditions, the state traces a trajectory,
$\Omega=\\{\bm{x}(t);t_{0}\leq t\leq t_{f}\\}$, governed by the equations of
motion $\dot{\bm{x}}=\bm{F}[\bm{x}(t)]$ in an $n=2dN$-dimensional phase space,
$\mathcal{M}$, Fig. 1(a).

Infinitesimal phase space volumes, $\delta\mathcal{V}$, surrounding the phase
point have an intrinsic rate at which they deform – the divergence of the
flow,
$\Lambda:=\delta\dot{\mathcal{V}}/\delta\mathcal{V}=\gradient\cdot\dot{\bm{x}}=\sum^{2dN}_{i}\partial
F_{i}/\partial x_{i}$. The dynamics are dissipative if the time averaged
$\Lambda$ is negative in the long-time limit [32]. The instantaneous
contraction rate $\Lambda$ is readily calculated for biological and engineered
dynamical systems [33] as well as myriad (bio)molecular systems simulated with
molecular dynamics [34]. When these systems interact with a deterministic
thermostat [35], the phase space contraction rate gives the thermodynamic
entropy production rate and transport coefficients [14, 15].

When exchanging energy with its surroundings, the deformation of phase space
volumes also determines the time evolution of the probability density over
phase space $\rho(\bm{x},t)$. The density evolves according to Liouville’s
equation $d_{t}\rho=\partial_{t}\rho+\dot{\bm{x}}\nabla\rho=-\rho\Lambda$.
However, if the number of members of the ensemble is constant, then the
equation of motion becomes $d_{t}\delta\mathcal{V}=+\Lambda\delta\mathcal{V}$
with the solution
$\delta\mathcal{V}(t)=\delta\mathcal{V}(t_{0})e^{+\int_{t_{0}}^{t}\Lambda(t^{\prime})dt^{\prime}}$.
As a trajectory evolves from an initial state ($\bm{x}(t_{0})$,
$\delta\mathcal{V}(t_{0})$) to a final state ($\bm{x}(t)$,
$\delta\mathcal{V}(t)$), the infinitesimal volumes are measure of their
distinguishability through the compression factor,

$\displaystyle{\Gamma(t,t_{0})}:=e^{2\int_{t_{0}}^{t}\Lambda[\bm{x}(t^{\prime})]dt^{\prime}}.$
(1)

We use the square volumes because they are numerically computable as the
determinant of classical density matrices [36] (SM I, II). While the rate
$\Lambda$ at given phase point does not generally have a definite sign, there
are finite time intervals $\Delta t=t-t_{0}$ in which the net contraction
$\int_{t_{0}}^{t}\,\Lambda(t^{\prime})dt^{\prime}$ is negative if the system
is dissipative.

Classical speed limit from phase space contraction.— Choosing a value of the
compression factor defines a criterion that we can use for a speed limit on
the contraction of phase space. Say the dynamics over a time interval $\Delta
t$ meet the condition $\Gamma\equiv\eta$ by contracting local volume,
$\int_{t_{0}}^{t}\Lambda dt^{\prime}<0$. Applying a bilinear approximation
111The bilinear approximation of $e^{2x}$ leads to the inequality
$e^{2x}\geq(1+x)/(1-x)$ for $x\leq 0$. to Eq. 1 gives

$\displaystyle\Delta
t_{\eta}\bar{\Lambda}\geq\frac{1-\eta}{1+\eta}=:f(\eta)\quad\quad(\bar{\Lambda}>0,0<\eta\leq
1),$ (2)

a trade-off between time and the time average phase space contraction rate,
$\bar{\Lambda}=-\Delta t^{-1}\int_{t_{0}}^{t}\Lambda dt$ 222A similar
inequality exists for net expanding phase volumes, $\bar{\Lambda}\leq 0$ (SM
LABEL:sm-speed-limit).. This inequality applies to any dissipative, continuous
time dynamical system. To illustrate, Fig. 1(b) shows $\Lambda$ as a function
of time for a singly thermostatted oscillator for $\eta=1/3$ with intervals
marking time windows that satisfy this speed limit.

This inequality has three immediate interpretations:

(1) Rearranging to $\bar{\Lambda}/f(\eta)\geq\Delta t_{\eta}^{-1}$, we see
that the dissipation rate $\bar{\Lambda}$ is the maximum speed $\Delta
t_{\eta}^{-1}$ at which the system can transition between the distinct initial
and final states. It is also the maximum number of distinguishable states the
system can pass through per unit time. That is, it sets a speed limit on
classical dynamical systems transitioning between any two
$\eta$-distinguishable states. Alternatively, if the dissipate rate
$\bar{\Lambda}$ is unknown, the time interval can be seen as a lower bound.

(2) This inequality is a classical analogue of recent extensions to the
Margolus-Levitin speed limit on quantum dynamics 333A weaker form of the
inequality in Eq. 2 follows from the linear approximation, $e^{x}\geq 1+x$,
which gives $\Delta t\bar{\Lambda}/n\geq 1/2$ and a factor of two that one
might expect comparing to the Margolus-Levitin speed limit. For a systems with
a time-independent $\Lambda$ and constant rate of phase space contraction,
such as the damped harmonic oscillator and the Lorenz-Fetter model, the
inequality also simplifies to $-\Delta t\Lambda\geq 1$.. Quantum mechanical
systems with a finite average energy cannot evolve between two states
distinguishable by a fidelity $\epsilon$ in a time shorter than $\Delta
t_{\perp}\geq h\alpha(\epsilon)/4\langle E\rangle$ [40], where $\alpha$ is the
fidelity factor [41, 42]. Here, a classical dynamical system with a finite
average dissipation rate $\bar{\Lambda}$ cannot evolve between two states
distinguishable by $\eta$ on a time shorter than $\Delta t_{\eta}\geq
f(\eta)/\bar{\Lambda}$. If $\Delta t_{\nicefrac{{1}}{{3}}}$ is time for a
given trajectory to reach the condition $\eta=\nicefrac{{1}}{{3}}$, then
$\Delta t_{\nicefrac{{1}}{{3}}}\geq 1/2\bar{\Lambda}$, which is mathematically
similar to the original Margolus-Levitin bound.

(3) Lastly, in the form $\bar{\Lambda}\Delta t_{\eta}\geq f(\eta)$, the
inequality sets a finite lower bound on the accessible phase space volume (or
the count of the distinct physical states) the system visits over the time
interval.

Regardless of the form we choose, there are limiting cases of interest for the
maximum speed set by this, and other measures of, dissipation. The bound
saturates if the net contraction rate is sufficiently small for higher powers
of $\int dt\Lambda$ (cubic and beyond) to be negligible. This condition is a
result of using Padé approximation, which ensures the power series of
$f(\eta)$ matches the power series of the exponential up to second order:
$e^{x}\approx\frac{2+x}{2-x}\approx 1+x+x^{2}/2$. Another limiting case is
when the system is not dissipative $\bar{\Lambda}=0$. Over the course of a
trajectory, the minimum time would then diverge because the states are
indistinguishable, $\eta\to 1$. However, at the other extreme, if the phase
space contraction rate diverges or when the two states are perfectly
distinguishable, $\eta\to 0$, then the minimum time is zero. As confirmed by
our numerical simulations, we expect actual dissipative dynamics systems to
fall between these extremes where the trade-off between time and phase space
contraction is significant.

This first speed limit on phase space contraction is well connected to
existing results in nonequilibrium statistical mechanics. As in quantum
dynamics, the distinguishability criterion $\Gamma=\eta$ determines the
numerical value for the minimum time or maximum speed. Physical motivation for
the criterion in Eq. 1 comes from the intuition that phase space contracts
when a system is dissipative but also from the appearance of this criterion in
deterministic fluctuation theorems. The compression factor $\Gamma$ in Eq. 1
is the probability that phase space contracts by an amount $A$ relative to the
probability it expands by $-A$ over a time interval $\Delta t$ up to a two.
According to the transient fluctuation theorem [43, 13]

$\displaystyle\frac{p_{\Delta t}(\bar{\Lambda}=+A)}{p_{\Delta
t}(\bar{\Lambda}=-A)}=e^{+\Delta tA}=:+\sqrt{\eta},$ (3)

increasing the system size or the time interval $\Delta t$ suppresses second
law “violations” 444The derivation of this fluctuation theorem assumes that
the initial phase space distribution is time symmetric but otherwise
arbitrary.. From this perspective, the speed limit is the minimum time for the
relative probability to reach a value $\eta$. Like the fluctuation theorem,
the speed limit applies to finite-size systems away from equilibrium.
Recognizing this direct connection to fluctuation theorems makes it possible
to derive a family of speed limits with physical observables, including the
entropy production, heat, and transport coefficients.

Speed limit from the entropy production rate.— The first step to deriving
statistical mechanical speed limits is to instead consider an ensemble of
systems. Again, the probability density $\rho(\bm{x},t)$ of the ensemble
evolves between two arbitrary distributions $\rho(\bm{x},t_{0})$ and
$\rho(\bm{x},t)$ over a time interval $\Delta t=t-t_{0}$ according to
Liouville’s equation. At any moment of time in this interval, the Gibbs
entropy rate is the ensemble average of the phase space contraction rate [45,
46, 47]:

$\dot{S}_{G}/k_{B}=+\int
d\bm{x}\rho(\bm{x},t)\Lambda(\bm{x},t)=+\langle\Lambda\rangle.$ (4)

This rate is negative for strongly chaotic and dissipative dynamical systems
with no vanishing Lyapunov exponents [46, 47], so a convention is to add an
overall negative sign and interpret this as the entropy production rate in the
environment [48].

Another fluctuation theorem applies to these statistical ensembles of a system
producing entropy away from equilibrium. The integrated fluctuation theorem
measures probability of second law violations of a small system in a time
interval $\Delta t$ [13]. However, as our first speed limit here suggests, we
can also use fluctuation theorems as a condition for the statistical
distinguishability of states and derive corresponding speed limits. Over a
finite time interval, a system produces more entropy than it consumes by a
factor $P_{\Delta t}(\bar{\Lambda}<0)/P_{\Delta t}(\bar{\Lambda}>0)=\langle
e^{-\Delta t\bar{\Lambda}}\rangle_{\bar{\Lambda}>0}$.

The probabilities here are for any negative (positive) values of
$\bar{\Lambda}$ in the time interval $\Delta t$ across the phase space.
Ensemble averaging Eq. 1,

$\displaystyle\langle\eta\rangle:=\langle\Gamma(t,t_{0})\rangle=\left\langle
e^{-2\Delta t\bar{\Lambda}}\right\rangle_{\bar{\Lambda}>0},$ (5)

the value of $\langle\eta\rangle$ defines a statistical distinguishability
criterion based upon how much more likely a set of trajectory segments is to
produce than consume entropy.

From the exponential in Eq. 5, Jensen’s inequality gives:
$\langle\eta\rangle\geq e^{-2\Delta t\overline{S}_{e}/k_{B}}$

where $\overline{S}_{e}/k_{B}=-\Delta t_{\eta}^{-1}\int
dt\langle\Lambda\rangle\geq 0$ is the total entropy produced in the
environment. Applying the bilinear approximation as before and defining the
function
$f(\langle\eta\rangle):=(1-\langle\eta\rangle)/(1+\langle\eta\rangle)$, we get

$\displaystyle\Delta t_{\langle\eta\rangle}\overline{S}_{e}/k_{B}\geq
f(\langle\eta\rangle),$ (6)

a trade-off between the Gibbs entropy rate and the time to transition between
states.

This speed limit has three interpretations like the speed limit in Eq. 2. It
also has the advantage that the ensemble average makes the bound independent
of the initial condition for the underlying deterministic dynamics. The bound
instead applies to an ensemble transitioning between equilibrium or
nonequilibrium states described by a probability density. For
deterministically thermostatted systems, the Gibbs entropy rate is the
thermodynamic entropy production rate [49, 50], so this inequality is a
tradeoff between thermodynamic dissipation and time for any system
transitioning between two statistical states distinguishable by
$\langle\eta\rangle$. There is a corresponding bound for intensive
observables, SM III.

Stepping back, these two speed limits are closely connected and part of a much
larger family of speed limits involving observables. One immediate bound
applies to mixing dynamical systems: The rate $\Lambda(t)$ is the sum over all
instantaneous Lyapunov exponents, $\Lambda(t)=\sum_{i}^{n}\lambda_{i}(t)$, but
we can simply restrict the sum to positive exponents to derive a speed limit
set by the finite-time Kolmogorov Sinai entropy [25], another signature of
deterministic chaos $\Delta t_{\eta}h_{\scriptscriptstyle\textrm{KS}}\geq 1$.
However, looking at the many known fluctuation theorems [13] suggests the
approach here gives bounds on many physical quantities. These follow from
recognizing Equations 2 and 6 as two fundamental realizations of the
dissipation function, $\bar{\Lambda}$ and $\overline{\langle\Lambda\rangle}$.
Following the steps here for other expressions of the dissipation function
then leads directly to speed limits on other observables [13], including the
dissipative flux (SM IV), heat rate, and transport coefficients. From a given
fluctuation theorem, the minimum transition time (a quantity that is difficult
to predict a priori) can be predicted by a physically measurable quantity. To
illustrate, we take a couple of important examples: the energy dissipated as
heat and the electrical conductivity.

Speed limit from energy dissipated as heat.—

Another immediate physical observable subject to a speed limit is the energy
exchanged as heat between a many-particle system and its environment. For
physical systems interacting with deterministic, time-reversible thermostats
(e.g., Gaussian, Nosé-Hoover), the phase space contraction rate and the Gibbs
entropy rate are a link between the properties of microscopic dynamics [49,
51] and nonequilibrium states [34]. Not only is the (negative) Gibbs entropy
rate the thermodynamic entropy production rate [52, 30, 32, 51, 53]

it is also directly related to the energy is dissipated as heat at a mean rate
$\dot{Q}$, ${\dot{S}_{e}/k_{B}=\dot{S}_{G}/k_{B}=\dot{Q}/k_{B}T}$ 555We use
the convention that the heat rate is negative when the system dissipates
energy to the surroundings..

Using this identity, Eq. 6 becomes

$\displaystyle\Delta t_{\eta}\bar{Q}\geq k_{B}T,$ (7)

an energy-time tradeoff with the temperature setting a lower bound on the
minimum time for the system to exchange an amount of energy as heat with its
environment, $\bar{Q}/k_{B}T=-\Delta t_{\eta}^{-1}\int\langle\Lambda\rangle
dt$. A higher heat exchange rate relative to the typical energy fluctuations,
$k_{B}T$, will lower the minimum transition time between two statistical
states.

Figure 2:  Increasing the entropy rate suppresses the transition time in a
thermally-conducting oscillator.— (a) The time elapsed for the net dissipation
as measured by $\langle\Gamma\rangle$ to reach the threshold $0<\eta<1$ from
numerical simulations of the 0532 model, $(\alpha,\beta)=(0.05,0.32)$. Data
are for three representative values of the temperature gradient parameter
$\epsilon\in\\{0.1,0.3,0.5\\}$ (solid blue, green, and pink lines). The
minimum time for the ensemble to reach those $\eta$ thresholds (dashed lines)
set by the entropy produced in the environment
$f(\langle\eta\rangle)k_{B}/\bar{S}_{e}$ with
$f(\langle\eta\rangle)=(1-\langle\eta\rangle)/(1+\langle\eta\rangle)$.
Averages are over $10^{5}$ initial conditions $(q_{0},p_{0},\xi_{0}=0)$
randomly chosen from the uniform distribution in the interval $\\{-1.5,1.5\\}$
for each value of $\eta$. (b) In the $\zeta=0$ phase plane $10^{6}$ initial
conditions show dissipative structures: the time $t_{\eta}$ to reach
$\eta=\nicefrac{{1}}{{3}}$ is longer for trajectories initiated from points on
logarithmic spirals emanating from the fixed point at $(0,0)$. An orbit of the
harmonic oscillator is shown for comparison (white circle).

Speed limit from transport coefficients.—

As another example, Eq. 6 also extends to transport coefficients. For systems
interacting with Gaussian and Nosé-Hoover thermostats, the electrical
conductivity, diffusion, and viscosity are functions of
$\langle\Lambda\rangle$ [55, 48]. Consider the Lorentz gas in which charged
particles (neglecting Coulomb interactions) move in a two-dimensional array of
fixed hard disk scatterers. Particles are also driven by an external electric
field $\mathbf{E}$ and thermostatted by a reservoir at temperature $T$ [56].
The electrical conductivity is
$\sigma(t)=T/\mathbf{|E|^{2}}\overline{\langle\Lambda(t)\rangle}$. From Eq. 6,
the electric field and temperature set the lower bound: $\Delta
t_{\eta}\sigma(t)\geq T/2\mathbf{|E|^{2}}$ for $\eta=1/3$. Intuitively,
stronger electric fields or higher conductivity can expedite the transition
between states, while higher temperatures will delay the transition. Similar
intuition comes from the simpler case of the harmonic oscillator with friction
coefficient $\kappa$ and mass $m$: for $\eta=1/3$, the speed limit $\Delta
t_{\eta}\kappa/m\geq 1/2$ suggests that decreasing the friction or increasing
the mass will lower the minimum the time to dissipate a given amount of
energy.

Numerical simulations of a thermally-conducting oscillator.— As a numerical
test case, we analyze the dynamics of a one-dimensional oscillator interacting
with a Nosé Hoover thermostat [31], SM V. The position and momentum $(q,p)$ of
the oscillator are thermostatted at temperature $T$ by the linear friction
force $-\zeta p$. For any nonzero with friction coefficient $\zeta$, the phase
space contraction rate is $\Lambda(t)=-\zeta\left(\alpha+3\beta
p^{2}/T\right)$. A negative (positive) $\Lambda$ indicates heat lost to
(gained from) the thermostat. To simulate the “0532 model”, we choose
$(\alpha,\beta)=(0.05,0.32)$, which is ergodic with a Gibbsian phase space
distribution [30] at thermal equilibrium. More importantly here, adding a
local temperature gradient with a profile $T=1+\epsilon\tanh(q)$ and parameter
$0<\epsilon<1$ drives the system away from equilibrium. To compute the speed
limit from the entropy production, we chose three representative values of
$\epsilon\in\\{0.1,0.3,0.5\\}$.

Figure 2 shows the transition time and the lower bound
$\overline{\langle\Lambda\rangle}$. These data are averaged over $10^{5}$
initial conditions $(q_{0},p_{0},0.0)$ randomly chosen in the interval
$[-1.5,1.5]$ on the $\zeta=0$ plane for 50 equally spaced values of $\eta$
between 0 and 1. They show that as the ensemble relaxes to the nonequilibrium
steady state attractor, the time intervals in which trajectories first satisfy
the $\eta$ threshold is lower bounded by $f(\eta)/\bar{S}_{e}$. Increasing the
nonequilibrium drive from $\epsilon=0.1$ to $0.5$ increases the entropy
produced and decreases the transition time.

Across the $(q,p)$ phase space, the speed limit resolves structures that are
dissipative on a particular timescale. For $\eta=\nicefrac{{1}}{{3}}$, the
$\zeta=0$ plane in the second panel of Fig. 2. In this case, this speed limit
is: $\Delta t\bar{s}_{i}\geq\frac{1}{2}$. For improved resolution, we randomly
sample $10^{6}$ initial points and color them according to their $\Delta t$ at
which the corresponding trajectories take to attain the
$\eta=\nicefrac{{1}}{{3}}$. Those regions on the $\zeta=0$ plane with higher
values of $\Delta t_{\eta}$ are correlated with the smaller values of
$\bar{\Lambda}$, satisfying our speed limit on dissipation in Eq. 2 in each
case. These phase points dissipating less appear as logarithmic spirals
approaching the fixed point of the system at the origin. These structures show
that there are some regions of phase space that are slower and others that are
faster to transition between states, up to but not exceeding the maximum speed
set by phase space contraction rate.

The speed limits on dissipation here are statistical and mechanical, applying
to many physical observables, and quantifying the natural tension between
dissipation and the time to produce it during irreversible processes. They
results bridge the deterministic, potentially chaotic, dynamics of physical
systems, avoiding any hypotheses about chaos and instead making measures of
chaos explicit through the contraction of phase space. These bounds are a
quantitative statement of how the production of entropy, expulsion of energy
as heat, or transport of matter will decrease the minimum time to evolve
between two distinguishable physical states.

Acknowledgments.— This material is based upon work supported by the National
Science Foundation under Grant No. 2124510.

## References

  * Callen [1985] H. B. Callen,  _Thermodynamics and an introduction to thermostatistics; 2nd ed._ (Wiley, New York, NY, 1985).
  * Goel and Vogel [2008] A. Goel and V. Vogel, Harnessing biological motors to engineer systems for nanoscale transport and assembly, Nat. Nanotechnol 3, 465 (2008).
  * Bustamante  _et al._ [2005] C. Bustamante, J. Liphardt, and F. Ritort, The nonequilibrium thermodynamics of small systems, Phys. Today 58, 43 (2005).
  * Bryant and Machta [2020] S. J. Bryant and B. B. Machta, Energy dissipation bounds for autonomous thermodynamic cycles, Proc. Natl. Acad. Sci. U.S.A. 117, 3478 (2020).
  * Skinner and Dunkel [2021] D. J. Skinner and J. Dunkel, Improved bounds on entropy production in living systems, Proc. Natl. Acad. Sci. U.S.A. 118, e2024300118 (2021).
  * Shiraishi  _et al._ [2018] N. Shiraishi, K. Funo, and K. Saito, Speed limit for classical stochastic processes, Phys. Rev. Lett. 121, 070601 (2018).
  * Hasegawa and Van Vu [2019] Y. Hasegawa and T. Van Vu, Uncertainty relations in stochastic processes: An information inequality approach, Phys. Rev. E 99, 062126 (2019).
  * Ito and Dechant [2020] S. Ito and A. Dechant, Stochastic time evolution, information geometry, and the Cramér-Rao bound, Phys. Rev. X 10, 021056 (2020).
  * Nicholson  _et al._ [2018] S. B. Nicholson, A. del Campo, and J. R. Green, Nonequilibrium uncertainty principle from information geometry, Phys. Rev. E 98, 032106 (2018).
  * Falasco and Esposito [2020] G. Falasco and M. Esposito, Dissipation-time uncertainty relation, Phys. Rev. Lett. 125, 120604 (2020).
  * Nicholson  _et al._ [2020] S. B. Nicholson, L. P. García-Pintos, A. del Campo, and J. R. Green, Time–information uncertainty relations in thermodynamics, Nat. Phys. 16, 1211 (2020).
  * Deffner and Campbell [2017] S. Deffner and S. Campbell, Quantum speed limits: from Heisenberg’s uncertainty principle to optimal quantum control, J. Phys. A Math. Theor. 50, 453001 (2017).
  * Evans and Searles [2002] D. J. Evans and D. J. Searles, The fluctuation theorem, Adv. Phys. 51, 1529 (2002).
  * Gaspard [1998] P. Gaspard,  _Chaos, Scattering and Statistical Mechanics_, Cambridge Nonlinear Science Series (Cambridge University Press, Cambridge, 1998).
  * Dorfman [1999] J. R. Dorfman,  _An Introduction to Chaos in Nonequilibrium Statistical Mechanics_, Cambridge Lecture Notes in Physics (Cambridge University Press, Cambridge, 1999).
  * Shanahan  _et al._ [2018] B. Shanahan, A. Chenu, N. Margolus, and A. del Campo, Quantum speed limits across the quantum-to-classical transition, Phys. Rev. Lett. 120, 070401 (2018).
  * Okuyama and Ohzeki [2018] M. Okuyama and M. Ohzeki, Quantum speed limit is not quantum, Phys. Rev. Lett. 120, 070402 (2018).
  * Poggi  _et al._ [2021] P. M. Poggi, S. Campbell, and S. Deffner, Diverging quantum speed limits: A herald of classicality, PRX Quantum 2, 040349 (2021).
  * Das and Green [2023] S. Das and J. R. Green, Speed limits on deterministic chaos and dissipation, Phys. Rev. Res. 5, L012016 (2023).
  * Evans and Searles [1994] D. J. Evans and D. J. Searles, Equilibrium microstates which generate second law violating steady states, Phys. Rev. E 50, 1645 (1994).
  * Gallavotti and Cohen [1995] G. Gallavotti and E. G. D. Cohen, Dynamical ensembles in nonequilibrium statistical mechanics, Phys. Rev. Lett. 74, 2694 (1995).
  * Cohen and Gallavotti [1999] E. G. D. Cohen and G. Gallavotti, Note on two theorems in nonequilibrium statistical mechanics, J. Stat. Phys. 96, 1343 (1999).
  * Tailleur and Kurchan [2007] J. Tailleur and J. Kurchan, Probing rare physical trajectories with Lyapunov weighted dynamics, Nat. Phys. 3, 1745 (2007).
  * Banigan  _et al._ [2013] E. J. Banigan, M. K. Illich, D. J. Stace-Naughton, and D. A. Egolf, The chaotic dynamics of jamming, Nat. Phys. 9, 288 (2013).
  * Green  _et al._ [2013] J. R. Green, A. B. Costa, B. A. Grzybowski, and I. Szleifer, Relationship between dynamical entropy and energy dissipation far from thermodynamic equilibrium, Proc. Natl. Acad. Sci. U.S.A. 110, 16339 (2013).
  * Evans and Morris [1990] D. Evans and G. Morris,  _Statistical Mechanics of Nonequilibrium Liquids_ (New York, 1990).
  * Bosetti and Posch [2014] H. Bosetti and H. A. Posch, What does dynamical systems theory teach us about fluids?, Commun. Theor. Phys. 62, 451 (2014).
  * Das and Green [2017] M. Das and J. R. Green, Self-averaging fluctuations in the chaoticity of simple fluids, Phy. Rev. Lett. 119, 115502 (2017).
  * Das and Green [2019] M. Das and J. R. Green, Critical fluctuations and slowing down of chaos, Nat. Commun. 10, 2155 (2019).
  * Hoover  _et al._ [2016a] W. G. Hoover, J. C. Sprott, and C. G. Hoover, Ergodicity of a singly-thermostated harmonic oscillator, Commun. Nonlinear Sci. Numer. Simul. 32, 234 (2016a).
  * Patra  _et al._ [2016] P. K. Patra, W. G. Hoover, C. G. Hoover, and J. C. Sprott, The equivalence of dissipation from Gibbs’ Entropy Production with Phase-Volume loss in Ergodic Heat-Conducting Oscillators, Int. J. Bifurc. Chaos Appl. Sci. Eng. 26, 1650089 (2016).
  * Ramshaw [2017] J. D. Ramshaw, Entropy production and volume contraction in thermostated Hamiltonian dynamics, Phys. Rev. E 96, 052122 (2017).
  * Pikovsky and Politi [2016] A. Pikovsky and A. Politi,  _Lyapunov Exponents_ (Cambridge University Press, 2016).
  * Evans and Morriss [2008] D. J. Evans and G. Morriss,  _Statistical Mechanics of Nonequilibrium Liquids_, 2nd ed. (Cambridge University Press, 2008).
  * Hoover [1999] W. G. Hoover,  _Time Reversibility, Computer Simulation, and Chaos_ (World Scientific, 1999).
  * Das and Green [2022] S. Das and J. R. Green, Density matrix formulation of dynamical systems, Phys. Rev. E 106, 054135 (2022).
  * Note [1] The bilinear approximation of $e^{2x}$ leads to the inequality $e^{2x}\geq(1+x)/(1-x)$ for $x\leq 0$.
  * Note [2] A similar inequality exists for net expanding phase volumes, $\bar{\Lambda}\leq 0$ (SM LABEL:sm-speed-limit).
  * Note [3] A weaker form of the inequality in Eq. 2 follows from the linear approximation, $e^{x}\geq 1+x$, which gives $\Delta t\bar{\Lambda}/n\geq 1/2$ and a factor of two that one might expect comparing to the Margolus-Levitin speed limit. For a systems with a time-independent $\Lambda$ and constant rate of phase space contraction, such as the damped harmonic oscillator and the Lorenz-Fetter model, the inequality also simplifies to $-\Delta t\Lambda\geq 1$.
  * Margolus and Levitin [1998] N. Margolus and L. B. Levitin, The maximum speed of dynamical evolution, Phys. D: Nonlinear Phenom. 120, 188 (1998), proceedings of the Fourth Workshop on Physics and Consumption.
  * Giovannetti  _et al._ [2003] V. Giovannetti, S. Lloyd, and L. Maccone, Quantum limits to dynamical evolution, Physical Review A 67, 052109 (2003).
  * Frey [2016] M. R. Frey, Quantum speed limits—primer, perspectives, and potential future directions, Quantum Information Processing 15, 3919 (2016).
  * Evans  _et al._ [1993] D. J. Evans, E. G. D. Cohen, and G. P. Morriss, Probability of second law violations in shearing steady states, Phys. Rev. Lett. 71, 2401 (1993).
  * Note [4] The derivation of this fluctuation theorem assumes that the initial phase space distribution is time symmetric but otherwise arbitrary.
  * Andrey [1985] L. Andrey, The rate of entropy change in non-Hamiltonian systems, Phys. Lett. A 111, 45 (1985).
  * Ruelle [1996] D. Ruelle, Positivity of entropy production in nonequilibrium statistical mechanics, J. Stat. Phys. 85, 1 (1996).
  * Ruelle [1997] D. Ruelle, Entropy production in nonequilibrium statistical mechanics, Commun. Math. Phys. 189, 365 (1997).
  * Klages [2007] R. Klages,  _Microscopic Chaos, Fractals and Transport in Nonequilibrium Statistical Mechanics_ (World Scientific, 2007).
  * Jepps and Rondoni [2010] O. G. Jepps and L. Rondoni, Deterministic thermostats, theories of nonequilibrium systems and parallels with the ergodic condition, J. Phys. A Math. Theor. 43, 133001 (2010).
  * Harish and Patra [2021] M. S. Harish and P. K. Patra, Temperature and its control in molecular dynamics simulations, Mol. Simul. 47, 701 (2021).
  * Hoover and Hoover [2020] W. G. Hoover and C. G. Hoover, From hard spheres and cubes to nonequilibrium maps with thirty-some years of thermostatted molecular dynamics, J. Chem. Phys. 153, 070901 (2020).
  * Hoover  _et al._ [2016b] W. G. Hoover, C. G. Hoover, and J. C. Sprott, Nonequilibrium systems: hard disks and harmonic oscillators near and far from equilibrium, Mol. Simul. 42, 1300 (2016b).
  * Ramshaw [2020] J. D. Ramshaw, Thermodynamic vs statistical entropy production in thermostatted Hamiltonian dynamics, J. Phys. A Math. Theor. 53, 495001 (2020).
  * Note [5] We use the convention that the heat rate is negative when the system dissipates energy to the surroundings.
  * Klages  _et al._ [2000] R. Klages, K. Rateitschak, and G. Nicolis, Thermostating by deterministic scattering: Construction of nonequilibrium steady states, Phys. Rev. Lett. 84, 4268 (2000).
  * van Beijeren and Dorfman [1995] H. van Beijeren and J. R. Dorfman, Lyapunov exponents and Kolmogorov-Sinai entropy for the Lorentz gas at low densities, Phys. Rev. Lett. 74, 4412 (1995).

Supplemental Material: Maximum speed of dissipation

##  I Equation of motion for the trace of the density matrix

A perturbation vector in the tangent space of a classical trajectory evolves
according to:

$d_{t}\ket{\delta\bm{x}_{i}(t)}={\bm{A}}[\bm{x}(t)]\ket{\delta\bm{x}_{i}(t)}.$
(1)

For many-body systems, these vectors are perturbations to positions and
momenta $\delta\bm{x}_{i}=(\delta\bm{q}^{i},\delta\bm{p}_{i})$ . The stability
matrix, ${\bm{A}}:={\bm{A}}[\bm{x}(t)]=\gradient\bm{F}$ generating these
dynamics has the elements $({\bm{A}})^{i}_{j}=\partial\dot{x}^{i}(t)/\partial
x^{j}(t)$. Dirac’s notation here represents a finite-dimensional vector (its
transpose) with the ket (bra): $\ket{\delta\bm{x}(t)}:=[\delta x^{1}(t),\delta
x^{2}(t),\ldots,\delta x^{n}(t)]^{\top}\in T\mathcal{M}$. Solving the equation
of motion [33]:

$\displaystyle\ket{\delta\bm{x}_{i}(t)}$
$\displaystyle=\mathcal{T}_{+}e^{\int_{t_{0}}^{t}{\bm{A}}(t^{\prime})\,dt^{\prime}}\ket{\delta\bm{x}_{i}(t_{0})},$
(2)

gives the perturbation vector at a time $t$ (with the time ordering operator
$\mathcal{T}_{+}$ because stability matrices at two times do not generally
commute).

To span the volume $\delta\mathcal{V}$, we can choose a complete set of $n$
linear independent unit vectors $\\{\ket{\delta{\bm{\phi}}_{i}}$} and
normalize them $\delta{\bm{\phi}}_{i}=\delta\bm{x}_{i}/\|\delta\bm{x}_{i}\|$
so that the outer product defines a classical density matrix,
${\bm{\xi}}(t)=\sum_{i=1}^{n}\outerproduct{\delta{\bm{\phi}}_{i}}{\delta{\bm{\phi}}_{i}}$.
As the basis vectors form complete set of independent vectors, the determinant
of ${\bm{\xi}}$ gives the volume $\delta\mathcal{V}$ span by them in the phase
space.

In our previous work, see Ref. [36], we have shown that the determinant of
${\bm{\xi}}$ denoted by $|{\bm{\xi}}|$ has the following equation of motion:

$\displaystyle d_{t}\ln|{\bm{\xi}}|=2\Tr{\bm{A}}=2\Lambda.$ (3)

This gives rate of change of phase space volume under any non-conservative
dynamics. The evolution of its trace $\Tr{\bm{\xi}}$ is governed by:

$\displaystyle d_{t}\ln\Tr{\bm{\xi}}=2\Lambda/n.$ (4)

For a single tangent vector $\ket{\delta\bm{x}_{i}(t)}$, the trace
$\Tr{\bm{\xi}}_{i}(t)=\bra{\delta\bm{x}_{i}(t)}\ket{\delta\bm{x}_{i}(t)}$
gives the rate at which the perturbation represented by the vector
$\ket{\delta\bm{x}_{i}(t)}$ evolves in time.

The solution to Eq. 4 reads:

$\displaystyle\Tr{\bm{\xi}}(t)=\Tr{\bm{\xi}}(t_{0})e^{2n^{-1}\int_{t_{0}}^{t}\Lambda(t^{\prime})dt^{\prime}}.$
(5)

We use this solution to define the compression factor
$\Gamma(t,t_{0})=\Tr{\bm{\xi}}(t_{0})/\Tr{\bm{\xi}}(t_{0})$ which leads to the
dissipation-time bound in the main text.

##  II Geometric interpretation of the density matrix

$\ket{\delta{\bm{\phi}}_{1}}$$\ket{\delta{\bm{\phi}}_{2}}$$\ket{\delta{\bm{\phi}}_{3}}$$\propto(\Tr{\bm{\xi}})^{1/2}$$\delta\mathcal{V}=|{\bm{\xi}}|^{1/2}$
Figure SM1: A compressed phase space volume at any time $t$ can be mapped onto
a hypersphere in the phase space using the classical density matrix framework.
This schematic shows a contracted phase volume spanned by a complete set of
basis
$\\{\ket{\delta{\bm{\phi}}_{1}},\ket{\delta{\bm{\phi}}_{2}},\ket{\delta{\bm{\phi}}_{3}}\\}$
in a three dimensional phase space with an arbitrary shape. Its volume
$\delta\mathcal{V}$ is given by the square root of determinant of the density
matrix $|{\bm{\xi}}(t)|^{1/2}$, where
${\bm{\xi}}(t)=\sum_{i=1}^{3}\outerproduct{\delta{\bm{\phi}}_{i}}{\delta{\bm{\phi}}_{i}(t)}$.
This contracted phase space is equivalent to a sphere of the same volume
$\delta\mathcal{V}$ and radius proportional to $\Tr{\bm{\xi}}(t)^{1/2}$ at any
time $t$. This geometric picture holds for an $n$-dimensional phase space.

The dynamic of a unit vector in the tangent space can be normalized at every
time. The equation of motion of unit tangent vector
$\ket{\delta{\bm{u}}_{i}}:=\ket{\delta\bm{x}_{i}(t)}/\|\delta\bm{x}_{i}(t)\|$
is given by:

$\displaystyle
d_{t}\ket{\delta{\bm{u}}_{i}(t)}:={\bm{A}}\ket{\delta{\bm{u}}_{i}(t)}-r_{i}(t)\ket{\delta{\bm{u}}_{i}},$
(6)

where
$r_{i}(t)=\bra{\delta{\bm{u}}_{i}(t)}{\bm{A}}(\bm{x})\ket{\delta{\bm{u}}_{i}(t)}$
can be identified as the instantaneous Lyapunov exponent.

For this normalized dynamics, we can construct a density matrix
${\bm{\varrho}}_{i}\ \ket{\delta{\bm{u}}_{i}(t)}\bra{\delta{\bm{u}}_{i}(t)}$
which accounts for the normalized evolution of $\ket{\delta{\bm{u}}_{i}}(t)$
with time. The trace of ${\bm{\varrho}}_{i}$ remains 1 for all times.

As before, we consider the complete set of linearly independent unit vectors
${\ket{\delta{\bm{\phi}}_{i}}}$ to use their outer product to construct
density matrix
${\bm{\varrho}}=n^{-1}\sum_{i=1}^{n}\outerproduct{\delta{\bm{\phi}}_{i}}{\delta{\bm{\phi}}_{i}}$.
Due to the normalized evolution of $\ket{\delta{\bm{\phi}}_{i}}$ for all
times, $\Tr{\bm{\varrho}}=1$ also remains time invariant. The matrices
${\bm{\varrho}}$ and ${\bm{\xi}}$ and their determinants are related as
follows:

$\displaystyle{\bm{\varrho}}$ $\displaystyle=\frac{{\bm{\xi}}}{\Tr{\bm{\xi}}}$
(7) $\displaystyle\implies|{\bm{\varrho}}|$
$\displaystyle=\frac{|{\bm{\xi}}|}{(\Tr{\bm{\xi}})^{n}},$ (8)

at any time $t$ in an $n$-dimensional phase space.

The determinants $|{\bm{\xi}}(t)|$ and $|{\bm{\varrho}}(t)|$ give different
phase space volumes at time $t$. But, since ${\bm{\varrho}}$ generates
normalized evolution, the volume it represents is a constant of motion. This
observation leads to the generalized Liouville’s theorem and equation we have
proposed in Ref.[36].

This shows that $|{\bm{\xi}}(t)|\propto(\Tr{\bm{\xi}}(t))^{n}$, see Fig. SM1.
For a geometric interpretation, imagine a hypersphere in the phase at time $t$
attached to the given trajectory. The volume of this hypersphere is given by
$|{\bm{\xi}}(t)|^{1/2}$ and radius $(\Tr{\bm{\xi}}(t))^{1/2}$. This $n-$sphere
represents a contracted phase space volume in a dissipative system at time
$t$. We have thus used the density matrix formalism to map the compressed
phase space volume of an arbitrary shape at time $t$ to a hypersphere.

In the next section, we use the observation
$|{\bm{\xi}}(t)|\propto(\Tr{\bm{\xi}}(t))^{n}$ to obtain intensive bounds.
These are another set of speed limit on dissipation

##  III Derivation of the speed limit on dissipation (System-size intensive)

Consider the compression factor,

$\displaystyle\Gamma(t,t_{0})=e^{2\int_{t_{i}}^{t_{f}}dt\Lambda/n}.$ (9)

Let $\Delta t_{\eta}$ is the time interval in which $\Gamma$ reduces by a
factor of $0\leq\eta<1$. This leads to

$\displaystyle\eta=e^{2\int_{t_{i}}^{t_{f}}dt\Lambda/n}.$ (10)

Applying a bilinear approximation of $e^{2x}$ for $x\leq 0$, we obtain

$\displaystyle\eta\geq\frac{1-\Delta t_{\eta}\bar{\Lambda}/n}{1+\Delta
t_{\eta}\bar{\Lambda}/n},$ (11)

where $\bar{\Lambda}=-\Delta t_{\eta}^{-1}\int_{t_{i}}^{t_{f}}dt\Lambda$.
Rearranging the inequality,

$\displaystyle\Delta
t_{\eta}\bar{\Lambda}/n\geq\frac{1-\eta}{1+\eta};\quad(\bar{\Lambda}>0,0<\eta\leq
1).$ (12)

This inequality is the system-size intensive version of Eq. 2 (main text). We
also give the entropy speed limit, which is obtained following the steps which
yield Eq. 6 (main text):

$\displaystyle\Delta
t_{\langle\eta\rangle}\overline{S}_{e}/nk_{B}\geq\frac{1-\langle\eta\rangle}{1+\langle\eta\rangle};\quad(0<\langle\eta\rangle\leq
1),$ (13)

where $\overline{S}_{e}/k_{B}=-\Delta t_{\eta}^{-1}\int
dt\langle\Lambda\rangle\geq 0$ is the total entropy produced in the
environment.

A similar set of relations can be found if we consider the expansion of phase
volume $\Lambda\geq 0$ i.e. $\bar{\Lambda}\leq 0$. In dissipative systems,
there may be time intervals in which $\Gamma$ increases implying that
$\eta\geq 1$. The bound in Eq. 2 (main text) then takes the following form:

$\displaystyle\quad-\Delta
t_{\eta}\bar{\Lambda}/n\geq\frac{\eta-1}{\eta+1};\quad(\bar{\Lambda}<0,\eta\geq
1).$ (14)

This inequality holds for small net expansion $0\leq-\Delta
t\bar{\Lambda}=\int dt\Lambda<1$. The ensemble average version of the bound in
Eq. 14 follows as before.

##  IV Speed limit from dissipative flux

In this section, we will show a possible speed limit from the dissipative flux
$J$ in isothermal systems. A typical example is a set of $N$-particle in the
presence of an external dissipative field $F_{e}$. Following Ref. [13], we
give the transient fluctuation theorem for isokinetic dynamics:

$\displaystyle\frac{p(\bar{J}=A)}{p(\bar{J}=-A)}=e^{AtF_{e}V/(k_{B}T)},$ (15)

where $\bar{J}=-\Delta t^{-1}\int dt^{\prime}J$ defines the time average of
the dissipative flux $J$ and $V$ represents the volume of the system of
interest. This theorem holds for all times for the isokinetic ensemble when
all initial phases are sampled from an equilibrium isokinetic ensemble.

We define the compression factor $\eta=e^{2AtF_{e}V/(k_{B}T)}$ in accordance
with our formulation described in the main text and follow Eq. 2 (main text)
to write the following bound :

$\displaystyle\Delta t_{\eta}\bar{J}F_{e}V\geq f(\eta)k_{B}T.$ (16)

The speed limit from the rate of the ensemble average $\langle J\rangle$ has a
form similar to Eq 6 (main text):

$\displaystyle\Delta t_{\eta}\overline{\langle J\rangle}F_{e}V\geq
f(\langle\eta\rangle)k_{B}T.$ (17)

We have derived these bounds with assumption that $F_{e}$ is a constant.
Similar speed limits can be derived from all physical observables which follow
the fluctuation theorem. For a list, see table 4.1 in Ref. [13].

##  V Singly-thermostated harmonic oscillator: The 0532 model

This recently introduced system belongs to the family of Nosé-Hoover
thermostats. It is one of the simplest example to investigate thermodynamic
dissipation in the thermostatted dynamics of the one-dimensional harmonic
oscillator. The system is given by oscillator coordinate $q$, velocity $p$,
and friction coefficient $\zeta$ at a thermostat temperature $T$:

$\displaystyle\dot{q}$ $\displaystyle=p,\quad\dot{p}=-q-\zeta\left(\alpha
p+\frac{\beta p^{3}}{T}\right),$ (18) $\displaystyle\dot{\zeta}$
$\displaystyle=\alpha\left(\frac{p^{2}}{T}-1\right)+\beta\left(\frac{p^{4}}{T^{2}}-\frac{3p^{2}}{T}\right),$

with parameters, $\alpha$ and $\beta$, set to $0.05$ and $0.32$. These
equations impose a weak control of the second and fourth moments of $p$. The
control variable $\zeta$ is responsible for stabilizing the kinetic energy of
the system.

The system is deterministic, time-reversible, and has been shown to be ergodic
through extensive numerical tests. At equilibrium [31, 30] the model
reproduces the Gibbs distribution $f(q,p,\zeta)\propto
e^{-q^{2}/2T}e^{-p^{2}/2T}e^{-\zeta^{2}/2T}$ in units with $k_{B}=1$.
Moreover, in these ergodic thermostats, the time-averaged heat current
$\overline{p^{2}}/2$ measured directly from the model has been found to be in
excellent agreement with the Green-Kubo linear-response theory [31].

To study thermodynamic entropy production, a localized temperature gradient
can be introduced using a smooth temperature profile $T=1+\epsilon\tanh(q)$
with $0\leq\epsilon<1$. The localized temperature within the system thus
varies between $1-\epsilon$ and $1+\epsilon$:

$\displaystyle T(-\infty)=1-\epsilon\leq T(q)\leq 1+\epsilon=T(+\infty).$ (19)

For a nonzero value of $\epsilon$, the system is out of equilibrium. The
temperature gradient $dT/dq=\epsilon/\cosh^{2}(q)$ generates a heat transfer
from (to) the system to (from) the thermostat making the dynamics dissipative
and irreversible.

Figure SM2:  The NESS value of the thermodynamic entropy production in the
0532 thermostat with increasing temperature gradient parameter $\epsilon$.
Averages are over $10^{5}$ initial conditions $(q_{0},p_{0},\xi_{0}=0)$
randomly chosen from the uniform distribution in the interval $\\{-1.5,1.5\\}$
for 12 different values of $\eta$ in the interval $[0.0,0.5]$. For time
averaging, the $\Lambda$ values were taken after $10^{6}$ iterations for each
trajectory. The long time average of $\Lambda$ followed by its large ensemble
average gives the NESS value of $\langle\Lambda\rangle$ for each $\epsilon$.

The phase space volume contraction rate $\Lambda$ is given by:

$\displaystyle\Lambda=\frac{\partial\dot{q}}{\partial
q}+\frac{\partial\dot{p}}{\partial
p}+\frac{\partial\dot{\zeta}}{\partial\zeta}=-\zeta\left(\alpha+\frac{3\beta
p^{2}}{T}\right).$ (20)

At any time $t$, the instantaneous rate of heat loss/gain $\dot{h}$ by the
system can be obtained as follows:

$\displaystyle\dot{h}=q\dot{q}+p\dot{p}=-\zeta\left(\alpha p^{2}+\frac{\beta
p^{4}}{T}\right).$ (21)

To fix the kinetic energy, the time average of the square of thermostatted
variable, $\bar{\zeta}^{2}$ needs to remain a constant. This means that
$\dot{\zeta}$ must vanish at all times and we can write:

$\displaystyle\overline{\zeta\dot{\zeta}}=\overline{\zeta\alpha\left(\frac{p^{2}}{T}-1\right)+\zeta\beta\left(\frac{p^{4}}{T^{2}}-\frac{3p^{2}}{T}\right)}$
$\displaystyle=0.$

We can rearrange this equation to find:

$\displaystyle\overline{-\zeta\left(\alpha+\frac{3\beta p^{2}}{T}\right)}$
$\displaystyle=\overline{-\frac{\zeta}{T}\left(\alpha p^{2}+\frac{\beta
p^{4}}{T}\right)}$ $\displaystyle\implies\overline{\Lambda}$
$\displaystyle=\overline{\frac{\dot{h}}{T}},$ (22)

from Eqs. 20 and 21. We now take ensemble average both sides of Eq. 22 and
identify $\langle\Lambda\rangle$ as the Gibbs entropy rate $\dot{S}_{G}$ to
find the following identities:

$\displaystyle\frac{\dot{S}_{e}}{k_{B}}=\frac{\dot{S}_{G}}{k_{B}}=-\langle\Lambda\rangle=\frac{\dot{\mathcal{Q}}}{T},$
(23)

where $\dot{\mathcal{Q}}=\langle\dot{h}\rangle$ is the rate of average
exchange of heat of the system with its environment. These ensemble average
can be replaced by time averages as the system is ergodic. The Gibbs entropy
rate gives the thermodynamic entropy production rate in thermostatted
Hamiltonian systems of Nosé-Hoover and Gaussian types.

The non-equilibrium steady state (NESS) value of $\dot{S}_{e}/k_{B}$ increases
with $\epsilon$ in the interval $[0.0,0.5]$ as shown in Fig. SM2.