# Screening of electric field and nuclear EDM in non-stationary states of
atoms and molecules

V.V. Flambaum  School of Physics, University of New South Wales, Sydney 2052,
Australia

###### Abstract

According to the Schiff theorem, an external electric field vanishes at atomic
nucleus in a neutral atom in a stationary state, i.e. it is completely
shielded by electrons. This makes a nuclear electric dipole moment (EDM)
unobservable. We show that if atom or molecule is not in a stationary state
(e.g. in a superposition of two stationary states), electric field on the
nucleus is not zero and interaction with nuclear EDM does not vanish. In
molecules this effect is enhanced by the ratio of nuclear mass to to electron
mass, $M_{n}/m_{e}$ , since nuclei in a molecule are slow (compare to
electrons) and do not provide efficient screening in a non-stationary
environment. Electric field on the nucleus may also affect nuclear reactions.

##  I Introduction

Existence of electric dipole moments (EDM) of elementary particles, nuclei,
atoms and molecules in a state with a definite angular momentum violates time
reversal invariance (T) and parity (P). EDM also violates CP invariance if the
CPT invariance holds. A very extensive experimental and theoretical activity
related to EDM is motivated by the need to test unification theories
predicting T, P, and CP violation.

However, there is a problem here. A homogeneous static electric field does not
accelerate neutral atom. This means that the total electric field ${\bf E}$
acting on the atomic nucleus is zero since otherwise the charged nucleus would
be accelerating, i.e. the external field is completely shielded by atomic
electrons. The absence of the electric field means that the nuclear EDM $d$ is
unobservable, ${\bf d\cdot E}=0$. One may also present this result
differently: total atomic EDM is zero even if the nucleus has EDM, i.e. EDM of
electron cloud in atom is exactly opposite in sign to nuclear EDM.

A quantum-mechanical derivation of this result for an arbitrary non-
relativistic system of point-like charged particles with EDMs has been done by
Schiff Schiff . Schiff also mentioned that his theorem is violated by the
finite nuclear size. The effect of the finite nuclear size was implemented as
the nuclear Schiff moment which was introduced in Refs. Sandars ; Hinds ; SFK
; FKS1986 . An electrostatic interaction between the nuclear Schiff moment and
electrons produces atomic and molecular EDM. Refs. Sandars ; Hinds  calculated
the finite nuclear size effect of the proton EDM. Refs. SFK ; FKS1986
calculated (and named) the nuclear Schiff moment produced by the P,T-odd
nuclear forces. It was shown in SFK  that the contribution of the P,T-odd
forces to the nuclear EDM and Schiff moment is $\sim 40$ times larger than the
contribution of the nucleon EDM. An additional 2-3 orders of magnitude
enhancement appears in nuclei with the octupole deformation Auerbach .

The suppression factor for the atomic EDM relative to the nuclear EDM,
proportional to a very small ratio ($\sim 10^{-9}$) of the squared nuclear
radius to the squared atomic radius, is partly compensated by the factor
$Z^{2}R_{S}$, where $Z$ is the nuclear charge and $R_{S}$ is the relativistic
factor SFK .

The Schiff theorem is also violated by the magnetic interaction Schiff ;
Khriplovich . Corresponding atomic EDMs produced by the nuclear EDM and
electron-nucleus magnetic interaction have been calculated in Ref. FGP . In
light atoms this mechanism of atomic EDM dominates but in heavy atoms it is
smaller than the effect of the finite nuclear size since the latter very
rapidly increases with the nuclear charge, as $Z^{2}R_{S}$, while the magnetic
effect increases slower, as $ZR_{M}$ where $R_{M}$ is the relativistic factor
for the magnetic effect FGP  .

There is no complete shielding in ions. For example, in a molecular ion the
shielding factor for the nuclear EDM is $(Z_{i}/Z)(M_{n}/M_{m})$, where
$Z_{i}$ is the ion charge, $Z$ is the nuclear charge, $M_{n}$ is the nuclear
mass and $M_{m}$ is the total molecular mass FlambaumKozlov .

Screening of time-dependent electric field is incomplete. External electric
field can even be enhanced if the frequency of electric filed oscillations is
in resonance with atomic or molecular transition. Screening of oscillating
field has been investigated in Refs. DFSS1986 ; Flambaum2018 ; DBGF2018 ;
FlambaumSamsonov ; Bao .

Penetration of electric field to atomic nucleus may affect nuclear reactions -
see e.g. Refs. LaserNucl ; Ho , where neutron capture to a nucleus, enhanced
by a laser field, is discussed. For example, electric field may admix large
s-wave neutron capture amplitude to a kinematically suppressed p-wave
amplitude in a p-wave resonance and produce effects such as asymmetry in the
photon distribution correlated with the direction of the electric field. This
effect is somewhat similar to the effects of parity violating interaction
which mixes s-wave and p-wave resonances. Such kinematic enhancement factor
(ratio of s-wave and p-wave amplitudes) is $\sim 10^{3}$ for slow neutrons.
Another enhancement factor is due to a very small energy interval between the
energy levels in a compound nucleus. These two enhancement factors lead to a
$10^{6}$ enhancement of the parity violating effects in neutron reactions
predicted in Refs. SF1980 ; SF1982 ; Nucl1 ; Nucl2  and confirmed in numerous
experiments involving a hundred of p-wave resonances in many nuclei - see
reviews Gribakin ; Bowman ). The same enhancement may appear if the s and p
resonances mixing is produced by the electric field.

The aim of this paper is to show that for atoms and molecules in a non-
stationary state the electric filed on the nucleus does not vanish and to
derive formulas for this electric field. This electric field may interact with
nuclear EDM and affect nuclear reactions.

##  II  Shielding theory for non-stationary atomic states

The Hamiltonian of a neutral atom in an external electric field along the
$z$-axis $E^{ext}=E^{ext}_{z}$ may be presented as

$\displaystyle H_{E}=H_{0}-E^{ext}_{z}D_{z}\,,$ (1) $\displaystyle
D_{z}=-e\sum_{k=1}^{Z}z_{k}\,,$ (2)

where $H_{0}$ is the Schrodinger or the Dirac Hamiltonian for the atomic
electrons in the absence of the external field $E^{ext}_{z}$, $Z$ is the
number of the electrons and protons, $-e$ is the electron charge (i.e. $e$ is
the proton charge), $z_{k}$ is the $z$-axis projection of the electron
position relative to the nucleus. We assume that the nuclear mass is infinite
and neglect very small effects of the Breit and magnetic interactions. The
total electric field on the nucleus may be presented as ${\bf E}^{t}={\bf
E}^{ext}+<{\bf E}^{e}>$, where the $z$ component of the electron electric
field on the nucleus is

$E^{e}_{z}=e\sum_{k=1}^{Z}\frac{z_{k}}{r_{k}^{3}}=\frac{i}{Ze\hbar}[P_{z},H_{0}]\,,$
(3)

where $P_{z}=\sum_{k=1}^{Z}p_{z,k}$ is the total momentum of the atomic
electrons. The second equality follows from the differentiation of the nuclear
Coulomb potential in the Dirac or Schrodinger Hamiltonian $H_{0}$ since the
total electron momentum $P_{z}$ commutes with the electron kinetic energy and
the electron-electron interaction. Similarly, the $z$ component of the total
electric field on the nucleus may be presented as

$E^{t}_{z}=\frac{i}{Ze\hbar}[P_{z},H_{E}]\,.$ (4)

In agreement with the Schiff theorem, in a stationary state $|n>$ expectation
value of the total electric field on the nucleus vanishes,
$<n|E^{t}_{z}|n>=0$, since
$<n|[P_{z},H_{E}]|n>=(\epsilon_{n}-\epsilon_{n})<n|P_{z}|n>=0$.

A non-stationary state may be presented as a sum over stationary states. For
brevity, we include two states in the sum:

$\psi=c_{a}\psi_{a}\exp{(-\frac{i}{\hbar}\epsilon_{a}t)}+c_{b}\psi_{b}\exp{(-\frac{i}{\hbar}\epsilon_{b}t)}$
(5)

In such state the $z$ component of the total electric field on the nucleus is

$\displaystyle<E^{t}_{z}>=-\frac{i(\epsilon_{a}-\epsilon_{b})}{Ze\hbar}[c_{a}^{*}c_{b}<a|P_{z}|b>\exp{(\frac{i(\epsilon_{a}-\epsilon_{b})t}{\hbar})}$
$\displaystyle-
c_{a}c_{b}^{*}<b|P_{z}|a>\exp{(-\frac{i(\epsilon_{a}-\epsilon_{b})t}{\hbar})}]\,.\,$
(6)

It is also instructive to present $<E^{t}_{z}>$ using a substitution of the
nonrelativistic expression for the momentum,
$P_{z}=-\frac{im}{e\hbar}[H_{E},D_{z}]$:

$\displaystyle<E^{t}_{z}>=-\frac{(\epsilon_{a}-\epsilon_{b})^{2}m}{Ze^{2}\hbar^{2}}\times$
$\displaystyle[c_{a}^{*}c_{b}<a|D_{z}|b>\exp{(\frac{i(\epsilon_{a}-\epsilon_{b})t}{\hbar})}$
$\displaystyle+c_{a}c_{b}^{*}<b|D_{z}|a>\exp{(-\frac{i(\epsilon_{a}-\epsilon_{b})t}{\hbar})}],$
(7)

where $m$ in this expression is electron mass. Note that
$\psi_{a}=\psi^{(0)}_{a}+\delta\psi_{a}$ and
$\psi_{b}=\psi^{(0)}_{b}+\delta\psi_{b}$ are eigenfunctions of the Hamiltonian
$H_{E}$ including interaction with the external electric field $E^{ext}$. Here
$\psi^{(0)}_{a}$ and $\psi^{(0)}_{b}$ are eigenfunctions of the unperturbed
Hamiltonian $H_{0}$.

We have two different cases here. If the matrix element between unperturbed
wave functions satisfy selection rules for the electric dipole matrix element,
i.e. $<\psi^{(0)}_{a}|D_{z}|\psi^{(0)}_{b}>$ is not equal to zero, we have
oscillating electric filed on the nucleus even in the absence of external
electric field $E^{ext}$. Indeed, atom in such state Eq. (5) has oscillating
electric dipole moment $<D_{z}>$ which produces electric field on the nucleus.
For real $c_{a}$, $c_{b}$ and $<a|D_{z}|b>$ we obtain

$<D_{z}>=2c_{a}c_{b}<a|D_{z}|b>\cos{(\frac{(\epsilon_{a}-\epsilon_{b})t}{\hbar})}$
(8)

and

$<E^{t}_{z}>=-\frac{2(\epsilon_{a}-\epsilon_{b})^{2}m}{Ze^{2}\hbar^{2}}c_{a}c_{b}<a|D_{z}|b>\cos{(\frac{(\epsilon_{a}-\epsilon_{b})t}{\hbar})}.$
(9)

Numerical estimate for the amplitude of this oscillating field is

$|E^{t}_{z}|\sim\frac{(\epsilon_{a}-\epsilon_{b})^{2}}{Z(\text{eV})^{2}}10^{7}\text{V}/\text{cm}.$
(10)

We assumed $c_{a}c_{b}\sim 1$. For $Z\sim 1$ and $(\epsilon_{a}-\epsilon_{b})$
equal to few eV, this field may exceed by three orders of magnitude electric
fields, which have been used to measure neutron and atomic EDM. However, this
is a very rapidly oscillating electric field. In the case of a small
oscillation frequency, the electric field is strongly suppressed by the factor
$(\epsilon_{a}-\epsilon_{b})^{2}$.

When $<\psi^{(0)}_{a}|D_{z}|\psi^{(0)}_{b}>=0$, we should consider effect
produced by the external electric field $E^{ext}$. Substitution of the
perturbation theory expression for $\delta\psi$ gives, for real $c_{a}$,
$c_{b}$ and matrix elements of electric dipole moment operator $<a|D_{z}|n>$
and $<b|D_{z}|n>$,

$<E^{t}_{z}>=\frac{2c_{a}c_{b}(\epsilon_{a}-\epsilon_{b})^{2}m}{Ze^{2}\hbar^{2}}\alpha_{a,b}E^{ext}_{z}\cos{(\frac{(\epsilon_{a}-\epsilon_{b})t}{\hbar})},$
(11)

where

$\displaystyle\alpha_{a,b}=\sum_{n}\frac{<a|D_{z}|n><n|D_{z}|b>}{\epsilon_{a}-\epsilon_{n}}+$
$\displaystyle\frac{<a|D_{z}|n><n|D_{z}|b>}{\epsilon_{b}-\epsilon_{n}}$ (12)

is the Stark amplitude between the states $a$ and $b$, and $|n>$ are
intermediate states in the perturbation theory sum for $\delta\psi$. We see
that constant external electric field $E^{ext}_{z}$ is transformed into
oscillating electric field on the nucleus. Numerical estimate for the
amplitude of the field on the nucleus for $c_{a}c_{b}\sim 1$ is

$|E^{t}_{z}|\sim\frac{(\epsilon_{a}-\epsilon_{b})^{2}}{Z(27\text{eV})^{2}}E^{ext}_{z}.$
(13)

