Geometric Analysis of a System with Chemical Interactions

In this paper, we present some initial results aimed at defining a framework for the analysis of thermodynamic systems with additional restrictions imposed on the intensive parameters. Specifically, for the case of chemical reactions, we considered the states of constant affinity that form isoffine submanifolds of the thermodynamic phase space. Wer discuss the problem of extending the previously obtained stability conditions to the considered class of systems.


Introduction
The study of thermodynamic stability is one of the main and traditional tasks of general and chemical thermodynamics. The general stability criteria were proposed by Gibbs [1] at the end of 19th century, but further investigations showed the usefulness and significance of other various forms of thermodynamic stability conditions, especially for the systems subject to complex external influences and specific restrictions.
With reference to Gibbs, the deviations from equilibrium are typically assumed to be virtual, although in his fundamental work [1], Gibbs discussed the possibility of real deviations or, specifically, fluctuations. In Prigogine's later work [2], the analysis of stability conditions was carried out taking into account fluctuations, but the inclusion of such elements really requires a significant expansion of the theory, as was noted by Tisza [3]. The analysis of fluctuations should include finite deviations, and this opens up the possibility to consider the difference between stable and metastable states. On the other hand, such a consideration requires knowledge of the complete state diagram of the system state, and general conclusions cannot be obtained, cf. Münster [4]. Already Gibbs noted that an extension of his stability conditions to finite perturbations would require defining the boundaries within which they will be valid: ". . . it must be possible to assign limits within which it shall hold true of finite differences" [1], (p. 106).
Gibbs also proposed a variant of the stability condition that transforms the stability matrix into the product of second derivatives of internal energy w.r.t. the extensive variables at different fixings of the sets of intensive and extensive parameters [1], (p. 114). A thorough analysis of these conditions was performed in the papers [5,6] based upon a contact geometric formulation. The geometric framework has proved to be rather helpful for analyzing complex thermodynamic relations. Specifically, several additional results have been obtained within this framework-in particular, those related to the Le Chatelier-Brown principle. It was also indicated how these conditions can be used to conclude about the stability of a thermodynamic system.
In this paper, we aim at extending the developed approach to the thermodynamic systems with additional restrictions imposed on the intensive parameters. The most important, and by far the most practically relevant problem statement is presented by chemical reactions. In this case, the condition of constant affinity defines the isoaffine submanifolds of the thermodynamic phase manifold, thus resulting in a natural geometric formulation of the problem.
It is worth noting that the idea of using a geometric approach to analyzing thermodynamic relations was also suggested by Gibbs at the end of 19th century. Later, Hermann [7] recognized a deep connection between the equilibrium energy manifold and a Legendre manifold of the specific 1-form, later referred to as the thermodynamic 1-form. This contact geometry-based approach had been recognized and further developed in a number of papers [8][9][10]; see also more recent papers [11][12][13][14][15] and references therein. Almost concurrently with the first papers on contact geometric description of thermodynamic systems, the first substantial work on applying the geometric approach to chemical reaction dynamics was carried out by Oster, Perelson and Katchalsky in the 1970's; see, e.g., Oster and Perelson [16]. Later, van der Schaft and Maschke [17] extended these results and formulated the dynamics of reaction networks within the port-Hamiltonian framework.
A somewhat different approach consists in studying the structure of chemical reactions as determined by the stoichiometric matrix. This matrix possesses a number of nontrivial properties and has been subject to a thorough study for decades; see [18,19] for an overview of recent results. Quite remarkably, within the framework developed in [5,6], the matrix theoretic analysis of chemical reactions turns out to be well connected to the geometric study of the isoaffine submanifolds. This connection will be demonstrated below.

Notation
Throughout the paper, we use upright capital letters to denote matrices, slanted capital vectors for thermodynamic variables and functions, and small letters for auxiliary parameters and constants. Vectors are denoted either by upright capital letters, e.g., Y or X, or by lower-case bold letters, e.g., x, depending on the context. We write 0 and 1 to denote the vector (matrix) all of whose elements are equal to 0, resp. 1. It is assumed that the dimensions of all matrices are clear from the context. Otherwise, the dimensions are written as lower subscripts, e.g., 0 [2×k] .
Let U(X), U : R n → R, be a sufficiently smooth function. We let ∇U and ∇ 2 U denote the row vector of partial derivatives of U(X) (the gradient), and the [n × n] matrix of second-order partial derivatives of U(X) (the Hessian). Let M ∈ R (n×m) be a rectangular matrix; we then define the null and the range spaces as N (M) = {X ∈ R m |MX = 0} and R(M) = {Y ∈ R n |Y = MX, X ∈ R m }, respectively. Respectively, the co-range and the co-null spaces are coR(M) = N (M ), and coN (M) = R(M ). Note that for a symmetric matrix, the range and co-null spaces coincide, as well as the null and co-range spaces. Furthermore, we write M[I, J] to denote a submatrix of the [n × n] matrix M, obtained by removing from M the rows with indices i ∈ I ⊂ {1, . . . , n} and the columns with indices j ∈ J ⊂ {1, . . . , n}. The notation M[I, ∅] (resp., M[∅, J]) means that only rows (resp., columns) are removed.

