Toward Formal Analysis of Thermodynamic Stability: Le Chatelier—Brown Principle

In this contribution, we carry on with the research program initiated in J. Math. Chem., 58(6), 2020. Using the methods from geometric thermodynamics, we formally derive and analyze different conditions for thermodynamic stability and determine the limits of their use. In particular, we study, in detail, several versions of the Le Chatelier—Brown principle and demonstrate their application to the analysis of thermodynamic stability.


Introduction
The conditions of thermodynamic stability belong to the basic principles of general and chemical thermodynamics and are of the utmost importance for various fundamental and applied problems. Recall, for instance, the quotation from Herbert Callen in [1]: "Considerations of stability lead to some of the most interesting and significant predictions of thermodynamics". Thus it is not strange that there are different variants of the stability conditions, both in classical monographs, e.g., [2][3][4] and also in research papers, e.g., [5][6][7][8][9][10][11]. The main difference between these formulations is the choice of the initial thermodynamic potentials, as well as the set of additional conditions imposed on the thermodynamic system. While most criteria are formulated in terms of the stability matrix, that is, the matrix formed by the second derivatives of thermodynamic potentials, there are several different approaches based upon expressing stability conditions in terms of inequalities involving partial derivatives of different thermodynamic potentials under different fixations of variables. Such approaches have a long history and date back to the works of Josiah Willard Gibbs [12] and Paul Ehrenfest [13]. In these papers, some initial formulations were proposed that were later extended and generalized. However, no formal mathematical analysis of these conditions and the limits of their applicability has been undertaken to date.
The idea of geometric interpretation of thermodynamics, originally suggested by Gibbs at the end of the 19th century, later formulated within the context of contact geometry by Robert Hermann in 1973 in [14], and further developed in a number of seminal works [15][16][17][18], forms a natural framework for a rigorous study of (equilibrium) thermodynamic systems. There have been a series of works extending this geometric approach to the study of dynamical evolution of thermodynamic systems, see, e.g., [19][20][21][22][23][24] and references therein. See also [25] for a comparison of different approaches.
In our research, we take a different route and apply the developed geometric approach to the rigorous mathematical analysis of stability conditions for an equilibrium thermodynamic system. In doing so, we further the research program initiated in [26]. Specifically, using the contact geometric formulation we describe the equilibrium energy manifold as an n-dimensional Legendre manifold of a particular thermodynamic differential 1-form, [27], and apply the tools from linear algebra and multivariable calculus, in particular the implicit function theorem and its ramifications, [28,29]. Along these lines, we not only formally derive some classical relations such as the Le Chatelier-Brown principle, but also carry out a rigorous mathematical analysis of these relations and determine the limits of their applicability. It is worth noting that the formal analysis of the studied relations require the use of rather non-trivial results from the matrix theory, some of which are not presented in the classical text books, e.g., [30], and hence were proved in the paper.
The paper is organized as follows: in Section 2, we introduce basic thermodynamic axioms and use them to derive some classical results from equilibrium thermodynamics. We also give a formal mathematical definition of generic thermodynamic stability and present initial results aimed at characterizing this property. In Section 3, the implicit function theorem is formulated and its application is illustrated. In particular, we give a previously obtained formulation of the Gibbs stability criterion and carry out a detailed analysis of possible situations that can occur when this criterion does not hold. Section 4 is devoted to the thorough analysis of the Le Chatelier-Brown principle. It is shown that while the respective inequalities can be used to indicate instability, they cannot guarantee that the system is stable. The paper is concluded with a discussion and two Appendices. Appendix A containing necessary results from the matrix analysis, while Appendix B contains proofs of some theorems that were removed from the main body of the paper to make it more straightforward.

Notation
We use straight capital letters to denote matrices and vectors, slanted capital vectors for thermodynamic variables and functions and small letters for auxiliary parameters and also for specific values of extensive thermodynamic parameters, as defined in Section 2. . , x n ) ∈ R n |x i > 0, i = 1, . . . , n}.

