Reaction Kinetics Path Based on Entropy Production Rate and Its Relevance to Low-Dimensional Manifolds

Abstract

The equation that approximately traces the trajectory in the concentration phase space of chemical kinetics is derived based on the rate of entropy production. The equation coincides with the true chemical kinetics equation to first order in a variable that characterizes the degree of quasi-equilibrium for each reaction, and the equation approximates the trajectory along at least final part of one-dimensional (1-D) manifold of true chemical kinetics that reaches equilibrium in concentration phase space. Besides the 1-D manifold, each higher dimensional manifold of the trajectories given by the equation is an approximation to that of true chemical kinetics when the contour of the entropy production rate in the concentration phase space is not highly distorted, because the Jacobian and its eigenvectors for the equation are exactly the same as those of true chemical kinetics at equilibrium; however, the path or trajectory itself is not necessarily an approximation to that of true chemical kinetics in manifolds higher than 1-D. The equation is for the path of steepest descent that sufficiently accounts for the constraints inherent in chemical kinetics such as element conservation, whereas the simple steepest-descent-path formulation whose Jacobian is the Hessian of the entropy production rate cannot even approximately reproduce any part of the 1-D manifold of true chemical kinetics except for the special case where the eigenvector of the Hessian is nearly identical to that of the Jacobian of chemical kinetics.

1. Introduction

One of the classical problems in describing the stability of dynamic systems is to find a Lyapunov function for a process in which a chemically reacting system converges to chemical equilibrium. For example, after a pioneering study by Wei [1], who characterized chemically reacting systems by some axioms, Shear [2] and Higgins [3] found and discussed a Lyapunov function for detailed-balanced systems. Krambeck [4] also characterized chemically reacting systems using some axioms with thermodynamic relations instead of a Lyapunov function. Feinberg, Horn, Jackson, and others [5–7] developed “chemical reaction network theory” to give attention to mass-action kinetics; those studies describe the stable dynamic behavior of such systems independently of the values of rate parameters [8–10]. Although the results from [5–7] are summarized in [10] from a modern viewpoint to review recent developments in the theory of chemical reaction networks ([9,11–13] for examples), an important issue in these classical papers is that complex balanced systems (including detailed-balanced systems) have a unique positive equilibrium concentration that is locally asymptotically stable, at least for closed systems with ideal mixing-rules [14]; this mathematically supports that the second law of thermodynamics holds for such systems. Horn and Jackson [7] have shown this using a Lyapunov function similar to that of Shear [2] (though there is another history about this problem, which is briefly reviewed by Powers and Paolucci [15]), whereas Johnston [10] recently showed it by linearizing the mass-action kinetics equation near the equilibrium state. Typically, this illustrates that there are two ways to confirm the existence of an equilibrium state: one is the use of a Lyapunov function and the other is the rather direct use of kinetics equations. However, use of a Lyapunov function provides no trajectories to arrive at equilibrium in the concentration phase space of a chemical reaction system.

In non-equilibrium thermodynamics, Ziegler [16] has postulated that the rate of entropy production is maximized with certain prescribed dissipative forces. The principal results obtained from this hypothesis were known for linear non-equilibrium systems and the second law [17]. However, for chemical reaction systems, there is no example of deriving the equation that describes a chemical change directly based on entropy production or on a Lyapunov function that does not use a chemical kinetics equation. Gorban et al. [18–20] and Lebiedz [21–23] utilized the second law of thermodynamics to determine the so-called slow invariant manifold (SIM) or the so-called low-dimensional manifold (LDM). Gorban et al. [18–20] utilized the basis orthonormal with respect to their “entropic” scalar product that uses the Hessian of a Lyapunov function (the free energy of a perfect gas in a constant volume at constant temperature) because their almost orthogonal “projector” that defines SIM helps convergence to determine SIM [20]. Lebiedz [21–23] proposed a variational principle to identify or approximate SIM as a geodesic to minimize the distance from equilibrium in concentration phase space [22,23], and they [21–23] proposed and compared various types of Riemannian metrics for concentration phase space. Two of these metrics are related to thermodynamics: one [22] is the Hessian of a Lyapunov function, which is the same as used by Gorban et al. [18–20], and the other [23] is the metric of a diagonal matrix whose diagonal components correspond to each species and is defined as the sum of the products of the corresponding stoichiometric coefficients and entropy production rates of respective reactions. The success of these methods in approximating SIM suggests that there might be some thermodynamic equation that determines trajectories for arriving at equilibrium in concentration phase space, although it is uncertain how the thermodynamic methods used by Gorban et al. [18–20] and Lebiedz [21–23] are related to a chemical kinetics equation. Furthermore, Al-Khateeb et al. [24] have concluded that a system’s dynamics cannot be deduced from the Gibbs free energy or entropy production rate (thermodynamic Lyapunov functions) even near the equilibrium state, and they stated that equilibrium thermodynamic potentials do not alone determine the dynamics of a reactive system during the approach to physical equilibrium.

Based on the rate of entropy production in chemical kinetics, this study derives the equation that approximately traces the trajectory in the concentration phase space of chemical kinetics, and argues that the above statement of Al-Khateeb et al. [24] seems to need careful understanding and interpretation. In the derivation, a type of quasi-equilibrium approximation is used in which the difference between the forward and backward rates of each reversible reaction is assumed to be sufficiently small with respect to the backward rate. The equation coincides with the true chemical kinetics equation to first order in a defined variable that characterizes the degree of quasi-equilibrium of each reaction. We call this new equation for the linear non-equilibrium path (ELNEP) because the equation is justified as an approximation to true chemical kinetics for some linear non-equilibrium domain and identifies the trajectory that arrives at equilibrium in concentration phase space for such a domain, although such an approximation is effective only for detailed-balanced quasi-equilibrium and detailed-balanced completely equilibrium state, i.e., complex balancing is not considered below in this paper. (The mathematical difference between detailed balancing and complex balancing has been discussed in [25] based on the theory of chemical reaction networks).

For trajectories to arrive at equilibrium, there is a hierarchy of manifolds of lower dimension that tends to equilibrium [18–23,26–28], i.e., a bundling of trajectories near “manifolds of slow motion” of successively lower dimension (the so-called SIM or LDM already mentioned) as time progresses [21]. The ELNEP reproduces such a bundling into a one-dimensional (1-D) manifold and even higher dimensional manifolds that are almost the same as those of true chemical kinetics when the contour of the entropy production rate in the concentration phase space is not highly distorted.

This paper is organized as follows: Section 2 presents the mathematical foundation for closed spatially homogeneous reaction systems and the ELNEP (including its derivation), together with equations for two artificial paths (the steepest descent and Newton direction paths) in addition to the true chemical kinetics path. Section 3 presents and compares all four types of paths for simple examples of reacting systems, including a non-isothermal case and a practical example used by Al-Khateeb et al. [24]. Section 4 gives a theoretical proof (although not in a mathematically formal fashion) that ELNEP reproduces almost the same LDMs of true chemical kinetics. Section 5 briefly comments on the relevance of ELNEP to Ziegler’s hypothesis [16], and Section 6 summarizes this study.

2. Methodology (Mathematical Basis)

The pure reaction contribution to the entropy production rate of a thermally perfect gas in a closed system is as follows [29]:

$$\frac{dS}{dt}=\widehat{R}V\sigma$$

