A Hessian Geometric Structure of Chemical Thermodynamic Systems with Stoichiometric Constraints

We establish a Hessian geometric structure in chemical thermodynamics which describes chemical reaction networks (CRNs) with equilibrium states. In our setup, the ideal gas assumption and mass action kinetics are not required. The existence and uniqueness condition of the equilibrium state is derived by using the Legendre duality inherent to the Hessian structure. The entropy production during a relaxation to the equilibrium state can be evaluated by the Bregman divergence. Furthermore, the equilibrium state is characterized by four distinct minimization problems of the divergence, which are obtained from the generalized Pythagorean theorem originating in the dual flatness. For the ideal gas case, we confirm that our existence and uniqueness condition implies Birch's theorem, and that the entropy production represented by the divergence coincides with the generalized Kullback-Leibler divergence. In addition, under mass action kinetics, our general framework reproduces the local detailed balance condition.

Historically, the conventional equilibrium chemical thermodynamics was established by the seminal papers by Gibbs [22], in which the chemical equilibrium state is variationally and globally characterized as the state to minimize the free energy. In the same period, the chemical kinetic theory was also being developed in parallel with Gibbs' equilibrium chemical thermodynamics.
By combining the law of mass action by Guldberg and Waage in 1864 [23] with Boltzmann's characterization of the equilibrium state by detailed balancing, Wegscheider clarified the condition which the rate constants of chemical reactions must satisfy to have equilibrium states [24]. The characterization of an equilibrium state by the detailed balancing of the reaction fluxes is kinetic and local but consistent with the global free-energy characterization of the equilibrium state under the ideal-gas or dilutesolution assumption [25].
Since then, theories of chemical reaction systems and their thermodynamics have been developed, mainly based on the detailed balancing characterization of equilibrium states and mass action systems. For example, in the 1970s, the chemical reaction network (CRN) theory emerged [26]. Here, Horn and Jackson formalized the complex-balanced state of a mass action system, which extended the uniqueness and stability of the equilibrium state [27]. In relation to the stability of the equilibrium state, it was found that the Gibbs free energy difference is identical to the generalized Kullback-Leibler divergence (also known as the pseudo-Helmholtz potential) and behaves as a Lyapunov function of mass action systems [27][28][29][30]. This result could be regarded as a chemical version of Boltzmann's H theorem and was used for proving the convergence of a mass action system to the unique equilibrium state characterized by detailed balancing [31,32]. Also, around 1970, Hill and Schnakenberg extended the theory to stochastic linear reaction systems [2,33,34].
The applicability of chemical thermodynamics has recently been extended further in various ways. The generalized mass action kinetics was proposed in the field of applied mathematics as a broader class of kinetics in which the properties of the equilibrium state can be conserved [18-20, 35, 36]. The authors of Refs. [14][15][16][17]37] established a theory for open CRNs and derived the conditions under which an open system has an equilibrium state .
However, most of recent developments are not purely thermodynamic, because they are based on the characterization of the equilibrium state by detailed balancing, assuming a specific kinetics analogous to the mass action. As a result, it is unclear whether the results are obtained just by mathematical extensions that crucially depend on the specific kinetics of the mass action and its variants, or whether they are truly consistent with the general framework of thermodynamics. If they are consistent, the same results must be derived from a purely thermodynamic argument in the line of Gibbs without assuming kinetics and detailed balancing. Once the consistencies are confirmed, one could apply the previously-established results to a much broader class of non-ideal and non-mass action chemical systems, because thermodynamics can describe the properties of systems independent of the details of their kinetics. However, there have been few attempts to establish the link [36]. In this work, we reveal the link in chemical thermodynamics of general CRNs.
Before outlining our main results, we recall the general framework of thermodynamics, which should be formulated as follows. The space of extensive variables is endowed with a concave function, called entropy [1,3]. According to the second law of thermodynamics, a system should evolve with time such that the total entropy function increases under any given constraints imposed on the extensive variables. If the constraints are trivial, the system converges to the maximum of the total entropy function, namely, the equilibrium state. Furthermore, there is a well-established procedure to evaluate the dissipation during a relaxation to the equilibrium state.
However, the constraints can have non-trivial impacts when such a general theory is applied to CRNs, especially to those with complex stoichiometry of chemical reactions. In a chemical thermodynamic system, the extensive variables include the numbers of molecules constituting the system. These numbers cannot change independently as they are algebraically constrained by the stoichiometry. It is these constraints which provide significant geometric structure into the problem, and also yield several properties obtained from detailed balancing and mass action kinetics.
[Outline of Main Results] In this paper, we develop a thermodynamic theory for chemical reaction systems with complex constraints. With this theory, we obtain the following four main results. Theorem 1. Necessary and Sufficient Condition for Existence of Equilibrium States for Open CRNs: This is a generalization of the Wegscheider condition and the equilibrium condition of open CRNs obtained recently in Refs. [14,15].
Theorem 2. Uniqueness Condition of the Equilibrium State: This is a generalization of the uniqueness condition of the equilibrium state obtained so far under the assumption of mass action kinetics.
Theorem 3. The Difference of Total Entropy between the Equilibrium State and Any State is Evaluated by Bregman Divergence: This is a generalization of the fact that the free energy difference is identical to the generalized Kullback-Leibler divergence. In particular, we clarify that a convex function characterizing the Bregman divergence corresponds to the thermodynamic potential of the system. Theorem 4. Four Variational Characterizations of the Equilibrium State: One of them is the generalization of the variational characterization of the equilibrium state as the minimizer of the free energy (the generalized Kullback-Leibler divergence). The other three are newly obtained as a result of our Hessian geometric formulation.
We emphasize again that these results are derived purely thermodynamically, without using any specific kinetics such as the mass action laws and the local characterization of equilibrium states such as detailed bal-ancing. In particular, we derive these generalizations by identifying and employing the Hessian geometric structure [38] in constrained thermodynamic systems. The Hessian geometry of thermodynamic systems plays an essential role, when the constraints between the variables become non-trivial and complex. This paper is organized as follows. We devote Sec. II to review the conventional thermodynamics of CRNs in light of the entropy maximization problem. Also, we recapitulate how different kinds of thermodynamic potentials are linked to each other. In Sec. III, we derive the existence and uniqueness condition for the equilibrium state. To this purpose, we introduce two spaces which are connected to each other by Legendre duality. This pair of the spaces and their duality is the basis of the Hessian geometric structure. The equilibrium state is then uniquely determined by the intersection of two submanifolds (see Theorem 2). Sec. IV reformulates the second law from the geometric point of view. We show that the dissipation can be evaluated by the Bregman divergence (see Theorem 3). Furthermore, we obtain four distinct characterizations of the equilibrium state (see Theorem 4). In Sec. V, we demonstrate our geometric structure in ideal gas cases to rederive the previously known results: our theorems for the existence and uniqueness condition reduce to Birch's theorem (see Theorem 5), and the entropy production can be represented by the generalized Kullback-Leibler divergence. In addition, by assuming the law of mass action, we reproduce the local detailed balance condition. Finally, we summarize our results with further discussions in Sec. VI.
While we derive all the results without the assumptions of mass action systems and detailed balancing, the assumptions are familiar to researchers working on stochastic thermodynamics and CRNs. In our accompanying paper [39], we reproduce our results for the special case starting from mass action systems and detailed balancing.

