# Profiling Visual Dynamic Complexity Using a Bio-Robotic Approach

Qinbing Fu1,2, Tian Liu2, Xuelong Sun2, Huatian Wang2, Jigen Peng1,3, Shigang
Yue1,2, Cheng Hu1,2 1Machine Life and Intelligence Research Centre, Guangzhou
University, Guangzhou, 510006, China 2School of Computer Science, University
of Lincoln, Lincoln, LN6 7TS, United Kingdom 3School of Mathematics and
Information Science, Guangzhou University, Guangzhou, 510006, China Q. Fu, T.
Liu, and X. Sun are joint first authors. Correspondence:
{chu,syue}<EMAIL_ADDRESS>
###### Abstract

Visual dynamic complexity is a ubiquitous, hidden attribute of the visual
world that every dynamic vision system is faced with. However, it is implicit
and intractable which has never been quantitatively described due to the
difficulty in defending temporal features correlated to spatial image
complexity. To fill this vacancy, we propose a novel bio-robotic approach to
profile visual dynamic complexity which can be used as a new metric. Here we
apply a state-of-the-art brain-inspired motion detection neural network model
to explicitly profile such complexity associated with spatial-temporal
frequency (SF-TF) of visual scene. This model is for the first time
implemented in an autonomous micro-mobile robot which navigates freely in an
arena with visual walls displaying moving sine-wave grating or cluttered
natural scene. The neural dynamic response can make reasonable prediction on
surrounding complexity since it can be mapped monotonically to varying SF-TF
of visual scene. The experiments show this approach is flexible to different
visual scenes for profiling the dynamic complexity. We also use this metric as
a predictor to investigate the boundary of another collision detection visual
system performing in changing environment with increasing dynamic complexity.
This research demonstrates a new paradigm of using biologically plausible
visual processing scheme to estimate dynamic complexity of visual scene from
both spatial and temporal perspectives, which could be beneficial to
predicting input complexity when evaluating dynamic vision systems.

##  I Introduction

When assessing dynamic vision systems, it would be more convincing to carry
out the research progressively by using a range of input visual stimuli
signals from simple, normally synthesised with clean backgrounds to
handcrafted samples with more complex backgrounds like natural cluttered
visual scene, finally to unstructured indoor and outdoor real world scenarios.
The moving object and background are all controllable in handcrafted testing
visual scene; in this case, the visual system models can be systematically
investigated as desired. The real world visual scenes are nevertheless
intractable; the models usually perform more randomly.

Many visual system modelling studies insisted on working effectively and
robustly in complex and dynamic environments. However, none of those have
quantitatively estimated the complexity of dynamic visual scene leaving the
said “dynamic-complex” implicit and fully empirical. Accordingly, it could be
hard to define or find the boundary of visual system models challenged by
unstructured, real physical scene. Under constrained parameters, the models
could also be overestimated without quantitative analysis on dynamic
complexity of visual scene.

In this regard, we reflected on how to estimate the visual dynamic complexity
with quantitative description. Image processing researchers have presented
some prominent methods on evaluating the spatial image complexity. Yu and
Winkler reviewed relevant methods including compression-based complexity,
spatial information indicated by edge energy [1]. Some methods have been
requesting observers to rate the perceived complexity of images [2, 3, 4].
There have been complexity-based similarity measures used to understand
content-based image recognition problems [5], cluster and classify images [6],
aesthetic pictures [7]. Briefly speaking, the definition of image complexity
is still in debate, as most methods are based on subjective measures that are
inconsistent. On the aspect of objective measures, Corchs et al. proposed a
new complexity measurement by linearly combining different image features
including spatial, frequency and colour properties [8]. More recently, there
were also methods focusing on quantifying the complexity of black-and-white
images based on the detection of typical length scales [9], as well as
modelling visual complexity via image features fusion of multiple kernels
[10].

However, to the best of our knowledge, none of these methods have analysed
visual complexity from a temporal perspective. A potential challenge is how to
defend correlated features in both space and time scales. To this end,
biological visual neural systems are capable of measuring then adapting to
changing environment very quickly and robustly. Moreover, the neural circuits
balance well performance and power consumption where some neural activities
could be flexible to rise up in relatively more complex dynamic environment
whilst resting in easy conditions. The underlying mechanisms on measuring
visual dynamic complexity for such self-adaptation have not yet been
understood.

Neuroscientists have found the spatial-temporal frequency (SF-TF) of movements
are important attributes of dynamic visual scene [11]. Visual neurons show
selectivity to different frequencies corresponding to distinct responsive
activities [12, 13]. Besides that, a vast majority of physiological works
apply sine-wave gratings as input stimuli to neural circuits and pathways
[14], as the gratings can precisely reflect the attributes of visual scene in
both SF-TF domains. Accordingly, visual sensory systems are likely encoding
such neural dynamic characteristics for complexity measurement.

Recently, a brain-inspired motion detection neural network model was proposed
which can decode spatiotemporal attributes of visual scene into interpretable
angular velocity, largely independent of spatial frequency and contrast
information [15]. This fits well with our objective to further develop a
metric to involve also temporal features coupled with spatial image
complexity. Firstly, we investigate the relation between model response and
visual dynamic complexity associated with both SF-TF of visual scene.
Importantly, we have found a mathematically monotonic relation indicating that
higher complexity (SF or TF) of visual scene brings about stronger model
dynamic response. This can be established as a new metric to estimate visual
complexity from both spatial and temporal perspectives. After that, we
implement this approach into embedded vision of an autonomous micro-mobile
robot. During visual navigation in an arena encompassed by visually dynamic
walls displaying either grating or drifting natural scene at varying
frequencies, such complexity can be profiled by on-line visual processing,
effectively and efficiently, in different situations. Finally, we demonstrate
a practical use of this new metric as a predictor to find the boundary of a
collision detection visual system in changing environment with increasing
complexity. The bio-robotic experiments demonstrate this biologically
plausible metric is feasible and computing efficient to profile visual dynamic
complexity which could be beneficial to evaluating other dynamic vision
systems as a predictor of input complexity. This approach could also be
scalable to a variety of robotic visual tasks to measure changing environment.

The rest of this paper is organised as follows. Section II presents the
methods and materials. Section III shows the metric validation through off-
line simulation and on-line robot tests. Section IV exhibits the profiled
complexity through robot visual navigation. Section V demonstrates a valid use
of this new metric. Section VI concludes this paper.

##  II Methods and Materials

###  II-A Model

This proposed metric is established by a biologically plausible visual
processing scheme implemented in a vision-based micro-mobile robot. First of
all, the dynamic vision system is a recently published brain-inspired motion
detection neural network model called “angular velocity decoding model” (AVDM)
mimicking honeybee’s motion vision neural pathways to transform spatiotemporal
image features into interpretable visual angular velocity [15], as shown in
Fig. 1. Previously, the AVDM was only tested by handcrafted visual stimuli. In
this research, it is for the first time implemented in robot vision. The
emphasis herein is laid on further exploring its efficacy in robot visual
navigation, and potential on profiling dynamic complexity. There have been
many similar visual system models that work effectively to encode diverse
motion cues [16]. The motivation to select the AVDM to estimate visual dynamic
complexity is its large independence on spatial frequency and image contrast
compared to the current repository.

To be more specific, the AVDM mainly consists of two separate visual neural
pathways. The texture pathway measures image complexity in space only, which
works effectively to deal with SF and image contrast. The motion pathway
focuses on TF to correlate features in time. Here we do not repeat the
algorithms in the two visual pathways, but emphasise the last output “decoding
layer” which is crucial to estimate the dynamic complexity associated with SF-
TF of visual scene. Given the outputs from texture ($\lambda$) and motion
($R$) pathways, the AVDM’s initial idea is to decode the angular velocity
($\omega$) from the response $R(\omega,\lambda)$ using an approximation
method. Consequently, the actual angular velocity can be approximated
combining a power function of texture pathway output, and a square root
function of motion pathway response, as follows:

$\omega\approx\hat{a}\lambda^{\hat{b}}\frac{1+\hat{C}}{2\hat{C}}\sqrt{R},$ (1)

where $\hat{C}$ is the estimated spatial contrast. Parameters $\hat{a}$ and
$\hat{b}$ can be learnt by minimizing the difference from ground truth,
iteratively as:

$(\hat{a},\hat{b})=\text{arg}\underset{a,b}{\text{min}}(\omega-a\lambda^{b}\frac{1+C}{2C}\sqrt{R(\omega,\lambda)}).$
(2)

Figure 1:  Illustration of the AVDM model adapted from [15] for profiling
visual dynamic complexity in this research. The neural network model mainly
consists of two pathways, the motion and the texture pathways to encode
spatiotemporal features. At the output layer, the coupled signals are
interpreted to indicate visual angular velocity of movements.  Figure 2:
Illustration of the robot arena with visual walls displaying different dynamic
visual scenes including grating and cluttered natural scenes. The micro-robot
is used for profiling the visual dynamic complexity.

(a) moving grating stimuli in time

(b) SF=1upc, dynamic TF

(c) TF=1Hz, dynamic SF

(d) SF=10upc, dynamic TF

(e) TF=5Hz, dynamic SF

(f) SF=50upc, dynamic TF

(g) TF=50Hz, dynamic SF

Figure 3:  Dynamic response of AVDM challenged by different SF-TF gratings.
Here the SF metric is short for “units per circle”.

(a) three camera snapshots during approaching

(b) grating wall TF=0Hz

(c) grating wall TF=1Hz

(d) grating wall TF=2Hz

(e) grating wall TF=3Hz

Figure 4:  Dynamic response of AVDM during approaching the grating wall: each
process is repeated ten times (variance in grey shadow).

(a) profile grating visual scene

(b) profile cluttered natural scene

Figure 5:  Arena settings (from top-down view) for profiling visual dynamic
complexity in (a) grating, (b) cluttered natural scenes.

###  II-B Configuration