Thus, for $\epsilon_{a}-\epsilon_{b}$ smaller than atomic unit of energy $27$
eV, electric field on the nucleus is smaller than external electric field.

##  III Shielding in non-stationary molecular states

In molecules in a stationary rotational state the screening of external
electric field is produced by both electrons and nuclei. However, in a non-
stationary state electric field on the nucleus is proportional to mass of the
particles which produce this screening - see Eq. (II). Mass of nuclei is from
three to six orders of magnitude bigger than mass of electron. Therefore,
electric field on the nucleus may be significantly bigger in molecules compare
to atoms (for equal values of $\epsilon_{a}-\epsilon_{b}$). Indeed, nuclei in
molecules are slow, they are not as efficient screeners as electrons in the
case when (electron) electric field varies. We observed similar enhancement
when considered screening of oscillating external electric field in molecules
Bao .

Let us consider electric field on the nucleus 1 in a diamagnetic diatomic
molecule which is given by the following expression:

$\displaystyle{\bf E^{(1)}}=-\frac{i}{Z_{1}e\hbar}[{\bf P^{(1)}},H_{E}]\,,$
(14)

where ${\bf P^{(1)}}$ is the momentum of the nucleus 1, $Z_{1}$ is its charge
and the Hamiltonian $H_{E}$ includes both nuclei and electrons. We may
subtract from the momentum ${\bf P^{(1)}}$ the contribution of the center of
mass motion with velocity ${\bf v}={\bf P}_{t}/M_{t}$, where ${\bf P}_{t}$ is
the total momentum and $M_{t}$ is the total mass of the molecule:

${\bf\Pi}^{1}={\bf P}^{(1)}-M^{(1)}{\bf P}_{t}/M_{t}\,.$ (15)

Center of mass momentum commutes with the Hamiltonian, therefore, we may re-
write Eq. (14) using commutator with ${\bf\Pi}^{1}$:

$\displaystyle{\bf E^{(1)}}=-\frac{i}{Z_{1}e\hbar}[{\bf\Pi}^{1},H_{E}]\,,$
(16)

The expectation value of the total electric field on the nucleus 1 in the
state (5) is

$\displaystyle<{\bf
E^{(1)}}>=-\frac{i(\epsilon_{a}-\epsilon_{b})}{Z_{1}e\hbar}\times$
$\displaystyle[c_{a}^{*}c_{b}<a|{\bf\Pi}^{1}|b>\exp{(\frac{i(\epsilon_{a}-\epsilon_{b})t}{\hbar})}$
$\displaystyle-
c_{a}c_{b}^{*}<b|{\bf\Pi}^{1}|a>\exp{(-\frac{i(\epsilon_{a}-\epsilon_{b})t}{\hbar})}].\,\,\,\,\,$
(17)

Now we can use the following relation from Ref. Bao , where we have neglected
terms proportional to electron mass in comparison with the terms proportional
to the nuclear masses:

${\bf\Pi}^{1}=i\mu[{\bf R},H^{\prime}_{E}]\,,$ (18)

where $\mu=M_{1}M_{2}/(M_{1}+M_{2})$ is the reduced nuclear mass, ${{\bf
R}=\bf R^{(1)}-R^{(2)}}$ is the relative coordinate for the first and second
nucleus and $H^{\prime}_{E}$ is the Hamiltonian with subtracted contribution
of the center of mass motion. Using this relation we obtain

$\displaystyle<{\bf
E^{(1)}}>=-\frac{(\epsilon_{a}-\epsilon_{b})^{2}\mu}{Z_{1}e\hbar^{2}}\times$
$\displaystyle[c_{a}^{*}c_{b}<a|{\bf
R}|b>\exp{(\frac{i(\epsilon_{a}-\epsilon_{b})t}{\hbar})}$
$\displaystyle+c_{a}c_{b}^{*}<b|{\bf
R}|a>\exp{(-\frac{i(\epsilon_{a}-\epsilon_{b})t}{\hbar})}].$ (19)

We consider rotational molecular states $|a>=Y_{00}(\theta,\phi)\Psi$ and
$|b>=Y_{10}(\theta,\phi)\Psi$ in a diamagnetic diatomic polar molecule, where
$Y_{00}$ and $Y_{10}$ describe rotational $L=0$ and $L=1$ molecular states and
$\Psi$ is an internal molecular wave function, describing electron state and
nuclear vibrational state, which are the same for $<a|$ and $<b|$. Assuming
real $c_{a}$ and $c_{b}$ we obtain

$\displaystyle<E^{(1)}_{z}>=-\frac{2(\epsilon_{a}-\epsilon_{b})^{2}\mu
R}{\sqrt{3}Z_{1}e\hbar^{2}}c_{a}c_{b}\cos{(\frac{(\epsilon_{a}-\epsilon_{b})t}{\hbar})},$

where $R$ is the distance between the nuclei. Note that the second nucleus
experiences electric field of opposite sign (${\bf R}\to-{\bf R}$). The
electric field is inversely proportional to the nuclear charge, therefore, the
electric forces are equal in magnitude and have opposite sign, so there is no
acceleration of the center of mass.

The oscillation frequency of the electric field on the nucleus may be many
orders of magnitude smaller than frequency of optical transitions in atoms.
Indeed, the interval between molecular rotational levels is $\sim\mu/m_{e}$
times smaller than the interval between the electron levels. Moreover,
molecules may have doublets of levels with the energy interval much smaller
than the rotational energy. Note however that in this case there is strong
suppression of the field by the small factor
$(\epsilon_{a}-\epsilon_{b})^{2}$.

##  IV Conclusion

The Schiff’s shielding theorem about total screening of electric field on
atomic nucleus does not apply to atoms and molecules in a non-stationary
state, i.e. electric field on a nucleus in a neutral system in a non-
stationary state may be not equal to zero. The field on the nucleus is much
bigger in molecules than in atoms (for equal values of
$\epsilon_{a}-\epsilon_{b}$) since nuclei in molecules are slow and do not
produce such efficient screening in a non-stationary case as electrons do. The
field enhancement factor is the ratio of nuclear mass to electron mass. In
principle, electric field penetrating to the nucleus makes nuclear EDM
observable and may also affect nuclear reactions.

Schiff theorem has a simple classical explanation. Application of a
homogeneous electric field to a neutral system does not lead to its motion. In
this case the nucleus does not move, despite the fact that it has non-zero
charge. This means that the total electric field on the nucleus is zero. The
case of a neutral system in a non-stationary state is different. Nuclei may
move, keeping center of mass of the system at rest.

For brevity, we presented non-stationary state as a sum of two basis
stationary states. In general case the number of terms in the sum may be
arbitrary large and we have to add summation over basis states ($\sum_{b>a}$)
to the formulas for electric field on the nucleus.

Acknowledgements — The work was supported by the Australian Research Council
Grants No. DP230101058 and DP200100150.

## References

  * (1) L.I. Schiff, Measurability of Nuclear Electric Dipole Moments, Phys. Rev. 132, 2194 (1963).
  * (2) P. G. H. Sandars, Measurability of the Proton Electric Dipole Moment, Phys. Rev. Lett. 19, 1396 (1967).
  * (3) E. A. Hinds and P. G. H. Sandars, Electric dipole hyperfine structure of TIF, Phys. Rev. A 21, 471 (1980).
  * (4) O.P. Sushkov, V.V. Flambaum, I.B. Khriplovich, On the possibility of investigation of P- and T-odd nuclear forces in atomic and molecular experiments, Zh. Exp. Teor. Fiz. 87, 1521 (1984) [ Sov. JETP 60, 873 (1984)].
  * (5) V.V. Flambaum, I.B. Khriplovich, O.P. Sushkov, On the P- and T-nonconserving nuclear moments, Nucl. Phys. A 449, 750, 1986.
  * (6) I.B. Khriplovich, Parity Nonconservation in Atomic Phenomena (Gordon and Breach, Amsterdam, 1991).
  * (7) N. Auerbach, V. V. Flambaum, and V. Spevak, Collective $T$\- and $P$-odd electromagnetic moments in nuclei with octupole deformations, Phys. Rev. Lett. 76, 4316 (1996).
  * (8) S.G. Porsev, J.S.M. Ginges, and V.V. Flambaum, Atomic electric dipole moment induced by the nuclear electric dipole moment: the magnetic moment effect, Phys. Rev. A 83,042507 (2011).
  * (9) V. V. Flambaum and A. Kozlov, Extension of the Schiff theorem to ions and molecules, Phys. Rev. A 85, 022505 (2012).
  * (10) V.A. Dzuba, V.V. Flambaum, P.G. Silvestrov, O.P. Sushkov, Shielding of external electric field in atoms, Phys. Lett. A 118, 177,1986.
  * (11) V. V. Flambaum, Shielding of an external oscillating electric field inside atoms, Phys. Rev. A 98, 043408 (2018).
  * (12) V. A. Dzuba , J.C. Berengut, J.S.M. Ginges, V. V. Flambaum, Screening of oscillating external electric field in atoms, Phys. Rev. A 98, 043411 (2018).
  * (13) V.V. Flambaum, I. B. Samsonov, Resonant enhancement of an oscillating electric field in atom, Phys. Rev. A 98, 053437 (2018),
  * (14) H. B. Tran Tan, V. V. Flambaum, I. B. Samsonov, Screening and enhancement of oscillating electric field in molecules, Phys. Rev. A 99, 013430 (2019).
  * (15) Y.K. Ho, F.C. Khanna, M.A. Lone. Laser-induced neutron radiative capture reaction, Nucl. Phys. A 440, 189 (1985).
  * (16) Yu-Kun Ho, Zheng-Ying Pan, Laser-induced two-step population of nuclear resonances near neutron binding energies, Phys. Rev. C 39, 94 (1989).
  * (17) O.P. Sushkov, V.V. Flambaum, On the possibility to observe parity nonconservation in neutron optics, Pis’ma Zh.Exp.Teor.Fiz. 32, 377, 1980. [JETP Lett. 32, 352 (1980)]. Preprint INP 80-148 (1980).
  * (18) O.P. Sushkov, V.V. Flambaum, Parity violation in interaction of neutrons with heavy nuclei, Usp. Fiz. Nauk. 136, 3, 1982. [Sov. Phys. Usp. 25, 1 (1982).]
  * (19) V.V. Flambaum, O.P. Sushkov, Theory of neutron reactions and weak interaction in heavy nuclei, Nucl. Phys. A412, 13-33, 1984.
  * (20) V.V.Flambaum, O.P. Sushkov, Angular and polarizatron correlations in $(n,\gamma)$ reaction, Nucl. Phys. A 435, 352, 1985.
  * (21) V. V. Flambaum and G. F. Gribakin, Enhancement of parity and time invariance violating effects in compound nuclei, Progress in Particle and Nuclear Physics 35, 423 (1995).
  * (22) G. E. Mitchell, J. D. Bowman, and H. A. Weidenmuller, Parity violation in the compound nucleus, Rev. Mod. Phys. 71, 445 (1999).


# New universal Lyapunov functions
for non-linear reaction networks

Alexander N. Gorban<EMAIL_ADDRESS>Department of Mathematics, University
of Leicester, Leicester, LE1 7RH, UK
and Lobachevsky University, Nizhny Novgorod, Russia

###### Abstract

In 1961, Rényi discovered a rich family of non-classical Lyapunov functions
for kinetics of the Markov chains, or, what is the same, for the linear
kinetic equations. This family was parameterized by convex functions on the
positive semi-axis. After works of Csiszár and Morimoto, these functions
became widely known as $f$-divergences or the Csiszár–Morimoto divergences.
These Lyapunov functions are universal in the following sense: they depend
only on the state of equilibrium, not on the kinetic parameters themselves.

Despite many years of research, no such wide family of universal Lyapunov
functions has been found for nonlinear reaction networks. For general non-
linear networks with detailed or complex balance, the classical thermodynamics
potentials remain the only universal Lyapunov functions.

We constructed a rich family of new universal Lyapunov functions for any non-
linear reaction network with detailed or complex balance. These functions are
parameterized by compact subsets of the projective space. They are universal
in the same sense: they depend only on the state of equilibrium and on the
network structure, but not on the kinetic parameters themselves.

The main elements and operations in the construction of the new Lyapunov
functions are partial equilibria of reactions and convex envelopes of families
of functions.

###### keywords:

reaction networks; non-linear kinetics; Lyapunov function; partial
equilibrium; detailed balance

##  1 Introduction

The classical Lyapunov functions in kinetics are closely related to the
concepts of entropy and free energy. The first example was provided by
Boltzmann in 1872 [1]. He proved that in an isolated system the functional

$H(f)=\int f(x,v)\ln(f(x,v)){\mathrm{d}}^{3}v{\mathrm{d}}^{3}x$ (1)

does not increase in time due to kinetic equation (the Boltzmann equation).
Here, $f(x,v)$ is the distribution density of particles in the phase space,
$x$ is position in space, $v$ is velocity of a particle. Boltzmann’s proof [1]
used the principle of detailed balance. Several years later he found more
general conditions for $H$-theorem [2] and invented what we call now
semidetailed balance or cyclic balance or complex balance (for modern theory
of chemical reaction networks we refer to tutorial [3]).