Internal Energy Function
We denote by Y = [Y 1 , . . . , Y n+1 ] ∈ R n+1 + the n + 1-dimensional column vector of (strictly positive) extensive parameters, e.g., S, V, m i (entropy, volume, and quantities of substances, respectively) and define Y 0 = U(Y 1 , . . . , Y n+1 ) to be a sufficiently smooth function that corresponds to the internal energy of the system. The intensive parameters X i , which correspond to T (temperature), −p (pressure), µ i (chemical potential of compound i), etc., are defined as The function U is postulated to be positively homogeneous of degree 1, i.e., we have By Euler's homogeneous function theorem we have Using rather straightforward manipulations, we recover the equation which is known as the Gibbs-Duhem equation. The first consequence of (3) is that the Hessian matrix of U, denoted by D 2 U, is singular: D 2 U · Y = 0. Finally, we will often make use of the Maxwell reciprocal

Stability Condition
Rather informally, we can say that a thermodynamic system is stable at some stateȲ if a sufficiently small perturbation applied to the system will be attenuated with time and will not lead to any change in the system state. We give a formal definition of local (generic) stability in terms of the Hessian of the internal energy. It should be noted that the proposed definition is formulated in terms of second-order partial derivatives and can be applied only when the respective matrix of second derivatives is well defined. Later on, we will consider the critical, i.e., boundary cases and point out possible approaches to analyzing them.
Beyond that, the suggested approach provides a characterization of stability only within a sufficiently small neighborhood of the point of interest. To infer the system behavior beyond that limit, higher order derivatives are to be analyzed. Keeping these reservations in mind, we proceed to the definition. This specific structure can be characterized in a number of ways. Below, we present a result that will be used later on. Theorem 1. Let M be an [k × k] positive semi-definite matrix. The following two conditions are equivalent:

1.
M has a single zero eigenvalue and v = v 1 , v 2 , . . . v k ∈ R k + is the corresponding right eigenvector; 2.

Proof. See Appendix B.
Note that the property of the energy Hessian with one zero eigenvalue is a structural invariant as follows from the Gibbs-Duhem equations. It is therefore of interest to determine the points where a second zero eigenvalue appears as such points correspond to the true tipping points, where the system loses its stability.
We can rewrite (4) to express the relation between u(y 1 , . . . , y n ) and U(Y 1 , . . . , Y n+1 ) as Applying the chain rule to (5), one can write the relationship between the first-order derivatives of both energy functions: Noting that This result agrees with the fact that the first-order derivative of a 1st order homogeneous function is a 0th order homogeneous function and hence, the first-order derivatives of both energy functions are equal. Following the same logic, one can obtain the relationship between the second-order derivatives: This result implies that the Hessian of the specific energy function normalized w.r.t. the variable Y i is equal to the principal submatrix D 2 U[i] taken with the factor Y i . Since Y i > 0, this implies the following result.

Lemma 1.
LetȲ ∈ R n + be a vector of extensive state variables. The following two statements are equivalent: The Hessians of specific energy function formulated with respect to the variables y i are positive definite at The fundamental problem of equilibrium thermodynamics consists of determining the conditions under which the system passes from stability to instability. This corresponds to the situation when the energy Hessian has at least two zero eigenvalues. We recall that the first zero eigenvalue is structurally determined, and hence cannot vary. Therefore, the loss of stability corresponds to the appearance of a second zero eigenvalue. The following theorem gives the required characterization.

Theorem 2.
Let M be an [k × k] positive semi-definite symmetric matrix. The following two conditions are equivalent:

1.
M has a zero eigenvalue of algebraic multiplicity 2 and one of the right eigenvectors corresponding to the zero eigenvalue has all nonzero components.
Proof. See Appendix B.
We will not consider the case when the transition from stability to instability corresponds to the simultaneous appearance of more than two zero eigenvalues, as such a situation is non-generic. To give an idea, we note that a space of all [3 × 3] real symmetric matrices has dimension 6 and the space of all such matrices with a double zero eigenvalue has dimension 3, while the set of [3 × 3] real symmetric matrices with 3 zero eigenvalues has dimension 0, as the only matrix satisfying this condition is the zero matrix.
Theorem 2 tells us that when studying the stability of the system with the energy function U(Y 1 , . . . , Y n+1 ), it is sufficient to consider any order n principal submatrix of the energy Hessian. When a zero eigenvalue appears, we conclude that the system state is located on the stability border. Now, recall that any such submatrix enjoys the same properties as the Hessians of the specific energy function taken with respect to some extensive variable. Therefore, it appears to be convenient to work with a Hessian of the specific energy function. Let us assume, for certainty, that this function is obtained by fixing the last extensive variable Y n+1 . This assumption will be the starting point for our subsequent analysis.
From now on, we adopt the following convention: the specific energy function and the specific extensive variables are denoted by lower-case letters, while the intensive variables are still denoted by capital letters. That is to say, we will write u(y 1 , . . . , y n ) and X i = ∂u/∂y i .

Implicit Description of the Equilibrium Energy Manifold and the Implicit Function Theorem
The energy equilibrium manifold can be considered as an n-dimensional manifold embedded in the n + 1-dimensional space and implicitly defined by Ψ(y 0 , y 1 , . . . , y n ) = y 0 − u(y 1 , . . . , y n ) = 0.
Alternatively, one can embed this manifold into a higher-dimensional space of dimension 2n + 1 with coordinates (y 0 , y 1 , . . . , y n , X 1 , . . . , X n ), which we will refer to as the thermodynamic state space. Following this approach, the energy equilibrium manifold U is defined by the following system of equations: We write (8) as Φ(y 0 , y, X) = 0, where Φ(y 0 , y, X) is a smooth vector-valued function, Φ : R 2n+1 → R n+1 . Obviously, (y 1 , . . . , y n ) form a set of independent variables, i.e., one can express y 0 and X i , i = 1, . . . , n in terms of y j , j = 1, . . . , n. However, other variables can be chosen to form a (local) coordinate system on U . These variables are characterized using the implicit function theorem. Below, we will briefly outline the mathematical apparatus of implicit function theorem and its application to deriving thermodynamic quantities of interest. For more details, the interested reader is referred to the paper [26].
The Jacobian matrix of Φ(y 0 , y, X) writes as where u i and u ij denote partial derivatives: u i = ∂u(y 1 ,...,y n ) ∂y i and u ij = ∂ 2 u(y 1 ,...,y n ) The outer upper row in (9) indicates the variables to which we differentiate. Let the point Ξ = (ȳ 0 ,ȳ 1 , . . . ,ȳ n ,X 1 , . . . ,X n ) belong to the equilibrium energy manifold U , i.e., it holds that Φ(Ξ) = 0. The implicit function theorem states that, in a sufficiently small neighborhood of Ξ, one can express n + 1 variables as functions of the remaining n variables if the matrix formed by the columns corresponding to these n + 1 variables is non-singular at Ξ.
The introduced framework allows for the efficient computation of partial derivatives of thermodynamic variables with respect to other variables if such partials exist. Suppose that we have determined that the set of (n + 1) variables Z can be expressed through a complementary set of n variables W, where Z and W form a disjoint partition of the set of all thermodynamic variables, both extensive and extensive, i.e., Z ∩ W = ∅, Z ∪ W = {y 0 , y 1 , . . . , y n , X 1 , . . . , X n }. The partial derivative ∂z i ∂w j (W), where z i ∈ Z and w j ∈ W, can be written as where C z i = col z i DΦ is a shorthand for the column of (9) indexed by the variable z i ∈ Z. Thus, the matrices in the preceding expression are formed by aligning the respective columns. Note that the choice of dependent and independent variables is crucial here. Sometimes, students dealing with thermodynamics get a feeling that "everything can be differentiated with respect to everything". Obviously, this is not quite true, as the subsequent analysis will show. However, there is a large number of possible variants and the particular choice of variables can prove to be decisive in analyzing the properties of a thermodynamic system.