Next, we introduce the configuration of model, robot, and arena. The
parameters setting of the AVDM on embedded vision are consistent with [15]. As
exhibited in Fig. 2, the micro-mobile robot is called “Colias-IV”, which is an
autonomous, vision-based, multi-sensor mobile platform [17]. It has been
successfully employed for research into swarm robotics and social insects
[18], and widely utilised in bio-inspired vision systems study [19]. The
detailed robot configuration can be found in [17]. Note that the robot motion
and control strategy is out of the scope of this study. In all types of tests,
we set the forward linear speed constant at approximately 0.08m/s. The turning
angle is 80$\sim$100 degrees randomly to right or left. The sampling frequency
of robot for image processing is maintained at around 33Hz.

The arena used in this research has been upgraded with controllable, LED-
displaying walls in comparison with our previous arena experimentation with
static, patterned walls [20]. As shown in Fig. 2, with a top-down facing
camera system localising the ID-specific pattern on top of robot in real time,
we can obtain on-line data consisting of time, position. Moreover, these can
be aligned with the AVDM response to be sent to hosting PC with a newly
developed communication system. All of these components form the pre-requisite
of this research.

(a) TF=0Hz, duration=30min

(b) TF=0Hz, AVDM$>$100

(c) TF=1Hz, duration=30min

(d) TF=1Hz, AVDM$>$100

(e) TF=3Hz, duration=30min

(f) TF=3Hz, AVDM$>$100

Figure 6:  The profiled dynamic complexity of grating visual scene by robot
overtime navigation. The TF of grating scene was is at 0, 1, 3Hz,
respectively. The larger circles indicate the stronger model dynamic response
with respect to time, i.e., the higher dynamic complexity of visual scene. In
the distribution histograms, the AVDM responses over 100 are counted.

(a) TF=0Hz, duration=30min

(b) TF=0Hz, AVDM$>$100

(c) TF=3Hz, duration=30min

(d) TF=3Hz, AVDM$>$100

Figure 7:  The profiled dynamic complexity of cluttered natural scene by robot
overtime navigation. The TF of natural scene is set at 0 and 3Hz.

##  III Metric Validation

The first step is to verify our hypothesis that such a biologically plausible
visual processing scheme can reflect dynamic complexity in a quantitative
manner. The complexity here is associated with SF-TF of visual scene. We carry
out off-line simulations with input signal stimuli covering a wide range of
SF-TF frequency sine-wave gratings (see examples in Fig. 3a). Based on the
Nyquist’s law (sampling theorem), our sampling frequency of off-line video
sequence is set at 300Hz; the investigated temporal frequency in this type of
experiments is set below 150Hz. The resolution of videos is set at
100$\times$100 in pixels.

(a) TF=0Hz, time=1$\sim$30min

(b) TF=3Hz, time=31$\sim$60min

(c) TF=3Hz, time=61$\sim$90min, model adjusted

(d) time=1$\sim$30, 229 avoidances, SR$>$90%

(e) time=31$\sim$60, 432 avoidances, SR$>$90%

(f) time=61$\sim$90, 233 avoidances, SR$>$90%

Figure 8:  Collision detection-and-avoidance density maps generated by robot
visual navigation. In (a),(b),(c), green circles indicate the avoidance
positions whilst red circles denote the crashes on wall by failing detection.
Note that the proposed metric works as a predictor of surrounding visual
dynamic complexity. After perceiving such complexity, the tested model can be
adjusted to leverage performance and energy consumption.

Fig. 3 compares the results of AVDM response by changing the ST or TF in a
fairly wide range. The AVDM’s dynamic response shows perfect monotonic
increasing when the SF grows at each certain TF. When fixing the SF and
increasing the TF, such a monotonic mapping generally aligns the former case
though the dynamic response shrinks at some tested frequencies. As a result, a
new metric herein can be established to describe the dynamic complexity in the
context, that is, higher-frequency (either SF or TF) movements likely elicit
stronger model dynamic responses. This can be used as a reasonable indicator
or predictor for such an implicit attribute of visual scene.

To verify its scalability and efficacy, we design open-loop robot tests in
real physical scene. To this end, the robot is initiated to move toward the
visual wall displaying moving gratings with certain TF at 0, 1, 2, 3Hz,
respectively, which resemble the circumstances in Fig. 3c, 3e, and 3g. It can
be clearly seen from Fig. 4 the model dynamic response also generally
represents the monotonic mapping to the contextual complexity despite
fluctuations challenged by real physical grating wall. The robot evaluation
demonstrates the proposed metric works consistently in real world visual
scene.

##  IV Profiling Dynamic Complexity

In this section, we profile the visual dynamic complexity. We carry on closed-
loop experiments via robot visual navigation. Two types of on-line tests have
been involved. Fig. 2 and 5 show the arena settings including grating and
cluttered natural visual scenes. The robot motion is bound within the field of
white-light ring displayed by the TV at bottom, near the surround visual
walls. Technically speaking, this strategy is used for invoking the
aforementioned “turning response”. The light ring can be picked up by the
robot’s light sensors assembled at bottom to turn around, in order to avoid
any human interventions in visual navigation.

In the first type of grating scene tests, Fig. 6 exhibits the profiled
complexity distributed in the arena, against different TF of grating wall.
Intuitively, Fig. 6a, 6c, and 6e all articulate that the model dynamic
response is more intensely aroused near the edges of arena indicating higher
dynamic complexity in the context. The results are consistent with the
previous simulation and robot open-loop tests. Furthermore, Fig. 6b, 6d, and
6f compares the magnitude of model dynamic responses. Interestingly, the
profiled complexity becomes fairly higher by increasing the TF of drifting
grating scene, as the neural activity is stronger as illustrated in the
corresponding distribution.

In the second type of natural scene tests, the complexity of more realistic
visual scene is explicitly profiled. Fig. 7 compares the nature scene
complexity against two TF of drifting visual scene. The results can be aligned
with the former grating scene tests which demonstrate the proposed method
could be scalable to more challenging, realistic environment. Accordingly, the
visual dynamic complexity has been quantitatively estimated and explicitly
described.

##  V Using This Metric

In this section, we demonstrate a case study on using this metric. As
introduced in Section I, dynamic complexity is a hidden attribute of visual
scene that is very difficult to be described explicitly, and could affect
accurate assessment of dynamic vision systems especially in real world scenes.
The current literature lacks effective methods from both spatial and temporal
perspectives. Therefore, the visual system models are likely overestimated
with limited testing data described ambiguously as “complex dynamic visual
scene”.

Through this research, we have put forward a new metric which can alleviate
this situation. For example, the boundary of visual system models can be
plainly found by quantitatively increasing such complexity. Using this metric,
here we investigate a collision detection visual system recently proposed by
Fu et al. [20]. The robot in this study is applied to implement this model
interacting with surrounding visual walls for collision detection and
avoidance in navigation. Note that the motion control is not the focus of this
paper. Accordingly, the robot avoidance behaviour is set simply to turn
randomly to right or left at 80 $\sim$ 100 degrees. With respect to previous
research experience, two classic metrics on assessing collision detection
visual systems have been considered: (1) The success rate (SR) of collision
detection-and-avoidance indicates the system robustness which should be
maintained at an acceptable level ($>$90% with regard to recent studies using
the micro-robot [20]). (2) The distance to collision (DTC) means how far the
potential risk is detected and avoided which indicates the system sensitivity;
an appropriate DTC can reduce energy consumption in confined space that is
indicated by the amount of avoidances in an identical time window in this case
study.

This experiment involves three procedures. Fig. 8 illustrates the robot crash
and avoidance events with density maps in the arena. Initially, the TF of
grating wall is set at 0Hz, most of the successful avoidances happened near
the peripheries of arena (Fig. 8a and 8d). The DTC is relatively low and the
SR is maintained above the satisfactory level with 229 avoidances in total.
Under the same set of model parameters, we increase the complexity of visual
scene by setting the TF at 3Hz. Interestingly, Fig. 8b and 8e show obviously
that the avoidance event density is much more centralised toward the centre
field of arena, in contrast with the performance in relatively simpler visual
scene. The tested model becomes more sensitive. In this case, the DTC
significantly extends with more energy consumption. Accordingly, a boundary of
this tested model has been found in visual scene with increased dynamic
complexity.

Such a metric works effectively as a predictor of surrounding input complexity
to visual systems which can explain, and be coupled with the different
behaviours. After perceiving such complexity, a next interesting attempt can
be finding an either learning or experience based adjustment to make the
tested visual system adapt to the new environment with performance and power
consumption leveraged. Fig. 8c and 8f together show a satisfactory balance by
simply increasing the threshold for collision detection. To sum up, this case
study reveals a promising way of using the proposed metric to inform the agent
surrounding dynamic complexity to make sensible adjustment or choice.

##  VI Conclusion

This paper has presented a novel bio-robotic approach to estimate and profile
visual dynamic complexity that is essentially a concealed, intractable
attribute of visual scene. This approach has integrated biologically plausible
visual processing scheme and robotic techniques to set up a metric describing
such complexity explicitly and quantitatively, from both spatial and temporal
perspectives. We have progressed this research with several steps. Firstly, we
have shown a method to fill the gap of measuring complexity associated with
spatial and temporal frequency attributes of visual scene. To establish this
metric, we have then validated it through off-line simulation and on-line bio-
robot tests. Importantly, we have found a monotonic mapping from the neural
dynamic response to a wide range of frequencies of drifting visual scene.
Next, we have applied an autonomous micro-mobile robot to profile the dynamic
complexity in visual navigation to describe surrounding input complexity
within an arena with visual walls displaying drifting visual scene. This
approach has also been verified to be flexible and effective in more
realistic, natural visual scene. Finally, we have given a case study of using
this metric as a predictor of surrounding visual dynamic complexity to help
find the boundary of a collision detection visual system in more complex
dynamic environment.

Our future works are 1) consolidating this metric with improvement in dynamic
visual processing formula to cope with various realistic natural scenes; 2)
coupling this predictor of input complexity with learning methods to benefit
other dynamic vision systems; 3) introducing this work into relevant robot
visual tasks to help the robot measure, and adapt to changing environment.

## Acknowledgement

This research has received funding from the National Natural Science
Foundation of China under the Grant No 12031003, the China Postdoctoral
Science Foundation Grant 2019M662837, 2020M682651, and the European Union’s
Horizon 2020 research and innovation programme under the Marie Sklodowska-
Curie Grant Agreement No 778602 ULTRACEPT.