For perfect chemical mixtures with components $A_{1},\ldots,A_{n}$ in
isothermal isochoric conditions (fixed volume) the analogue of Boltzmann’s
$H$-function is:

$H=\sum_{i=1}^{n}c_{i}\left(\ln\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)-1\right),$ (2)

where $c_{i}\geq 0$ is the concentration of $A_{i}$ and $c_{i}^{\rm eq}>0$ is
an equilibrium concentration of $A_{i}$ (under the standard convention that
$x\ln x=0$ for $x=0$). The idealization of perfect mixtures is applicable to
rarefied gases or to reactions of small admixtures in solutions.

The $H$-function (2) is a Lyapunov function for all mass action law systems
with detailed or complex balanced equilibrium (see e.g. [3])). We call this
property ‘universality’.

In information theory, the function $H$ appears as a measure of relative
information (in the distribution $c_{i}$ with respect to the distribution
$c_{i}^{\rm eq}$) and analogue of the $H$-theorem states that random
manipulations with data decrease the relative information with respect to the
equilibrium that does not change under manipulations [4, 5, 6].

It is not much surprising that the $H$ function (2) is essentially the only
universal Lyapunov function for all imaginable perfect kinetic systems with
detailed balance. Nevertheless, if we restrict the choice of the reaction
mechanism then the class of Lyapunov functions, which are conditionally
independent of reaction rate constants for a given detailed balanced or
complex balanced equilibrium, can be extended. We call such Lyapunov functions
conditionally universal (for a given reaction mechanism).

In 1961, Rényi discovered a class of conditionally universal Lyapunov
functions for Markov chains [7]. After works [8, 9] these functions were
studied by many authors under the name $f$-divergences or Csiszár–Morimoto
divergences. It is known that any universal Lyapunov functions for Markov
chains has the form of $f$-divergence [10, 11, 12] or is a monotonic function
of such a divergence. The continuous time Markov kinetic equation coincide
with the kinetic equations for linear (monomolecular) reactions of perfect
systems, and $f$-divergences can be considered as a direct generalization of
(2):

$H_{f}(c)=\sum_{i}c_{i}^{\rm eq}f\left(\frac{c_{i}}{c_{i}^{\rm eq}}\right),$
(3)

where $f$ is a convex function on the positive semi-axis.

For $f(x)=x\ln x$ and after adding a constant term proportional to
$\sum_{i}c_{i}$ we get the classical formula (2). (Recall that $\sum_{i}c_{i}$
does not change in linear kinetics.)

Existence of a very rich family of conditionally universal Lyapunov functions
for the linear reaction mechanisms makes us guess that there should be many
conditionally universal Lyapunov functions for any given nonlinear reaction
mechanism as well. In this paper, we construct new conditionally universal
Lyapunov functions for any given reaction mechanism using partial equilibria
of all single reactions and detailed or complex balance conditions.

In the next Sec. 2, we give the necessary formal definitions and introduce
notations for mass action law systems. The necessary and sufficient conditions
that a convex function is a Lyapunov function for all reaction networks with
given reaction mechanism and equilibrium point under detailed or complex
balance assumption are proven in Sec.3. The construction of a new family of
conditionally universal Lyapunov functions for any reaction network is
presented and the main result, Theorem 2, is proven for mass action law
systems in Sec. 4. In Sec. 5 we outline the possible generalizations and
applications of the results. In Conclusion the main results of the work are
summarized and an open question is formulated.

##  2 Prerequisites: mass action law and classical Lyapunov functions

In this section, we formally introduce mass action law and equations of
chemical kinetics. For more detailed introduction, including thermodynamical
backgrounds, detailed kinetics, applied kinetics, and mathematical aspect of
kinetics, we refer to the modern book [13]. Tutorial [3] gives the
mathematical introduction in dynamics of chemical reaction networks. Formalism
of chemical kinetics with special attention to heterogeneous catalysis is
discussed in detail in the monograph [14].

###  2.1 Mass action law

Consider a closed system with $n$ chemical species $A_{1},\ldots,A_{n}$,
participating in a complex reaction network. The reaction network is
represented in the form of the system of stoichiometric equations of
elementary reactions (called also reaction mechanism):

$\sum_{i=1}^{n}\alpha_{ri}A_{i}\to\sum_{j=1}^{n}\beta_{rj}A_{j}\;\;(r=1,\ldots,m)\,,$
(4)

where $\alpha_{ri}\geq 0$, $\beta_{rj}\geq 0$ are the stoichiometric
coefficients, $r=1,\ldots,m$, $i,j=1,\ldots,n$, $m$ is the number of
elementary reactions, $n$ is the number of components. In this representation,
the direct and reverse elementary reactions are considered separately.

The stoichiometric vector $\gamma_{r}$ of the elementary reaction is
$\gamma_{r}=(\gamma_{ri})$, $\gamma_{ri}=\beta_{ri}-\alpha_{ri}$ (‘gain minus
loss’). The gain vector is $\alpha_{r}=(\alpha_{ri})$ and the loss vector is
$\beta_{r}=\beta_{ri}$.

Elementary reactions of the form $A_{i}\to A_{j}$ are called linear or
monomolecular reactions.

According to the mass action law, the reaction rate for the elementary
reactions (4) are

$w_{r}=k_{r}\prod_{i=1}^{n}c_{i}^{\alpha_{ri}},$ (5)

where $k_{r}\geq 0$ is the reaction rate constant and the standard convention
is used: for any $x\geq 0$, $x^{0}=1$.

The kinetic equations for a perfect system in isochoric isothermal conditions
have the form

$\frac{{\mathrm{d}}c}{{\mathrm{d}}t}=\sum_{r=1}^{m}\gamma_{r}w_{r}$ (6)

The stoichiometric subspace is $Span\\{\gamma_{r}|r=1,\ldots,m\\}$. For any
set of values $k_{r}\geq 0$ and $c_{i}\geq 0$ the time derivative of $c$
belongs to the stoichiometric subspace

$\frac{{\mathrm{d}}c}{{\mathrm{d}}t}\in Span\\{\gamma_{r}|r=1,\ldots,m\\}.$

Therefore, each linear functional $b(c)$ that annuls
$Span\\{\gamma_{r}|r=1,\ldots,m\\}$ (i.e. $b(\gamma_{r})=0$ for all $r$ is the
conservation law: ${\mathrm{d}}b(c)/{\mathrm{d}}t=0$ according to (6). Such
functionals are called stoichiometric conservation law.

We always assume that there exists a strictly positive stoichiometric
conservation law that is such a vector $b=(b_{i})$, $b_{i}>0$ that
$\sum_{i}b_{i}\gamma_{ri}=0$ for all $r$. In chemical kinetics, this may be
the conservation of mass or of total number of atoms, for example. Due to this
assumption, every stoichiometric vector $\gamma$ has both positive and
negative components. This assumption has many important consequences. For
example, if

$\mbox{For all }i,r\mbox{ either }\alpha_{ri}=0\mbox{ or }\alpha_{ri}\geq 1,$
(7)

then the right hand side of (6) is a Lipschitz function in any reaction
polyhedron

$(Span\\{\gamma_{r}|r=1,\ldots,m\\}+c)\cap\mathbb{R}^{n}_{\geq 0},$ (8)

where $c$ is an arbitrary vector with non-negative coordinates and
$\mathbb{R}^{n}_{\geq 0}$ is the cone (orthant) of such vectors. Therefore,
solution of (6) exists and is unique for every initial state with non-negative
coordinates. Moreover, it remains non-negative and, therefore, any reaction
polyhedron (8) is positively invariant: solutions do not leave it in positive
time. Algorithms for detailed analysis of the structure of reaction polyhedra
are developed in [15]. It involves the classical double-description problem
[16, 17]. Reaction polyhedra for various examples of chemical reactions are
presented in [18, 19, 20]. The boundedness of positive semi-trajectories has
strong stability implications for reaction networks [21].

In many applications, the stoichiometric coefficients are non-negative
integers. For global existence, uniqueness of solutions, and preservation of
positivity the condition (7) and existence of positive stoichiometric
conservation law are sufficient. Moreover, in this work we do not need
uniqueness and global existence of solutions and analyze just the local
conditions ${\mathrm{d}}G/{\mathrm{d}}t\leq 0$ for various functions $G$ at
strictly positive points $c$. Therefore, a weaker condition is needed: instead
of (7) we assume below that $\alpha_{ri}\geq 0$, $\beta_{ri}\geq 0$.

###  2.2 Detailed balance and classical Lyapunov function

First proof of Boltzmann’s $H$-theorem [1] and its analogues for chemical
kinetics were based on detailed balance . Boltzmann’s argument were analyzed
by Tolman [22]. In general situation, this principle means that at
equilibrium, each elementary process is equilibrated with its reverse process.

To formulate this principle for chemical kinetics, we combine direct and
reverse reactions together and rewrite the reaction mechanism in the form

$\sum_{i=1}^{n}\alpha_{ri}A_{i}\rightleftharpoons\sum_{j=1}^{n}\beta_{rj}A_{j}\;\;(r=1,\ldots,m)\,,$
(9)

(formally, the transformation from (4) to (9) is always possible if we allow
some of $k_{r}^{-}=0$ when there is no reverse reaction in the initial
mechanism).

We use notation $k_{r}^{+}$ and $k_{r}^{-}$ for reaction rate constants of
direct and reverse reactions, and $w_{r}^{+}$, $w_{r}^{-}$ for reaction rates
of these reactions:

$w_{r}^{+}=k_{r}^{+}\prod_{i=1}^{n}c_{i}^{\alpha_{ri}},\;\;w_{r}^{-}=k_{r}^{-}\prod_{i=1}^{n}c_{i}^{\beta_{ri}},$
(10)

For the total rate of the couple of direct and reverse reaction we use
$w_{r}=w_{r}^{+}-w_{r}^{-}$. In these notations the kinetic equations are the
same (6) (with different value of $m$).

The principle of detailed balance for the mass action law is: there exists a
positive equilibrium $c_{i}^{\rm eq}>0$ with detailed balance,

$w_{r}^{+}(c_{i}^{\rm eq})=w_{r}^{-}(c_{i}^{\rm eq})\mbox{ for all }r.$ (11)

If the set of the stoichiometric vectors $\\{\gamma_{r}\\}$ is linearly
dependent then this condition implies algebraic relations between reaction
rate constants. Each elementary reaction is equilibrated at the point of
detailed balance, $c^{\rm eq}$. For systems that obey the mass action law,
this means that

$k_{r}^{+}\prod_{i=1}^{n}(c_{i}^{\rm
eq})^{\alpha_{ri}}=k_{r}^{-}\prod_{i=1}^{n}(c_{i}^{\rm
eq})^{\beta_{ri}}=w_{r}^{\rm eq}.$ (12)

It is convenient to use the detailed balance relations (12) and introduce a
set of independent parameters instead of the reaction rate constants: the
equilibrium fluxes $w_{r}^{\rm eq}$ and the equilibrium concentrations
$c_{i}^{\rm eq}$. The reaction rate constants have a simple and explicit
expression through the equilibrium flows and equilibrium concentrations:

$k_{r}^{+}=w_{r}^{\rm eq}\left(\prod_{i=1}^{n}(c_{i}^{\rm
eq})^{\alpha_{ri}}\right)^{-1},\;k_{r}^{-}=w_{r}^{\rm
eq}\left(\prod_{i=1}^{n}(c_{i}^{\rm eq})^{\beta_{ri}}\right)^{-1}.$ (13)

In this parameterization, the mass action law (10) with detailed balance
condition takes the form:

$w_{r}^{+}=w_{r}^{\rm eq}\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\alpha_{ri}}\,w_{r}^{-}=w_{r}^{\rm
eq}\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm eq}}\right)^{\beta_{ri}},$
(14)

In this parameterization, the kinetic equations for the mass action law with
detailed balance are

$\frac{{\mathrm{d}}c}{{\mathrm{d}}t}=\sum_{r=1}^{m}\gamma_{r}w_{r}^{\rm
eq}\left[\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\alpha_{ri}}-\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\beta_{ri}}\right].$ (15)

The classical $H$-theorem for isochoric isothermal perfect chemical systems
with detailed balance can be produced now by simple straightforward
calculations. Notice that

$\frac{\partial H}{\partial c_{i}}=\ln\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)$ (16)

and calculate the time derivative of $H$ according to the system of kinetic
equations with detailed balance condition (15):

$\frac{{\mathrm{d}}H}{{\mathrm{d}}t}=-\sum_{r=1}^{m}(w_{r}^{+}-w_{r}^{-})(\ln
w_{r}^{+}-\ln w_{r}^{-})\leq 0,$ (17)

and ${{\mathrm{d}}H}/{{\mathrm{d}}t}=0$ for points of detailed balance an only
for such points. Therefore, existence of a positive point of detailed balance
for the mass action law kinetic equation implies that $H$ (2) does not
increase in time and all the non-negative steady states of kinetic equations
are detailed balance equilibria.