Remark 1.
Note that if y 0 is considered to be a dependent variable, it does not alter the respective determinant, i.e., Thus, we can drop the first equation in (8) where we write DΦ(y 0 , y, X) [1] to denote the principal minor, obtained by deleting the first row and the first column.

Gibbs Stability Condition
In this subsection, we formally state the Gibbs stability criterion and expand upon its physical implications.

Theorem 3.
Let Ξ = (ȳ 0 ,ȳ,X) be a point on the equilibrium energy manifold U , such that det D 2 u i 1 (Ξ) = 0 for all i = 1, . . . , n − 1. Then, the determinant of the Hessian of u evaluated at Ξ is equal to the product of first-order partial derivatives where W i are defined in (12).
The straightforward corollary of Theorem 3 is that the thermodynamic system is locally stable at Ξ if the product of the first derivatives of thermodynamic potentials obtained for different choices of dependent and independent parameters is non-zero (positive). However, since we are interested in determining the conditions under which the system is located on the stability boundary, we might wish to consider, in more detail, the case when the product (13) is equal to zero. In particular, we would like to know if it is possible to focus on considering a single product term instead of the whole chain.
We recall that the product (13) can be written in terms of determinants of respective matrices as The numerator of the last term corresponds to the Hessian D 2 u and, hence, the last term turns equal to zero if there is a zero eigenvalue, that is the system is on the boundary of stability. On the other hand, the expressions in the denominators of the product terms correspond to the principal minors of the Hessian. If at least one minor equals zero, the whole expression turns out to be badly defined, as some of the fraction will take on an infinite value. It is thus of interest to understand the conditions under which the above situation may occur. We will approach this issue from a geometrical viewpoint. As said above, we will consider the generic case when there is a single zero eigenvalue. Let v be the eigenvector corresponding to the zero eigenvalue. Following Theorem 1, we conclude that all principal submatrices of D 2 u(Ξ) are non-degenerate if, and only if, the vector v does not contain zero components (note that the sign of the respective components is not relevant here).
It is a known fact that the eigenvectors of the energy Hessian define the principal directions of the paraboloid (respectively, its level sets) that locally approximates the function u(y 1 , . . . , y n ) at Ξ. At a regular point, this paraboloid is non-degenerate elliptic, i.e., all its sections are either ellipses or ellipsoids. This paraboloid gets degenerate at the point where its Hessian has a zero eigenvalue. The eigenvector corresponding to the zero eigenvalue is particularly important as it defines the direction in which the paraboloid "gets flat", that is, it becomes locally linear (see Figure 1 for an illustration of the difference). Whether this flat direction indicates an inflection point or just a local flattening depends on the higher order derivatives of u computed along this specific direction. Note that in Figure 1 and in the subsequent figures, we draw a supporting hyperplane (see, e.g., [31] (Ch. 5) for a formal definition) to illustrate the convexity property of a surface or the lack thereof. Intuitively, a function is locally convex at the point Ξ if the graph of the function in a neighborhood U(Ξ) of this point is located on the side of the supporting hyperplane and touches the hyperplane at the single point Ξ. In this sense, Figure 1b depicts a semi-convex function that touches the supporting hyperplane along a line. Furthermore, Figure 2 illustrates the difference between two cases that occurs when the Hessian turns out to have a zero eigenvalue. In this case, the decision about the stability or instability of the system at point Ξ should be made upon further analysis of the higher order derivatives. Figure 2a illustrates the case of a flat, but still convex, function, while Figure 2b depicts a non-convex function, with an inflection point determined by the third-order derivatives.
To illustrate the scenario where there is at least one zero component of the eigenvector corresponding to the zero eigenvalue, we go back to the original representation of the energy function. Letȳ = (ȳ 1 , . . . ,ȳ n ) be a state such that there exists v ∈ R n , satisfying D 2 u(ȳ)v = 0, i.e., v is the eigenvector corresponding to the zero eigenvalue. Furthermore, assume that the ith component of v is equal to zero. Then, the energy function U will be approximately flat for all points where ε is a sufficiently small number. We can immediately conclude that degenerate principal submatrices can appear in the case where the direction in which the energy surfaces flattens is such that one or more variables have to be constant. The preceding analysis indicates the following scenarios of the transition to instability. When the system state crosses the stability border, both the determinant det D 2 u(Ξ) and the respective product ∂X i ∂y i (W i ) Ξ turn into zero, while the principal minors remain positive. As the system states moves beyond the stability border, the determinant turns negative and the principal minors subsequently become equal to zero. In this case, the expression ∏ n i=1 ∂X i ∂y i (W i ) Ξ is not well defined any longer.
Alternatively, when the system is on the border of stability, both the determinant det D 2 u(Ξ) and some of the principal minors det D 2 u(Ξ)[i] become equal to 0. In this case, the product of the partial derivatives, as shown above, may not be well defined. We note that the latter situation is very specific and occurs only under the condition that the eigenvector, corresponding to the zero eigenvalue, has at least one zero element. Such a condition is not generic as it can be destroyed by an arbitrary perturbation of the system parameters. However, it may be of importance in the subsequent analysis of the higher order derivatives of the energy function. Finally, note that the transition to instability is always accompanied by zeroing of the determinant, det D 2 u(Ξ) = 0, as follows from the Cauchy interlacing theorem: if any principal minor of a symmetric positive semi-definite matrix M is equal to zero, then the matrix M has a zero eigenvalue.