Geometric Representation of the Energy Manifold
Let U(Y) be a sufficiently smooth first-order positive homogeneous function of n variables Y ∈ R n + . We define n conjugate variables X as X i = ∂U(Y) ∂Y i . Within the thermodynamic context, the function U(Y) is referred to as the internal energy function, the vector Y = [S, V, N 1 . . . , N n−2 ] contains the extensive parameters S, V, N i (entropy, volume, and the molar numbers of substances, respectively), while the vector X = [T, −P, µ 1 . . . , µ n−2 ] contains n intensive parameters T, −P, µ i (temperature, pressure, and chemical potentials of respective components). In the following, we will mostly use the generic notations Y i and X i for the extensive and intensive variables, respectively, and will use the specific notation S, T etc. to emphasize the role of particular variables. We also note that our notation implies that there are n − 2 different components (substances); hence we assume that n ≥ 4 to ensure that the problem is not degenerate.
The evolution of the system can be represented by the trajectory of a point on a co-dimension 1 manifold in the space of n + 1 variables (Y 0 , Y), where Y 0 = U(Y). The position of this point is uniquely specified by the values of n + 1 extensive variables (Y 0 , Y) (note that Y 0 = U is an extensive variable as well because it is determined by a first-order homogeneous function). However, for the purpose of formal analysis, it is desirable to have a formulation that explicitly encompasses not only the extensive, but also the intensive variables X.
In [5,6], the authors proposed a systematic approach to the analysis of thermodynamic systems based upon the apparatus of contact geometry. Below, we will present a brief overview thereof.
The suggested approach consists in embedding the n-dimensional energy manifold into a (2n + 1)-dimensional space with coordinates (Y 0 , Y, X), which is referred to as the thermodynamic phase space. The energy equilibrium manifold U (which can be formally seen as a Legendre submanifold of a particular thermodynamic 1-form, see Gromov and Caines [20] for details) is thus defined by the following system of equations: System (1) can be compactly written as (1) along with the implicit function theorem can be used to derive many useful thermodynamic relations, as was demonstrated in [5]. The application of the implicit function theorem is based upon the computation of the Jacobian matrix of Φ(Y 0 , Y, X). We write it concisely as where the symbols ∇ and ∇ 2 are used to denote the gradient (the row-vector of partial derivatives) and the Hessian (the matrix of second-order partial derivatives) of the energy function U, respectively. The matrix DΦ has a full row rank, i.e., rank DΦ = n + 1. Note that the matrix ∇ 2 U has a rank equal to n − 1 as follows from the Gibbs-Duhem equation (see [5] for the formal derivation). Let the point Ξ = (Ȳ 0 ,Ȳ,X) belong to the equilibrium energy manifold U , i.e., it holds that Φ(Ξ) = 0. Thus, according to the implicit function theorem, in a neighborhood of Ξ one can express n + 1 variables as functions of the remaining n variables. However, the application of the implicit function theorem goes far beyond that. It allows one to compute the partial derivatives of some dependent parameters with respect to certain independent parameters for different choices of independent parameters, thus providing a detailed characterization of the energy manifold curvature. In [6], such calculations were used to derive equilibrium conditions for a thermodynamic system not undergoing a chemical reaction.