II. CONVENTIONAL THERMODYNAMICS FOR CHEMICAL REACTION SYSTEMS
In this section, we recall conventional thermodynamics for chemical reaction systems. Readers who are familiar with the topic can skip to Eqs. (23), (24) and (25).
Consider a thermodynamic chemical reaction system surrounded by a reservoir. We assume that the system is always in a local equilibrium state, i.e., a well-mixed state, and therefore we can completely describe it by extensive variables (Ω, E, N, X). Here, Ω and E represent the volume and the internal energy; N = {N m } denotes a vector, each component of which is the number of the corresponding open chemical. The open chemicals can diffuse across the boundary between the system and the reservoir. By contrast, X = X i is the numbers of chemicals confined in the system; the indices m and i run from m = 1 to N N and from i = 1 to N X , respectively, where N N and N X represent the numbers of species of the open and confined chemicals. Since we only discuss isochoric cases (i.e., Ω = const.) for theoretical simplicity [40], we employ the density variables (ǫ, n, x) = (E/Ω, N/Ω, X/Ω). In thermodynamics, a concave, smooth and homogeneous function Σ, which is called the entropy, is defined on (Ω, E, N, X). Owing to the homogeneity of the entropy function, without loss of generality, we can write it as where σ [ǫ, n, x] represents the entropy density. In this work, we additionally assume that σ [ǫ, n, x] is strictly concave, which implies a situation without phase transitions from the physical point of view. The reservoir is characterized by intensive variables (T ,μ), whereT is temperature andμ = {μ m } are chemical potentials corresponding to the open chemicals; also we denote the corresponding extensive variables by (Ẽ,Ñ ). We denote the entropy function for the reservoir byΣT ,μ [Ẽ,Ñ ], and therefore the total entropy can be expressed by where we use the additivity of the entropy. Next, we define a dynamics as  (70)). The index r runs from r = 1 to N R , where N R is the number of reactions. In this paper, we employ Einstein's summation convention for notational simplicity.
Since, in most cases, the time scale of reactions is much slower than the others (that is, i B (t) , k B (t) ≫ j (t)), we can analyze the dynamics, Eq. (3), by separating it into the fast and slow scales. By employing a scaling: we obtain the fast scale effective dynamics as where we use the fact that the diverging bare flux densities, i B and k B , converge to the finite effective flux densities, i and k, in the scaling limit. The formal solution T ,μ ẼÑ (ǫ, n, x) of Eq. (4) with an initial condition (ǫ 0 , n 0 ,Ẽ 0 ,Ñ 0 ) can be represented as where ι (τ ) and κ (τ ) are the integrals of i (τ ) and k (τ ) with the initial condition ι (0) = κ (0) = 0: Here, we note that the densities of confined chemicals, x, can be regarded as a constant in the fast dynamics. By substituting the solution, Eq. (5), into Eq. (2), we have the time evolution of total entropy as where we use properties of the reservoir, Ωι (t) ≪ E 0 , Ωκ (t) ≪Ñ 0 , and the Taylor expansion forΣT ,μ ; we also employ the thermodynamic relations: ∂ΣT ,μ /∂Ẽ = 1/T and ∂ΣT ,μ /∂Ñ m = −μ m /T . Although the constant term is explicitly given asΣT ,μ [Ẽ 0 ,Ñ 0 ], we abbreviate it to "const.", because it never affects the theoretical framework.
To introduce thermodynamics into our dynamics, we briefly summarize its significant statements. According to the first law, a heat dissipation Q 0→τ from the system to the reservoir during a time interval [0, τ ] is given by the entropy increment in the reservoir: where Ωι (τ ) represents the internal energy gain of the system and −Ωμ m κ m (τ ) is the work done by the system through the injection of chemicals into the reservoir.
The second law states that, for spontaneous changes, the flux density functions, i (τ ) and k (τ ), must be chosen such that Σ tot (τ ) becomes an increasing function with respect to time τ [41]. In other words, the system climbs up the landscape defined by the concave function Σ tot (ι, κ, x) with respect to ι and κ in the time evolution, and finally converges to the maximum if it exists. If we write (ǫ (τ ) , n (τ )) → (ǫ QEQ , n QEQ ) for τ → ∞, candidates of the converged state (ǫ QEQ , n QEQ ) can be evaluated by a variational form: and (ǫ QEQ , n QEQ ) = (ǫ 0 + ι QEQ , n 0 + κ QEQ ). In thermodynamics, the states maximizing the total entropy are called equilibrium states; therefore, (ǫ QEQ , n QEQ ) is an equilibrium state in the fast dynamics. However, we call it a quasi-equilibrium state, because we will treat the slow dynamics later. By using the argument shift, ǫ = ǫ 0 + ι, n = n 0 + κ, we can rewrite the variational form as and we directly obtain candidates of the quasiequilibrium state.
Since, in thermodynamics, the function maximized in Eq. (10) is bounded above [42], the quasi-equilibrium state always exists. Furthermore, since we have assumed the strict concavity for σ in this work, we can conclude that the quasi-equilibrium state is uniquely determined by Eq. (10); and, for an arbitrary initial condition (ǫ 0 , n 0 ), the system converges to the unique quasi-equilibrium state (ǫ QEQ , n QEQ ). The above conclusion, which is the existence and uniqueness of the quasiequilibrium state, originates from the simplicity of the fast dynamics, Eq. (4). In other words, the maximization is easily conducted, because ǫ and n can be varied independently. As shown later, the conclusion no longer holds for the slow reaction dynamics, because of complex stoichiometric constraints. Also, the total entropy at the quasi-equilibrium state can be represented as (11) Employing the above results for the fast dynamics, we analyze the slow dynamics, which is the chemical reaction dynamics. Owing to the variational form, Eq. (10), the time evolutions of the densities of the internal energy ǫ (t) and of the open chemicals n (t) in the slow dynamics are already solved. By using the time evolution of the confined chemicals x (t), we have (12) Substituting these equations into Eq. (3), we obtain the effective slow dynamics as The formal solution of Eq. (13) with the initial condition x 0 is represented as where ξ (t) = {ξ r (t)} is the integral of j (t) with the initial condition ξ (0) = 0. The vector ξ (t) is the density of the extent of reaction. Also, the initial conditions of the reservoir for the slow dynamics,Ẽ (0) andÑ (0), can be calculated from the fast dynamics as The substitution of Eqs. (12) and (14) into Eq. (2) leads to the time evolution of the total entropy in the reaction dynamics: where we use the Taylor expansion forΣT ,μ again. If we use the quasi-equilibrium entropy function Σ tot QEQ = Σ tot QEQ (T ,μ; x (t)) given by Eq. (11), we can rewrite Eq.
(17) From the second law, an equilibrium state in the reaction dynamics is evaluated by a variational form: ξ EQ = arg max ξ Σ tot (ξ), and x EQ = x 0 + Sξ EQ ; also, the equilibrium total entropy is Σ tot EQ = max ξ Σ tot (ξ). Furthermore, by following the same argument as in Eq. (8), the heat dissipation of this dynamics is given by denotes the system entropy density at the quasi-equilibrium state with the confined chemicals x (t).
The representation of the total entropy, Eq. (17), may be unfamiliar to the reader, therefore we rewrite it by employing thermodynamic potentials. First, consider the maximization in Eq. (11) with respect to ǫ: which is called the Massieu potential density. By using this potential, the Helmholtz free-energy density is defined as which coincides with the partial grand potential density. Owing to the definition of ϕ[T ,μ; x], the quasiequilibrium entropy function Σ tot QEQ can be represented as and therefore Eq. (17) can be rewritten in a familiar form: . Hence, the heat dissipation, Eq. (18), can be expressed as Since all important thermodynamic quantities for the reaction dynamics can be calculated from the potential ϕ[T ,μ; x], we will use it, instead of the entropy density σ [ǫ, n, x], hereafter. Before closing this section, we consider the equilibrium state of the slow reaction dynamics. Owing to the second law, candidates of the equilibrium state are given by the variational form: and x EQ = x 0 + Sξ EQ . However, differently from the case in the fast dynamics (see Eq. (10)), the existence and uniqueness of the equilibrium state are not guaranteed in this case because of S and O. In the following sections, we will analyze the equilibrium state from a geometric point of view.