A Classical Formulation of the Le Chatelier-Brown Principle
In this section, we consider a different concept, which is, however, somewhat closely related to the notion of thermodynamic stability. Namely, the Le Chatelier-Brown principle says that "when a system at equilibrium experiences external impacts, the latter cause counter-action forces to appear in the system that tend to decrease the effect of these impacts". This formulation resembles the definition of stability used in the theory of dynamical systems: "when disturbed from a condition of equilibrium, the system develops forces or moments that restore the original condition". The Le Chatelier-Brown principle remains valid for all regular states where the energy function has a strict minimum. However, it may not hold any longer if the energy manifold is not convex. This forms a connection to the results presented above.
We start by considering the variant of the Le Chatelier-Brown principle proposed by Paul Ehrenfest in 1909. It is formulated as an inequality expressed in terms of first-order partial derivatives of the same variables, but with a different choice of dependent variables: where ∂y i ∂X i y k = ∂y i ∂X i (y 1 , . . . , y k−1 , y k+1 , . . . , y n , X k ). Using the expression (10), we can compute the respective partial derivatives Note that the first partial derivative exists only under the condition that the Hessian D 2 u is non-singular, as otherwise the denominator would turn to zero. Using the Jacobi identity (A1), we get , that is, this version of the Le Chatelier-Brown principle boils down to saying that the ith diagonal element of the energy Hessian matrix is larger or equal to the inverse of the respective diagonal element of the original Hessian: At this point, Theorem A3 can be used to show that the inequalities (15), resp. (14) hold for all i = 1, . . . , n if the energy Hessian is positive definite, i.e., the system is stable. Recall that for the positive semi-definite Hessian D 2 u, the expression (14) is not defined any longer, as the first partial derivative may turn to infinity. Note that the result, converse to that stated in Theorem A3, is not true even under certain additional conditions such as, e.g., the positiveness of the diagonal elements of the matrix M. Thus, one cannot conclude that D 2 u is positive definite if (14) holds for all i. However, this condition can be used to reject specific cases. Namely, we can conclude the D 2 u is not positive definite if (14) does not hold for at least one i ∈ {1, . . . , n}.