Geometry of an Isoaffine Submanifold
We wish to expand the developed framework by considering the setup, where there are chemical interactions between the components of the system. Suppose that the individual components are involved in k reactions, k < n − 2. Each reaction is described by the balance equation that for the ith reaction reads as follows: where R j are the chemical reagents (species), and ν ij are the stoichiometric coefficients, which are used to balance the reaction (positive values for reaction products and negative for initial reactants). The entropy production σ for a unit volume V due to this nonequilibrium process is presented as follows [2]: where ξ i is the extent of the ith reaction, t is time, and A i is the affinity of the ith reaction, which is defined as: Affinity A i is associated with the thermodynamic force, while the velocity of the ith reaction, dξ i dt , is interpreted as the flow. Accordingly, the progression of the ith reaction can be characterized both by ξ i and A i . However, the extent of the ith reaction is not a thermodynamic variable [4]: it describes only the deviation from the equilibrium. In contrast to the extent of reaction, chemical affinity, defined by (3), is a thermodynamic variable (a function of chemical potentials). Formally, affinity of a reaction can be determined at any stage of the reaction by introducing an inhibitor (negative catalyst), that is, by "freezing" the reaction [4]. This implies that one can consider the states of constant affinity, i.e., isoaffine submanifolds of the thermodynamic phase space, which are defined by the equations A i = const.
We are specifically interested in studying the properties of the system when the chemical reactions are in equilibrium, i.e., when affinity is zero. Alternatively, this can be formulated as: ∑ n−2 j=1 ν ij µ j = 0, which can be concisely written in vector-matrix notation as Here, µ = [µ 1 , . . . , µ n−2 ] is the vector of chemical potentials, and Σ = [ν ij ], i = 1, . . . , k, j = 1, . . . , n − 2 is the matrix of stoichiometric coefficients. We assume that there are k < n − 2 independent reactions, which implies that the matrix Σ has full row rank, i.e., rank Σ = k. Note that the assumption of zero affinity is not restrictive in any way, as all of the subsequent analysis will be focused on the structure of the stoichiometric matrix, which remains unchanged. The conditions (4) determine an isoaffine submanifold U A of the energy manifold U , which we will later refer to as the isoaffine energy manifold. To analyze the dimensionality of this manifold one has to check the rank of the Jacobian matrix of the extended system of constitutive equations formed by (1) and (4). We denote the extended system by Φ and write the Jacobian as: where Σ = [0 [k×2] , Σ], i.e., we pad Σ with two zero columns. Obviously, both Σ and Σ have the same row rank. The co-dimension of the manifold U A at a point Ξ = (Y 0 , Y, X) is equal to the rank of D Φ(Y 0 , Y, X) at that point. One can easily see that the rank of D Φ is determined by the rank of the principal submatrix obtained by removing the first row and the column of D Φ.
The following theorem provides the required characterization.