Definition of $H$ and $c_{i}^{\rm eq}$ requires additional comments.
Equilibrium concentrations of chemical mixture depend on the conserved
quantities those do not change in the course of chemical reactions. For
example, they depend on the atomic balances. Without fixing these values, the
vector of positive equilibria $c_{i}^{\rm eq}$ ($i=1,\ldots,n$) is defined
ambiguously. Any vector of positive equilibria $c_{i}^{\rm eq}>0$ can be used
in the definition of $H$ (2). After that, the $H$-function (2) can be used for
all values of $c_{i}\geq 0$ for all possible values of conserved quantities.
It is a simple exercise to show that for the mass action law kinetics with
detailed balance the difference between $H$-functions with different choices
of equilibria $c_{i}^{\rm eq}$ does not change in time.

The $H$-function (2) was utilised by Zeldovich in his proof of uniqueness of
positive equilibrium for given values of conserved quantities (1938, reprinted
in 1996 [23]). It was recognised as the main instrument for analysis of
stability of perfect kinetic systems in 1960s-1970s [24, 25].

###  2.3 Complex balance and Lyapunov functions

In 1887, Lorentz stated that the collisions of polyatomic molecules are
irreversible and, therefore, Boltzmann’s $H$-theorem is not applicable to the
polyatomic media [26]. Boltzmann found the solution immediately and invented
what we call now semidetailed balance or cyclic balance or complex balance
[2].

Now, it is proven that the Lorentz objections were wrong and the detailed
balance conditions hold for polyatomic molecules [27]. Nevertheless, this
discussion was seminal. The complex balance is a popular assumption in
chemical kinetics beyond the detailed balance [28, 29]. The comparative
analysis of detailed and complex balance assumption in practice of modeling of
chemical reaction networks was presented in work [30]. It is demonstrated how
the generalized mass action law with complex balance appears as a macroscopic
limit of the microscopic Markov kinetics [31]. The formal structures of
complex balance are also useful for analysis of systems with time delays [32].

Return to the ‘irreversible’ representation of the reaction mechanism (4),
where direct and reverse reaction are considered separately. Let $c^{\rm eq}$
be a positive vector of concentrations (it will be a point of complex
balance). Similarly to (13) introduce new parameters

$k_{r}=w_{r}^{\rm eq}\left(\prod_{i=1}^{n}(c_{i}^{\rm
eq})^{\alpha_{ri}}\right)^{-1}$ (18)

$w_{r}^{\rm eq}$ is the rate of the $r$th reaction at the point $c^{\rm eq}$.
The kinetic equations with these parameters are

$\frac{{\mathrm{d}}c}{{\mathrm{d}}t}=\sum_{r=1}^{m}\gamma_{r}w_{r}^{\rm
eq}\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm eq}}\right)^{\alpha_{ri}}.$
(19)

Let us calculate the time derivative of the function $H$ (2) by virtue of the
system (19).

$\frac{{\mathrm{d}}H}{{\mathrm{d}}t}=\sum_{i}\frac{\partial H}{\partial
c_{i}}\frac{{\mathrm{d}}c_{i}}{{\mathrm{d}}t}=\sum_{r}(\gamma_{r},\nabla
H)w_{r}^{\rm eq}\exp(\alpha_{r},\nabla H),$ (20)

where $\nabla H$ is given by (16)

An auxiliary function $\theta(\lambda)$ of one variable $\lambda\in[0,1]$ is
convenient for analysis of ${\mathrm{d}}H/{\mathrm{d}}t$ (see [15, 31]):

$\theta(\lambda)=\sum_{\rho}w_{r}^{\rm eq}\exp[(\nabla
H,(\lambda\alpha_{r}+(1-\lambda)\beta_{r}))]$ (21)

With this function, ${{\mathrm{d}}H}/{{\mathrm{d}}t}$ defined by (20) has a
very simple form:

$\frac{{\mathrm{d}}H}{{\mathrm{d}}t}=-\left.\frac{{\mathrm{d}}\theta(\lambda)}{{\mathrm{d}}\lambda}\right|_{\lambda=1}$
(22)

The condition $\theta(1)=\theta(0)$ is sufficient for the inequality
$\theta^{\prime}(1)\leq 0$, because $\theta(\lambda)$ is a convex function.
Hence, if $\theta(1)=\theta(0)$ then ${{\mathrm{d}}H}/{{\mathrm{d}}t}\leq 0$.
The explicit form of the condition $\theta(1)=\theta(0)$ is

$\sum_{r}w_{r}^{\rm eq}\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\alpha_{ri}}=\sum_{r}w_{r}^{\rm
eq}\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm eq}}\right)^{\beta_{ri}}.$
(23)

Let us consider the family of all gain and loss vectors
$\\{\alpha_{r},\beta_{r}\\}$ ($r=1,\ldots,m$). Usually, some of these $2m$
vectors coincide. Assume that there are $q$ different vectors among them. Let
$y_{1},\ldots,y_{q}$ be these vectors. For each $j=1,\ldots,q$ we take

$R_{j}^{+}=\\{r|\,\alpha_{r}=y_{j}\\}\,,\;R_{j}^{-}=\\{r|\,\beta_{r}=y_{j}\\}$
(24)

We can rewrite (23) in the form

$\sum_{j=1}^{q}\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{y_{ji}}\left[\sum_{r\in R_{j}^{+}}w_{r}^{\rm eq}-\sum_{r\in
R_{j}^{-}}w_{r}^{\rm eq}\right]=0$ (25)

The monomials

$\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm eq}}\right)^{y_{ji}}$

are linearly independent functions of $c$ for every finite set of vectors
$y_{j}$. Therefore for any $j=1,\ldots,q$

$\sum_{r\in R_{j}^{+}}w_{r}^{\rm eq}-\sum_{r\in R_{j}^{-}}w_{r}^{\rm eq}=0$
(26)

This complex balance condition is equivalent to the condition
$\theta(0)=\theta(1)$ and is sufficient for the inequality
${{\mathrm{d}}H}/{{\mathrm{d}}t}\leq 0$. Therefore, $H$ is a Lyapunov function
for systems with complex balance.

##  3 General $H$ theorem for perfect systems with detailed and complex
balance

The systems with detailed or complex balance have a classical Lyapunov
function (2) but we are interested in construction of new Lyapunov functions.
In this section, we find necessary and sufficient conditions that a convex
function is a conditionally universal Lyapunov function for the reaction
network with given reaction mechanism and equilibrium point under detailed or
complex balance assumption.

Let us start from the systems with detailed balance. By definition, a function
$G(c,c^{\rm eq})$ is a conditionally universal Lyapunov function for this
reaction network if it is a Lyapunov function for system (15) for every set of
non-negative values of the equilibrium fluxes $w_{r}^{\rm eq}$.

The following lemma is a simple consequence of the form of the kinetic
equations (15). Let $G(c,c^{\rm eq})$ be continuous function and a convex
function of $c$ for all $c^{\rm eq}$.

######  Lemma 1.

Time derivative of a function $G(c,c^{\rm eq})$ by virtue of system (15) is
non-negative, ${\mathrm{d}}G(c,c^{\rm eq})/{\mathrm{d}}t\leq 0$, for all
positive values of $w_{r}^{\rm eq}$, $c$, and $c^{\rm eq}$ if and only if the
derivatives of this function by virtue of the following systems are non-
negative for all $r$, all positive $c$ and given $c^{\rm eq}$:

$\frac{{\mathrm{d}}c}{{\mathrm{d}}t}=\gamma_{r}\left[\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\alpha_{ri}}-\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\beta_{ri}}\right].$ (27)

###### Proof.

For smooth $G(c,c^{\rm eq})$, calculate the time derivative of this function
by virtue of the system (15):

$\begin{split}\frac{{\mathrm{d}}G(c,c^{\rm
eq})}{{\mathrm{d}}t}&=\sum_{r=1}^{m}w_{r}^{\rm eq}(\nabla_{c}G(c,c^{\rm
eq}),\gamma_{r})\left[\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\alpha_{ri}}-\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\beta_{ri}}\right]\\\ &=\sum_{r=1}^{m}w_{r}^{\rm
eq}D_{r}G(c,c^{\rm eq}),\end{split}$ (28)

where $(\nabla_{c}G(c,c^{\rm eq}),\gamma_{r})$ is the standard inner product
(derivative of $G(c,c^{\rm eq})$ in the direction $\gamma_{r}$) and
$D_{r}G(c,c^{\rm eq})$ is the time derivative of $G(c,c^{\rm eq})$ by virtue
of the system (27).

The coefficients $w_{r}^{\rm eq}$ are independent non-negative variables.
Therefore, the time derivative of $G(c,c^{\rm eq})$ by virtue of the original
system (15) is a conical combination of $D_{r}G(c,c^{\rm eq})$. Non-negativity
of all conical combinations of $D_{r}G(c,c^{\rm eq})$ means that each term is
non-negative. Thus, non-negativity of ${{\mathrm{d}}G(c,c^{\rm
eq})}/{{\mathrm{d}}t}$ for the reaction network with any values of equilibrium
fluxes is equivalent to non-negativity of ${{\mathrm{d}}G(c,c^{\rm
eq})}/{{\mathrm{d}}t}$ for each single reaction subsystem of the network (i.e.
to the inequality $D_{r}G(c,c^{\rm eq})\geq 0$ for all $r$). For continuous
convex functions $G(c,c^{\rm eq})$ we have to use in (32) the subgradients of
$G(c,c^{\rm eq})$ instead of the gradients. (Recall that according to A.D.
Alexandrov theorem continuous convex functions are almost everywhere twice
differentiable and, therefore, the criterion (32) with classical gradients is
valid for them almost everywhere.) The time derivative in (28) becomes an
interval of non-positive numbers and the proof remains the same. For more
detail about continuous but non-smooth Lyapunov functions we refer to [33, 34]
∎

Lemma 1 allows us to reduce a complex validation of inequalities with $m+2n$
variables $w_{r}^{\rm eq}$, $c_{i}^{\rm eq}$, and $c_{i}$ to a series of $m$
simpler inequalities with $2n$ variables. Moreover, non-negativity of
$D_{r}G(c,c^{\rm eq})$ means that the function $G(c,c^{\rm eq})$ does not
increase with time along solutions of one-dimensional kinetic equations (27):
in this system, ${{\mathrm{d}}c}/{{\mathrm{d}}t}$ is proportional to vector
$\gamma_{r}$ and for any positive solution $c(t)$ ($t>0$) the difference
$c(t_{1})-c(t_{0})$ is always proportional to $\gamma_{r}$,
$c(t_{1})-c(t_{0})=\xi\gamma_{r}$ with some scalar multiplier $\xi$
($t_{0,1}>0$). This one-dimensional dynamics admits a strongly convex Lyapunov
functions $H$ (2) with time derivative $-(w_{r}^{+}-w_{r}^{-})(\ln
w_{r}^{+}-\ln w_{r}^{-})$ (17).

For each positive concentration vector $c$, the interval

$I=(c+\mathbb{R}\gamma_{r})\cap\mathbb{R}_{>0}^{n}$ (29)

is positively invariant with respect to (27). Restriction of $H$ (2) on this
interval is a strongly convex function. The derivative of this function has
the logarithmic singularity at the ends of the interval. $H$ has the unique
minimizer on $I$. It is the unique positive equilibrium point of (27) on $I$
(and, by definition, a partial equilibrium of the complete system (15)). These
observations allow us to formulate the following criterion for the
conditionally universal convex Lyapunov functions of the reaction kinetics
(15) for reaction networks with detailed balance.

Let $G(c,c^{\rm eq})$ be continuous function and a convex function of
$c\in\mathbb{R}_{>0}^{n}$ for all $c^{\rm eq}$.

######  Theorem 1.

$G(c,c^{\rm eq})$ is a monotonically non-increasing function of time on the
positive solutions of the kinetic equations (15) for all non-negative values
of equilibrium fluxes if and only if for every positive concentration vector
$c$ and every stoichiometric vector $\gamma_{r}$ of the reaction mechanism the
minimizer of $H(c,c^{\rm eq})$ on the interval $I$ (29) is, at the same time,
a minimizer of $G(c,c^{\rm eq})$ on this interval:

$\underset{{c+\gamma_{r}x\in\mathbb{R}_{>0}^{n}}}{\operatorname{argmin}}H(c+\gamma_{r}x,c^{\rm
eq})\subset\underset{{c+\gamma_{r}x\in\mathbb{R}_{>0}^{n}}}{\operatorname{argmin}}G(c+\gamma_{r}x,c^{\rm
eq})$ (30)

###### Proof.

According to Lemma 1, it is necessary and sufficient to prove this theorem for
the one-step reaction ($m=1$). Consider restriction of the one-step kinetic
equation (27) on the interval $I$ (29). On this interval, the system has one
equilibrium (the partial equilibrium $c_{\gamma_{r}}^{*}$). It is stable, and
the restriction of $H$ on this interval is the Lyapunov function of the
system. The point $c_{\gamma_{r}}^{*}$ is the minimizer of $H$ on $I$.