III. A GEOMETRIC REPRESENTATION OF EQUILIBRIUM STATES
In this section, we consider a geometric interpretation of the variational form, Eq. (25). As a result, we reveal the existence and uniqueness condition for the equilibrium state.
The geometry we use here is Hessian geometry [38], which is based on a pair of linearly dual spaces. These spaces are endowed with a second dual structure resulting from Legendre transformation with a given convex function. The two dualities yield a generalized orthogonality relation between affine subspaces in the two spaces. Also, the convex function induces the Bregman divergence, which works as an asymmetric distance on the dual spaces.
As we demonstrate, Hessian geometry quite naturally captures the duality between chemical densities and chemical potentials linked by the thermodynamic convex function, and disentangle the algebraic constraints imposed by the stoichiometry of CRNs.

A. Preparation for geometry
We write x ∈ X = R NX >0 for the density space of the confined chemicals, where N X is the number of species. Also, we define its dual space: y ∈ Y = R NX , which is the corresponding chemical potential space. Consider a map from X to Y by using the convex function ϕ (x) := ϕ[T ,μ; x] as where, to focus on x, we omit the argumentsT andμ in ϕ, and the convexity of ϕ (x) is guaranteed by the definitions of thermodynamic potentials, Eqs. (19), (20) and (21). In physical interpretation, the map ∂ϕ gives the value of the chemical potential of a state x. Since we have assumed strict concavity for σ, which implies strict convexity of ϕ (x), the map ∂ϕ is injective. Furthermore, in the ordinary setting of chemical reaction systems, the range of ∂ϕ is R NX ; thus ∂ϕ is bijective (see also [42]).
To construct the inverse map of ∂ϕ, we define the strictly convex function ϕ * (y) on the dual space Y by the Legendre transformation: Employing ϕ * (y), we can represent the inverse map as The diagrammatic summary of these spaces and maps is shown in FIG. 2. With the above setup, we analyze the equilibrium state given by Eq. (25). The critical equation of the variational form, Eq. (25), is represented as where we define the affinity A(T ,μ; x) [44]. This measures how far a state x is from the equilibrium state [14,15]. The solutions of Eq. (29) with respect to ξ give candidates of the equilibrium extent of reaction, ξ EQ .
Since it is difficult to directly analyze Eq. (29), we introduce its geometric representation. Define the following two submanifolds (subsets) in the density space X (see FIG. 2). One is the equilibrium manifold: which represents a set of candidates of the equilibrium state. The other is the stoichiometric manifold: which describes an affine subspace in X and expresses the domain in which the system can evolve by the reaction dynamics with an initial condition x 0 [18][19][20]27]. The important points here are that the equilibrium manifold V X EQ (T ,μ) is determined by the reservoir condition (T ,μ), whereas the stoichiometric manifold P X (x 0 ) is given by an initial condition x 0 .
By using these two submanifolds, we can identify candidates of the equilibrium state x EQ = x 0 + Sξ EQ with the intersection between them (see FIG. 2): If the intersection consists of precisely one point, the equilibrium state is uniquely determined by the variational form, Eq. (25), under a given initial condition x 0 The left and right spaces represent the density space X and the chemical potential space Y, which are mapped each other by ∂ϕ and ∂ϕ * . The red manifold represents the equilibrium manifold, which is curved in X and is flat in Y. By contrast, the blue manifold denotes the stoichiometric manifold, which is flat in X and is curved in Y. The intersection between these two submanifolds gives the equilibrium state.
[46]. By contrast, if the intersection is empty, the equilibrium state does not exist. In the next subsection, we will derive the existence condition for the equilibrium state and prove its uniqueness.