An Extended Version of the Le Chatelier-Brown Principle
We have considered the formulation of the Le Chatelier-Brown principle in which we compared the cases where the set of independent variables consisted of all X's and of all Y's with a single X i . This setup can be extended by allowing Y-variables to substitutte the respective X-variables one by one. The respective sets of dependent (Z i ) and independent (W i ) variables would be as follows: . . , X n , y i } W n = {y 1 , . . . , y i−1 , y i+1 , . . . , y n , X i }.
Therefore, we can write down the partial derivatives for different sets of dependent and independent variables (16) as follows: . . . . . .
Now, we can formulate an extended version of the Le Chatelier-Brown theorem that was initially formulated in [32] and [33] (as well as in some other form, using second derivatives of thermodynamic potentials in [7,9]), albeit in a less formal way.
Theorem 4. If the system is stable at Ξ = (ȳ 0 ,ȳ,X) then for any i = 1, . . . , n the following chain of inequalities holds at Ξ and for the sets W i defined in (16): Proof. The proof is based upon the observation that any inequality can be written as for some choice if indices I. This inequality holds as follows from Theorem A4.
We first note that the formulation of the Le Chatelier-Brown principle as formulated in [13] (equation (14)) follows from Theorem 4. Similarly to the mentioned condition, Theorem 4 does not provide a sufficient condition for the stability of a system. However, similarly to the condition (14), Theorem 4 can be used to determine the cases where the system loses stability.

Discussion
In this paper, we presented further developments of the research program aimed at carrying out a rigorous mathematical analysis of classical relations and obtaining novel thermodynamic conditions in a systematic way. In particular, we analyzed the classical Gibbs stability condition and indicated several possible scenarios leading to the violation of these conditions, both accompanied and unaccompanied by a transition to instability. The obtained results form a basis for the analysis of high-order stability conditions that will be presented in the subsequent papers. We also formally derived several versions of the Le Chatelier-Brown principle and indicated the conditions under which these conditions can be used to conclude about the stability of a thermodynamic system. All presented results are based on the use of certain rather non-trivial results from matrix analysis that are formulated and proved in the paper. It is worth noting that the simplest cases following from the resulting equations are quite expected and conversant, e.g., it follows from Theorem 4 that the rate of changes in extensive properties increases with a decrease in the number of fixed intensive parameters. The discussion of some other more complicated cases would be a subject of our following works.
The presented results can be extended and generalized along several directions. One possible extension would consist in considering an irreversible version of the Le Chatelier-Brown principle following the line of research that was initiated in the works by I. Gyarmati [34] and later developed by M. Grmela and collaborators, see, e.g., [35]. Another possibility consists of making a bridge from the developed thermodynamic stability criteria to the second-law-based proof of Lyapunov stability in homogeneous nonequilibrium thermodynamic systems, see, e.g., [36]. These directions will be the subject of our future research.

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

Appendix A. Necessary Results from Matrix Analysis
In this paper, we will mostly use classical results from matrix analysis. The interested reader can refer to [30] for more details. Below, we will present some results that are either not available from classical textbooks or require some additional elaboration. Assume that rank B C > rank B and choose the vectors x and y of appropriate dimensions such that x B = 0, x C = 0, and x Cy = 0 (this is always possible because the matrix B C has less linearly dependent rows than B). Now, we can write the respective quadratic form x y B C C D x y = 2x Cy + y Dy, where y Dy ≥ 0 (note that D is a principal submatrix of M w.r.t. the indices N \ I). For any fixed y, one can choose x such that the sum 2x Cy + y Dy turns negative. This leads to a contradiction, hence we conclude that There is a particularly useful relationship between the minors of a matrix and its inverse, that is referred to as the Jacobi identity in [30], albeit without a proof. Here, we prove a special version of it that we will need when deriving results on the Le Chatelier-Brown principle.
Theorem A2. Let M ∈ R n×n be a non-singular matrix and let I ⊂ N = {1, . . . , n} and I c = N \ I be a set of indices and its complement. Then the following identity holds:  The following theorem presents an extension of the Jacobi identity. However, the authors are not aware whether it was mentioned previously.