## Appendix

We have a supplementary video demo with this paper for the description of
methods.

## References

  * [1] H. Yu and S. Winkler, “Image complexity and spatial information,” in _2013 Fifth International Workshop on Quality of Multimedia Experience_ , 2013, pp. 12–17.
  * [2] H. Wu, M. Claypool, and R. Kinicki, “A study of video motion and scene complexity,” in _Tech. Rep. WPI-CS-TR06-19_ , 2006.
  * [3] V. Chikhman, V. Bondarko, M. Danilova, A. Goluzina, and Y. Shelepin, “Complexity of images: Experimental and computational estimates compared,” _Perception_ , vol. 41, no. 6, pp. 631–647, 2012.
  * [4] M. I. Chacon-Murguia, A. D. Corral-Saenz, and R. Sandoval-Rodriguez, “Image complexity measure: A human criterion free approach,” in _Proc. Annual Meeting of the North American Fuzzy Information Processing Society (NAFIPS)_ , 2005, pp. 241–246.
  * [5] J. Perkio and A. Hyvarinen, “Modelling image complexity by independent component analysis, with application to content-based image retrieval,” in _Proc. 19th International Conference on Artificial Neural Networks (ICANN)_ , 2009, pp. 704–714.
  * [6] T. Guha and R. K. Ward, “On image similarity, sparse representation and kolmogorov complexity,” _arXiv Computing Research Repository_ , vol. abs/1206.2627, 2012.
  * [7] J. Romero, P. Machado, A. Carballal, and A. Santos, “Using complexity estimates in aesthetic image classification,” _J. Mathematics and the Arts_ , vol. 6, no. 2-3, pp. 125–136, 2012.
  * [8] S. E. Corchs, G. Ciocca, E. Bricolo, and F. Gasparini, “Predicting complexity perception of real world images,” _PLOS ONE_ , vol. 11, no. 6, 2016.
  * [9] D. H. Zanette, “Quantifying the complexity of black-and-white images,” _PLOS ONE_ , vol. 13, no. 11, 2018.
  * [10] C. Fernandez-Lozano, A. Carballal, P. Machado, A. Santos, and J. Romero, “Visual complexity modelling based on image features fusion of multiple kernels,” _PeerJ_ , 2019.
  * [11] D. Shoham, M. Hubener, S. Schulze, A. Grinvald, and T. Bonhoeffer, “Spatio-temporal frequency domains and their relation to cytochrome oxidase staining in cat visual cortex,” _Nature_ , vol. 385, no. 6, pp. 529–533, 1997.
  * [12] K. H. Foster, J. P. Gaska, M. Nagler, and D. A. Polien, “Spatial and temporal frequency selectivity of neurones in visual cortical areas V1 and V2 of the macaque monkey,” _Journal of Physiology_ , vol. 365, pp. 331–363, 1985\.
  * [13] N. J. Priebe, S. G. Lisberger, and J. A. Movshon, “Tuning for spatiotemporal frequency and speed in directionally selective neurons of macaque striate cortex,” _Journal of Neuroscience_ , vol. 26, no. 11, pp. 2941–2950, 2006\.
  * [14] A. Borst and T. Euler, “Seeing things in motion: Models, circuits, and mechanisms,” _Neuron_ , vol. 71, no. 6, pp. 974–994, 2011.
  * [15] H. Wang, Q. Fu, H. Wang, P. Baxter, J. Peng, and S. Yue, “A bioinspired angular velocity decoding neural network model for visually guided flights,” _Neural Networks_ , vol. 136, pp. 180–193, 2020.
  * [16] Q. Fu, H. Wang, C. Hu, and S. Yue, “Towards computational models and applications of insect visual systems for motion perception: A review,” _Artificial Life_ , vol. 25, no. 3, pp. 263–311, 2019.
  * [17] C. Hu, Q. Fu, and S. Yue, “Colias IV: The affordable micro robot platform with bio-inspired vision,” in _Towards autonomous robotic systems conference_ , 2018, pp. 197–208.
  * [18] X. Sun, T. Liu, C. Hu, Q. Fu, and S. Yue, “ColCOS $\phi$: A multiple pheromone communication system for swarm robotics and social insects research,” in _Proceedings of the 2019 IEEE 4th International Conference on Advanced Robotics and Mechatronics (ICARM)_. IEEE, 2019, pp. 59–66.
  * [19] Q. Fu, C. Hu, P. Liu, and S. Yue, “Towards computational models of insect motion detectors for robot vision,” in _Towards autonomous robotic systems conference_ , 2018, pp. 465–467.
  * [20] Q. Fu, H. Wang, J. Peng, and S. Yue, “Improved collision perception neuronal system model with adaptive inhibition mechanism and evolutionary learning,” _IEEE Access_ , vol. 8, pp. 108 896–108 912, 2020.


# A Nernst heat theorem for nonequilibrium jump processes

Faezeh Khodabandehlou  Instituut voor Theoretische Fysica, KU Leuven, Belgium
<EMAIL_ADDRESS>Christian Maes  Instituut voor Theoretische
Fysica, KU Leuven, Belgium  Karel Netočný  Institute of Physics, Czech Academy
of Sciences, Prague, Czech Republic

###### Abstract

We discuss via general arguments and examples when and why the steady
nonequilibrium heat capacity vanishes with temperature. The framework is the
one of Markov jump processes on finite connected graphs where the condition of
local detailed balance allows to identify the heat fluxes, and where the
discreteness more easily enables sufficient nondegeneracy of the stationary
distribution at absolute zero, as under equilibrium. However, for the
nonequilibrium extension of the Third Law, a dynamic condition is needed as
well: the low-temperature dynamical activity and accessibility of the dominant
state must remain sufficiently high so that relaxation times do not start to
dramatically differ between different initial states. It suffices in fact that
the relaxation times do not exceed the dissipation time.

nonequilibrium jump process, excess heat, Nernst postulate

##  I Introduction

The Third Law of Thermodynamics is a cornerstone of physical and chemical
thermodynamics cal ; fer . It sets and constrains the scene of low-temperature
phenomena, at least where it concerns thermal properties of equilibrium
systems. To witness the importance, the quantum revolution was much inspired
by the difficulty to account in classical physics for the low-temperature
behavior of specific heats kle . The Nernst postulate deals explicitly with
these questions, arguing that, at zero temperature, the equilibrium entropy
becomes a constant independent of any external parameter, ner . Combined with
Boltzmann’s interpretation of entropy, Planck’s version states that this
entropy constant is zero, pla . It provides a fixed reference point for the
absolute entropy of any substance at any temperature. Further conditions imply
that the heat capacity must vanish with temperature cal ; mas  and the state
of absolute zero is physically unattainable via a finite number of equilibrium
steps. However, in contrast with the First and Second Law, the Third Law is
not always obeyed by substances and corresponding Nernst theorems, therefore,
require nontrivial conditions, lie1 ; lie2 .

Since a long time, it has been considered important to develop calorimetry and
to find a similar general structure in the study of thermodynamic processes
connecting steady states of nonequilibrium systems, pri ; oon ; kom2 . The
present work elaborates on such a general and nonperturbative extension of the
Third Law, valid at an arbitrary distance from equilibrium and operationally
formulated using a physically significant quasipotential from which heat
capacities and latent heats can be defined. We explain a dynamical condition
for those heats to vanish at absolute zero, requiring residual dynamical
activity, which, when violated, enables the possibility of a zero–temperature
phase transition.
As is well-known, general results on thermal properties of nonequilibrium
systems (steady or transient for that matter) are very few, with the
positivity of the entropy production, steady and transient fluctuation
theorems, thermodynamic and kinetic uncertainty relations, and fluctuation-
response relations being the most important ones; even there however the zero-
temperature limit is seldom discussed (see sor  for an interesting exception).

The point is not at all to formally extend the definition of entropy scr , and
we do not need a Clausius heat theorem out-of-equilibrium. Rather, we wish to
obtain quantitative measures of the heat released or absorbed during a
quasistatic relaxation to a (new) nonequilibrium condition because of (small)
changes in control parameters. That also applies to small systems and need not
be restricted to the context of linear irreversible thermodynamics or to
close-to-equilibrium processes. Far-from-equilibrium examples within
electronics are micro- and nano-scale devices as digital logic circuits with a
rich collection of nonlinear dissipative effects due to quantum tunneling,
current leakage, etc. pop . Other examples include driven hopping dynamics on
a lattice in contact with an electron- or photon-bath, where the driving can
be electromagnetic or optical, dls ; qa . As a third class of examples and
motivation, we can turn to thermochemistry. Under detailed balance, the energy
exchanges during chemical reactions are measured with the methods of
calorimetry, and molar heat capacities inform about the landscape of
thermodynamic potentials. Thus, at (very) low ambient temperatures, cryogenic
chemistry mostly focuses on the ability of molecules to lower the total
energy. We remember that Nernst got the 1920 Nobel Prize in chemistry “in
recognition of his work in thermochemistry.” However, under nonequilibrium
conditions where detailed balance is violated, kinetic aspects will matter
greatly and it is essential to understand the constraints set by an extended
Third Law. While, in the above examples, the driving can be time-dependent
(like when applying an AC-field or periodic feeding), here we imagine time-
scales where the external field can be taken almost constant. The extension to
truly time-dependent external driving is conceptually not different though,
and we give an example of random switching in Section III.3.