B. Existence and uniqueness condition for the equilibrium state
Let us begin with the derivation of the existence condition, which is composed of two steps: (1) finding the condition for V X EQ (T ,μ) = ∅ and (2) proving V X EQ (T ,μ) = ∅ ⇒ V X EQ (T ,μ)∩P X (x 0 ) = ∅. To obtain the condition for V X EQ (T ,μ) = ∅, we introduce the equilibrium manifold in the chemical potential space Y by using the map ∂ϕ:   [14,15]; Eq. (34) argues that the system does not feel external chemical gradients, and therefore the existence of the equilibrium state is expected. If and only if the condition, Eq. (34), is satisfied, we obtain V X EQ (T ,μ) = ∅ because the inverse map ∂ϕ * exists and V X EQ (T ,μ) = ∂ϕ * (V Y EQ ). We proceed to the second step: V X EQ (T ,μ) = ∅ ⇒ V X EQ (T ,μ) ∩ P X (x 0 ) = ∅. If the consistency condition, Eq. (34), holds, the simultaneous equations y i S i r + µ m O m r = 0 have solutions. Denoting a particular solution by y P = y P i , we getμ m O m r = −y P i S i r . The substitution of it into Eq. (25) leads to By using the argument change x = x 0 + Sξ, we obtain Here, we note that the function maximized in Eq. (36) is bounded above on the density space X [45], which means that the function is also bounded above on the affine subspace P X (x 0 ). Thus, the equilibrium state x EQ must exist, that is, V X EQ (T ,μ) ∩ P X (x 0 ) = ∅. Combining the above two steps, we obtain the following theorem: Theorem 1 Equilibrium states exist, if and only if the consistency condition Eq. (34) is satisfied. In that case, the intersection between the equilibrium and stoichiometric manifolds (Eq. (32) or Eq. (25)) is not empty.
An analogous theorem was originally stated by Wegscheider [24,35] and has been recently reported in Refs. [14,15], under the ideal gas assumption and mass action kinetics. Therefore, our theorem is a generalization of their statement because we use neither ideal gas nor mass action kinetics assumptions. Also, if the stoichiometric matrices S and O satisfy OV = 0 (i.e., Im[O T ] ⊂ Im[S T ]), the system is a so-called unconditionally equilibrium system [14,15]. This means that, for any choice of reservoir condition (T ,μ), the system must converge to an equilibrium state.
Next, we show uniqueness of the equilibrium state under a given initial condition x 0 . Since the function maximized in Eq. (36) is strictly concave and also bounded above on the affine subspace P X (x 0 ), the point x EQ is uniquely determined by Eq. (36). Hence, we obtain the following theorem: Theorem 2 If the consistency condition, Eq. (34), is satisfied, under a given initial condition x 0 , the system converges to the equilibrium state that is uniquely determined by Eq. (36), that is, the intersection, Eq. (32), consists of precisely one point.
This theorem is a generalization of the Horn-Jackson theory for detailed-balanced CRNs [13][14][15]27], which was more recently rephrased as Birch's theorem in the language of algebraic geometry [18][19][20]. As will be shown in Sec. V, if we assume ideal gas conditions, this theorem reduces to Birch's theorem.
In the derivation of the theorems, one may be concerned with the arbitrariness in choosing a particular solution. However, even if we choose another particular solution, the derived equilibrium state x EQ is unchanged because a particular solution is just a reference point for V Y EQ (μ). That is, the choice of a particular solution amounts to fixing a "gauge" in the theory.
Finally, we comment on the equilibrium state in the chemical potential space Y. On the one hand, by denoting a particular solution of the simultaneous equations y i S i r +μ m O m r = 0 by y P , the general solution can be represented as where η = {η l } represents coordinates on Ker[S T ]; Here, U is a basis matrix: U l i := U 1 , U 2 , ... T whose rows U l form a basis of Ker S T (i.e. U S = 0). Thus, we get a parameter representation of the equilibrium manifold On the other hand, by using ∂ϕ, we can map the stoichiometric manifold P X (x 0 ) into Y: where y 0 = y 0 i = ∂ϕ (x 0 ) represent the chemical potential at the initial state x 0 [47]. Here, we note that this manifold no longer describes an affine subspace in Y, but a curved one in general (see FIG. 2). By employing the above two submanifolds, V Y EQ (μ) and P Y y 0 , the equilibrium state can be characterized in Y as