(1)

where t is time, S is the entropy of the system, $\widehat{R}$ is the universal gas constant, V is the volume of the system, and σ is the total reaction contribution to the rate of entropy production.

The latter is minus one times the time derivative of the Lyapunov function found by Shear [2] and used by Gorban et al. [18–20], Lebiedz [22], and many others. Therefore, σ can be interpreted as another Lyapunov function of the system in the linear non-equilibrium regime [30], although this interpretation is not necessarily required for the derivation of ELNEP:

$$\sigma ={\displaystyle \sum _{l=1}^{N_r}{\sigma }_l}$$

(2)

where:

$${\sigma }_l=\left(R_{fl}-R_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}$$

(3)

Equation (3) gives the contribution of the l th reaction that has a forward reaction rate R_fl and reverse reaction rate R_bl. Chemical reactions proceed in accordance with the following ordinary differential equation of mass-action kinetics:

$$\frac{d\overrightarrow{y}}{dt}={\displaystyle \sum _{l=1}^{N_r}\left({\overrightarrow{\beta }}_l-{\overrightarrow{\alpha }}_l\right)\left(R_{fl}-R_{bl}\right)}$$

(4)

where:

$$R_{fl}=k_{fl}{\displaystyle \prod _{k=1}^N{y}_k^{{\alpha }_{kl}}}$$

(5a)

$$R_{bl}=k_{bl}{\displaystyle \prod _{k=1}^N{y}_k^{{\beta }_{kl}}}$$

(5b)

$$\overrightarrow{y}={\left[\begin{array}{llll}{y}_1\hfill & y_2\hfill & \cdots \hfill & y_N\hfill \end{array}\right]}^T$$

(6)

$${\overrightarrow{\alpha }}_l={\left[\begin{array}{llll}{\alpha }_{1l}\hfill & {\alpha }_{2l}\hfill & \cdots \hfill & {\alpha }_{Nl}\hfill \end{array}\right]}^T$$

(7a)

$${\overrightarrow{\beta }}_l={\left[\begin{array}{llll}{\beta }_{1l}\hfill & {\beta }_{2l}\hfill & \cdots \hfill & {\beta }_{Nl}\hfill \end{array}\right]}^T$$

(7b)

Here k_fl and k_bl are the reaction rate constants of the lth forward and reverse reactions, respectively; y_k is the concentration of the kth species; and α_kl and β_kl are the stoichiometric coefficients of the kth species in the lth forward and reverse reactions, respectively. In Section 2.1, the ELNEP is derived; however, before the derivation, two types of artificial paths that are apparently based on the rate of entropy production in Equation (2) are discussed briefly. These paths can be used to compare with the ELNEP.

Steepest descent path. Ziegler [16] proposed an orthogonality principle for irreversible processes, which is equivalent to a certain extremal principle (i.e., the maximal rate of entropy production). However, the thermodynamic and statistical bases of the principle are not clear, although the second law of thermodynamics can be obtained as a corollary of the principle [17]. In the present study, as a candidate for the uncertain thermodynamic equation that determines trajectories to arrive at equilibrium in concentration phase space, the principle is attempted to apply by assuming that reactions proceed in the direction of steepest descent of the entropy production rate σ given in Equation (1), i.e., in the direction of the maximal rate of entropy production:

$$\frac{d{\overrightarrow{y}}_{ind}}{ds}=-\nabla \sigma$$

(8)

where S is the pseudo-time (i.e., an artificial time-marching parameter whose unit is not necessarily that of physical time) and ${\overrightarrow{y}}_{ind}$ is the vector that consists of the species concentrations chosen as the set of independent variables consistent with element conservation (see Section 2.1). Note that the Jacobian of Equation (8) is the Hessian of σ, and that the Hessian of the Lyapunov function is used by Gorban et al. [18–20] and Lebiedz et al. [22] in their methods for identifying SIM, as described in Section 1. In addition, the left-hand side of Equation (8) would have a different unit from that of the right-hand side if the unit of S were that of physical time. Hereafter, the steepest descent path is termed to be the one that is in accordance with Equation (8).

Newton direction path. Equation (8) for the steepest descent path can be interpreted as a version of a well-known optimization method (i.e., the method of steepest descent or ascent). The standard algorithm uses a line search along a definite length in the search direction, whereas Equation (8) uses the pseudo-time and a differential. Another well-known optimization method is Newton iteration in which the function to be optimized is approximated by a second-order Taylor series expansion in the direction of the line search, the so-called Newton direction [31]. The Newton direction path can be calculated in accordance with:

$$H\frac{d{\overrightarrow{y}}_{ind}}{ds}=-\nabla \sigma$$

(9)

where:

$$H\equiv \nabla \otimes \left(\nabla \sigma \right)=\frac{\partial \nabla \sigma }{\partial {\overrightarrow{y}}_{ind}}$$

(10)

is the Hessian (matrix of second derivatives) of σ with respect to each component of ${\overrightarrow{y}}_{ind}$. Equation (9) might be another candidate for the uncertain thermodynamic equation, although it is rejected, as obviously shown below in this study. In spite of that, the left-hand side of Equation (9) has the same unit as that of right-hand side when the unit of S is that of physical time.

2.1. Derivation of ELNEP (Proof of First-order Coincidence)

First, the form of the ELNEP is introduced, and its derivation is described. The ELNEP is the following:

$$Q\frac{d{\overrightarrow{y}}_{ind}}{dt}=-\frac{1}{2}\nabla \sigma$$

(11)

where Q can be decomposed into the two parts (as the following Equation 12): One appears due to the constraints inherent in chemical kinetics such as element conservation and another is the part independent of the constraints. When the constraints do not exist, the first part disappears and consequently ELNEP becomes essentially equivalent to Equation (8) as explained later in this section:

(12)

Note that t is not the pseudo-time, and that the matrix Q is symmetric; hence, Q is always invertible because it is positive definite, as proved in Appendix E. The matrix B₃ is defined by:

$$B_3\equiv B_2^{-1}{B}_1$$

(13)

where B₂ [32] and B₁ are the coefficient matrices of the constraints:

$$B_1{\overrightarrow{y}}_{ind}+B_2{\overrightarrow{y}}_{dep}=\overrightarrow{c}$$

(14)

and (B₃)_i,j is the component in the ith row and jth column of matrix B₃. Here:

$${\overrightarrow{y}}_{ind}={\left[\begin{array}{llll}{y}_1\hfill & y_2\hfill & \cdots \hfill & y_K\hfill \end{array}\right]}^T$$

(15a)

$${\overrightarrow{y}}_{dep}={\left[\begin{array}{llll}{y}_{K+1}\hfill & y_{K+2}\hfill & \cdots \hfill & y_N\hfill \end{array}\right]}^T\text{and}$$

(15b)

$$\overrightarrow{c}={\left[\begin{array}{llll}{c}_1\hfill & c_2\hfill & \cdots \hfill & c_M\hfill \end{array}\right]}^T$$

(16)

where K is the number of independent species concentrations (the components of the vector ${\overrightarrow{y}}_{ind}$), M is the number of constraints (such as the number of elements) and also the number of dependent species concentrations (the components of the vector ${\overrightarrow{y}}_{dep}$), N = K + M; and c_i for i = 1, ⋯, M is the constants (such as the total number density of the i th element in the system). In principle, the choice of independent and dependent species is arbitrary if the constraints are satisfied. However, the choice may affect the speed or ease of numerical calculation, which is not investigated at present.