In driven systems, the energy exchange always runs against a dissipation
background. In the stationary situation, work done by nonconservative forces
is constantly dissipated as (Joule) heat; see the cartoon in Fig.1 (left). Or,
when in contact with different chemical reservoirs such as in chemostats and
bioreactors, particle currents are maintained between external sources and
sinks at some temperature. Changing an external parameter, such as that
ambient temperature, creates excess heat in the relaxation to the new steady
condition. Then, our interest may go to that excess part coming from changes
in the temperature and/or other parameters defining the system dynamics HS ;
GP ; cla . However, how the total heat exchange is split into the “background”
and the “excess” is not unique and various schemes exist exhibiting different
features. One such a scheme is known as the Hatano-Sasa approach HS , where
the total heat is decomposed into the “house-keeping” and the “excess”
components in such a way that the excess heat exactly satisfies the (extended)
Clausius heat theorem. Other variants of this approach can be found in ber ;
cla ; GP ; st . Yet, measuring that excess heat is not obvious. Another
scheme, which goes back to the original proposal by Oono and Paniconi oon ,
considers the excess heat defined by subtracting the time-integrated steady
state heat flux kom2 , sometimes called the “dynamical component” of heat wang
. While such an excess heat does not satisfy the Clausius relation unless
close to equilibrium, it enjoys the important properties of being both
“geometrical” (= dependent only on the thermodynamic trajectory) and
“reversible” (= antisymmetric with respect to protocol-reversal) in the
quasistatic limit. The mathematical analysis of that limit is obtained in
arxiv . Here we take the result in terms of the so-called quasipotential and
since that is also operationally accessible (e.g., by an AC-measurement with
low-frequency temperature variations), see e.g. the cartoon of Fig. 1 (right),
that becomes a preferential option for the calorimetric analysis of the
present paper. It allows a natural definition of the nonequilibrium heat
capacity, eu ; jir ; calo . The main result then is that for a large and
physically motivated class of nonequilibrium systems modeled as Markov jump
processes, that heat capacity vanishes at absolute zero. We explain the
context, significance, and main ideas in the present paper, also adding simple
examples.

In the next section, we review the main elements of nonequilibrium
calorimetry. Section III describes the framework of Markov jump processes and
how calorimetry is set up in the absence of (global) detailed balance.
Section IV is devoted to stating and explaining the extended Nernst postulate.
The simplest sufficient condition for the validity of the extended Third Law
is that relaxation times remain smaller than the dissipation time at low
temperatures. Various illustrations are added in Section V, showing that our
conditions are natural and not just sufficient. Sections VI and VII provide
extra background and discussion for better interpretation and significance of
the results. The conclusion section VIII gives a short summary.

We note that the present paper focuses on the underlying ideas and intuition.
We not only skip the detailed analysis of the quasistatic limit, but we also
present the result only for the heat capacity. In all it makes an essential
complement to the mathematical proof which can be found in arxiv .

##  II Nonequilibrium calorimetry

We collect in the present Section a review of the general thermal physics that
leads to the notion of a steady nonequilibrium heat capacity. As such and
independent of a specific mathematical setup or model, it has not appeared
before but the main ingredients are also contained in the papers eu ; jir ;
calo ; ner . We note that other definitions of nonequilibrium heat capacity
have been proposed, see e.g. man ; subas ; dls , such as the temperature-
derivative of the mean energy on which we comment briefly in Section II.3.

###  II.1 Thermal response

Consider an open system that is in a steady nonequilibrium condition,
maintaining stationary currents and a strictly positive mean dissipative
power. A simple example is depicted as a cartoon in the left part of Fig. 1: a
current-carrying resistor creates Joule heating of the surrounding liquid. It
stands of course for many other examples, including systems in nonlinear
regimes where the resistance depends in complicated ways on the bath
temperature $T$ or where the voltage is a nonlinear function of the current.
At any rate, if the bath is sufficiently large and the coupling with the
system is weak, the temperature $T$ can be considered constant and the bath is
and remains in thermal equilibrium. Putting $\beta=T^{-1}$, we often move
between the two notations for (inverse) temperature.

Figure 1: left: Changing the temperature $T$ of the bath changes the
dissipated power in Joule heating. right: In equilibrium, the mean dissipated
power ${\dot{q}}_{\beta}={\dot{q}}_{\beta+\textrm{d}\beta}=0$ vanishes, but
under driving, they may take different positive values. The dissipated power
is variable and fluctuates around its steady value ${\dot{q}}(\beta)$, as
indicated for negative times. At time zero, a condition $x$ is obtained for
the relaxation of the expected heat flux to its stationary value, indicated
for positive times. The shaded area gives the meaning of the quasipotential.

Before time zero, there is mean dissipated power $\dot{q}_{\beta}\geq 0$.
Imagine next that very slowly (in the so-called quasistatic limit), the
temperature of the bath is changed to a new value. Such a procedure requires
small steps in temperature, at a rate that is much slower than the inverse
relaxation times. To keep it simple, we consider in Fig. 1 (right) only one
step in temperature: at time $t=0$ the inverse temperature is changed
$\beta+\textrm{d}\beta$ and kept fixed at $\beta+\textrm{d}\beta$. Suppose the
state is $x$ at time zero. After a possibly long time, the system reaches a
new steady nonequilibrium condition, now with power
$\dot{q}_{\beta+\textrm{d}\beta}\geq 0$, independent of $x$. The shaded area
in Fig. 1 is the integral of a difference in power, and is given by the so-
called quasipotential

$V(x,\beta+\textrm{d}\beta)=\int_{0}^{\infty}\textrm{d}t\,\left[{\dot{q}}_{t}(x,\beta+\textrm{d}\beta)-{\dot{q}}_{\beta+\textrm{d}\beta}\right]$
(II.1)

where ${\dot{q}}_{t}(x,\beta+\textrm{d}\beta)$ is the expected dissipated
power at time $t>0$ when starting (at time zero) from state $x$. Note that we
assume here that, at least for all finite $\beta$ and relevant conditions $x$,
${\dot{q}}_{t}(x,\beta+\textrm{d}\beta)\rightarrow{\dot{q}}_{\beta+\textrm{d}\beta}$
sufficiently fast as $t\uparrow\infty$ so that the integral in (II.1) is
converging.

The quasipotential (II.1) measures an excess heat. There is obviously
important literature on such excesses; we already mentioned oon ; kom2 ; wang
, but there is also the notion of path-dependent heat and work, as used by
many people studying nonequilibrium fluctuations, including e.g. ken ; 2law .
Remember now that $x$ in (II.1) and in Fig. 1 (right) is a state which is
statistically selected from the nonequilibrium ensemble which describes the
static fluctuations at time zero when the system is still coupled with the
thermal bath at temperature $T$. The shaded area (II.1) has an average (over
$x$) $\langle V(x,\beta+\textrm{d}\beta)\rangle_{\beta}=O(\textrm{d}\beta)$ of
order $\textrm{d}\beta$, since clearly $\langle V(x,\beta)\rangle_{\beta}=0$
from (II.1). The (nonequilibrium) heat capacity is therefore mathematically
and physically well-defined as

$C(T)=\beta^{2}\,\frac{\langle
V(x,\beta+\textrm{d}\beta)\rangle_{\beta}}{\textrm{d}\beta}=\beta^{2}\,\big{\langle}\,\frac{\textrm{d}V}{\textrm{d}\beta}(x,\beta)\,\big{\rangle}_{\beta}$
(II.2)

measuring the extra heat toward the system per temperature change while we fix
possible other external controls.

###  II.2 Fluctuation expression

The formula (II.2) for the heat capacity has an average over a steady
nonequilibrium ensemble. It is determined by the time-invariant probabilities,
denoted by $\rho(x)=\rho_{\beta}(x)$ where we often ignore the temperature-
dependence in the notation, to find the system in state $x$ when the
equilibrium bath (to which it is weakly coupled) is kept at inverse
temperature $\beta$. Since $\langle
V(x,\beta)\rangle_{\beta}=\sum_{x}\rho_{\beta}(x)V(x,\beta)=0$, we can
continue from (II.2) to write

$C(\beta)=\beta^{2}\sum_{x}\rho_{\beta}(x)\frac{\textrm{d}V}{\textrm{d}\beta}(x,\beta)=-\beta^{2}\sum_{x}\frac{\textrm{d}\rho_{\beta}}{\textrm{d}\beta}(x)\,V(x,\beta)$
(II.3)

where for convenience we assumed a finite number of possible states $x$.
(Otherwise, we would write an integral with $\rho(x)$ as probability density,
or even more general expressions.) As a consequence, after a last rewriting of
(II.3),

$C(\beta)=-\beta^{2}\big{\langle}\frac{\textrm{d}\log\rho_{\beta}(x)}{\textrm{d}\beta}\,\,V(x,\beta)\big{\rangle}_{\beta}$
(II.4)

which is a correlation function between the quasipotential $V(x,\beta)$,
related to excess heat, and the temperature-dependence
$\frac{\textrm{d}}{\textrm{d}\beta}\log\rho(x)$ of the stationary state-
probabilities.

###  II.3 Equilibrium case

The basic idea to define heat capacity $C(\beta)$ as in (II.2) is unchanged
when moving between an equilibrium and a nonequilibrium situation, but
obviously, equilibrium comes with a simplification.
To evaluate (II.1) in thermal equilibrium, we first note that then excess heat
just means heat and $\dot{q}(\beta)=\dot{q}(\beta+\textrm{d}\beta)=0$.
Furthermore, for reversible transformations when no work is done on the
system, the First Law reads
${\dot{q}}_{t}(x,\beta+\textrm{d}\beta)=-\frac{\textrm{d}}{\textrm{d}t}\langle
E(x_{t})\rangle_{\beta+\textrm{d}\beta}$ for $x_{0}=x$, in terms of the energy
function $E$ of the system. Hence, when there is no external (nonequilibrium)
driving, (II.1) becomes
$V(x,\beta+\textrm{d}\beta)=E(x)-U(\beta+\textrm{d}\beta)$, where
$U(\beta)=\langle E(x)\rangle_{\beta}$ is the (average) equilibrium energy. As
a result, (II.2) then gives the well-known expression for the heat capacity at
constant volume:
$C(\beta)=-\beta^{2}\frac{\textrm{d}}{\textrm{d}\beta}U(\beta)$. That formula
no longer holds in steady nonequilibrium where, even at constant volume, there
is the possibility of irreversible work.
For comparison with the fluctuation formula (II.4), we recall that when the
system is in canonical equilibrium with partition function $Z$,

$\displaystyle\frac{\textrm{d}\log\rho}{\textrm{d}\beta}(x)$ $\displaystyle=$
$\displaystyle\beta^{2}\,E(x)+\beta^{2}\frac{\textrm{d}\log
Z}{\textrm{d}\beta}$ $\displaystyle V(x,\beta)$ $\displaystyle=$
$\displaystyle E(x)-U(\beta)$ (II.5)