IV. THE SECOND LAW AS MINIMIZATION OF DIVERGENCE
If the consistency condition, Eq. (34), is satisfied, the time evolution of the total entropy, Eq. (23), on the stoichiometric manifold P X (x 0 ) can be written as follows. By usingμ m O m r = −y P i S i r , we get where y P ∈ V Y EQ (μ) and x (t) ∈ P X (x 0 ). We also note that the form of Eq. (41) does not depend on the choice of a particular solution y P (see details in [48]). In this section, we give a geometric representation of Eq. (41) through the Bregman divergence [38,49,50]. Moreover, we reformulate the second law from the viewpoint of Hessian geometry. As a result, we obtain four distinct characterizations of the equilibrium state; one of them is equivalent to Eq. (36).
A. Entropy production during a relaxation to the equilibrium state The Bregman divergence on X is defined by (42) It measures the deviation at a point x between the convex function ϕ (x) and the hyperplane tangent to it at a point x ′ (see FIG. 3). This divergence has the following The line denotes the hyperplane tangent to ϕ(·) at the point x ′ . The divergence is given by the deviation between them at the point x, which is shown in red.
property: D X [x||x ′ ] ≥ 0, and equality holds if and only if x = x ′ , i.e., it acts as an asymmetric distance from x ′ to x. By employing the divergence, we can calculate the production (increment) of the total entropy, Eq. (41), during a time interval [t ′ , t] as where x P := ∂ϕ * y P ∈ V X EQ (T ,μ). Here, we also used the fact that both x(t) and x(t ′ ) are on the stoichiometric manifold P X (x 0 ) to cancel out the term depending on the initial condition x 0 in Eq. (41). Using this representation, we can evaluate the heat dissipation, Eq. (24), involving the divergence as For a relaxation to the equilibrium state, the production of the total entropy, Eq. (43), is computed as where we choose the equilibrium state x EQ as a particular state x P ∈ V X EQ (T ,μ) in the second equality. Thus, the heat dissipation during the relaxation can be represented as The above result can be summarized as follows: Theorem 3 If a CRN relaxes to the equilibrium state (i.e., the consistency condition, Eq. (34), is satisfied), then the total entropy production during a relaxation from an initial state x 0 to the corresponding equilibrium state x EQ can be evaluated by the Bregman divergence given by Eq. (45). Furthermore, the heat dissipation during the relaxation is calculated by Eq. (46).
This theorem represents a generalization of the result by Rao and Esposito [15], which was also reported in the context of mass action systems in Refs. [27][28][29]. As shown in Sec. V, if we assume ideal gas conditions, the Bregman divergence reduces to the generalized Kullback-Leibler divergence, and our statement corresponds to their result.