Assume that $c_{\gamma_{r}}^{*}$ is a minimizer of $G(c,c^{\rm eq})$ for $c\in
I$. Then convexity of $G(c)$ implies that $G$ is monotonically non-increasing
function of time on the solution $c(t)$ of (27) on $I$. ($G(c(t))$ decreases
monotonically to the minimal value when the point $c(t)$ approaches its
minimizer $c_{\gamma_{r}}^{*}$.)

Assume now that $G(c(t))$ does not increase in time due to dynamic of (27) on
$I$. This dynamics lead $c(t)$ to the unique equilibrium $c_{\gamma_{r}}^{*}$.
This equilibrium should be a minimizer of $G$ on $I$. Indeed, if at some point
$y\in I$ the function $G$ takes smaller value then in $c_{\gamma_{r}}^{*}$,
then in the motion from $y$ to $c_{\gamma_{r}}^{*}$ the value of $G$ should
increase, which contradicts the assumption. ∎

Thus, to check that a convex function is a conditionally universal function
for the reaction network with detailed balance, it is sufficient to check that
its minimizers in the direction of the stoichiometric vectors of the reaction
mechanism include the minimizers of $H$ (i.e. one-step partial equilibria). Of
course, this is a much simpler task than analysis of the signs of
${\mathrm{d}}G/{\mathrm{d}}t$ for all states and all values of parameters.
Nevertheless, this simple check gives necessary and sufficient conditions for
a function to be a conditionally universal Lyapunov function for the kinetic
equations with a given reaction mechanism.

A cone of possible velocities is a convenient tool for analysis of
conditionally universal Lyapunov functions for reaction networks. Consider all
systems with detailed balance, a given reaction mechanism and a positive
equilibrium $c^{\rm eq}$. According to the mass action law kinetic equations
with detailed balance (15) the time derivatives ${\mathrm{d}}c/{\mathrm{d}}t$
at given point $c$ form a cone $\mathbb{Q}_{\rm DB}(c)$:

$\mathbb{Q}_{\rm
DB}(c)=\left.cone\left\\{\gamma_{r}sign\left[\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\alpha_{ri}}-\prod_{i=1}^{n}\left(\frac{c_{i}}{c_{i}^{\rm
eq}}\right)^{\beta_{ri}}\right]\;\right|r=1,\ldots,m\right\\},$ (31)

where $cone$ stands for the conic hull and $sign(x)=1$ if $x>0$, $-1$ if $x<0$
and 0 if $x=0$.

######  Remark 1 (Differential criterion of conditionally universal Lyapunov
functions).

A smooth function $G(c,c^{\rm eq})$ is a conditionally universal Lyapunov
functions for a given reaction mechanism and positive equilibrium $c^{\rm eq}$
if

$(\nabla_{c}G(c,c^{\rm eq}),x)\leq 0\mbox{ for all non-negative }c\mbox{ and
all }x\in\mathbb{Q}_{\rm DB}(c).$ (32)

For continuous convex functions we have to use in (32) the subgradients of
$G(c,c^{\rm eq})$ instead of gradients (see also the proof of Lemma 1).

The systems of complex balance are defined by linear relationships (26)
between non-negative values $w_{r}^{\rm eq}$. Therefore, the time derivatives
${\mathrm{d}}c/{\mathrm{d}}t$ at given point $c$ for all systems with complex
balanced equilibrium $c^{\rm eq}$ also form a cone, $\mathbb{Q}_{\rm CB}(c)$.

It is obvious that $\mathbb{Q}_{\rm DB}\subseteq\mathbb{Q}_{\rm CB}$.
Surprisingly, these cones coincide [35, Theorem 2 (Local equivalence of
detailed and complex balance)]. This means that for every mass action law
system with the complex balanced positive equilibrium and any concentration
vector $c$ there exists a mass action law system with detailed balance and the
same positive equilibrium such that the velocity vectors
${\mathrm{d}}c/{\mathrm{d}}t$ at point $c$ for these systems coincide. If some
reactions of the complex balance system are irreversible, a reverse reaction
should be added. Its rate constant will be zero for the original complex
balance system and non-zero for the detailed balance system. This theorem was
proven in [35] even for more general kinetic law, the generalized mass action
law. A bit earlier, such theorem was proven for Markov processes [36].
Continuous-time Markov kinetics with a finite number of states and a given
positive equilibrium was studied. This class of systems is significantly wider
than the systems with detailed balance. Nevertheless, for an arbitrary
probability distribution $P$ and a general system there exists a system with
detailed balance and the same equilibrium that has the same velocity
${\mathrm{d}}P/{\mathrm{d}}t$ at point $P$. The results are extended to
nonlinear systems with the generalized mass action law.

######  Remark 2 (Coincidence of conditionally universal Lyapunov functions
for detailed and complex balance systems).