so that (II.4) reproduces then the equilibrium expression

$C(\beta)=\beta^{2}\,\big{\langle}\left(E-U(\beta)\right)^{2}\big{\rangle}_{\beta}\geq
0$

in terms of the energy variance. However. the nonequilibrium heat capacity and
hence the correlation (II.4) can get negative, as we will illustrate via
Example III.3.

##  III Markov jump processes

To understand the low-temperature properties of the heat capacity defined
above, we need to add a class of models that imitate some of the main aspects
of the quantum world. The present section presents the setup and Section VI
gives extra explanation and interpretation. In particular, the (open) systems
that we model should be thought of as consisting of discrete degrees of
freedom, where the zero-temperature condition is well-separated from its
excitations.
Furthermore, because we explore the nonequilibrium regime, we take a dynamics
that violates the condition of detailed balance.

###  III.1 Setup

When reactions occur far from equilibrium, modeling with Markov jump processes
becomes very relevant to understand their calorimetry, dm ; kkm .
Additionally, there is a large literature on chemical reaction networks with
papers starting from exactly the same setup, lux1 ; lux2 ; lux3 ; lux4 ; viv ;
sch . Such Markov dynamics with discrete states can be obtained as the weak-
coupling limit of a quantum system, which is in incoherent contact with one or
several heat baths; see more in Section VI. We refer to kamp  for more context
on stochastic processes in physics and chemistry.

For a formal description, we consider Markov jump processes on an (arbitrary)
finite, simple, and connected graph $\mathcal{G}$. Vertices $x,y,z,\ldots$
denote the physical states or conditions on a mesoscopic scale of description.
We emphasize that the states $x,y,\dots$ may represent many-body
configurations such as indicating the quantity of a certain substance. The
edges indicate the possible chemomechanical transitions between
configurations, including discrete changes in position or energy. We denote
the transition rates by $k(x,y)>0$ for a jump $x\rightarrow y$, and $k(y,x)>0$
is the rate for a jump $y\rightarrow x$. See Fig. 2 for a summary.

Figure 2: left: graph with vertices $x$ and $y$, connected via an edge.
middle: the transition rates for the directed jumps over the edge. right:
indication of heat given or taken during the transition between the states
$x\rightarrow y$; see (III.4).

The transition rates determine a unique stationary distribution $\rho(x)>0$,
solution of the stationary Master equation: for all states $x$,

$\sum_{y}\big{[}k(y,x)\rho(y)-k(x,y)\rho(x)\big{]}=0$ (III.1)

We also introduce the backward generator $L$ of the Markov process, satisfying

$\displaystyle Lh\,(x)$ $\displaystyle=$
$\displaystyle\sum_{y}k(x,y)[h(y)-h(x)]$ (III.2) $\displaystyle e^{tL}h\,(x)$
$\displaystyle=$ $\displaystyle\langle h(X_{t})\,|\,X_{0}=x\rangle$ (III.3)

for the expected value at time $t\geq 0$ of a function $h$, when starting the
process at $X_{0}=x.$ We refer to textbooks like kamp ; grim  for more
introduction to Markov jump processes.

The rates $k(x,y)$ contain kinetic and thermodynamic information summarizing
the weak coupling of the system to a heat bath at inverse temperature $\beta$.
The rates may depend on external parameters but, for simplicity, we
concentrate here on the dependency of the rates on temperature (often as a
subscript in $k_{\beta}(x,y)$).

###  III.2 Quasipotential and heat capacity

The crucial link between the mathematics of jump processes and their
thermodynamic interpretation is contained in the condition of local detailed
balance ldb ; Kat ; les ,

$\log\frac{k_{\beta}(x,y)}{k_{\beta}(y,x)}=\beta\,q(x,y)$ (III.4)

where $q(x,y)=-q(y,x)$ is required to be the heat dissipated in the heat bath
at inverse temperature $\beta$; see ldb  and references therein (and Fig. 2
for a graphical reminder).
The local detailed balance has nothing to do with the notion of local
equilibrium. Rather, it must be contrasted with global detailed balance, where
all currents vanish. When there is global detailed balance (which we take as
defining equilibrium dynamics), there is an energy state function $E$ for
which that heat in (III.4) equals $q_{\text{eq}}(x,y)=E(x)-E(y)$, for all
edges $(x,y)$.
We also note that local detailed balance should be formulated for each
dissipative transition channel separately, and it is easily broken by taken
further coarse-graining or combining rates; see e.g. ast .

Recall (II.1) and how the quasipotential $V(x,\beta)$ appears in the
expression (II.2) of the nonequilibrium heat capacity. From now on, we write
$V(x)=V(x,\beta)$ if no confusion with temperatures can arise. Then, using the
above formalism we have

$V(x)=\int_{0}^{\infty}\Big{[}\big{\langle}\dot{q}(X_{t})\,|\,X_{0}=x\big{\rangle}-\dot{q}_{\beta}\Big{]}$
(III.5)

where

$\displaystyle\dot{q}(x)$ $\displaystyle=$
$\displaystyle\sum_{y}k_{\beta}(x,y)\,q(x,y)$ (III.6)
$\displaystyle\dot{q}_{\beta}$ $\displaystyle=$
$\displaystyle\sum_{x}\rho_{\beta}(x)\,\dot{q}(x)$ (III.7)

and $\rho_{\beta}(x)$ is the stationary distribution, solution of (III.1). The
above formulæ realize the previous (II.1). In other words, the quasipotential
(III.5) must be used in the heat capacity (II.2),

$C(\beta)=\beta^{2}\,\big{\langle}\frac{\textrm{d}V}{\textrm{d}\beta}\big{\rangle}_{\beta}$
(III.8)

and, similarly, in the fluctuation formula (II.4) for the correlation function

$C(\beta)=-\beta^{2}\,\big{\langle}\frac{\textrm{d}\log\rho}{\textrm{d}\beta}(x)\,V(x)\big{\rangle}_{\beta}$
(III.9)

Note here that the stationary distribution $\rho$, solution of (III.1), is in
general not explicitly known.

###  III.3 Example: molecular switch

When a molecule can reversibly switch between two stable states, we have a
molecular two-level switch. It has various realizations, see e.g. ms ; ms2 ,
and can be engineered according to various protocols of switching. Here we
take a toy-version, without being specific about the molecule, but focusing on
the essential nonequilibrium character and its thermal response.

Consider the reaction scheme $A\rightleftarrows B$ and
$A^{\prime}\rightleftarrows B^{\prime}$, where each time we imagine two
molecular levels separated by energy $\varepsilon$, but where, by a random
environmental stimulus, shifts $A\leftrightarrow A^{\prime}$ and
$B\leftrightarrow B^{\prime}$ occur at a certain rate $\alpha$; see Fig. 3.

Figure 3: Molecular switch as described in (III.10). left: The possible
transitions between the four states. States $A-B$ and states
$A^{\prime}-B^{\prime}$ are separated by an energy $\varepsilon$. The
horizontal transitions correspond to the switches. right: A sample trajectory
with the switching indicated by the change in color. The heat is suggested as
in Fig. 2(right) when the system jumps between energy levels.

(a)

(b)

Figure 4: Heat capacity for the two-level switch defined in (III.10), (a)
without kinetic barrier, (b) with a barrier at $\Delta=2$.

We thus distinguish four states $x\in\\{A,A^{\prime},B,B^{\prime}\\}$. The
transitions $A\leftrightarrow A^{\prime}$ and $B\leftrightarrow B^{\prime}$
occur at fixed rate
$k(A,A^{\prime})=k(A^{\prime},A)=\alpha=k(B^{\prime},B)=k(B,B^{\prime})$ for
$\alpha\geq 0$. The other four transitions are thermal at inverse temperature
$\beta$ with rates

$\displaystyle k(A,B)=k(B^{\prime},A^{\prime})$ $\displaystyle=$
$\displaystyle\nu\,e^{\beta\varepsilon/2}$ (III.10) $\displaystyle
k(A^{\prime},B^{\prime})=k(B,A)$ $\displaystyle=$
$\displaystyle\nu\,e^{-\beta\varepsilon/2}$

The kinetic factor $\nu>0$ is a characteristic frequency for transitions
between both levels. Detailed balance is broken for $\alpha\varepsilon\neq 0$
since there is a difference between the product of rates along two opposite
orientations of a loop: the orientation $A\rightarrow B\rightarrow
B^{\prime}\rightarrow A^{\prime}\rightarrow A$ has
$\alpha^{2}\nu^{2}e^{\beta\varepsilon}$ as product of its rates, whereas the
reversed order makes the product $\alpha^{2}\nu^{2}e^{-\beta\varepsilon}$.

The heat capacity (III.8) can be calculated exactly along the lines of the
previous subsection. Fig. 4(left) gives the heat capacity as a function of
temperature for $\varepsilon=\nu=1$. Observe that whenever $\alpha\neq 0$ the
heat capacity becomes massively enhanced at (very) low temperatures. In Fig. 4
(right), we take $\nu=e^{-\beta\Delta}$ which strongly depends on temperature
with kinetic barrier $\Delta>0$. In Fig. 4 (right), we use $\Delta=2$.
Interestingly, the heat capacity takes both positive and negative values. The
low-temperature behavior is very different both from the (nonequilibrium)
heat-bath dynamics in Fig. 4(left) and from the equilibrium case, while still
$C(T)\rightarrow 0$ as the temperature reaches absolute zero.

##  IV Extended Third Law for Nonequilibria

Even for equilibrium system, the Third Law of Thermodynamics comes with a
condition: the ground state must be nondegenerate to avoid residual entropy.
For the extension of the present paper we need a similar (and first) static
condition, to be discussed in the following subsection. However, for our
nonequilibrium purpose we also have a (second) dynamic condition which is the
subject of Section IV.2. The main issue then for an extended Nernst Postulate
is to explain when and how the nonequilibrium heat capacity (III.8) vanishes
with temperature. More specifically, for (III.8)–(III.9), to have