B. Characterizations of the equilibrium state
Next, we characterize the equilibrium state by four distinct variational forms based on the divergence. For any three points, x, x ′ and x ′′ in X , the following equality holds: Since we have assumed that the consistency condition, Eq. (34), holds, the equilibrium manifold is not empty, V X EQ (T ,μ) = ∅, and the unique equilibrium state x EQ ∈ V X EQ (T ,μ) ∩ P X (x 0 ) exists. If we choose x ∈ P X (x 0 ) , x ′ = x EQ and x ′′ = x P ∈ V X EQ (T ,μ), the last term in the right hand side of Eq. (47) vanishes, because where we use the facts that x−x EQ ∈ Im [S] , y EQ −y P ∈ Ker S T and Im [S] ⊥ Ker S T . This represents the orthogonality between P X (x 0 ) and V X EQ (T ,μ) at x EQ [38,49]. Thus, we get the generalized Pythagorean theorem (see FIG. 4): From this equality, we can derive two distinct variational forms to characterize the equilibrium state x EQ . First, we minimize Eq. (49) with respect to x in the stoichiometric manifold P X (x 0 ). Then, we obtain where we use x EQ = arg min x∈P X (x0) D X [x||x EQ ]. Taking Eqs. (41) and (43) into account, we find that this variational form coincides with Eq. (36); that is, Eq. (50) implies the conventional characterization of the equilibrium state by the second law. Second, if we minimize Eq. (49) with respect to x P in the equilibrium manifold V X EQ (T ,μ) and set x = x 0 , we get another non-trivial variational form: where we use x EQ = arg min xP ∈V X EQ (T ,μ) D X [x EQ ||x P ]. In addition, from Eq. (45), the total entropy production for a relaxation can be evaluated as (52) Owing to Eq. (52), we can evaluate the heat dissipation during the relaxation by using Eq. (46).
The above framework constructed in the density space X can be mapped to the chemical potential space Y. We define the Bregman divergence on Y by holds, we get the Pythagorean theorem in Y as where y P ∈ V Y EQ (μ), y ∈ P Y y 0 and y EQ = ∂ϕ (x EQ ); y 0 = ∂ϕ (x 0 ). By employing the same discussion as for the density space X , the equality, Eq. (54), yields the other two variational forms in Y to characterize the equilibrium state y EQ . One is given by the minimization of Eq. (54) with respect to y in P Y y 0 : which corresponds to Eq. (50). The other is obtained by the minimization of Eq. (54) with respect to y P in V Y EQ (μ): which correspond to Eqs. (51) and (52). Of course, from these variational forms, Eqs. (56) and (57), we can evaluate the heat dissipation during the relaxation as where we write all arguments of ϕ * (y) as ϕ * [T ,μ, y].
The above four characterizations of the equilibrium state are the main results of this work, which is summarized as follows: Theorem 4 Consider a CRN such that the stoichiometric matrices S and O satisfy the consistency condition, Eq. (34) (i.e. the CRN relaxes to the equilibrium state). Define the Bregman divergences in the density space X and the chemical potential spaces Y by Eqs. (42) and (53), respectively; the convex function ϕ(x) represents the partial grand potential density given by Eq. (21), and ϕ * (y) is its Legendre dual function as in Eq. (27). Then, in the density space X , the equilibrium state x EQ for a given initial state x 0 is characterized by the two distinct variational forms, Eqs. (50) and (51). Also, the total entropy production during a relaxation to x EQ is evaluated by Eq. (52). Furthermore, in the chemical potential space Y, the equilibrium state y EQ for a given initial state y 0 is determined by the other two distinct variational forms, Eqs. (55) and (56); the total entropy production during a relaxation to y EQ is computed by Eq. (57).
In particular, the variational forms, Eqs. (56) and (57), lead us to the following simple prescription to identify the equilibrium state: (5) If one wants to know the equilibrium density of confined chemicals, x EQ , it is given by using the inverse map ∂ϕ * . Also, the heat dissipation during the relaxation is computed by Eqs. (57) and (58). 5. The initial state in Y is denoted by y 0 . The total entropy production is given by the minimized divergence from y 0 to the equilibrium manifold V Y EQ (μ), which is an affine subspace in Y. The orthogonal projection of y 0 to V Y EQ (μ) represents the equilibrium state y EQ .