It can be proved that Equation (11) coincides with Equation (4) (true chemical kinetics equation) to first order in the variable ε_l (the difference is in the second order of ε_l), which is defined by:

$$R_{fl}=\left(1+{\epsilon }_l\right)R_{bl}$$

(17)

Here, it should be noted that sufficiently small |ε_l| means partial equilibrium or quasi-equilibrium of the lth reaction, and detailed-balanced equilibrium means ε_l = 0 for every lth reaction in the system. Accordingly, ELNEP is supposed to be applicable to the relaxation process leading to detailed- balanced equilibrium, i.e., a series of successive processes in which the number of elementary reactions that participate in quasi-equilibrium of the system increases (although more careful description of the process is given in Section 4 and Note [33]). In fact, the paths determined by ELNEP trace the low-dimensional manifolds of true chemical kinetics, which is shown in Section 3 and theoretically proved in Section 4; although ELNEP cannot be applied to the complex balancing that does not guarantee ε_l = 0 and even ε_l ≈ 0. In brief, the existence of detailed-balanced equilibrium guarantees that each ε_l approaches zero sooner or later, and the magnitude of each |ε_l| is a rough indicator how closely ELNEP approximates the path of true chemical kinetics, i.e., ELNEP does not work best in the domain far from equilibrium (large |ε_l|s). Now, Equation (11) is called as ELNEP, which is derived as follows and its first-order coincidence with Equation (4) is proved.

To begin with, it shall be considered how to properly use −∇σ as a type of thermodynamic “force” to generate paths that are close to the true chemical kinetics path. In such an interpretation, σ would be a thermodynamic “potential” determined by the “position” vector $\overrightarrow{y}$. Such forces are usually employed in optimization techniques, for example, the method of steepest descent or ascent and the Newton iteration, as mentioned earlier for introducing two artificial paths. With respect to this thinking, basic Equations (2) and (3) are calculated to obtain:

$$\frac{\partial \sigma }{\partial y_i}={\displaystyle \sum _{l=1}^{N_r}\frac{\partial {\sigma }_l}{\partial y_i}={\displaystyle \sum _{l=1}^{N_r}\left\{\left(\frac{\partial R_{fl}}{\partial y_i}-\frac{\partial R_{bl}}{\partial y_i}\right)\ln \frac{R_{fl}}{R_{bl}}+\left(R_{fl}-R_{bl}\right)\left(\frac{1}{R_{fl}}\frac{\partial R_{fl}}{\partial y_i}-\frac{1}{R_{bl}}\frac{\partial R_{bl}}{\partial y_i}\right)\right\}}}$$

(18)

Substituting Equation (5) into Equation (18) gives:

$$\frac{\partial \sigma }{\partial y_i}=\frac{1}{y_i}\left\{{\displaystyle \sum _{l=1}^{N_r}\left({\alpha }_{il}{R}_{fl}-{\beta }_{il}{R}_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}+{\displaystyle \sum _{l=1}^{N_r}\left(R_{fl}-R_{bl}\right)\left({\alpha }_{il}-{\beta }_{il}\right)}}\right\}.$$

(19)

At this stage, the constraint must be constituted as follows: As described in Appendix A, the mathematical identity:

$$\frac{\partial f}{\partial x_{ind}}={\left(\frac{\partial f}{\partial x_{ind}}\right)}_{\begin{array}{l}\text{without}\\ \text{constraints}\end{array}}+{\displaystyle \sum _{k=dep}{\left(\frac{\partial f}{\partial x_k}\right)}_{\begin{array}{l}\text{without}\\ \text{constraints}\end{array}}}\frac{\partial x_k}{\partial x_{ind}}$$

(A4)

holds for an arbitrary function f, whose arguments (x_ind and x_dep) are not independent of each other because of constraints between the x_ind and x_dep variables. Equation (A4) is applied to Equation (19) by treating σ and y in Equation (19) as f and x in Equation (A4), respectively. The result is:

(20)

where y_ind stands for each component of the vector ${\overrightarrow{y}}_{ind}$. Since Equation (4) holds, Equation (20) is evidently equivalent to:

$$\begin{array}{l}\frac{\partial \sigma }{\partial y_{ind}}=\frac{1}{y_{ind}}\left\{{\displaystyle \sum _{l=1}^{N_r}\left({\alpha }_{ind,l}{R}_{fl}-{\beta }_{ind,l}{R}_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}-\frac{d{y}_{ind}}{dt}}\right\}\\ +{\displaystyle \sum _{k=dep}\frac{1}{y_k}\left\{{\displaystyle \sum _{l=1}^{N_r}\left({\alpha }_{kl}{R}_{fl}-{\beta }_{kl}{R}_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}-\frac{d{y}_k}{dt}}\right\}\frac{\partial y_k}{\partial y_{ind}}}\end{array}$$

(21)

In this paper, Equation (14) provides a relation between the two concentration vectors ${\overrightarrow{y}}_{ind}$ and ${\overrightarrow{y}}_{dep}$ whose components are the y_ind and y_dep variables in Equation (21), respectively. Accordingly, $\frac{\partial y_{dep}}{\partial y_{ind}}$ corresponds to ${\left(-B_3\right)}_{dep,ind}$ defined by Equation (13), since the derivative of Equation (14):

$$B_{1}d{\overrightarrow{y}}_{ind}+B_{2}d{\overrightarrow{y}}_{dep}=\overrightarrow{0}$$

(22)

holds. Thus, Equation (21) becomes:

$$\begin{array}{l}\frac{\partial \sigma }{\partial y_{ind}}=\frac{1}{y_{ind}}\left\{{\displaystyle \sum _{l=1}^{N_r}\left({\alpha }_{ind,l}{R}_{fl}-{\beta }_{ind,l}{R}_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}-\frac{d{y}_{ind}}{dt}}\right\}\\ -{\displaystyle \sum _{k=dep}\frac{1}{y_k}\left\{{\displaystyle \sum _{l=1}^{N_r}\left({\alpha }_{kl}{R}_{fl}-{\beta }_{kl}{R}_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}-\frac{d{y}_k}{dt}}\right\}}{\left(B_3\right)}_{k,ind}\end{array}$$

(23)

Furthermore, $\frac{d{y}_k}{dt}$ in Equation (23) can be replaced as follows. From Equation (22), $d{y}_k=-{\displaystyle \sum _{j=ind}{\left(B_3\right)}_{k,j}d{y}_j}$ and:

$$\frac{d{y}_k}{dt}=-{\displaystyle \sum _{j=ind}{\left(B_3\right)}_{k,j}\frac{d{y}_j}{dt}}$$

(24)

Substituting Equation (24) into Equation (23), the following is obtained:

$$\begin{array}{l}\frac{\partial \sigma }{\partial y_{ind}}=\frac{1}{y_{ind}}\left\{{\displaystyle \sum _{l=1}^{N_r}\left({\alpha }_{ind,l}{R}_{fl}-{\beta }_{ind,l}{R}_{bl}\right)}\ln \frac{R_{fl}}{R_{bl}}-\frac{d{y}_{ind}}{dt}\right\}\\ -{\displaystyle \sum _{k=dep}\frac{1}{y_k}\left\{{\displaystyle \sum _{l=1}^{N_r}\left({\alpha }_{kl}{R}_{fl}-{\beta }_{kl}{R}_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}+{\displaystyle \sum _{j=ind}{\left(B_3\right)}_{k,j}\frac{d{y}_j}{dt}}}\right\}{\left(B_3\right)}_{k,ind}}\end{array}$$