$C(\beta)\rightarrow 0,\;\text{ with }\,\beta\uparrow\infty$ (IV.1)

The arguments make the contents of the following two subsections. We focus
here on some of the essential ingredients; mathematical details are in arxiv .
All considerations about thermal properties of low-temperature processes must
be compatible with (IV.1) (when the static and dynamic criterion are
satisfied). The illustrations of Section V show how the conditions are exactly
on target.

###  IV.1 Low–temperature asymptotics of static fluctuations

From Section III.1, there is a unique stationary probability distribution
$\rho_{\beta}>0$ satisfying (III.1). At low temperatures, there exists an
approximate expression for $\rho_{\beta}(x)$. Not too surprisingly, it is
determined by the asymptotics of the transition rates,

$\phi(x,y)=\lim_{\beta}\frac{1}{\beta}\log k_{\beta}(x,y)$ (IV.2)

From an analysis of the Kirchhoff formula in heatb ; lowT ; intr , we have
found that for $\beta\uparrow\infty$,

$\rho_{\beta}(x)\simeq e^{-\beta[\phi^{*}-\phi(x)]}$ (IV.3)

with

$\phi^{*}=\max_{z}\phi(z),\quad\phi(z)=\max_{\mathcal{T}}\sum_{(x,y)\in\mathcal{T}_{z}}\phi(x,y)$
(IV.4)

where the maximum is over all spanning trees $\mathcal{T}$ on the graph, and
with $\mathcal{T}_{z}$ the oriented tree obtained from $\mathcal{T}$ by
directing it toward $z$. Any state $x^{*}$ satisfying $\phi^{*}=\phi(x^{*})$
is called dominant; they have the largest probability in the large-$\beta$
asymptotics. We think of dominant states as the nonequilibrium analogue of
equilibrium ground states, but note that they do not need to minimize energy.
The asymptotics (IV.3) has to be understood in the logarithmic sense, meaning
that there may be a prefactor, depending on $x$, which is not exponentially
small in $\beta$.

Actually, (IV.3) is a general result for the Markov jump processes we consider
here, but for the present purpose of establishing an extended Third Law, it
suffices that there is a unique state $x^{*}$ so that

$\rho_{\beta}(x)=\delta_{x,x^{*}}+O(\beta^{-2})$ (IV.5)

as our first criterion. It means that $\rho_{\beta}(x)\rightarrow 0$ for
$x\neq x^{*}$, decays to zero with temperature at least as fast as $T^{2}$. It
implies
$\frac{\textrm{d}}{\textrm{d}\beta}\,\rho_{\beta}^{\text{s}}\rightarrow 0$ as
$\beta\uparrow\infty$. Checking with (III.9), the behavior (IV.5) will send
the heat capacity to zero when the quasipotential $V_{\beta}$ remains bounded
in the limit $\beta\uparrow\infty$. Requiring the behavior (IV.5) is therefore
a first (static) condition ensuring the validity of the extended Third Law.

Note also that the decay rate $\rho_{\beta}(x)\rightarrow 0$ for $x\neq x^{*}$
reveals the support of the zero-temperature distribution, much as in
equilibrium. See (IV.3): if $x^{*}$ is the unique dominant state, i.e., the
unique maximizer of $\phi(z)$ in the sense that all other states are
suppressed by $O(e^{-\tilde{\varepsilon}\beta}$) for some
$\tilde{\varepsilon}>0$, then the convergence in (IV.5) is actually
exponentially fast. If there are multiple dominant states (i.e., maximizers of
$\phi(z)$), then the convergence $\textrm{d}\rho_{\beta}\rightarrow 0$ may be
polynomial in temperature. As a result, the heat capacity goes to zero
exponentially fast in $\beta\uparrow\infty$ in the case of a unique dominant
state and possibly much slower due to the presence of more than one dominant
state.
All that is shared with the equilibrium situation. There is however another
and second condition which is dynamic, and addresses the boundedness of
$V_{\beta}$ (which is never an issue in equilibrium; see Section II.3).

###  IV.2 Boundedness of the quasipotential

We come to the statement and explanation of a second and dynamical condition
for the nonequilibrium version of the Third Law: a sufficient criterion for
the quasipotential $V$ to be uniformly bounded at vanishing temperature.
As is visible from (III.5), boundedness of $V$ as a function of
$\beta\uparrow\infty$ strongly relates with low-temperature relaxation
properties and in particular with the accessibility of the states that support
the zero-temperature condition; we will see that explicitly in Section V.

For the actual argument, which is different from the graphical arguments in
arxiv , we start with a remark. If for all edges $(xy)$ in the graph, the
differences $|V(x,\beta)-V(y,\beta)|$ remain bounded, then also the
quasipotential $V$ must remain bounded, because, under our conditions,
$V(x^{*},\beta)$ cannot be diverging.
Let us, therefore, consider differences of (III.5),

$V(x,\beta)-V(y,\beta)=\int_{0}^{+\infty}\textrm{d}t\,[\,\langle{\dot{q}}(X_{t})\,|\,X_{0}=x\rangle-\langle{\dot{q}}(Y_{t})\,|\,X_{0}=y\rangle\,]$
(IV.6)

where we indicated with $X_{t}$, respectively $Y_{t}$, random walkers starting
from initial position $x$, respectively $y$. We couple those two walkers,
i.e., we consider the process of coupled random walkers $(X_{t},Y_{t})$ with
$(X_{0},Y_{0})=(x,y)$ on the double graph. In an optimal coupling, once they
meet for the first time, $X_{t}=Y_{t}$ at a random time
$t=t_{M}(x,x^{\prime};\beta)$, they keep on taking the same path. Surely, for
every finite $\beta$ there is such a first time $t_{M}(x,y;\beta)$, and
$X_{t}=Y_{t}$ for all $t\geq t_{M}(x,y;\beta)$. Let $\tau_{M}=\langle
t_{M}(x,y;\beta)\rangle$ be the average first meeting time. Then, an upper
bound on the integral (IV.6) is given by

$|V_{\beta}(x)-V_{\beta}(y)|\leq
R\,\tau_{M}\,\max_{x^{\prime},y^{\prime}}|{\dot{q}}(x^{\prime})-{\dot{q}}(y^{\prime})|$
(IV.7)

where $R$ is the linear size of the graph (maximal graph distance between any
two vertices). The time $\tau_{M}$ is bounded by the relaxation time $\tau$.
On the other hand, the expected difference in the power $\dot{q}$, while
fluctuating, can be estimated, at least far away from equilibrium, by a
quantity proportional to the dissipation rate $\tau_{d}^{-1}$. The latter
refers to the dissipation time $\tau_{d}^{-1}$, estimated as the mean power
divided by the heat. Therefore, to have (IV.7) bounded, it suffices that
relaxation times $\tau$ do not exceed the dissipation time $\tau_{d}$:

$\tau<\tau_{d}$ (IV.8)

Or, our second condition on the transition rates states that the mutual
differences in relaxation times do not exceed the dissipation time scale. In
the case of detailed balance, the dissipation time $\tau_{d}=\infty$ and
$V_{\beta}(x)=E(x)-\langle E\rangle_{\beta}$ is automatically bounded indeed.

For a different yet manageable mathematical approach, we observe that (III.5)
is of the form

$V(x,\beta)=\int_{0}^{\infty}\textrm{d}t\,(e^{tL}\dot{q}(x)-\dot{q}_{\beta})$

where we used (III.3) with the backward generator $L$. Doing the integral over
time $t$, we see that (III.5) is equivalent with

$LV\,(x,\beta)=\dot{q}_{\beta}-\dot{q}(x)$

Then, using (III.2) and substituting (III.6), we get the linear equations

$\sum_{y}k_{\beta}(x,y)[V(y,\beta)-V(x,\beta)+q(x,y)]=\dot{q}_{\beta}$ (IV.9)

to be satisfied for all states $x$ (while the right-hand side does not depend
on $x$). The solution is in terms of a graphical representation of the
$V(x,\beta)-V(y,\beta)$; see arxiv  where the main step is an application of
the matrix-forest theorem. To understand the issue, note that when, in the
limit $\beta\uparrow\infty$, the transition rates
$k_{\beta}(x,y)\sim\exp(-c\beta)$ for some $c>0$, get exponentially smaller
than the dissipated power $\dot{q}_{\beta}$, the differences
$|V(x,\beta)-V(y,\beta)|$ may grow exponentially with $\beta$ and yet satisfy
(IV.9). Since we want bounded $V$, that needs to be avoided. That will be
illustrated in Example V.2.

For getting more concrete, suppose that in (IV.2)

$\displaystyle\phi(x,y)$ $\displaystyle=$ $\displaystyle 0\quad\qquad\text{
when }q(x,y)\;\geq 0$ $\displaystyle=$ $\displaystyle q(x,y)\;\,\text{ when
}q(x,y)<0$

in terms of the heat $q(x,y)$ obtained in (III.4), so that $\phi^{*}\leq 0$.
The $\phi^{*}$ was introduced in (IV.4) and is dominated by the accessibility
of the most probable state $x^{*}$; see (IV.5). An extension of the Kirchhoff
formula applied to (IV.9) (the matrix-forest theorem as explained in arxiv ),
shows that asymptotically $V(x,\beta)=v_{\beta}(x)\,e^{-\beta\phi^{*}}$ where
$v_{\beta}(x)$ is uniformly bounded in $\beta\uparrow\infty$. Therefore, to
have $V(x,\beta)$ bounded, it suffices that $\phi^{*}=0$ in (IV.4). That is
the case when there is a spanning tree with edges directed towards $x^{*}$
along which all $q(x,y)\geq 0$. Hence, it certainly suffices that the digraph,
the directed subgraph obtained by only keeping the edges where $q(x,y)\geq 0$,
is sufficiently well-connected to allow reaching $x^{*}$ from everywhere in
the graph.
We conclude that for the extended Third Law to hold, in addition to (IV.5), it
suffices there is a rooted spanning tree in the graph along edges where the
dissipated heat $q(x,y)\geq 0$. Then, the heat capacity vanishes at absolute
zero, as announced in (IV.1).
As an extreme example, on the complete graph and for every choice of rates
realizing (IV.2), the extended Third Law always holds, because then every edge
has a direction where $q(x,y)\geq 0$, and we know that a tournament has a
Hamiltonian path; see book . The fact that for denser graphs the extended
Third Law gets satisfied easier (because of greater accessibility of $x^{*}$),
is of course compatible with having smaller relaxation times, as we used
before. It also indicates the role of zero-temperature tunneling: as it
increases the connectivity of the graph, it promotes the validity of the
extended Third Law.