V. CONNECTION TO PREVIOUS WORK
In the preceding sections, we have not imposed detailed functional forms on the thermodynamic potential or the flux density. We have only assumed for them that the potential is a convex (or a concave) function and the flux density satisfies the second law, which guarantee the increasing property of the total entropy function Σ tot (ξ (t)). In this section, we take the ideal gas potential and the mass action kinetics as specific forms of the thermodynamic potential ϕ[T ,μ; x] and the reaction flux density j (t). As a result, a connection to previous work is clarified.
Readers, who are familiar with kinetic modeling of CRNs, can refer to our accompanying paper [39]. There, we derive the results of this section starting from the mass action kinetics and detailed balancing.

A. Ideal gas
In this subsection, under the ideal gas assumption, we demonstrate the geometric structure of thermodynamics constructed in the preceding sections. The form of the Helmholtz free-energy density for the ideal gas is known as (59) where R represents the gas constant;  21), we can calculate the partial grand potential density as (61) This argues that, under a given constantμ, the density of the open chemicals is kept to be constant in the reaction dynamics for the ideal gas cases. This is a natural consequence, because the ideal gas does not have any interactions among chemicals. Furthermore, if the reservoir also consists of ideal gas, its chemical potentials can be represented as which gives the chemical potential for the confined chemicals at a state x. The dual convex function ϕ * (y) on Y is calculated by the Legendre transformation, Eq. (27), as Therefore, the map ∂ϕ * from Y to X is which is the inverse map of ∂ϕ.
If the consistency condition, Eq. (34), holds, by using the inverse map, Eq. (65), we get a parameter representation of the equilibrium manifold in X as This theorem is known as Birch's theorem [18][19][20], which is employed not only for chemical reaction systems but also for the maximum likelihood estimation in statistics [52].
Next, we calculate the Bregman divergences on X and Y, Eqs. (42) and (53): where we use Eqs. (60), (63), (64) and (65). Note that, for the ideal gas cases, the Bregman divergences on X reduces to the generalized Kullback-Leibler divergence [13,15,27]. Therefore, if the consistency condition, Eq.  The substitution of y EQ into Eq. (65) leads to the equilibrium density x EQ . Also, the heat dissipation during a relaxation is computed by Eq. (58).