(25)

Equation (25) is equivalent to Equation (4) (true chemical kinetics equation) because Equation (25) is just a variety of Equation (4) based on an identity Equation (20).

The ELNEP is obtained by substituting Equation (17) into Equation (25) as follows. Substituting Equation (17) into the term $\left({\alpha }_{•,l}{R}_{fl}-{\beta }_{•,l}{R}_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}$ in Equation (25) gives:

$$\left({\alpha }_{•,l}{R}_{fl}-{\beta }_{•,l}{R}_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}=\left({\alpha }_{•,l}-{\beta }_{•,l}\right){\epsilon }_l{R}_{bl}+O\left({\epsilon }_l^2\right)=\left({\alpha }_{•,l}-{\beta }_{•,l}\right)\left(R_{fl}-R_{bl}\right)+O\left({\epsilon }_l^2\right)$$

(26)

Therefore, Equation (26) is an approximation that is accurate to first order in the variable ε_l. Then:

$${\displaystyle \sum _{l=1}^{N_r}\left({\alpha }_{•,1}{R}_{fl}-{\beta }_{•,l}{R}_{bl}\right)\ln \frac{R_{fl}}{R_{bl}}=}{\displaystyle \sum _{l=1}^{N_r}\left({\alpha }_{•,l}-{\beta }_{•,l}\right)\left(R_{fl}-R_{bl}\right)+{\displaystyle \sum _{l=1}^{N_r}O\left({\epsilon }_l^2\right)}=-\frac{d{y}_{•}}{dt}+{\displaystyle \sum _{l=1}^{N_r}O\left({\epsilon }_l^2\right)}}$$

(27)

Using Equations (24) and (27), Equation (25) becomes:

(28)

which is the same as Equation (11) except for the second order term of each ε_l. Accordingly, the ELNEP Equation (11) is equivalent to Equation (4) (true chemical kinetics equation) to first order in each variable ε_l.

Note that ∇σ on the right-hand side (rhs) of Equations (8), (9), and (11) must also account for the constraints Equation (14) to eliminate y_dep in the calculation; whereas Equation (4) intrinsically satisfies the constraints. If any constraint is not considered, Equation (28) can be transformed into the following equation:

$$\frac{\partial \sigma }{\partial y_i}=-\frac{2}{y_i}\frac{d{y}_i}{dt}+{\displaystyle \sum _{l=1}^{N_r}O\left({\epsilon }_l^2\right)}$$

(29)

which is essentially equivalent to the following Equation (8) for the steepest descent path to first order in each variables ε_l (as described in Appendix B):

$$\frac{d{\overrightarrow{y}}_{ind}}{ds}=-\nabla \sigma$$

(8)

Accordingly, the ELNEP is the equation of the steepest descent path that sufficiently accounts for element conservation. In addition, note that Equation (29) and ELNEP are described with the variable t whose unit is that of physical time, whereas Equation (8) is described with the variable s, pseudo-time whose unit is not that of physical time.

2.2. Relevance to Linear Non-Equilibrium

Since the following Equation (30) also holds to second order in ε_l, the path that is in accordance with Equation (11) can be termed the path of linear non-equilibrium:

$$\sigma ={\displaystyle \sum _{l=1}^{N_r}\left({\epsilon }_l^2{R}_{bl}+O\left({\epsilon }_l^3\right)\right)}={\displaystyle \sum _{l=1}^{N_r}{X}_l{J}_l}+{\displaystyle \sum _{l=1}^{N_r}O\left({\epsilon }_l^3\right)},\text{with}$$

(30)

$$J_l\equiv {\epsilon }_l{R}_{bl}=R_{bl}{X}_l,\text{and}$$

(31)

$$X_l\equiv {\epsilon }_l$$

(32)

where J_l and X_l are, respectively, interpreted as the thermodynamic flux and thermodynamic force due to the lth reaction.

3. Results: Sample Calculations

In addition to the path of linear non-equilibrium (the path according to ELNEP, Equation 11), the chemical kinetics path (Equation 4) and the paths given by Equations (8) and (9) were calculated for the following examples of simple reaction systems, all of which are detailed-balanced systems with mass-action kinetics because ELNEP is intrinsically based on the reversibility of each reaction and the mass-action law. All calculations in this paper are for a constant volume with fixed or variable (non-isothermal) reaction rate constants, and are dimensionless except for the last example (Section 3.4). In the following, “one-dimensional” is abbreviated as “1-D,” “two-dimensional” is abbreviated as “2-D,” etc.

3.1. A ⇆ 2B ⇆ C: 1-D Manifold

Figures 1 and 2 show the two groups of paths, respectively, for the reactions A⇆2B⇆C projected onto the concentration (y_A, y_B) phase plane with a contour plot of σ. Figure 1 shows the chemical kinetics and linear non-equilibrium paths, and Figure 2 shows the steepest descent and Newton direction paths. Variable y_C is determined by element conservation:

$$y_A+0.5y_B+y_C=c$$

(33)

These calculations were performed with the fixed rate constants given in

Table 1. Dimensionless values [37] of the rate constants employed for the calculations shown in Figures 1 and 2.

Reaction	A→2B		2B→C
Figures	Forward	Reverse	Forward	Reverse
1 and 2	1	1	0.5	0.5

for the same equilibrium composition but with different initial compositions (twelve in each figure [34]), and c = 3. The dimensionless equilibrium composition is (y_A, y_B, y_C) = (1.223, 1.106, 1.223). Note that all paths reside in the so-called reaction simplex [35] (kinetics subspace or stoichiometric subspace, both of which coincide [10,36] for the present system) that is determined by c and the stoichiometric coefficients (or the so-called reaction vector [36]) of each reaction. The present reaction simplex is a 2-D plane in the 3-D concentration (y_A, y_B, y_C) phase space. The Following Observations are made from Figures 1 and 2:

(1)

The paths of linear non-equilibrium (magenta lines in Figure 1) most closely resemble that of chemical kinetics. They have the same 1-D manifold as that of chemical kinetics (one of the eigenvectors of the Jacobian of Equation (4) at equilibrium) expected from the path coinciding with the chemical kinetics path to first order in each variable ε_l. (The two eigenvectors of the Jacobian of Equation (4) at equilibrium are (1, 0.1544) and (1, −2.928) in the coordinates (y_A, y_B). The former is the 1-D manifold.)

(2)