According to [35, Theorem 2], a function $G(c,c^{\rm eq}$ is a Lyapunov
function for all systems (6) with given reaction mechanism and complex
balanced equilibrium $c^{\rm eq}$ if and only if it is the universal Lapunov
functions for all systems (15) with detailed balance, the same equilibrium and
the same mechanism (supplemented by the reverse reactions, if necessary).

For some reaction mechanisms there exist Lyapunov functions without any
relation to detailed balance or complex balance. For example, assume that all
the elementary reactions have the form

$\alpha_{ri}A_{i}\to\sum_{j=1}^{n}\beta_{rj}A_{j}$ (33)

(only one $\alpha_{ri}$ can be non-zero; direct and reverse reactions are
considered separately and some reactions can be irreversible). If there exists
a positive balance $\sum_{i}m_{i}c_{i}=M=const$ then for any two solutions of
the kinetic equations $c^{1}(t)$, $c^{2}(t)$ with the same value of $M$ the
weighted $l_{1}$ distance between them
$\sum_{i}m_{i}|c^{1}_{i}(t)-c^{2}_{i}(t)|$ monotonically decreases [37].

Convergent dynamics with quadratic Lyapunov norms $\|x\|^{2}=(x,Px)$, where
$P$ is a symmetric positive definite matrix, was studied by Demidovich in
1960s and widely used [38]. Systems (33) give us example of a class of non-
linear convergent systems in weighted $l_{1}$ norm. Some other examples of
reaction mechanisms with such convergence property were produced in [39, 40]
on the basis of monotonicity idea. All these selected mechanisms are rather
simple. They have convergence property for any values of reaction rate
constants. The monotonicity idea is useful for selection of reaction kinetic
equations with stable dynamics without detailed or complex balance properties
[39, 41].

On the contrary, in this paper we consider reaction networks with an arbitrary
(presumably, nonlinear) reaction mechanism but with specific restrictions on
the reaction rate constants. They should obey the principle of detailed
balance. (Later on we explain why the same results are valid for systems with
the so-called complex balance.) For such systems, there exist thermodynamic
Lyapunov functions. For perfect systems under isothermal isochoric conditions
the explicit form of this function is presented by (2). The situation with
linear kinetics was similar when Rényi revealed $f$-divergences for Markov
chains. The decrease of relative entropy (information) in time was well-known
but there were no other Lyapunov functions until Rényi work [7]. Below we
construct a wide family of additional Lyapunov functions for any nonlinear
reaction network, obeying the mass action law and the principle of detailed
balance.

##  4 Partial equilibria and new Lyapunov functions for mass action law

Let vector $\gamma$ have both positive and negative components. For every
vector of concentrations $c$ we define the corresponding partial equilibrium
in direction $\gamma$ as

$c_{\gamma}^{*}(c)=\underset{{c+\gamma
x\in\mathbb{R}_{>0}^{n}}}{\operatorname{argmin}}H(c+\gamma x).$ (34)

This partial equilibrium $c_{\gamma}^{*}(c)$ is the minimizer of $H$ on the
interval

$(c+\mathbb{R}\gamma)\cap\mathbb{R}_{>0}^{n}.$

This interval is bounded. For a positive point $c$ the minimizer
$c_{\gamma}^{*}(c)$ is also positive. This is an elementary consequence of the
logarithmic singularity of $(c\ln c)^{\prime}$ at zero. Here and below, argmin
is the set of points where the function gets its minimum. The functions $H(c)$
is strongly convex on each bounded set because its Hessian has the form

$\frac{\partial^{2}H(c)}{\partial c_{i}\partial
c_{j}}=\frac{1}{c_{i}}\delta_{ij},$

where $\delta_{ij}$ is the Kronecker delta. Therefore, each the argmin set in
(34) consists of one point.

For monomolecular and for bimolecular reactions there are simple analytic
expression for partial equilibria. Consider a monomolecuar reaction
$A_{i}\rightleftharpoons A_{j}$. The non-zero components of the stoichiometric
vector $\gamma$ are: $\gamma_{i}=-1$, $\gamma_{j}=1$. For a given vector $c$,
the partial equilibrium is given by the equation
$k^{+}c_{i}^{*}(c)=k^{-}c_{j}^{*}(c)$ under condition that
$c^{*}(c)=c+x\gamma$. Simple algebra gives:

$\begin{split}&c_{i}^{*}(c)=\frac{k^{-}}{k^{+}+k^{-}}\left(c_{i}+c_{j}\right);\\\
&c_{j}^{*}(c)=\frac{k^{+}}{k^{+}+k^{-}}\left(c_{i}+c_{j}\right).\end{split}$
(35)

Other components of $c^{*}(c)$ coincide with those of $c$. The sum
$c_{i}+c_{j}=b$ does not change in the reaction $A_{i}\rightleftharpoons
A_{j}$. Rewrite (35) using this ‘partial balance’ $b$:

$c_{i}^{*}(c)=\frac{k^{-}b}{k^{+}+k^{-}};\;\;c_{j}^{*}(c)=\frac{k^{+}b}{k^{+}+k^{-}}.$
(36)

The ‘rate constants’ in (35) and (36) can be defined through a positive
equilibrium point $c^{\rm eq}$: for the linear reaction,
$A_{i}\rightleftharpoons A_{j}$, $k^{+}c_{i}^{\rm eq}=k^{-}c_{j}^{\rm eq}$ and
we can take, for example,

$k^{+}=\frac{c_{j}^{\rm eq}}{c_{i}^{\rm eq}+c_{j}^{\rm
eq}};\;k^{-}=\frac{c_{i}^{\rm eq}}{c_{i}^{\rm eq}+c_{j}^{\rm eq}}$

(we use the normalization condition $k^{+}+k^{-}=1$ to select one solution
from the continuum of proportional sets of constants). The expression for the
partial equilibrium for the linear reaction (36) is

$c_{i}^{*}(c)=\frac{c_{i}^{\rm eq}b}{c_{i}^{\rm eq}+c_{j}^{\rm
eq}};\;\;c_{j}^{*}(c)=\frac{c_{j}^{\rm eq}b}{c_{i}^{\rm eq}+c_{j}^{\rm eq}}.$
(37)

For a bimolecular reaction $A_{i}+A_{j}\rightleftharpoons A_{k}$ the non-zero
components of the stoichiometric vector $\gamma$ are:
$\gamma_{i}=\gamma_{j}=-1$, $\gamma_{k}=1$. Two independent ‘partial balances’
that do not change in the reaction are:

$b_{1}=c_{i}+c_{j}+2c_{k},\;\;b_{2}=c_{i}-c_{j}.$

The partial equilibrium $c^{*}(c)$ is the positive solution of the equation
$k^{+}c_{i}^{*}(c)c_{j}^{*}(c)=k^{-}c_{k}^{*}(c)$ under condition that
$c^{*}(c)=c+x\gamma$. After solving of quadratic equation for $x$ we get:

$\begin{split}&c_{i}^{*}(c)=\frac{b_{2}}{2}-\frac{k^{-}}{2k^{+}}+\sqrt{\frac{b_{2}^{2}}{4}+\frac{k^{-}b_{1}}{2k^{+}}+\left(\frac{k^{-}}{2k^{+}}\right)^{2}};\\\
&c_{j}^{*}(c)=-\frac{b_{2}}{2}-\frac{k^{-}}{2k^{+}}+\sqrt{\frac{b_{2}^{2}}{4}+\frac{k^{-}b_{1}}{2k^{+}}+\left(\frac{k^{-}}{2k^{+}}\right)^{2}};\\\
&c_{k}^{*}(c)=\frac{b_{1}}{2}+\frac{k^{-}}{2k^{+}}-\sqrt{\frac{b_{2}^{2}}{4}+\frac{k^{-}b_{1}}{2k^{+}}+\left(\frac{k^{-}}{2k^{+}}\right)^{2}}.\end{split}$
(38)

The signs in front of square root are selected to provide positivity of
$c^{*}(c)$.

We can rewrite (38) using any positive equilibrium $c^{\rm eq}$. Indeed, at
this point, the detailed balance gives $k^{+}c^{\rm eq}_{i}c^{\rm
eq}_{j}=k^{-}c^{\rm eq}_{k}$ and

$\frac{k^{-}}{k^{+}}=\frac{c^{\rm eq}_{i}c^{\rm eq}_{j}}{c^{\rm eq}_{k}}.$

Therefore,

$\begin{split}&c_{i}^{*}(c)=\frac{b_{2}}{2}-\frac{c^{\rm eq}_{i}c^{\rm
eq}_{j}}{2c^{\rm eq}_{k}}+\sqrt{\frac{b_{2}^{2}}{4}+b_{1}\frac{c^{\rm
eq}_{i}c^{\rm eq}_{j}}{2c^{\rm eq}_{k}}+\left(\frac{c^{\rm eq}_{i}c^{\rm
eq}_{j}}{2c^{\rm eq}_{k}}\right)^{2}};\\\
&c_{j}^{*}(c)=-\frac{b_{2}}{2}-\frac{c^{\rm eq}_{i}c^{\rm eq}_{j}}{2c^{\rm
eq}_{k}}+\sqrt{\frac{b_{2}^{2}}{4}+b_{1}\frac{c^{\rm eq}_{i}c^{\rm
eq}_{j}}{2c^{\rm eq}_{k}}+\left(\frac{c^{\rm eq}_{i}c^{\rm eq}_{j}}{2c^{\rm
eq}_{k}}\right)^{2}};\\\ &c_{k}^{*}(c)=\frac{b_{1}}{2}+\frac{c^{\rm
eq}_{i}c^{\rm eq}_{j}}{2c^{\rm
eq}_{k}}-\sqrt{\frac{b_{2}^{2}}{4}+b_{1}\frac{c^{\rm eq}_{i}c^{\rm
eq}_{j}}{2c^{\rm eq}_{k}}+\left(\frac{c^{\rm eq}_{i}c^{\rm eq}_{j}}{2c^{\rm
eq}_{k}}\right)^{2}}.\end{split}$ (39)

Similar formulas can be easily obtained for the bimolecular reactions
$A_{i}+A_{j}\rightleftharpoons A_{k}+A_{l}$ and $A_{i}+A_{j}\rightleftharpoons
2A_{k}$. They require nothing more than the quadratic formula and the detailed
balance condition.

Partial equilibria appear in non-equilibrium thermodynamics from the very
beginning. Already Jaynes considered conditional maximization of entropy as a
basic method of equilibrium and non-equilibrium statistical physics [42].
According to Grmela, the equilibrium and nonequilibrium thermodynamics as well
as the equilibrium and nonequilibrium statistical mechanics should be
considered as particular representations of the Dynamical Maximum Entropy
Principle [43]. Conditional maximization of entropy gives thermodynamic basis
to thermodynamic of driven systems [44]. Partial equilibria (or
quasiequilibria - conditional maximizers of entropy for more general
conditions) are often considered as constrained equilibria as a result of
introducing external or internal constraints. For example, the models of
extended irreversible thermodynamics can be produced by the conditional
maximization of entropy subject to various dynamical constrains and hypotheses
about slow and fast variables [45]. Methods of invariant manifolds produces
dynamical correction to these conditional maximum entropy models [46].

In this work, we use a very simple case of partial equilibria: these are
equilibria with constrains under which only one reversible elementary reaction
is on, while other reactions are frozen. To be more precise, we use the term
‘$H$-function’ instead of ‘entropy’.

The partial equilibrium $H$-function is associated with the partial
equilibria:

$H_{\gamma}^{*}(c)=H(c_{\gamma}^{*}(c))=\min_{c+\gamma
x\in\mathbb{R}_{>0}^{n}}H(c+\gamma x).$ (40)

We do not assume here any time separation and constrained dynamics and use
partial equilibria in a completely different way. Dynamics of $c(t)$ remains
undeformed and foll;owes the original kinetic equations. For each
stoichiometric vector $\gamma$ we consider the partial equilibrium
$c^{*}_{\gamma}(c(t))$ as a projection of the genuine dynamics on the
hypersurface of partial equilibria. This set of projections is a ‘shadow’ of
kinetic curves. The new Lyapunov function for the genuine kinetics is combined
from the entropies of these projections. For every compact set $\Gamma$ of
vectors $\gamma$ with both positive and negative components, we define

$H_{\Gamma}^{*}(c)=\max_{\gamma\in\Gamma}H_{\gamma}^{*}(c).$ (41)

This definition has a simple explanation: for a given initial positive
concentration vector $c$ and each stoichiometric vector $\gamma\in\Gamma$ we
find the partial equilibrium of the one-step system ($m=1$) with this
stoichiometric vector. This partial equilibrium, $c_{\gamma}^{*}(c)$, is the
projection of the initial vector $c$ parallel to the vector $\gamma$ onto the
hypersurface of partial equilibria defined by the equation

$(\nabla H,\gamma)=0.$

From all these projections ($\gamma\in\Gamma$) we select the most non-
equilibrium state, i.e., the state with the maximal value of
$H=H(c_{\gamma}^{*}(c))$ (40). This maximal value is the new function
$H_{\Gamma}^{*}(c)$ (41). It is defined by the set $\Gamma$ that should
include all the stoichiometric vectors of the reaction mechanism.

It is necessary to stress that for every reaction mechanism there exists
continuum of compact sets $\Gamma$, which include the stoichiometric vectors
of this mechanism. The partial equilibrium is the same for the vector $\gamma$
and $x\gamma$ for any $x\neq 0$. Therefore, rigorously speaking, the functions
$H_{\Gamma}^{*}(c)$ should be indexed by subsets of the projective space (the
space of one-dimensional subspaces), and not by sets of vectors.

The hypersurface of partial equilibria for a given vector $\gamma$ is the
hyperplane orthogonal to $\gamma$ in the entropic inner product. The level
sets of the partially equilibrium function $H_{\gamma}^{*}(c)$ are cylinders
with the spheric base and axis parallel to $\gamma$.

Fig. 1 represents a very simplified example of the partial equilibrium
$c_{\gamma}^{*}(c)$ and the Lyapunov function $H_{\Gamma}^{*}(c)$
construction. The system with three components of the same molecular weight is
presented in the triangle $c_{1}+c_{2}+c_{3}=const$ drawn in barycentric
coordinates. The reaction mechanism consists of three reactions
$A_{1}\rightleftharpoons A_{2}$,$A_{2}\rightleftharpoons A_{3}$, and
$2A_{1}\rightleftharpoons A_{2}+A_{3}$. The equilibrium $c^{\rm eq}$ is
assumed in the center of the triangle ($c^{\rm eq}_{1}=c^{\rm eq}_{2}=c^{\rm
eq}_{3}$). The partial equilibria of the first two reactions form the straight
lines in the triangle, the medians, while the partial equilibria of the non-
linear reaction form a parable $c_{2}c_{3}/c_{1}^{2}=1$. These three lines
intersect in the equilibrium due to detailed balance.

Explicit expressions for the partial equilibria $c^{*}_{\gamma}(c)$ are:

  * 1.

For the reaction $A_{1}\rightleftharpoons A_{2}$

$c^{*}_{1,\gamma_{1}}(c)=c^{*}_{2,\gamma_{1}}(c)=\frac{1}{2}(c_{1}+c_{2}),\;c^{*}_{3,\gamma_{1}}(c)=c_{3};$
(42)

  * 2.

For the reaction $A_{2}\rightleftharpoons A_{3}$

$c^{*}_{1,\gamma_{2}}(c)=c_{1},\;c^{*}_{2,\gamma_{2}}(c)=c^{*}_{3,\gamma_{2}}(c)=\frac{1}{2}(c_{2}+c_{3});$
(43)

  * 3.

For the reaction $2A_{1}\rightleftharpoons A_{2}+A_{3}$

$\begin{split}&c_{1,\gamma_{3}}^{*}(c)=\frac{1}{3}\left(-b_{1}+\sqrt{4b_{1}^{2}-3b_{2}^{2}}\right);\\\
&c_{2,\gamma_{3}}^{*}(c)=\frac{1}{6}\left(4b_{1}+3b_{2}-\sqrt{4b_{1}^{2}-3b_{2}^{2}}\right);\\\
&c_{3,\gamma_{3}}^{*}(c)=\frac{1}{6}\left(4b_{1}-3b_{2}-\sqrt{4b_{1}^{2}-3b_{2}^{2}}\right),\end{split}$
(44)

where $b_{1}=c_{1}+c_{2}+c_{3}$ and $b_{2}=c_{2}-c_{3}$ are the independent
‘partial balances’ for this reaction.

Fig. 1a shows the partial equilibria for an arbitrarily selected point $c$. In
Fig. 1b, one level set of $H_{\Gamma}^{*}$ is presented, where
$\Gamma=\\{\gamma_{1},\gamma_{2},\gamma_{3}\\}$:

$H_{\Gamma}^{*}(c)=\max\\{H(c^{*}_{\gamma_{1}}(c)),H(c^{*}_{\gamma_{2}}(c)),H(c^{*}_{\gamma_{3}}(c))\\},$

$c^{*}_{\gamma_{i}}$ $i=1,2,3$ are given by (42), (43), (44) and $H$-function
is the classical one (2).

The sublevel set of $H_{\Gamma}^{*}(c)$ is the intersection of strips with
sides parallel to the stoichiometric vectors $\gamma_{i}$. These strips are
sublevel sets for the partial equilibrium entropies $H_{\gamma_{i}}^{*}(c)$.
In higher dimensions, the level sets of partial equilibrium $H$-function
$H_{\gamma}^{*}(c)$ are cylindrical hypersurfaces with the generatrix parallel
to $\gamma$. The base (or the directrix) of this cylindric surface is the
level set of $H$ on the surface of partial equilibrium. In Fig. 1b, the
‘surfaces’ of partial equilibria are lines, the level sets of $H$ on these
lines are couples of points. These points are highlighted. Note that the
sublevel areas for the new function $H^{*}_{\Gamma}$ in Fig. 1 are convex
polygons, whereas for the classical $H$-function they have smooth border.

Figure 1: The stoichiometric vectors $\gamma_{1},\gamma_{2},\gamma_{3}$ and
the partial equilibria for the reaction mechanism $A_{1}\rightleftharpoons
A_{2}$, $A_{2}\rightleftharpoons A_{3}$, $2A_{1}\rightleftharpoons
A_{2}+A_{3}$. The concentration triangle $c_{1}+c_{2}+c_{3}=b$ is split by the
lines of partial equilibria into six compartments. In each compartment, the
dominated direction of each reaction (towards the partial equilibrium) is
defined unambiguously. a) Partial equilibria (highlighted points) for an
arbitrary positive concentration vector $c$. b) A level set of
$H_{\Gamma}^{*}$ for $\Gamma=\\{\gamma_{1},\gamma_{2},\gamma_{3}\\}$.

Assume that the principle of detailed balance holds and the vector of
equilibrium concentrations $c_{i}^{\rm eq}$ in the definition of $H$-function
(2) is the point of detailed balance.

######  Theorem 2.

Let all the stoichiometric vectors of the reaction mechanism (4) belong to
$\Gamma$ then $H_{\Gamma}^{*}(c)$ (40) is monotonically non-increasing
function along the solutions of the kinetic equations (15) for all positive
$c$ and all non-negative values of equilibrium fluxes $w^{\rm eq}$.

###### Proof.

To use Theorem 1 we have to prove two statements:

  1. 1.

The function $H_{\Gamma}^{*}(c)$ is convex in $\mathbb{R}_{>0}^{n}$;

  2. 2.

For each $\gamma\in\Gamma$ and a positive vector $c$ the minimizer of $H$ on
the interval $(c+\gamma\mathbb{R})\cap\mathbb{R}_{>0}^{n}$ is, at the same
time, a minimizer of $H_{\Gamma}^{*}$ on this interval.

We prove convexity of the function $H_{\Gamma}^{*}(c)$ in
$\mathbb{R}_{>0}^{n}$ in two steps.

  * 1.

Convexity of $H_{\gamma}^{*}(c)$ for one-element sets $\Gamma=\\{\gamma\\}$.

  * 2.

Convexity of $H_{\Gamma}^{*}(c)=\max_{\gamma\in\Gamma}H_{\gamma}^{*}(c)$.

Let us prove convexity of $H_{\Gamma}^{*}(c)$ for one-element sets $\Gamma$,
that is, we will prove convexity of the partial equilibrium $H$-function
$H_{\gamma}^{*}(c)$ (40). By definition, convexity of $H_{\gamma}^{*}(c)$ in
$\mathbb{R}_{>0}^{n}$ means that for each two positive concentration vectors
$c^{1}$ and $c^{2}$ and a number $\lambda\in[0,1]$ the inequality holds:

$H_{\gamma}^{*}(\lambda c^{1}+(1-\lambda)c^{2})\leq\lambda
H_{\gamma}^{*}(c^{1})+(1-\lambda)H_{\gamma}^{*}(c^{2}).$

First, notice that due to the convexity of $H$

$H(\lambda c_{\gamma}^{*}(c^{1})+(1-\lambda)c_{\gamma}^{*}(c^{2}))\leq\lambda
H(c_{\gamma}^{*}(c^{1}))+(1-\lambda)H(c_{\gamma}^{*}(c^{2})).$

Secondly, $H(c_{\gamma}^{*}(c^{1,2}))=H_{\gamma}^{*}(c^{1,2})$ by the
definition of the partial equilibrium $H$-function (40). Therefore, the
previous inequality can be rewritten as