B. Mass action kinetics
In this subsection, by using mass action kinetics as a specific form of the reaction flux density j (t), we discuss, in terms of the kinetics, the chemical reaction systems composed of ideal gas. As a result, we obtain the local detailed balance condition and find that the entropy production can be represented by the flux density.
For modeling reaction flux densities in ideal gas chemical reaction systems, we here employ mass action kinetics [2][3][4], which is defined as follows: Consider a set of chemical equations, the rth reaction of which is represented as i r state, which is determined by the intersection of equilibrium and stoichiometric manifolds. Also, the entropy production during a relaxation to the equilibrium state is evaluated by the Bregman divergence. Furthermore, the equilibrium state is characterized by four distinct minimization problems of the divergence, two of which are in the density space and the other two are in the chemical potential space. For the ideal gas cases, we have confirmed that our Theorem 2 reduces to Birch's theorem, and the entropy production represented by the divergence coincides with the generalized Kullback-Leibler divergence; the additional assumption of the mass action kinetics leads to the local detailed balance condition.
Although we have only treated the isochoric ideal gas cases in Sec. V, the application is straightforward to conventional CRNs appearing in isobaric ideal-dilutesolution situations in chemistry and biology. To move from an isochoric to an isobaric situation, we replace the Helmholtz free energy with the Gibbs one. At this step, one may be concerned that the volume of the system can change under a constant pressure. However, we can effectively identify the isobaric situation with the isochoric one, because the solvent dominates the volume and the amount of solvent is constant in the reaction dynamics. Thus, the Gibbs free energy is obtained just by modifying the standard chemical potentials in Eq. (59) [1,15]. This direct correspondence between Helmholtz and Gibbs free energies originates from the fact that we can regard the solvent as the background of the reaction dynamics.
However, there exist situations, such as cellular growth, which do not have the simple correspondence between isochoric and isobaric free energies. In this situation, the volume Ω is no longer constant with time, and thus the system may not have any conserved quantities. Due to that, the homogeneity of the entropy function (see Eq. (1)) gives a non-trivial impact to the structure of our theory, and a further extension is required [53].
Much work has been devoted to interpret thermodynamics with geometric frameworks [54][55][56][57][58]. They revealed the geometric dual structure by Legendre transformations in thermodynamics. However, if non-trivial constraints such as stoichiometric ones enter the problem, the constraints introduce important submanifolds (equilibrium and stoichiometric manifolds) into the Legendre dual spaces. We have clarified how the resulting Hessian geometric structure enables us to handle the complex constraints in CRNs. We have also demonstrated that the characteristic thermodynamic properties obtained for mass action systems with the local detailed balance condition emerge from this fundamental structure without assuming any of them. Not limited to CRNs, such geometric structure with the submanifolds can appear in a wide variety of systems with complex constraints, which implies general applicability of our theory.
In this work, we have only dealt with the cases that the system converges to the equilibrium state, that is, the reservoir satisfies the condition, Eq. (34). Otherwise, the equilibrium state does not exist, and the total entropy keeps increasing and finally diverging. Even in such cases, it is known that the system may converge to a certain stable state in a time evolution, which is called the nonequilibrium steady state (NESS) [13][14][15][16][17]. A typical example of NESS in CRNs is the complex-balanced state [13,15,18,27]. However, we can not characterize NESS solely by the entropy function, because the variational form based on the entropy maximization as in Eq. (25) can no longer be employed. The extension of our geometric structure to the cases of NESS is future work [59].