The steepest descent paths (blue lines in Figure 2) have a slightly different 1-D manifold from that of chemical kinetics because the Jacobian of Equation (8) is the Hessian H in Equation (10); this Jacobian differs from that of Equation (4). (The two eigenvectors of the Hessian at equilibrium are (1, 0.1193) and (1, −8.380) in the coordinates (y_A, y_B).

(3)

The Newton direction paths (orange lines in Figure 2) have no 1-D manifold or three quasi-1-D manifolds.

It should be noted again that the Jacobian of Equation (8) is the Hessian H, and the difference between ELNEP (Equation 11) and Equation (8) is that Equation (11) does and Equation (8) does not sufficiently account for the constraints such as element conservation, as described in Section 2.1.

Appendix C shows the specific form of the ELNEP (Equation 11) for this case.

3.2. A ⇆ B ⇆ C ⇆ D: 2-D and 1-D Manifolds

Figures 3 and 4 show the two groups of paths, respectively, for the reactions A ⇆ B ⇆ C ⇆ D projected into the concentration (y_A, y_B, y_C) phase space. Figure 3 shows the chemical kinetics and linear non-equilibrium paths, and Figure 4 shows the steepest descent and Newton direction paths. Variable y_D is determined by element conservation:

$$y_A+y_B+y_C+y_D=c$$

(34)

These calculations were performed with the fixed rate constants given in

Table 2. Dimensionless values [37] of the rate constants employed for the calculations shown in Figures 3–7.

Reaction	A→B		B→C		C→D
Figures	Forward	Reverse	Forward	Reverse	Forward	Reverse
3–7	1	1	0.3	0.3	0.09	0.09

for the same equilibrium composition but with different initial compositions (seven in each figure [34]), and c = 3. Each of the seven initial compositions gives all four types of paths plotted in Figures 3 and 7. The dimensionless equilibrium composition is (y_A, y_B, y_C, y_D) = (0.75, 0.75, 0.75, 0.75). Note that all paths reside in the reaction simplex determined by c and the stoichiometric coefficients of each reaction, as described in Section 3.1.

The present reaction simplex is a 3-D space in the 4-D concentration (y_A, y_B, y_C, y_D) phase space. Figures 5 and 7 for examining the existence of 2-D manifolds of each path are also shown. The three eigenvectors of the chemical kinetics Jacobian at equilibrium are given in

Table 3. Three eigenvectors of the chemical kinetics Jacobian at equilibrium shown in Figures 3–7.

In the coordinates (yA, yB, yC)	Remarks
(1, 0.8991, 0.2606)	1-D manifold (α-eigenvector) of true chemical kinetics
(1, 0.5111, −1.951)	β-eigenvector of true chemical kinetics
(1, −1.190, 0.1988)	γ-eigenvector of true chemical kinetics

, and those of the Hessian of σ at equilibrium are given in

Table 4. Three eigenvectors of the Hessian of σ at equilibrium shown in Figures 3–7.

In the coordinates (yA, yB, yC)	Remarks
(1, 0.8144, −0.2330)	1-D manifold of the steepest decent path
(1, 4.3424, 19.467)	Almost on the plane that consists of α- and β-eigenvectors of true chemical kinetics
(1, −1.168, 0.2092)	Almost the same as γ-eigenvector of true chemical kinetics

.

The Following Observations are made from Figures 3 and 7:

(1)

The paths of linear non-equilibrium (Figure 3) most closely resemble that of chemical kinetics.

(2)

The path of linear non-equilibrium (Figures 3 and 5) has the same 2-D manifold (the plane formed by the two eigenvectors attached to the equilibrium point, denoted as α and β in Figure 5) and the same 1-D manifold (the eigenvector of the chemical kinetics Jacobian at equilibrium, denoted as α in Figure 5) as those of chemical kinetics. However, the linear non-equilibrium path is considerably different from the chemical kinetics path in the noted 2-D manifold (Figure 5).

(3)

The steepest descent path (Figures 4 and 6) has a 2-D manifold, i.e., the plane formed by the two eigenvectors of the Hessian of σ (blue straight lines); the two eigenvectors are adjacent to the eigenvectors of the chemical kinetics Jacobian, which are denoted as α and β (green straight lines). The plane (in Figure 6) is almost the same as those of chemical kinetics and linear non-equilibrium paths (in Figure 5), but the 1-D manifold is notably different between the steepest descent path (Figure 6) and chemical kinetics (likewise, the path of linear non-equilibrium) (Figure 5). The latter is the eigenvector of the chemical kinetics Jacobian denoted as α (α-eigenvector), whereas the former is the eigenvector of the Hessian of σ that is adjacent to the α-eigenvector.

(4)

The Newton direction path (Figures 4 and 7) is unique and distinct compared to the other paths. A quasi-1-D manifold similar to those of the chemical kinetics path (likewise the path of linear non-equilibrium) and the steepest descent path was observed; however, it could not be described definitely. A quasi-2-D manifold parallel to the γ-eigenvector (Figure 7) was also observed, which is peculiar because this quasi-2-D manifold does not encompass a 1-D manifold.

From these observations and those described in Section 3.1, it is obvious that Equation (9), which governs the Newton direction path, is not the thermodynamic equation (or even its approximation) that is sought, and Equation (8), which governs the steepest descent path, cannot correctly reproduce the 1-D manifold of true chemical kinetics, whereas the ELNEP does. More importantly, point 2 (above) denotes that ELNEP pursues a course in the same LDMs as those of true chemical kinetics but appreciably deviates from the true kinetics path within the noted LDM. The generality of the observation that the ELNEP almost reproduces the higher dimensional manifolds of true chemical kinetics is proved in Section 4. Appendix C shows the specific form of the ELNEP for this case.

3.3. Non-isothermal Case

Observations similar to those previously described for the isothermal case can be expected even for the non-isothermal case because Equation (11) (ELNEP) is an approximation to Equation (4) (to first order in each variable ε_l); this holds irrespective of variations in temperature. Figure 8 shows the chemical kinetics and linear non-equilibrium paths for the non-isothermal reactions A ⇆ B ⇆ C projected into the phase space (y_A, y_B, T / 200), where T is the dimensionless temperature of the system. We determined y_C by element conservation:

$$y_A+y_B+y_C=c$$

(35)

The contour of σ can be plotted on intersections orthogonal to the axis of T in Figure 8, but these intersectional contours are different from each other because the corresponding temperature and reaction rate constants vary. These calculations were performed with the variable rate constants (that must obey the Arrhenius law and hence depends on the gas temperature T) given in

Table 5. Dimensionless values [37] of rate parameters and heats of reaction employed for the calculations shown in Figures 8 and 9. The rate constant is expressed as k = A exp (−E/T).

Reaction	A→B					B→C
Figures	Forward		Reverse		Heat of reaction	Forward		Reverse		Heat of reaction
8–10	A	E	A	E		A	E	A	E
	1	0	1	700	0	0.3	700	0.3	10−7	−1.5

and c = 3. The dimensionless equilibrium composition and temperature are (y_A, y_B, y_C, T) = (0.7014, 1.597, 0.7014, 850.7).

Figures 9a,b show the four paths (chemical kinetics, linear non-equilibrium, steepest descent, and Newton direction) projected onto the concentration (y_A, y_B) phase plane for the same calculation as for Figure 8a, together with the contour of σ at the final equilibrium temperature. Figures 9c,d also show the paths for the isothermal case, the temperature of which was set as the final equilibrium temperature of the non-isothermal case; therefore, the dimensionless equilibrium composition and temperature are the same as those of the non-isothermal case.

Since the system is considered to be isolated in a constant volume, internal energy is conserved. This constraint can be treated as an equation for just determining the temperature T from the species concentrations (as described in Appendix D), and T disappears in the formulation. Because ∇σ in the rhs of Equation (11) has no derivative with respect to T, the specific form of the ELNEP Equation (11) is the same as that for the isothermal case (Appendix C shows the form for this example).

All species (A, B, and C) are treated as calorically and thermally perfect with one (dimensionless) specific heat of 10⁻³ at constant volume; Table 5 gives the (dimensionless) heats of reactions. The reaction A→B is neither endothermic nor exothermic, and the reaction B→C is exothermic.

As of yet, the Shear Lyapunov function [2,3] is always applied to isothermal cases [18–20,22,35]; therefore, it is expected that, likewise, the ELNEP cannot be rationalized for non-isothermal cases. However, σ can be deduced from thermodynamic considerations that allow arbitrary rate constants [29], and the mathematical derivation of the ELNEP (in Section 2.1) does not need an interpretation of σ or even its thermodynamic origin. Furthermore, the proof [10] of the Shear Lyapunov function allows variable (i.e., non-isothermal) rate constants, even though the mathematical form of the Shear Lyapunov function coincides, for some reason, with a variant of the Helmholtz free energy of an ideal gas whose usage is physically limited to the isothermal and isochoric condition.

The following observations are made from Figures 8 and 9:

• Since the temperature T is determined from the composition through energy conservation, the number of degrees of freedom for this reaction system is two, the same as that for the corresponding isothermal system. Accordingly, the observed LDM is 1-D (Figure 8). Note that all paths reside in the reaction simplex determined by c and the stoichiometric coefficients of each reaction, as described in Sections 3.1 and 3.2. The present reaction simplex is a 2-D space in the 3-D concentration (y_A, y_B, y_C) phase space, although Figure 8 shows the projection into the phase space (y_A, y_B, T/200) that is a pile (along the axis of T/200) of instantaneous projections at respective instants and temperatures.

• Similar results are observed as those found for the examples described in Sections 3.1 and 3.2, but the following results are also noted.

• The paths of linear non-equilibrium with their initial compositions closer to equilibrium (Figure 8a) approximate more closely the paths of true chemical kinetics than initial compositions far from equilibrium (Figure 8b). This can be really observed for the previous examples given in Sections 3.1 and 3.2.

• The 1-D manifold for the non-isothermal case substantially differs from that of the isothermal case because the temperature dependence of the rate constants distorts the contour of σ for the non-isothermal case (Figure 8 and Figures 9a,b).

• In the isothermal case, the steepest descent paths (blue curved lines in Figure 9d) have no LDM, whereas those in the non-isothermal case (blue curved lines in Figure 9b) have a LDM. This indicates that the rate constants employed in this case (Table 5) render the isothermal (non-isothermal) cases insufficiently (sufficiently) to produce the gap in reaction rate between the two reversible reactions, respectively. Nevertheless, the paths of linear non-equilibrium (magenta lines in Figures 9a,c) have a 1-D manifold in each case and approximate the true chemical kinetics irrespective of the difference in the gap.

3.4. An Example of Reference [24]

This example is adopted from Section IV. A of Al-Khateeb et al. [24], where it serves as a simple but realistic reactive system for constructing a 1-D slow invariant manifold (SIM). Also, Al-Khateeb et al. [24] have demonstrated their conclusion (mentioned in Section 1 of this paper) mainly using this example.

The kinetic model of this example is the well-known Zel’dovich mechanism that consists of five species, two elements, and two reversible reactions: N + O2 ⇆ NO + O and N + NO ⇆ N2 + O, which has a 1-D manifold in two-dimensional reaction simplex owing to the following three constraints in its chemistry:

$$y_{\mathrm{NO}}+y_{\mathrm{O}}+2y_{{\mathrm{O}}_2}=c_{\mathrm{O}}$$

(36a)

$$y_{\mathrm{NO}}+y_{\mathrm{N}}+2y_{{\mathrm{N}}_2}=c_{\mathrm{N}}$$

(36b)

$$y_{\mathrm{NO}}+y_{\mathrm{N}}+y_{\mathrm{O}}+y_{{\mathrm{O}}_2}+y_{{\mathrm{N}}_2}=c_{\text{total}}$$

(36c)

where y_NO, y_N, y_O, $y_{{\mathrm{O}}_2}$, and $y_{{\mathrm{N}}_2}$ are concentrations of the species NO, N, O, O₂, and N₂, respectively. These constraints are the conservations of element oxygen, element nitrogen, and the total number of species in the system, respectively. The calculation is done with:

$${\overrightarrow{y}}_{ind}={\left[\begin{array}{ll}{y}_{NO}\hfill & y_N\hfill \end{array}\right]}^T\text{and}$$

(37a)

$${\overrightarrow{y}}_{dep}={\left[\begin{array}{lll}{y}_{\mathrm{O}}\hfill & y_{{\mathrm{O}}_2}\hfill & y_{{\mathrm{N}}_2}\hfill \end{array}\right]}^T$$

(37b)

and Appendix C shows the specific form of the ELNEP for this case.

All the data employed here (including kinetic and thermodynamic data) have been kindly provided by the author of Reference [24], and the calculation in this section is dimensional as in Reference [24]. Figure 10 shows the chemical kinetics, linear non-equilibrium, and the steepest descent paths for twelve initial compositions at constant temperature of 4,000 K in constant volume (the same condition and the same final equilibrium as in Reference [24] except for the initial compositions). The Newton direction path is not shown because of its unimportance indicated by the previous examples.

From this figure the following observations are made:

(1)

1-D manifold of true chemical kinetics is obviously different from those of linear non-equilibrium paths and the steepest descent paths.

(2)

1-D manifold of linear non-equilibrium paths alomost conincides with that of the steepest descent paths.

(3)

At equilibrium the eigenvectors of the Jacobian of true chemical kinetics and those of the Hessian of σ are almost the same in agreement with the description in Reference [24].

(4)

In agreement with the theory given in Section 2, 1-D manifold of linear non-equilibrium paths is tangential to the eigenvectors of the Jacobian of true chemical kinetics. Also, 1-D manifold of the steepest descent paths is tangential to the eigenvectors of the Hessian of σ, which is consistent with Points 2 and 3.

Points 1 and 2 are very striking features that are quite different from those of the previous examples (Sections 3.1–3.3), where 1-D manifold of linear non-equilibrium paths and that of the steepest descent paths are clearly different and 1-D manifold of linear non-equilibrium paths are approximately the same as that of true chemical kinetics. However, Points 2 and 3 seem to have no direct relation to Point 1 [38]. The most important is Point 1, which is directly relevant to the applicability of ELNEP. Since near equilibrium ELNEP is the first order approximation of true chemical kinetics equation, closer look into such domain is helpful to understand this example: Figure 11 shows the paths of true chemical kinetics and linear non-equilibrium for the same condition and initial composition as in Reference [24]. It can be seen from this figure that at the time of 10⁻⁵ s the path of linear non-equilibrium starts to nearly ride on the eigenvector of the Jacobian at equilibrium and from this riding point the path of linear non-equilibrium looks very close to the chemical kinetic path, which means that the final part of 1-D manifold of the kinetics is approximated by linear non-equilibrium path (the chemical kinetics path also nearly rides on the eigenvector at the same time 10⁻⁵ s, which suggests that the ride on the eigenvector may be one of the topological events that described later, but it needs more evidences). Therefore, it can be said that the linear non-equilibrium path in this example behaves in accordance with the theory, but it is clarified that ELNEP does not necessarily approximate the part of 1-D manifold of true chemical kinetics in some distance away from equilibrium. The reason for the latter is the highly distorted contour of σ whose eigenvectors largely vary from those at equilibrium, and consequently the eigenvector (that defines 1-D manifold) of the Jacobian of ELNEP (described in Appendixes F and I) considerably differs from that of the Jacobian of true chemical kinetics in that part. Also it has been numerically confirmed that for the two paths in Figure 11 the condition [39] that the method of ILDM is safely applicable (mentioned in Appendix F) does not hold at the part of 1-D manifold farther from equilibrium beyond the point of the time 10⁻⁵ s at least.

Now the conclusion of Reference [24] is recalled: Reference [24] states that a system’s dynamics cannot be deduced from the Gibbs free energy or entropy production rate even near the equilibrium state. However, ELNEP approximates at least the final part of 1-D manifold of true chemical kinetics as shown above, and ELNEP itself is the near-equilibrium approximation of true chemical kinetics equation, Equation (4) [40]. Accordingly, the statement of Reference [24] should be replaced here by “a system’s dynamics near the equilibrium state can be approximately deduced from the rate of entropy production. The approximation improves as each ε_l approaches zero.”

4. Relevance to LDMs

The phenomenon in which the concentration phase space of reaction kinetics is progressively reduced to a lower dimensional manifold can be interpreted as successive processes in which the number of quasi-equilibrium reactions increases. For example, Chiavazzo and Karlin [41] devised an algorithm to utilize the so-called quasi-equilibrium manifold as a first approximation to SIM. Regarding this hierarchy of manifolds of lower dimension [18–23,26–28], the observation that the paths determined by ELNEP trace at least final part of the 1-D and (possibly) 2-D manifolds of true chemical kinetics in the samples described in Section 3 is proved or, at least, theoretically justified to hold in more general cases where the method of intrinsic low-dimensional manifolds (ILDM) [26,39] is applicable. However, the proof below might be less formal from a mathematician’s viewpoint.

The coincidence of ELNEP with the true chemical kinetics equation to first order in each variable ε_l is obviously consistent with the ability of ELNEP to emulate a 0-D manifold, i.e., an equilibrium state, because the variable ε_l for every reaction vanishes at equilibrium. For higher dimensional manifolds, some elementary reactions participate in quasi-equilibrium (|ε_l|s are sufficiently small [33]), and the others do not; i.e., for each reduction in dimension of the reaction simplex (each step to descend the hierarchy), a group of elementary reactions constitutes quasi-equilibrium, and the size (number of constituent elementary reactions) of such groups increases as the system approaches its equilibrium. If the ascent on the hierarchy (in a time-reversed way from equilibrium) is observed, the size of such groups decreases as the system moves away from equilibrium because a split in a group occurs. The theoretical proof of the ability of ELNEP to trace every higher dimensional manifold when the method of ILDM is applicable to ELNEP is the following.

At first, ELNEP produces the trajectory that finally reaches the equilibrium state, i.e., a unique solution of ∇σ = 0; this is easily proved on the basis of the fact that the matrix Q of ELNEP is symmetric and positive definite, as described in Appendix E. Therefore, from any initial point in concentration phase space, paths determined by ELNEP reach the same equilibrium point as that of true chemical kinetics. From this starting point, the ability of ELNEP to emulate a 1-D manifold, 2-D manifold, etc. in a time-reversed way from equilibrium (ascent on the hierarchy) is proved inductively.

Base case

Along a 1-D manifold, all elementary reactions constitute one group in quasi-equilibrium and concertedly approach true equilibrium (all ε_ls approach zero). At this stage, the path determined by ELNEP obviously an approximation to the 1-D manifold of true chemical kinetics because otherwise the coincidence to first order in the variable ε_l for every lth reaction cannot hold between the two.

Inductive case

In the time-reversed way of thinking, at the entry (in ordinary time order) of 1-D manifold a group separates and its |ε_l|s grow larger. Now, consider the k -D manifold for m ≥ k > 1, where m is the dimension of the reaction simplex, and assume that the (k–1)-D manifold determined by ELNEP is an approximation to the (k–1)-D manifold of true chemical kinetics. Then, it can be proved that the k-D manifold of ELNEP is also an approximation to the k-D manifold of true chemical kinetics. The key to the proof is that the Jacobian and its eigenvectors of ELNEP are exactly the same as those of true chemical kinetics (in the form that explicitly accounts for the constraints inherent in chemical kinetics) at the equilibrium state, as described in Appendix I. Since each eigenvector of the Jacobian constitutes the local basis of the reaction simplex, the k-D manifold of ELNEP is created upon the (k–1)-D manifold when the k th group of reactions that appreciably contributes to a component of the basis (the eigenvector whose direction is the same as that of the most dominant component of acceleration in concentration phase space at this stage) separates (in a time-reversed way) from one of the previous groups that constitute the (k – 1)-D manifold. (The numbering of the groups of reactions also begins in a time-reversed way, i.e., the first group comprises all reactions and corresponds to a 1-D manifold.) The component of the basis is added as one degree of freedom to the trajectories that contract (in ordinary time order) into the (k – 1)-D manifold, and the k-D manifold of ELNEP is an approximation to that of true chemical kinetics because the eigenvector for ELNEP is approximately the same as that for true chemical kinetics even when the system is out of equilibrium, unless the contour of σ is highly distorted as the example in Section 3.4. The latter assumption holds when the method of ILDM is applicable to ELNEP (See Equation F20 and the passage above Equation F21). With this assumption, the respective eigenvectors for ELNEP and for true chemical kinetics are both approximately the same as those at equilibrium (but exactly the same for true chemical kinetics with a constant Jacobian) where the Jacobian and its eigenvectors of ELNEP are exactly the same as those of true chemical kinetics (Appendix I). In addition, since the eigenvalues for ELNEP are also approximately the same as those for true chemical kinetics, the order of the eigenvalues in magnitude and consequently the order of appearance of the manifolds are the same when the gaps among eigenvalues are sufficiently large. Thus, it is proved that the k -D manifold of ELNEP is an approximation to that of true chemical kinetics when the (k – 1)-D manifold does so (as long as the method of ILDM is applicable to ELNEP or the contour of σ is not highly distorted). Since the 1-D manifold (Base case) does so, the proof is inductively completed for all manifolds, although the path itself does not necessarily coincide even approximately with that of true chemical kinetics in manifolds higher than 1-D (as shown in Figures 3 and 5) because, for the k-D manifold (k > 1), the |ε_l|s of at least the k th group of reactions (that contributes to the basis of k -D manifold) are likely too large to approximate true chemical kinetics.

5. Relevance to Ziegler’s Hypothesis

The present results are nearly (not exactly) in accordance with Ziegler’s hypothesis [16] that the rate of entropy production is maximized with certain prescribed dissipative forces because (a) the path produced by ELNEP is an approximation to the true chemical kinetics path (though in the linear non-equilibrium domain or along at least a final part of the 1-D manifold) and (b) Equations (30)–(32) satisfy the so-called orthogonality conditions [17] that guarantees that the composition of the thermodynamic force X_l, each of which corresponds to the flux J_l, is orthogonal to the surface where σ is constant in the flux space [17].

6. Conclusions

The equation that approximately traces the trajectory in the concentration phase space of chemical kinetics is derived based on the rate of entropy production. In the derivation, a type of quasi-equilibrium approximation, which is automatically fulfilled in the relaxation process leading to detailed-balanced equilibrium, is used. The equation coincides with the true chemical kinetics equation to first order in a defined variable that characterizes the degree of quasi-equilibrium of each reaction. This new equation is named as ELNEP; the equation is justified as an approximation to true chemical kinetics along at least a final part of its 1-D manifold and identifies the trajectory that arrives at equilibrium in concentration phase space. The ELNEP is the equation of steepest descent path that sufficiently accounts for the constraints inherent in chemical kinetics. Although ELNEP is not a technique to extract low-dimensional manifolds, the trajectory of 1-D manifold in the concentration phase space of true chemical kinetics or at least its final part near equilibrium is approximately given by ELNEP because of the first-order coincidence with true chemical kinetics; whereas the simple steepest descent path whose Jacobian is the Hessian of σ does not sufficiently account for the constraints, and consequently cannot reproduce the 1-D manifold of true chemical kinetics except for the cases where the eigenvector of the Hessian of σ is nearly identical to that of the Jacobian of chemical kinetics.

In addition to the 1-D manifold, the paths determined by ELNEP approximately trace every higher dimensional manifold of true chemical kinetics when the contour of σ is not highly distorted or the method of ILDM is applicable to ELNEP, because the Jacobian and its eigenvectors of ELNEP are exactly the same as those of true chemical kinetics (in the form that explicitly accounts for the constraints inherent in chemical kinetics) at the equilibrium state. However, the path itself does not necessarily coincide even approximately with that of true chemical kinetics in manifolds higher than 1-D, as shown in the numerical calculations.

As for the time accuracy of ELNEP, each time of topological events such as the entry of lower dimensional manifold seems to be not badly approximated even in the domain far from equilibrium, but it must be further investigated.

Nomenclature

Here the symbols used only in Appendices D–J are omitted, where their own nomenclatures are consistently given for their respective equations.

B1, B2, B3	Matrixes defined by Equations (13) and (14) for element conservation
(B3)kj	Component of matrix B3 at the kth row and the jth column
c	Constant such as total number of a specified element per volume (a unit of concentration)
$\overrightarrow{c}$	Vector whose components are the constants such as the total number of respective elements per volume (a unit of concentration)
H	Hessian matrix of σ
Jl	Thermodynamic flux due to the lth reaction defined by Equation (31)
K	Number of independent species concentrations (the components of the vector ${\overrightarrow{y}}_{ind}$)
kbl	Backward reaction rate constant of the lth reaction
kfl	Forward reaction rate constant of the lth reaction
M	Number of constraints (such as the number of elements)
Nr	Number of reactions
Q	Matrix defined by Equation (12)
$\widehat{R}$	Universal gas constant
Rbl	Backward reaction rate of the lth reaction
Rfl	Forward reaction rate of the lth reaction
S	Entropy
s	Pseudo-time (an artificial time-marching parameter)
t	Physical time
V	Volume
X1	Thermodynamic force due to the lth reaction defined by Equation (32)
yA, yB, yC, yD	Concentration of species A, B, C, and D, respectively
yk	Concentration of the kth species
$\overrightarrow{y}$	Vector whose components are concentrations of respective species
${\overrightarrow{y}}_{dep}$	Vector whose components are concentrations of respective dependent species
${\overrightarrow{y}}_{ind}$	Vector whose components are concentrations of respective independent species
Greek symbols
αkl	Stoichiometric coefficient of the kth species as a reactant in the l th reaction
${\overrightarrow{\alpha }}_l$	Vector whose components are the stoichiometric coefficients of respective species as reactants in the lth reaction
βkl	Stoichiometric coefficient of the kth species as a product in the lth reaction
${\overrightarrow{\beta }}_l$	Vector whose components are the stoichiometric coefficients of respective species as products in the lth reaction
vkl	Difference of the two stoichiometric coefficients βkl − αkl
${\overrightarrow{v}}_{kl}$	Difference of the two vectors ${\overrightarrow{\beta }}_{kl}-{\overrightarrow{\alpha }}_{kl}$
εl	Quasi-equilibrium parameter of the lth reaction defined by Equation (17)
σ	Total reaction contribution to the rate of entropy production
σl	Contribution of the lth reaction to the rate of entropy production
Ω	Matrix defined by Equation (G2) in Appendix G
Suffixes
b	Backward reaction
dep	Dependent species
f	Forward reaction
ind	Independent species
l	lth reaction
k	kth species
•	An arbitrary species common in an equation where this symbol is used
(K×1)	Vector with K components for any K
(K×N)	Matrix with K rows and N columns for any K and any N

B1, B2, B3	Matrixes defined by Equations (13) and (14) for element conservation
(B3)kj	Component of matrix B3 at the kth row and the jth column
c	Constant such as total number of a specified element per volume (a unit of concentration)
$\overrightarrow{c}$	Vector whose components are the constants such as the total number of respective elements per volume (a unit of concentration)
H	Hessian matrix of σ
Jl	Thermodynamic flux due to the lth reaction defined by Equation (31)
K	Number of independent species concentrations (the components of the vector ${\overrightarrow{y}}_{ind}$)
kbl	Backward reaction rate constant of the lth reaction
kfl	Forward reaction rate constant of the lth reaction
M	Number of constraints (such as the number of elements)
Nr	Number of reactions
Q	Matrix defined by Equation (12)
$\widehat{R}$	Universal gas constant
Rbl	Backward reaction rate of the lth reaction
Rfl	Forward reaction rate of the lth reaction
S	Entropy
s	Pseudo-time (an artificial time-marching parameter)
t	Physical time
V	Volume
X1	Thermodynamic force due to the lth reaction defined by Equation (32)
yA, yB, yC, yD	Concentration of species A, B, C, and D, respectively
yk	Concentration of the kth species
$\overrightarrow{y}$	Vector whose components are concentrations of respective species
${\overrightarrow{y}}_{dep}$	Vector whose components are concentrations of respective dependent species
${\overrightarrow{y}}_{ind}$	Vector whose components are concentrations of respective independent species
Greek symbols
αkl	Stoichiometric coefficient of the kth species as a reactant in the l th reaction
${\overrightarrow{\alpha }}_l$	Vector whose components are the stoichiometric coefficients of respective species as reactants in the lth reaction
βkl	Stoichiometric coefficient of the kth species as a product in the lth reaction
${\overrightarrow{\beta }}_l$	Vector whose components are the stoichiometric coefficients of respective species as products in the lth reaction
vkl	Difference of the two stoichiometric coefficients βkl − αkl
${\overrightarrow{v}}_{kl}$	Difference of the two vectors ${\overrightarrow{\beta }}_{kl}-{\overrightarrow{\alpha }}_{kl}$
εl	Quasi-equilibrium parameter of the lth reaction defined by Equation (17)
σ	Total reaction contribution to the rate of entropy production
σl	Contribution of the lth reaction to the rate of entropy production
Ω	Matrix defined by Equation (G2) in Appendix G
Suffixes
b	Backward reaction
dep	Dependent species
f	Forward reaction
ind	Independent species
l	lth reaction
k	kth species
•	An arbitrary species common in an equation where this symbol is used
(K×1)	Vector with K components for any K
(K×N)	Matrix with K rows and N columns for any K and any N