$H(\lambda c_{\gamma}^{*}(N^{1})+(1-\lambda)c_{\gamma}^{*}(N^{2})\leq\lambda
H_{\gamma}^{*}(c^{1})+(1-\lambda)H_{\gamma}^{*}(c^{2}).$

Finally,

$c_{\gamma}^{*}(c^{1,2})\in(c^{1,2}+\gamma\mathbb{R})\cap\mathbb{R}_{>0}^{n},$

hence,

$\lambda c_{\gamma}^{*}(c^{1})+(1-\lambda)c_{\gamma}^{*}(c^{2})\in\lambda
c^{1}+(1-\lambda)c^{2}+\gamma\mathbb{R}$

and

$H(c_{\gamma}^{*}(c^{1})+(1-\lambda)c_{\gamma}^{*}(c^{2}))\geq
H_{\gamma}^{*}(\lambda c^{1}+(1-\lambda)c^{2})$

because the last value is the minimum of $H$ on the interval

$(\lambda c^{1}+(1-\lambda)c^{2}+\gamma\mathbb{R})\cap\mathbb{R}_{>0}^{n}.$

Convexity of the partial equilibrium $H$-function $H_{\gamma}^{*}(c)$ is
proven.

Convexity of $H_{\Gamma}^{*}(c)$ follows from the convexity of the partial
equilibrium $H$-function $H_{\gamma}^{*}(c)$, from the definition of
$H_{\Gamma}^{*}(c)$ as the maximum of $H_{\gamma}^{*}(c)$ ($\gamma\in\Gamma$),
and from the following fact from convex analysis: Maximum of a set of convex
functions is again convex. The shortest proof is based on the definition of a
convex function as a function with convex epigraph [47] and follows from the
observation that the epigraph of the maximum of a family of functions is the
intersection of their epigraphs.

Let us analyse the minimizers of $H_{\Gamma}^{*}$ on the interval
$I=(c+\gamma\mathbb{R})\cap\mathbb{R}_{>0}^{n}$ for a positive concentration
vector $c$. Select $\gamma\in\Gamma$. The minimizer of $H$ on the interval $I$
is the partial equilibrium $c_{\gamma}^{*}(c)$, by the definition (34).
Function $H_{\gamma}^{*}(c)$ is constant on the interval
$I=(c+\gamma\mathbb{R})\cap\mathbb{R}_{>0}^{n}$ and
$H_{\gamma}^{*}(c)=H(c_{\gamma}^{*}(c))$ on $I$. Therefore,
$c_{\gamma}^{*}(c)$ is a minimizer of $H_{\gamma}^{*}(c)$ on $I$ (trivially,
as all other points of $I$ do).

Notice, that for all $\gamma^{\prime}\in\Gamma$

$H_{\gamma}^{*}(c_{\gamma}^{*}(c))\geq
H_{\gamma^{\prime}}^{*}(c_{\gamma}^{*}(c))$ (45)

and the equality is strong if $c_{\gamma}^{*}(c)\neq
c_{\gamma^{\prime}}^{*}(c_{\gamma}^{*}(c))$. Indeed,
$c_{\gamma^{\prime}}^{*}(c_{\gamma}^{*}(c))$ is the unique minimizer of $H$ on
the interval
$(c_{\gamma}^{*}(c)+\gamma^{\prime}\mathbb{R})\cap\mathbb{R}_{>0}^{n}$. If
this minimizer does not coincide with $c_{\gamma}^{*}(c)$ then

$H_{\gamma^{\prime}}^{*}(c_{\gamma}^{*}(c))=H(c_{\gamma^{\prime}}^{*}(c_{\gamma}^{*}(c)))<H(c_{\gamma}^{*}(c))=H_{\gamma}^{*}(c_{\gamma}^{*}(c)).$

According to inequality (45),

$\max_{\gamma^{\prime}\in\Gamma}H_{\gamma^{\prime}}^{*}(c_{\gamma}^{*}(c))=H_{\gamma}^{*}(c_{\gamma}^{*}(c)).$
(46)

For any family of convex functions $\Phi$ on the interval $I$ the following
statement holds. Let $f\in\Phi$. If $y$ is a minimizer of $f(x)$ on $I$,
$F(x)=\max\\{\phi(x)|\phi\in\Phi\\}$, and $F(y)=f(y)$ then $y$ is a minimizer
of $F(x)$ on $I$. Indeed, for any $z\in I$ $f(z)\geq f(y)$ because $y$ is a
minimizer of $f$. At the same time, $F(z)\geq f(z)$ by the definition of $F$.
Hence, for any $z\in I$ $F(z)\geq F(y)$ and $y$ is a minimizer of $F(x)$ on
$I$.

Let us take $f=H_{\gamma}^{*}(c)$ and
$\Phi=\\{H_{\gamma^{\prime}}^{*}(c)|\gamma^{\prime}\in\Gamma\\}$. Then
$F=H_{\Gamma}^{*}(c)$. Select $y=c_{\gamma}^{*}(c)$. Notice that $y$ is a
minimizer of $f$ and $F(y)=f(y)$ according to (46). Therefore,
$c_{\gamma}^{*}(c)$ is a minimizer of $H_{\Gamma}^{*}(c)$ on $I$. By combining
this result with the proven convexity of $H_{\Gamma}^{*}(c)$ and applying
Theorem 1, we prove that $H_{\Gamma}^{*}(c)$ is a monotonically non-increasing
function on solutions of kinetic equations. ∎

######  Corollary 1.

Assume that the complex balance condition holds and the vector of equilibrium
concentrations $c_{i}^{\rm eq}$ in the definition of $H$-function (2) is the
point of complex balance.

Let all the stoichiometric vectors of the reaction mechanism (4) belong to
$\Gamma$ then the time derivative $H_{\Gamma}^{*}(c)$ (41) by virtue (15) is
non-positive for all positive $c$ and all non-negative values of equilibrium
fluxes $w^{\rm eq}$: ${\mathrm{d}}H_{\Gamma}^{*}(c)/{\mathrm{d}}t\leq 0$.

###### Proof.

This corollary follows immediately from Theorem 2 and from [35, Theorem 2
(Local equivalence of detailed and complex balance)] because $\mathbb{Q}_{\rm
DB}(c)=\mathbb{Q}_{\rm CB}(c)$ for all positive $c$. ∎

When $c(t)$ approaches an equilibrium, $c^{\rm eq}$, all the partial
equilibria also converge to $c^{\rm eq}$. Locally, near the equilibrium
$c^{\rm eq}$, this picture can be simplified. Let us use quadratic
approximation to $H(c)$. The entropic inner product is defined as $\langle
x|y\rangle=(x,(D^{2}H)_{c^{\rm eq}}y)$, where $(D^{2}H)_{c^{\rm eq}}$ is the
second differential of $H$ at the equilibrium (Hessian matrix) and $(\ ,\ )$
is the standard inner product.

######  Remark 3 (Universal Lyapunov functions for linearized kinetic
equations).

In the entropic inner product, the quadratic approximation to the classical
$H$-function near equilibrium $c^{\rm eq}$ is

$H(c)-H(c^{\rm eq})\approx\frac{1}{2}\langle c-c^{\rm eq}|c-c^{\rm
eq}\rangle.$

For this approximation, elementary linear algebra gives

$\begin{split}&c^{*}_{\gamma}(c)-c^{\rm eq}=c-c^{\rm eq}-\frac{\langle
c-c^{\rm eq}|\gamma\rangle\gamma}{\langle\gamma|\gamma\rangle};\\\
&H_{\gamma}^{*}(c)-H(c^{\rm eq})=\frac{1}{2}\langle c^{*}_{\gamma}(c)-c^{\rm
eq}|c^{*}_{\gamma}(c)-c^{\rm eq}\rangle=H(c)-H(c^{\rm
eq})-\frac{1}{2}\frac{\langle c-c^{\rm
eq}|\gamma\rangle^{2}}{\langle\gamma|\gamma\rangle};\\\
&H_{\Gamma}^{*}(c)-H(c^{\rm eq})=\frac{1}{2}\langle c-c^{\rm eq}|c-c^{\rm
eq}\rangle-\frac{1}{2}\min_{\gamma\in\Gamma}\left\\{\frac{\langle c-c^{\rm
eq}|\gamma\rangle^{2}}{\langle\gamma|\gamma\rangle}\right\\}.\end{split}$ (47)

Formula (47) gives universal Lyapunov functions for linearized kinetic
equations with the given reaction mechanism and detailed balance conditions.

Thus, the new family of universal Lyapunov functions for chemical kinetic
equations with detailed or complex balance is constructed.

##  5 Outline of possible generalizations and applications

In this section, we briefly outline the possible generalizations and
applications of the main results. This is basically a ‘to do’ section and
detailed analysis of all generalizations and applications is beyond the scope
of the paper.

#### Generalized mass action law

Analysis of the proven statements demonstrate that two properties of $H$ are
used: strict convexity in directions $\gamma$ in positive orthant
$\mathbb{R}_{>0}$ and logarithmic singularity of the $\partial H/\partial
c_{i}$ when $c_{i}\to 0$. One more requirement is the generalized mass action
laws for the reaction rates:

$w_{r}(c)=\phi_{r}\exp(\nabla H,\alpha_{r}),$

where $\exp(\nabla H,\alpha_{r})$ is the ‘Boltzmann factor’, and $\phi_{r}\geq
0$ is a non-negative quantity, ‘kinetic factor’ [31]. Detailed balance and
complex balance conditions for the generalized mass action law have the same
form as for perfect systems. For detailed balance such kinetics was introduced
by Feinberg [48] (Marselin–De Donder kinetics, see also [3]). Grmela studied
properties of geometry of nonlinear non-equilibrium thermodynamics for the
generalized mass action law [49]. Detailed balance in the limit of non-
equilibrium reactions was analyzed in [50]. General analysis of non-classical
entropies and their relations to the second law of thermodynamics was
presented in work [35]. Complex balance conditions for macroscopic kinetics
was proven for systems with Markov microscopic description in [31]. The basic
constructions of partial equilibria and functions $H^{*}_{\Gamma}$ can be
extended to the generalized mass action law.

#### Non-isochoric and non-isothermal conditions

The Lyapunov function (2) and kinetic equations (6) with mass action low (5)
are valid for isochoric isothermal conditions (constant volume and
temperature). For other conditions they should be modified. Such a
modification is rather simple and does not add any substantial change but the
equation and Lyapunov functions have different analytic form (see, for
example, [14, 51]). The analysis of Lyapunov functions for these conditions
can be provided in the formalism of the generalized mass action law. A general
form for the description of non-isothermal reactions in closed chemical
systems in terms of the Marcelin–De Donder kinetics and explicit form of the
$H$-function for the systems with detailed balance under various conditions is
presented in [52].

#### Spatially distributed systems

Transport in Boltzmann’s equation is conservative (the free flight) and
dissipative terms are local (the collision integral). Generalization of our
approach to the collision integral has to be done. In general, models of
complex transport processes can include both dissipative and conservative
terms. There are many attempts to create thermodynamic theory of such
processes. In the GENERIC approach, conservative and dissipative components
are explicitly separated with some commutativity conditions between them [53].
Generalization of proposed construction onto dissipative components of
transport processes seems to be a challenging task. If the dissipative part is
described in the language of quasi-chemical formalism with a finite number of
steps [54] then this generalization is more straightforward.

#### Possible applications

What can add the new Lyapunov functions to research tools? Of course, more
Lyapunov functions are better than less. The non-classical Lyapunov functions
for linear systems are widely used for various estimates and information
analysis in the situations, where linear Markov chains could serve as adequate
models of information transformation. Just for example of various
applications, we can refer to works [55, 56, 57]. Universal Lyapunov functions
are instruments for evaluation of possible dynamics when the reaction rate
constants are unknown or highly uncertain. Without any knowledge of the
reaction mechanism we use the thermodynamic potentials for evaluation of the
attainable sets of chemical reactions (the theory and algorithms are presented
in [15, 19, 58], some industrial applications are discussed in [18, 59, 60]).
Convexity allows us to transform the $n$-dimensional problems about
attainability and attainable sets into an analysis of one-dimensional continua
and discrete objects, thermodynamic trees [19]. When the reaction mechanism is
known, we can use this information about mechanism for sharper estimations. A
new class of estimates is needed and the new Lyapunov functions give a
collection of instruments for such estimates. The use of many Lyapunov
functions in the analysis of attainable sets makes these estimates more narrow
and close to reality.

##  6 Conclusion

###  6.1 New Lyapunov functions: the main result

A rich family of universal Lyapunov functions for any linear or non-linear
reaction network with detailed or complex balance is presented. Consider a
mass action law system with a given reaction mechanism, the set of
stoichiometric vectors $\Upsilon$ and detailed or complex balance. For
construction of the new Lyapunov functions for this system the following
operations are needed.

  1. 1.

Normalize vectors from $\Upsilon$. Let $\Upsilon_{0}$ be the set of normalized
stoichiometric vectors.

  2. 2.

Select a finite set of vectors $\Gamma$ such that $\Upsilon_{0}\subset\Gamma$
and each vector from $\Gamma$ has both positive and negative components.

  3. 3.

For each $\gamma\in\Gamma$ calculate