Theorem 1. Consider the [(n +
where Σ ∈ R (k×n) is a full row rank matrix, i.e., rank Σ = k < n − 2, and ∇ 2 U is a symmetric matrix with a rank deficiency equal to 1, i.e., rank ∇ 2 U = n − 1. Then, the rank of M is equal to n + k − 1 if N ( Σ) ⊂ R(∇ 2 U) and n + k otherwise.

Proof. See Appendix A.
To interpret the formulated theorem we recall that the null space N (∇ 2 U) changes as the system state moves along the manifold and coincides with αȲ, α = 0, at a point Ξ = (Y 0 ,Ȳ,X). Furthermore, note that R(∇ 2 U) ⊥ N (∇ 2 U) =Ȳ as ∇ 2 U is a symmetric matrix. Thus the condition N ( Σ) ⊂ R(∇ 2 U) is equivalent to saying that N ( Σ) ⊥Ȳ, which, in turn, implies thatȲ belongs to the span of the rows of matrix Σ. Note, however, that the first two elements of each row of the matrix Σ are zero, which implies that the manifold U A loses its rank at the points whereȲ 1 =Ȳ 2 = 0, which do not belong to the energy manifold U . Thus, we eventually conclude that the manifold U A has a (full row) constant rank equal to n + k + 1, where one dimension comes from the removed first row.

Dependent and Independent Thermodynamic Variables on an Isoaffine Energy Manifold
The analysis carried out in the previous subsection implies that on the manifold U A one can choose n + k + 1 dependent variables that can be expressed through the remaining n − k ones. To put it differently, in a small neighborhood of any point Ξ ∈ U A , one can choose exactly n − k independent coordinates.
To provide an intuition for this result, assume that we have chosen variables (S, V, µ) as independent ones and the remaining ones as dependent. If the variables µ are consistent with (4), one can express k chemical potentials in terms of the remaining (n − 2 − k) ones. Thus, the set of independent variables reduces to n − k: S, V and (n − k − 2) independent chemical potentials. However, the problem of determining which variables can serve as independent ones is in general rather intricate as will be discussed below.
In the following, we will assume that the internal energy U is considered as a dependent variable. We thus restrict ourselves to considering 2n thermodynamic variables Y and X. The problem of determining which variables can be taken as dependent boils down to finding n + k linearly independent columns of the [2n × 2n] matrix (6). The first n columns correspond to the extensive variables Y, while the remaining columns correspond to the intensive variables X.
One can immediately observe that-in contrast to the case without chemical reactionsthe intensive variables have a different status: while the variables T and −P are "free", the chemical potentials µ j are "bound" by (4). On the other hand, the chemical potentials are not in equal rights as well: some of them can be expressed through others, while some cannot. One should also take into account that a chemical potential, in our case the internal energy U, is determined up to an arbitrary constant; this is the property of any potential.
Note that the value of µ j is usually defined through the standard value for the pure component j at certain constant values of T and P: where µ mix j is the change in the chemical potential as a result of mixing the original pure substance with other substances at given T and P. Thus, the affinity (3) is also a function of T and P. Formally, an additional link between chemical potentials, temperature and pressure could also be imposed by assuming that, e.g., but such constructions go beyond the scope of realistic physical framework. The change in the standard value µ • j (T, P) can also be derived from the Gibbs-Helmholtz equation where G is the Gibbs energy and H is enthalpy. For pure substance, the molar Gibbs energy is equal to the chemical potential, but the determination of this value requires knowledge of the temperature dependence of enthalpy. As a result, no general conclusions can be obtained along this line. Before proceeding with the detailed analysis, observe that the whole set of extensive variables Y cannot be chosen to be dependent, as the columns of the matrix ∇ 2 U are linearly dependent. However, under the thermodynamic stability condition, cf. [5], any [(n − 1) × (n − 1)] principal submatrix of ∇ 2 U has full rank. This allows us to formulate the first result.

Proof. See Appendix A.
Note that the set of linearly independent columns of the matrix Σ is not uniquely determined. Thus, depending upon the structure of the matrix Σ, the set J can be chosen in different ways, which introduces additional degrees of freedom.
Theorem 2 can be extended as the following result shows.

Proof. See Appendix A.
Note that in all preceding theorems, the set of chemical potentials {µ j } j∈J was always chosen to be dependent. These variables can be chosen to be independent as well, but this choice would not be valid for the whole isoaffine submanifold. A characterization of the subsets of U A , where the respective columns of the Hessian matrix (6) turn out to be linearly dependent, is closely related to the problem of stability analysis of the respective system and will be reported in detail in a subsequent publication.

Discussion
In this paper, we presented some initial results aimed at defining a framework for the analysis of systems with additional restrictions imposed on the intensive parameters. Specifically, for the case of chemical reactions, we considered the equilibrium states or, alternatively, the states of constant affinity as the most important setup for physical chemistry. However, using the developed approach, one can consider the situations where some other links, e.g., between T and P, are imposed.
The presented results form a rigorous mathematical framework for the formulation and analysis of stability conditions in the course of chemical reactions, including chemical equilibrium. We emphasize that we are not discussing the stability of the chemical equilibrium itself, but rather the stability of the system as a whole with additional coupling equations (in this case, for the constancy of chemical affinity). These and similar problems have been considered in some works in the traditional thermodynamic way and for specific cases. Some examples may be the works of Gitterman for chemical equilibrium [21,22], Prigogine on configurational heat capacity (e.g., [2,23]), or the recent work by Toikka [24]. However, it is believed that the presented formulation would allow for a more systematic analysis of the problems of stability and equilibrium displacement using well-developed methods of vector algebra and matrix analysis.  Data Availability Statement: All necessary data are contained in the paper.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Proofs of the Theorems
Proof of Theorem 1. Note that rank M = 2n − dim(N (M)). Therefore, it suffices to compute the dimension of the null space of M. Consider the vectors x = x 1 x 2 ∈ R 2n , where x 1 , x 2 ∈ R n . A vector x belongs to N (M) if Mx = 0, which amounts to the following system: By definition, there are exactly (n − k) linearly independent vectors x 2 satisfying (A1b), {x i 2 }, i = 1, . . . , n − k. For each such vector, the system (A1a) has a solution if x 2 ∈ R(∇ 2 U). If this is the case, we have N (Σ) ⊂ R(∇ 2 U), that is, the null space of Σ forms a linear subspace of the range space of ∇ 2 U. Note that k < n − 2 and hence, this is always possible. Let us pick up a set of (n − k) orthogonal vectors x i 1 satisfying (A1a) for the corresponding x i 2 . The set of vectors x i , obtained by vertically stacking the vectors x i 1 and x i 2 , forms an (n − k)-dimensional subspace of R 2n . This set can be extended by the vector . This vector is linearly independent with the previously defined vectors as x ⊥ ⊥ coR(∇ 2 U) x i 1 . Noting that dim(N (∇ 2 U)) = 1, we conclude that there are no more linearly independent vectors that can be added to the =set. Therefore, dim(N (M)) = k + 1 and rank(M) = 2n − (n − k + 1) = n + k − 1. This proves the first part of the statement.
To prove the second part, we note that since rank ∇ 2 U = n − 1, N (Σ) ⊂ R(∇ 2 U) implies N (Σ) + R(∇ 2 U) = R n . Using the Grassmann theorem, we conclude that dim(N (Σ) ∩R(∇ 2 U)) = 1. The latter implies that there are (n − 1) linearly independent vectors x i 2 ∈ N (Σ) satisfying (A1a) for some (linearly independent) vectors x i 1 ∈ R(∇ 2 U). Extending this set by the vector x ⊥ , we obtain a set of (n − k) linearly independent vectors that form the basis of N (M). Computing the dimensions, we obtain the required result.