##  V Illustrations

The following examples explain in more concrete ways how to think of the
(extra with respect to equilibrium) dynamical condition for the nonequilibrium
Third Law. Counterexamples allow us to see in what sense that condition is
natural, hence necessary in some sense.

###  V.1 Driven reaction with source/sink

Consider the three-cycle $A\leftrightarrow B\leftrightarrow C\leftrightarrow
A$, while all three states are connected to $Z$.

Figure 5: right: Cartoon depicting a circuit of three states $A\leftrightarrow
B\leftrightarrow C\leftrightarrow A$, and possible transitions to a
sink/source state $Z$ in the center. left: the energy landscape modeling the
transitions between the states in the cycle and the source/sink.

The transition rates are

$\displaystyle k(A,B)=k(B,C)=k(C,A)$ $\displaystyle=$
$\displaystyle\frac{\nu}{1+e^{-\beta\zeta}}\;$ $\displaystyle
k(B,A)=k(C,B)=k(A,C)$ $\displaystyle=$
$\displaystyle\frac{\nu}{1+e^{\beta\zeta}}\;\text{ when
}\eta_{i+1}=1=1-\eta_{i}$

where $\zeta>0$ is the work done by a driving force, making $A\rightarrow
B\rightarrow C\rightarrow A$ a dissipative cycle.
Furthermore, there are transitions from and toward a sink/source state $Z$,
with

$\displaystyle k(A,Z)=k(B,Z)=k(C,Z)$ $\displaystyle=$
$\displaystyle\tilde{\nu}e^{-\beta\Delta}$ $\displaystyle
k(Z,A)=k(Z,B)=k(Z,C)$ $\displaystyle=$
$\displaystyle\tilde{\nu}\,e^{-\beta\Delta}e^{\beta\mu}$ (V.1)

where $\mu$ plays the role of a chemical potential, and $\Delta\geq 0$ is a
kinetic barrier for entering of leaving the dissipative cycle. for loading a
dot. Here, $\mu$ is the chemical potential at zero temperature (i.e., the
Fermi energy) and the energy $a\geq 0$ is a kinetic barrier.
Using the above scheme, the heat capacity (III.8) has been calculated (see
pnas  for a similar case) and was found to be,

$\,C(\beta)=\frac{\beta^{2}\mu^{2}e^{-\beta\mu}}{(1+e^{-\beta\mu})^{2}}-\beta^{2}\mu\zeta\frac{\nu}{\tilde{\nu}}\,\frac{(e^{\beta(a-3\mu)}-e^{\beta(a-2\mu)})}{(1+e^{-\beta\mu})^{4}}\tanh(\beta\zeta/2)\
$ (V.2)

The second term gathers the nonequilibrium (driving $\zeta$) and kinetic
($\Delta,\nu,\tilde{\nu}$) parameters.

We have Third Law behavior when $\mu<-\Delta$ or $\mu>\Delta/2$; then
$C(\beta\uparrow\infty)\to 0$. Yet, small $|\mu|/\Delta$ leads to violations:
$C(\beta\uparrow\infty)\to\infty$ diverges for $-\Delta<\mu<\Delta/2$ when
$\mu\,\Delta\,\zeta\,\nu\neq 0$, pnas . The previous section and condition
(IV.8) more specifically, make clear what happens: relaxation times $\tau$ get
too large compared to the dissipation time $\tau_{d}$. E.g., for $\mu<0$, the
state $Z$ has the largest stationary probability, but the cyclic reaction
typically runs for a time $\tau\propto\tilde{\nu}^{-1}\,e^{\beta\Delta}$
before the system enters state $Z$. On the other hand, the dissipation time
$\tau_{d}$ is given by the inverse of the current,
$\tau_{d}\propto\nu^{-1}\,e^{-\mu\beta}$. In other words, $\tau<\tau_{d}$ when
$\Delta<-\mu$. Similarly, for $\mu>0$, the dissipation time is
$\tau_{d}\propto\nu^{-1}\,e^{\mu\beta}$ and the relaxation time is
$\tau\propto\tilde{\nu}^{-1}\,e^{\beta(\Delta-\mu)}$ for the empty state (= to
get stationary). Hence, the extended Third Law holds when $\tau<\tau_{d}$ or
$\mu>\Delta/2$, but at a chemical potential $\mu=-\Delta$ and at
$\mu=\Delta/2$ zero-temperature transitions in thermal properties occur,
unseen in equilibrium $\zeta=0$.

###  V.2 Avoiding delay — staying connected

We illustrate the physics explained below (IV.9) with a more abstract example.
Consider the tree-loop of Fig. 6, where also the constant driving of strength
$\zeta$ and the energy landscape are indicated. The driving is counter-
clockwise for $\zeta>0$.

Figure 6:  A driven system with 5 states. Depending on $E(u)=U$ the
quasipotential $V_{\beta}$ remains uniformly bounded ($U\leq 6$).

The dissipative structure is the triangle $x\rightarrow z\rightarrow
y\rightarrow x$ carrying a nonzero steady current. It is easy to verify that
$z=x^{*}$ (see (IV.3)) is the unique dominant state and that the dissipative
power $\dot{q}_{\beta}$ is of order $\exp(-\beta)$, which is the escape rate
from $z$. The energy difference over the outer edge $(uw)$ is $E(w)-E(u)=5-U$,
which means that $k_{\beta}(w,u)\simeq\exp[-\beta(U-5)]$ for $U>5$. Hence, in
(IV.9), we can take a bounded edge difference $V(w)-V(u)$ as long as $U\leq
6$. There is a trap over the edge $(w,u)$ when $U$ is too large, with a too
long delay for relaxation for the quasipotential in $V(w)$ to remain bounded.
That stands for the general scenario. If we would add for example the
possibility of tunneling between the vertices $x\leftrightarrow w$, the
extended Third Law (IV.1) would be satisfied again.

##  VI Quantum context

The models we are using in the previous sections (Markov jump processes),
concern systems without quantum coherence or, equivalently, in the weak
coupling limit. Alternatively, we can consider it as the dynamics one obtains
in the Born-Oppenheimer approximation or via the Dirac-Fermi Golden Rule
beyond the decoherence time. In that sense, it is easiest to imagine a
microscopic system for which we assume, for simplicity, that the decoherence
time is much smaller than the relaxation time.
In that same spirit of weak coupling, the dynamics we consider are Markovian,
i.e. memoryless for the residence times. That is why the formalism can be
specific in its condition (III.4) of (local) detailed balance in terms of the
transition rates.
Note also that diffusion processes, as modeled by Langevin or Smoluchowski
equations, are no good for discussing very low temperatures in the same way as
classical (Newtonian) physics breaks down there. In particular, for exactly
the same reason as in equilibrium, we cannot expect an extended Third Law for
(nonequilibrium) diffusion processes.

Let us finally mention a point that we learn from the previous examples
illustrating the issue of obtaining the boundedness of the quasipotential.
Interestingly, the validity of the extended Third Law gets enhanced when the
low-temperature process does not get trapped or frozen in certain states that
do not contribute to the stationary dissipative condition. In other words,
dynamic activity helps. That then is a feature that is typically quantum as we
expect that due to quantum fluctuations, the dynamics remains sufficiently
active even at close-to-zero temperature, e.g. via tunneling events or, even
more generally, just because of the Heisenberg uncertainty relation.

##  VII Significance

Extending the Third Law to nonequilibrium processes, delivers a nontrivial
result on thermal properties of steady states. As a matter of fact, very few
such general results exist at all. The main operational significance of the
extended Third Law is that heat capacities become zero when moving to very low
temperatures, universally in all parameters. We emphasize that the operational
and physically-strongest meaning of the Third Law, both in equilibrium and in
nonequilibrium, does not need to proceed via the Second Law at all. Heat,
excess heat, and heat capacity are all perfectly well-defined and have
physical meaning without the existence of entropy. Obviously, the Second Law
is satisfied (independently) and there is strict positivity of the entropy
production, but it plays no role here and that is normal. Even in equilibrium,
entropies are measured via heat capacity, and the Third Law gets
experimentally validated by low-temperature calorimetry. As the standard
(equilibrium) Third Law but now extended to nonequilibria, that sets
constraints to all low-temperature transformations. Even though the total net
heat flux cannot vanish now, the heat capacity indeed does go to zero for a
large class of nonequilibrium systems.

While the extended Third Law may well be an important result on the thermal
susceptibility of nonequilibrium systems at low ambient temperature, measuring
nonequilibrium heat capacities does constitute a new challenge, dio ; cera ;
kkm . The extended Third Law begs for experimental confirmation. As explained
in calo ; dm ; ND ; pnas , an AC-calorimetric scheme should work perfectly
well and we hope the present work motivates further experimental work on low-
temperature calorimetry for out-of-equilibrium systems.

##  VIII Conclusion

Under a static “nondegeneracy” condition and a dynamic “no-trapping”
condition, we have argued that the steady nonequilibrium heat capacity tends
to zero with temperature. Both general theoretical arguments and specific
examples have been used to specify those conditions.

Concerning the static condition, exponential versus polynomial decay of the
heat capacity, when the ambient temperature tends to zero, informs about the
degeneracy structure of nonequilibrium steady systems. The dependency of the
low-temperature heat capacity on the behavior of the stationary condition
$\rho_{\beta}$ is an important instance of how heat, also in nonequilibrium
statistical mechanics, informs about low-temperature degeneracy, extending the
connection which exists in equilibrium between the Clausius heat and the
Boltzmann entropy, Balianbookvolume1 .
Interestingly, our nonequilibrium extension also necessitates an (extra)
dynamical condition, which asks for sufficient connectivity of the graph of
states, allowing low-temperature dynamical activity so that relaxation times
are smaller than the dissipation time.