$H^{*}_{\gamma}(c)=\min_{c+\gamma x\in\mathbb{R}_{>0}^{n}}H(c+\gamma x).$

  4. 4.

Find

$H^{*}_{\Gamma}(c)=\max_{\gamma\in\Gamma}H^{*}_{\gamma}(c).$

Thus, we have to perform several smooth one-dimensional convex minimizations
and then select the maximum of these minima.

Theoretically, $H^{*}_{\Gamma}$ is a Lyapunov function for any compact
$\Gamma$. Such a function can be approximated (from below) by the functions
calculated for finite $\varepsilon$-networks on $\Gamma$.

Functions $H^{*}_{\Gamma}$ are not presented as elementary functions even if
$H$ has explicit analytical form. If the reaction mechanism consists of mono-
and bimolecular reactions then the partial equilibria can be represented by
radicals and $H^{*}_{\gamma}(c)$ can be written as an elementary function but
in general case the less implicit form with min and max operations is
necessary. If we substitute $H$ by its quadratic approximation near
equilibrium then functions $H^{*}_{\Gamma}$ can be represented more explicitly
(47) and the difference between $H$ and $H^{*}_{\Gamma}$ becomes obvious.

The mystery about the fundamental difference between the rich family of
Lyapunov functions for linear networks and a very limited collection of
Lyapunov functions for non-linear networks in thermodynamic conditions is
resolved: there is no such crucial difference anymore.

###  6.2 An open problem

There remains an obvious difference between explicit analytic expression of
$f$-divergences (3) and not so obvious construction of Lyapunov functions for
general networks using partial equilibria of non-linear reactions. For linear
reactions, the partial equilibria have very simple analytic expression, for
bimolecular reactions they are given explicitly using quadratic formula, but
for trimolecular reactions the analytic formulas become too bulky.

Alt least one important question is still open. The new Lyapunov functions
$H_{\Gamma}^{*}$ are, at the same time, universal Lyapunov functions for
linear kinetics, if the stoichiometric vectors of the linear reaction
mechanism $A_{i}\rightleftharpoons A_{j}$ are included in $\Gamma$. Due to the
results of [10, 11], such a function should be, essentially, a $f$-divergence
(3) $H_{f}(c)$, or, more precisely, it should be a monotonic function of
$H_{f}(c)+\lambda_{i}c_{i}$ for some constant $\lambda$. Nevertheless, now we
know nothing about these $f$-divergences except their existence. Constructive
transformation of $H_{\Gamma}^{*}$ into $f$-divergence is desirable because an
explicit form (3) brings some benefits for analysis.

## Acknowledgments

The work was supported by the University of Leicester and the Ministry of
Science and Higher Education of the Russian Federation (Project No.
14.Y26.31.0022).

## References

  * [1] Boltzmann L. Weitere Studien über das Wärmegleichgewicht unter Gasmolekülen. Sitzungsber Keis Akad Wiss 1872;66:275–370.
  * [2] Boltzmann L. Neuer Beweis zweier Sätze über das Wärmegleichgewicht unter mehratomigen Gasmolekülen. Sitzungsber Keis Akad Wiss 1887;95(2):153–164.
  * [3] Angeli D. A tutorial on chemical network dynamics. Eur. J. Control, 2009;15:398–406.
  * [4] Shannon CE. A mathematical theory of communication. Bell Syst Tech J 1948;27:379–423,623–656.
  * [5] Cohen JE, Derriennic Y, Zbaganu GH. Majorization, monotonicity of relative entropy and stochastic matrices. Contemp Math 1993;149:251–59.
  * [6] Cohen JE, Iwasa Y, Rautu G, Ruskai MB, Seneta E, Zbaganu G. Relative entropy under mappings by stochastic matrices. Linear Algebra Appl 1993;179:211–35.
  * [7] Rényi A. On measures of entropy and information. In Proceedings of the 4th Berkeley Symposium on Mathematics, Statistics and Probability 1960, Volume 1; Berkeley, CA, USA: University of California Press; 1961; pp. 547–61.
  * [8] Csiszár I. Eine informationstheoretische Ungleichung und ihre Anwendung auf den Beweis der Ergodizität von Markoffschen Ketten. Magyar Tud Akad Mat Kutato Int Kozl 1963;8:85–108.
  * [9] Morimoto T. Markov processes and the $H$-theorem. J Phys Soc Jap 1963;12:328–331.
  * [10] Gorban PA. Monotonically equivalent entropies and solution of additivity equation. Physica A 2003;328:380–390.
  * [11] Amari, S.-I. Divergence, Optimization, Geometry. In Proceedings of the 16th International Conference on Neural Information Processing, Bangkok, Thailand, 1–5 December 2009; Leung, C.S., Lee, M., Chan, J.H., Eds.; Berlin, Germany: Springer; 2009; pp. 185–93.
  * [12] Gorban AN, Gorban PA, Judge G. Entropy: the Markov ordering approach. Entropy 2010;12(5):1145–1193.
  * [13] Marin GB, Yablonsky GS, Constales D. Kinetics of chemical reactions: Decoding complexity. Weinheim, Germany: Wiley-VCH; 2019.
  * [14] Yablonskii GS, Bykov VI, Gorban AN, Elokhin VI, Kinetic models of catalytic reactions. Amsterdam: Elsevier; 1991.
  * [15] Gorban AN. Equilibrium Encircling. Equations of Chemical Kinetics and Their Thermodynamic Analysis. Novosibirsk: Nauka; 1984.
  * [16] Chernikova NV. An algorithm for finding a general formula for nonnegative solutions of system of linear inequalities, USSR Comput. Math. Math. Phys., 5 (1965), pp. 228–233.
  * [17] Fukuda K, Prodon A. Double description method revisited, in Combinatorics and Computer Science, Lecture Notes in Comput. Sci. 1120, Springer-Verlag, Berlin, 1996, 91–111.
  * [18] Gorban AN, Kaganovich BM, Filippov SP, Keiko AV, Shamansky VA, Shirkalin IA. Thermodynamic Equilibria and Extrema: Analysis of Attainability Regions and Partial Equilibria. New York: Springer; 2006.
  * [19] Gorban AN. Thermodynamic Tree: The Space of Admissible Paths. SIAM J Appl Dyn Syst 2013;12(1):246–278.
  * [20] Pisarenko YA, Yakhyaev MA, Serna-Loaiza S, Cardona CA. Stoichiometric restrictions on operating modes in chemical technology. Chem. Eng. Sci. 2018;192:642–654.
  * [21] Hangos K, Szederkényi G. The effect of conservation on the dynamics of chemical reaction networks, In IFAC Workshop on Thermodynamic Foundations of Mathematical Systems Theory, July 13-16, Lyon, France. , pp. 30–35.
  * [22] Tolman RC. The Principles of Statistical Mechanics. London: Oxford University Press; 1938.
  * [23] Zeldovich YB. Proof of the uniqueness of the solution of the equations of the law of mass action. In: Selected Works of Yakov Borisovich Zeldovich; Volume 1, Ed. by JP Ostriker; Princeton, NJ, USA: Princeton Univ. Press; 1996; pp. 144–148.
  * [24] Shapiro NZ, Shapley LS. Mass action law and the Gibbs free energy function, SIAM J Appl Math 1965;16:353–375.
  * [25] Volpert AI, Khudjaev SI. Analysis in classes of discontinuous functions and equations of mathematical physics. Dordrecht, The Netherlands: Nijhoff; 1985. (Translation from the 1st Russian ed. Moscow: Nauka publ.; 1975.)
  * [26] Lorentz H-A. Über das Gleichgewicht der lebendigen Kraft unter Gasmolekülen. Sitzungsber Keis Akad Wiss 1887;95(2):115–152.
  * [27] Cercignani C, Lampis M. On the $H$-theorem for polyatomic gases. J Stat Phys 1981;26(4):795–801.
  * [28] Horn FJM. Necessary and sufficient conditions for complex balancing in chemical kinetics. Arch Ration Mech Anal 1972;49:172–186.
  * [29] Feinberg M. Complex balancing in general kinetic systems. Arch Rat Mechan Anal 1972;49:187–194.
  * [30] Szederkényi G, Hangos KM. Finding complex balanced and detailed balanced realizations of chemical reaction networks. J Math Chem 2011;49(6):1163–79.
  * [31] Gorban AN, Shahzad M. The Michaelis-Menten-Stueckelberg theorem. Entropy 2011;13:966–1019.
  * [32] Lipták G, Hangos KM, Pituk M, Szederkényi G. Semistability of complex balanced kinetic systems with arbitrary time delays. Syst Control Lett 2018;114:38–43.
  * [33] Sontag ED. Stability and stabilization: discontinuities and the effect of disturbances. In Nonlinear analysis, differential equations and control 1999 (pp. 551–598). Springer, Dordrecht.
  * [34] Clarke F. Nonsmooth analysis in control theory: a survey. European Journal of Control. 2001;7(2-3):145–59.
  * [35] Gorban AN. General H-theorem and entropies that violate the second law. Entropy 2014;16(5):2408–32.
  * [36] Gorban AN. Local equivalence of reversible and general Markov kinetics, Physica A 2013;392:1111–1121.
  * [37] Gorban AN, Bykov VI, Yablonskii GS. Thermodynamic function analogue for reactions proceeding without interaction of various substances, Chem Eng Sci 1986. Vol. 41, No. 11. P. 2739–2745.
  * [38] Pavlov A, Pogromsky A, van de Wouw N, Nijmeijer H. Convergent dynamics, a tribute to Boris Pavlovich Demidovich. Syst Control Lett 2004;52(3-4):257–261.
  * [39] De Leenheer P, Angeli D, Sontag ED. Monotone chemical reaction networks. J Math Chem 2007;41(3):295–314.
  * [40] Banaji M, Mierczyński J. Global convergence in systems of differential equations arising from chemical reaction networks. J Differ Equ 2013;254(3):1359–1374.
  * [41] Sontag ED. Monotone and near-monotone biochemical networks. Syst. Synth. Biol. 2007;1(2):59–87.
  * [42] E.T. Jaynes, Information Theory and Statistical Mechanics, I; II Phys. Rev. 1957 106, 620–63; 108, 171–190.
  * [43] Grmela, M. (2013). Role of thermodynamics in multiscale physics. Computers and Mathematics with Applications, 65(10), 1457–1470.
  * [44] Grmela, M. (2016). Externally driven macroscopic systems: Dynamics versus thermodynamics. Journal of Statistical Physics, 166(2), 282–316.
  * [45] Jou D, Casas-Vázquez J, Lebon G. Extended irreversible thermodynamics. Springer, Berlin, Heidelberg. 1996.
  * [46] Gorban AN, Karlin IV. Invariant manifolds for physical and chemical kinetics (Vol. 660). Springer, Berlin, Heidelberg. 2005.
  * [47] Rockafellar RT. Convex Analysis. Princeton, NJ, USA: Princeton University Press; 2015.
  * [48] Feinberg M. On chemical kinetics of a certain class. Arch Rat Mechan Anal 1972;46:1–41.
  * [49] Grmela M. Fluctuations in extended mass-action-law dynamics. Physica D 2012;241:976–986.
  * [50] Gorban AN, Mirkes EM, Yablonsky GS. Thermodynamics in the limit of irreversible reactions, Physica A 2013;392:1318–1335.
  * [51] Hangos KM. Engineering model reduction and entropy-based Lyapunov functions in chemical reaction kinetics. Entropy 2010;12:772–797.
  * [52] Bykov VI, Gorban AN, Yablonskii GS, Description of nonisothermal reactions in terms of Marcelin-De Donder kinetics and its generalizations, React. Kinet. Catal. Lett. 1982;20(3-4):261–265.
  * [53] Grmela M, Öttinger HC. Dynamics and thermodynamics of complex fluids. I. Development of a general formalism. Phys Rev E 1997;56:6620.
  * [54] Gorban AN, Sargsyan HP, Wahab HA. Quasichemical models of multicomponent nonlinear diffusion. Math Model Nat Phenom 2011;6(5):184–262.
  * [55] Taneja IJ, Kumar P. Relative information of type s, Csiszár’s f-divergence, and information inequalities. Inf. Sci. 2004;166(1-4):105–25.
  * [56] Rahmani H, Sahli N, Kamoun F. DDoS flooding attack detection scheme based on F-divergence. Comput. Commun. 2012;35(11):1380–91.
  * [57] Nowozin S, Cseke B, Tomioka R. f-gan: Training generative neural samplers using variational divergence minimization. In Proceedings of 30th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain, Advances in Neural Information Processing Systems 2016; pp. 271–279.
  * [58] Feinberg M, Hildebrandt D. Optimal reactor design from a geometric viewpoint–I. Universal properties of the attainable region. Chem. Eng. Sci. 1997;52:1637–1665.
  * [59] Hildebrandt D, Glasser D. The attainable region and optimal reactor structures. Chem Eng Sci 1990;45:2161–2168.
  * [60] Metzger MJ, Glasser D, Hausberger B, Hildebrandt D, Glasser BJ. Use of the attainable region analysis to optimize particle breakage in a ball mill. Chem. Eng. Sci. 2009;64:3766–3777