Acknowledgments: We are grateful to Pritha Dolai for discussions on the
example in Fig. 5, see pnas , and to Irene Maes for discussions on some
mathematical aspects, see arxiv .

Conflict of Interest: The authors have no conflicts to disclose.
Data Availability Statement: The data that supports the findings of this study
are available within the article .

## References

  * (1) H. B. Callen, Thermodynamics and an Introduction to Thermostatistics. Wiley; 2nd edition (1991).
  * (2) E. Fermi, Thermodynamics. Dover Publications; New Ed edition (June 1, 1956).
  * (3) M. J. Klein, Einstein, Specific Heats and the Early Quantum Theory. Science, New Series, American Association for the Advancement of Science 148, 3667, 173-180 ( April 9, 1965).
  * (4) W. Nernst, The New Heat Theorem. Methuen and Company, Ltd. (1926).
  * (5) M. Planck, Thermodynamik. 3rd edn (1911).
  * (6) L. Masanes and J. Oppenheim, A general derivation and quantification of the third law of thermodynamics. Nat. Commun. 8, 14538 (2017).
  * (7) E.H. Lieb, Residual Entropy of Square Ice. Phys. Rev. 162, 162 (1967).
  * (8) M. Aizenman and E. Lieb, The third law of thermodynamics and the degeneracy of the ground state for lattice systems. J. Stat. Phys. 24, 279–297 (1981).
  * (9) P. Glansdorff, G. Nicolis and I. Prigogine, The thermodynamic stability theory of non-equilibrium states. Proc. Nat. Acad. Sci. USA 71, 197–199 (1974).
  * (10) Y. Oono and M. Paniconi, Steady state thermodynamics. Prog. Theor. Phys. Suppl. 130, 29 (1998).
  * (11) T.S. Komatsu, N. Nakagawa, S.I. Sasa and H. Tasaki, Steady State Thermodynamics for Heat Conduction – Microscopic Derivation. Phys. Rev. Lett. 100, 230602 (2008).
—, J. Stat. Phys. 134, 401 (2009).

  * (12) S. Bhar, M. Swar, U. Satpathi, S. Sinha, R. D. Sorkin, S. Chaudhuri and S. Roy, Measurements and analysis of response function of cold atoms in optical molasses. Optics Continuum 1, 171–188 (2022).
  * (13) C. Maes, Nonequilibrium entropies. Physica Scripta 86, 058509 (2012).
  * (14) E. Pop, Energy Dissipation and Transport in Nanoscale Devices. Nano Res 3, 147 (2010).
  * (15) B. Schmittmann and R.K.P. Zia, Statistical mechanics of driven diffusive systems. Phase transitions and critical phenomena 17, Academic Press, London (1995).
  * (16) S. Chen, H. Xie, Y. Zhang, X. Cui and G. Chen, Quantum transport through an array of quantum dots. Nanoscale, 5, 169-173 (2013).
  * (17) T. Hatano and S.-I. Sasa, Steady-state thermodynamics of Langevin systems. Phys. Rev. Lett. 86, 3463 (2001).
  * (18) C. Maes and K. Netočný, Revisiting the Glansdorff-Prigogine criterion for stability within irreversible thermodynamics. J. Stat. Phys. 159, 1286–1299 (2015).
  * (19) C. Maes and K. Netočný, A nonequilibrium extension of the Clausius heat theorem. J. Stat. Phys. 54, 188–203 (2014).
  * (20) L. Bertini, D. Gabrielli, D. Jona-Lasinio, and C. Landim, Clausius inequality and optimality of quasistatic transformations for nonequilibrium stationary states. Phys. Rev. Lett. 110, 020601 (2013).
  * (21) K. Saito and H. Tasaki, Extended Clausius Relation and Entropy for Nonequilibrium Steady States in Heat Conducting Quantum Systems. J. Stat. Phys. 145, 1275-1290 (2011).
  * (22) Z. Wang, L. Wang, J. Chen, et al. Geometric heat pump: Controlling thermal transport with time-dependent modulations. Front. Phys. 17, 13201 (2022).
  * (23) F. Khodabandehlou, C. Maes, I. Maes and K. Netočný, The vanishing of excess heat for nonequilibrium processes reaching zero ambient temperature, arXiv:2210.09858v1 [cond-mat.stat-mech]
  * (24) E. Boksenbojm, C. Maes, K. Netočný, and J. Pešek, Heat capacity in nonequilibrium steady states. Europhys. Lett. 96, 40001 (2011).
  * (25) J. Pešek, E. Boksenbojm, and K. Netočný, Model study on steady heat capacity in driven stochastic systems. Cent. Eur. J. Phys. 10 (3), 692–701 (2012).
  * (26) C. Maes and K. Netočný, Nonequilibrium Calorimetry. J. Stat. Mech. 114004, (2019).
  * (27) D. Mandal, Nonequilibrium heat capacity. Phys. Rev. E 88, 062135 (2013).
  * (28) J.-T. Hsiang, C. H. Chou, Subaşı, Y. and B. L. Hu, Quantum thermodynamics from the nonequilibrium dynamics of open systems: Energy, heat capacity, and the third law, Phys. Rev. E, 97, 012135 (2018).
  * (29) K. Sekimoto, Stochastic energetics. Lecture Notes in Physics 799, (2010).
  * (30) M. Esposito and C. Van den Broeck, Three faces of the second law. I. Master equation formulation. Phys. Rev. E 82, 011143 (2010).
  * (31) P. Dolai, C. Maes and K. Netočný, Calorimetry for active systems. SciPost Physics (2023).
  * (32) F. Khodabandehlou, S. Krekels and I. Maes, Exact computation of heat capacities for active particles on a graph. Journal of Statistical Mechanics: Theory and Experiment, 123208, (2022).
  * (33) F. Avanzini, N. Freitas and M. Esposito, Circuit Theory for Chemical Reaction Networks. arXiv:2210.08035v1 [q-bio.MN].
  * (34) A. Wachtel, R. Rao and M. Esposito, Free-Energy Transduction in Chemical Reaction Networks: from Enzymes to Metabolism arXiv:2202.01316v2 [q-bio.MN].
  * (35) A. Lazarescu, T. Cossetto, G. Falasco and M. Esposito, Large deviations and dynamical phase transitions in stochastic chemical networks. arXiv:1902.08416v1 [cond-mat.stat-mech].
  * (36) M. Polettini and M. Esposito, Irreversible thermodynamics of open chemical networks. I. emergent cycles and broken conservation laws, The Journal of Chemical Physics 141, 024117 (2014).
  * (37) S. Dal Cengio, V. Lecomte and M. Polettini, Geometry of nonequilibrium reaction networks. arXiv:2208.01290v2 [cond-mat.stat-mech].
  * (38) J. Schnakenberg, Network theory of microscopic and macroscopic behavior of master equation systems, Rev. Mod. Phys. 48, 571 (1976).
  * (39) N. G. Van Kampen, Stochastic Processes in Physics and Chemistry. 3rd Edition (2007).
  * (40) G. Grimmett and D. Stirzaker, Probability and Random Processes. Oxford University Press, third edition (2001).
  * (41) C. Maes, Local detailed balance. SciPost Phys. Lect. Notes 32, 17pp (2021).
  * (42) S. Katz, J.L. Lebowitz, and H. Spohn, Phase transitions in stationary nonequilibrium states of model lattice systems. Phys. Rev. B 28, 1655(R) (1983)
– Nonequilibrium steady states of stochastic lattice gas models of fast ionic
conductors. J.Stat. Phys. 34, 497–537 (1984).

  * (43) J. L. Lebowitz and H. Spohn, A Gallavotti Cohen-Type Symmetry in the Large Deviation Functional for Stochastic Dynamics. J.Stat. Phys. 95 (1999).
  * (44) R.D. Astumian et al. Non-equilibrium kinetics and trajectory thermodynamics of synthetic molecular pumps. Materials Chemistry Frontiers 5 (2020).
R.D. Astumian. Kinetic asymmetry allows macromolecular catalysts to drive an
information ratchet. Nat. Commun. 10, 3837 (2019).

  * (45) P. A. Lewis, C. E. Inman, F. Maya, J. M. Tour, J. E. Hutchison and P. S. Weiss, Molecular Engineering of the Polarity and Interactions of Molecular Electronic Switches. J. Am. Chem. Soc. 127, 17421–17426 (2005).
  * (46) W. Liu, S. Yang, J. Li, G. Su and J. Chang Ren, One molecule, two states: Single molecular switch on metallic electrodes. WIREs Computational Molecular Science 11, 1759–0876 (2021).
  * (47) C. Maes and K. Netočný, Heat bounds and the blowtorch theorem. Ann. H. Poincaré 14, 1193–1202 (2013).
  * (48) C. Maes, K. Netočný and W. O’Kelly de Galway, Low temperature behavior of nonequilibrium multilevel systems. J. Phys. A: Math. Theor. 47, 035002 (2014).
  * (49) F. Khodabandehlou, C. Maes and K. Netočný, Trees and forests for nonequilibrium purposes: an introduction to graphical representations. J. Stat. Phys. 189, 41 (2022).
  * (50) D. B. West. Introduction to graph theory. Pearson Education, (2002).
  * (51) P. Dolai and C. Maes, Towards many-body nonequilibrium calorimetry: specific heat for a driven fermionic array. arXiv:2301.04524v1 [cond-mat.stat-mech], (2022).
  * (52) M. J. de Oliveira, Complex heat capacity and entropy production of temperature modulated systems. arXiv:1905.10306v1 [cond-mat.stat-mech] (2019).
  * (53) J. del Cerro and S.Ramos, Specific heat of latgs ferroelectric crystal under dissipative conditions. Ferroelectrics Letters, 16, 119 (1993).
  * (54) J. K. Nielsen and J. C. Dyre, Fluctuation-dissipation theorem for frequency-dependent specific heat. Phys. Rev. B 54, 15754 (1996).
  * (55) R. Balian, From Microphysics to Macrophysics: Methods and Applications of Statistical Physics, Volume 1. Springer Berlin, Heidelberg (2007).