Gradient and GENERIC Systems in the Space of Fluxes, Applied to Reacting Particle Systems.

In a previous work we devised a framework to derive generalised gradient systems for an evolution equation from the large deviations of an underlying microscopic system, in the spirit of the Onsager-Machlup relations. Of particular interest is the case where the microscopic system consists of random particles, and the macroscopic quantity is the empirical measure or concentration. In this work we take the particle flux as the macroscopic quantity, which is related to the concentration via a continuity equation. By a similar argument the large deviations can induce a generalised gradient or GENERIC system in the space of fluxes. In a general setting we study how flux gradient or GENERIC systems are related to gradient systems of concentrations. This shows that many gradient or GENERIC systems arise from an underlying gradient or GENERIC system where fluxes rather than densities are being driven by (free) energies. The arguments are explained by the example of reacting particle systems, which is later expanded to include spatial diffusion as well.

1. Introduction

By the Boltzmann–Einstein relation, the free energy of a system is inherently related to the fluctuations of an underlying microscopic particle system. In this sense, two systems with the same macroscopic behaviour can be driven by completely different free energies, if their corresponding microscopic systems are different. Therefore, one of the main objectives of (equilibrium) statistical mechanics is to derive the “physically correct” macroscopic free energy from fluctuations in microscopic systems. A similar principle can be applied to systems that evolve over time, where dynamic fluctuations may lead to a gradient flow, driven by the free energy. For stochastically reversible systems and close to equilibrium this is the classic Onsager–Machlup theory [1,2]. Such relations are known to hold for many reversible dynamics, not necessarily close to equilibrium [3]. More recently, it was shown that microscopic reversibility always implies the emergence of a macroscopic gradient flow [4], but in this generality one needs to allow for so-called generalised gradient flows. In brief, a generalised gradient structure (GGS) is defined by a possibly non-linear relation between velocities and affinities. Although there exist non-reversible models that lead to macroscopic gradient flows [5], these are considered non-typical; so in order to understand systems with non-reversible microscopic fluctuations, one needs to look for even further (thermodynamically consistent) generalisations of a gradient flow.

One such generalisation is a class of equations called GENERIC [6]. These equations can be seen as a coupling between a gradient flow of some non-increasing free energy (or non-decreasing entropy) and a Hamiltonian system of some conserved energy. One assumes that the two structures are in a sense orthogonal to each other, which then guarantees that also for the coupled evolution, the free energy is non-decreasing and the Hamiltonian energy is conserved. We explain this concept in more details in Section 3.3 and Section 4.2.

The first GENERIC structure that was derived from dynamical large deviations can be found in [7]. In order to pursue a similar procedure in a general setting, one again needs to allow for non-linear relations, just like generalised gradient flows. One is thus lead to study generalised GENERIC structures (GGEN) [8]. In a recent work, necessary and sufficient conditions were found under which microscopic fluctuations induce such a GGEN structure [9,10].

Naturally, many systems do not satisfy those conditions nor do they have a GGEN structure. Therefore there is a need for even more general structures that can still be given a meaningful thermodynamic interpretation. At this point we place ourselves in the context of Macroscopic Fluctuation Theory (MFT) [11]. Central to MFT is the idea that more thermodynamic properties of non-equilibrium fluctuations can be derived if, in addition to macroscopic state variables, the corresponding fluxes are taken into account. In general, fluxes hold more information than state variables due to the possible occurrence of “divergence-free” fluxes that do not alter the states.

One way in which this extra information can be exploited is to extract a generalisation of a GGS, where the affinity/driving force may no longer be the gradient of some free energy [12,13,14,15].

Another way to exploit the flux fluctuations, which is pursued in this paper, is to extract GGS/GGEN structures in a larger “flux space”. The heuristics behind this is that a GGS/GGEN, even in flux space, can be interpreted as a free energy balance. If there would be work done that results in a divergence-free flux, then one might expect an gap in the energy balance, so that such systems can not induce a GGS or GGEN. Could it be possible for such a system to have a GGS/GGEN structure in the flux space, when it fails to have a GGS/GGEN structure in the state space? The main point of this paper is that this is generally impossible. We will see in Theorem 3 that if the fluctuations induce a GGS or GGEN in the space of fluxes, this is, up to physical conditions, equivalent to the fluctuations inducing a GGS or GGEN in the state space.

Our leading example and main application will be that of a (non-spatial) isothermal chemical reaction network, as studied in [16]. In Section 2 we recall the main arguments from that paper, applied to concentrations undergoing reactions. In Section 3 we expand these ideas and show how they can be applied to reaction fluxes. It turns out that reactions that occur on a faster time scale may give rise to a GGEN Hamiltonian part. Based on this example, we then develop an abstract theory about induced GGSs and GGENs in flux space, in Section 4. To show the generality of these principles, we then show in Section 5 how the theory applies to transport fluxes in a diffusing particles model, and in Section 6 we combine the arguments from Section 3 and Section 6 to derive results for transport and reaction fluxes in a simple reaction-diffusion model.

2. Leading Example 1: Chemical Reactions

In this section we explain the main concepts, mostly by reiterating the arguments of [4] and [16]. In particular, we argue that large deviations/fluctuations provide a “physically correct” GGS for the evolution of concentrations undergoing chemical reactions.

Consider a network of isothermal chemical reactions, for example:

We denote the set of species by (in this example ) and the set of reactions by ; for each reaction r we consider a forward and backward reaction (so that here consists of one element). The stoichiometric coefficients are denoted by for the forward reactants (here ) and for the forward products (here ), which yields the state change matrix . The evolution of the concentrations is then described by the Reaction Rate Equation, with concentration-dependent reaction rates .

Typically (and often in this paper), the reaction rates will be of the form for some constants , using the notation . In that case the network is said to be of mass-action kinetics. For more background on chemical reaction networks we refer to the survey [17]. More details about the fluctuations and induced GGSs for chemical reaction networks can be found in [16]; for completeness we shall recall these results in this section.

To notationally stress the similarity and differences between different concepts throughout this paper, we shall always denote net quantities by a bar (), and we distinguish functionals on state space from functionals on flux space by a hat (). In this section we study concentrations only, which we consider to be states.

2.1. Reacting Particle System

A classical microscopic particle system underlying the evolution (1) is the following [18]. Let V be a large, well-mixed volume that contains, at time t, a total number of particles of species , . A reaction r occurs randomly with some propensities (jump rates) . Whenever a forward reaction r occurs, particles are removed and particles are created, and vice versa for a backward reaction. Hence each reaction requires a cumbersome relabelling of particles . It is therefore more practical to work directly with the (particles per volume) concentration . This quantity will also play the role of the macroscopic state variable. Whenever a forward or backward reaction r takes place, the concentration can now be simply updated by a jump . Then the -dimensional vector is a Markov jump process, which satisfies the master equation

It will be beneficial to work with the corresponding generator, which is the adjoint of the right-hand side of the master equation, with respect to the dual pairing with an arbitrary test function:

The propensities that are usually used in the so-called chemical master Equation [17] are derived from combinatoric considerations, and yield the mass-action kinetics in the limit [16,18]: using the notation and .

2.2. Equilibrium: Limit, Large Deviations and Free Energy

We will from now on (throughout this section) assume that the reaction network is of mass-action kinetics (2), and chemically detailed balanced, i.e., there exists a for which

Naturally, is an equilibrium under the deterministic evolution (1). It should be stressed that, given an initial concentration , both the deterministic evolution and the stochastic model is confined to the “stoichiometric compatibility class” . Therefore, it is not clear whether this equilibrium lies within the compatibility class that corresponds to the initial concentration. However, if there exists a detailed balanced concentration , then there exists a unique detailed balanced concentration within each such class, see [17] and the references therein. Without loss of generality, we can therefore implicitly assume that the detailed balanced equilibrium is unique, and lies within the correct compatibility class.

Under Assumption (5), the invariant distribution of the stochastic model is known to be [17]:

Letting , this invariant distribution concentrates on the equilibrium state:

One can then extract the free energy by considering the corresponding large deviations, i.e., the exponential rate which with converges to zero. Indeed, by Stirling’s formula,

Such limit is known as a large-deviation principle; the function on the right characterises the stochastic cost of microscopic fluctuations. If the reaction rates are related to an internal energy via Arrhenius’ law, then the expression is really the Helmholtz free energy, apart from a normalisation term and a constant scaling, as explained in more detail in (Section 2.2 & 2.3, [16]).

2.3. Dynamics: Limit and Large Deviations

Observe that in the microscopic model, the process speeds up as V increases with order V while the jump sizes are of size . Therefore by (4), as the generator converges to the limit generator,

Since this generator depends on the test function f through only, we can make the ansatz that the limit process is deterministic , for some curve . Plugging this into the definition of the generator yields:

As this relation holds for any test function f, we see that the ansatz was justified if the postulated curve satisfies the Reaction Rate Equation (1). Hence the stochastic process converges (pathwise in probability) to the deterministic solution of the Reaction Rate Equation.

Similar to the calculation of the fluctuations of the equilibrium (6), we can study the large deviations of the path probabilities ; this is known as a dynamic large-deviation principle. These dynamical fluctuations can be formally calculated with the framework of [19]. To this aim we study the non-linear generator:

As before, the limit depends on the test function through only, which is consistent with the fact that the limit is deterministic. We then define, by a slight abuse of notation,

The dynamic large-deviation principle now states that The rigorous definition of the large-deviation principle, the heuristics behind this method, and the rigorous proof of this statement is all beyond the scope of this paper. For the precise details we refer to [19], and for the rigorous proof for this particular system (by more classical methods) to [20,21,22]. For the sake of brevity, we assume that the randomness in the initial condition is sufficiently small (e.g., deterministic) so that we do not obtain initial fluctuations.

The function where

2.4. GGS, Energy Balance, and Relation with Fluctuations

The (naive) aim is to rewrite the macroscopic equation as a gradient flow of some free energy : where is some linear symmetric, positive definite (linear response) operator that maps thermodynamic forces to velocities. Mathematically, this operator can be interpreted as the inverse of the metric tensor of some manifold, so that the right-hand side is the gradient on this manifold.

Clearly, (10) is equivalent to requiring if we set and . Integrated over a time interval , this reads

The last two terms describe the free energy loss (or entropy production), and the first two terms describe the dissipation; as such this equation represents a free energy balance. (For a linear gradient flow (10), ; hence the dissipation can be seen as a kinetic energy. This is however no longer true for general GGSs).

Observe that this expression is always non-negative, and is 0 exactly on the gradient flow (10). Moreover, we see that this expression has the same dimension as , the free energy; it is indeed the free energy cost to deviate from the macroscopic dynamics. This interpretation shows that this cost should be equal to the cost of microscopic fluctuations (9).

This is in many cases, and particularly in this case of chemical reactions, impossible. Since the large-deviation function (8) is non-quadratic, one should allow for non-quadratic dissipation terms. We therefore replace the two squared norms by a pair of dual dissipation potentials: where, as in the quadratic case, the potentials are convex duals of each other, i.e., and . Moreover, we assume that and are both non-negative; from (12) we then see that the free energy is non-increasing on the flow.

The right-hand side of (12) is always non-negative, since by definition of the convex dual. Therefore, the right-hand side of (12) can only be 0 if is minimal in the right-hand side. By differentiation we then get that (12) implies:

We call such equation a generalised gradient flow, and the underlying structure a generalised gradient structure (GGS). Note that the generalisation with respect to (10) lies in the fact that we allow for a non-linear relation between forces and velocities. We moreover say that a gradient structure is induced by a cost function whenever .

In Section 4 we will recall the relation between fluctuation costs and GGSs, as described in [4]. Applied to the current setting of chemical reactions, we reiterate the following result from [16]. If we again assume mass-action kinetics (2) and chemical detailed balance (5), then there exists a unique GGS induced by the large-deviation cost from (8), where with . Since appears as a sum, the expression for becomes a so-called inf-convolution where all reactions are strongly intertwined. For more details on these inf-convolutions and the factor in front of the free energy, we again refer to (Section 3.4, [16]).

3. Leading Example 2: Fast-Slow Reaction Fluxes

As mentioned in the introduction, the motivation behind the current paper is to search for thermodynamically consistent structures for systems that are not detailed balanced. The idea is that we increase the space by taking fluxes into account. However, in order to see the connection between structures in flux space and structures in state space, we dedicate this section to an example that is “almost” detailed balanced. More general systems will be dealt with when discussing the general theory in Section 4.

3.1. Reacting Particle System and Macroscopic Equation

We consider a system of fast and slow chemical reactions, , where the slow reactions are assumed to be of mass-action kinetics and detailed balanced, as in the previous section. By contrast, we will not assume anything of the like for the fast reactions, neither shall we assume that each reaction consist of a forward and a backward reaction. However, in the microscopic model, we shall assume that the fast reactions happen on a faster time-scale, i.e., but their effect on the concentrations is smaller, i.e., whenever a reaction occurs at time t. To notationally distinguish the two time scales we use a tilde () to denote one-way fast quantities. Let us briefly mention that the additional fast reactions are not essential to see the relation between structures in flux and state space; they just lead to a richer example that is interesting in its own right.

The main idea is now to increase the state space by bookkeeping the events in the microscopic system, in this case, by counting the number of reactions that have occurred up to a given time. With the right scaling, this defines the following integrated reaction fluxes: and to shorten notation we sometimes write The term “integrated” signifies that these fluxes are cumulative over time. This simplifies the microscopic analysis since the corresponding process is Markovian; on a macroscopic scale only the time-derivatives will play a role. We also mention that we consider net rather than one-way slow fluxes, else the slow dynamics would not induce any force or gradient structure, see (Section 4.6, [15]).

We shall always assume that the initial condition is known (deterministically) a priori, so that there always holds:

In this sense the integrated fluxes encode more information than the concentrations. The equation above relates the change in concentration to the integrated fluxes, and can thus be interpreted as a continuity equation.

The integrated fluxes are again a Markov process, with generator (cf. Equation (3)):

3.2. Limit and Large Deviations

We now mimic the arguments of Section 2.3, but now in the space of fluxes. Let , and so . Then by the same argument as in Section 2.3, using the scalings (4) and (14), one finds that as , the random process converges (pathwise in probability) to the solution of the macroscopic equations: again using the notation (15). Indeed, combining these equations leads to the macroscopic Equation for the concentrations.

Similarly, we study the fluctuations through the non-linear generator:

As in Section 2.3, this limit depends on the gradient only, which is consistent with the deterministic limit (16). Again we set, by a slight abuse of notation, using the notation if and ∞ otherwise. Then, the dynamic large-deviation principle for the integrated fluxes states that (see [22] for a rigorous proof):

Comparing this dynamic large-deviations principle with (9), we see—not coincidentally—strong similarities. Let us assume that the limit fast fluxes do not influence the concentration, i.e.,

Naturally, this macroscopic condition entails a sort of decoupling between the slow and fast dynamics, e.g., when the species involved in the slow dynamics act as a catalyst for the fast dynamics. In that case is exactly the process with generator (3). The contraction principle of large deviations theory (Theorem 4.2.1, [23]) then states that the two large-deviation costs are related via:

In this setting, this infimum is only with respect to the second variable, because depends on w through only, which in turn arises naturally from the fact that the jump rates depend on the state and not on the integrated flux. Therefore the relation above simplifies further to: which is consistent with Remark 1. These relations will be the starting point of the general theory that we develop in Section 4. Let us only mention here that as a consequence we also have the relation for all , cf. (7) and (18).

3.3. Induced GENERIC Structure

We now investigate whether the large deviations (19) induces some structure in the space of fluxes. It turns out that this is indeed the case. As was found in [12,15] and further studied in [13], under the detailed balance assumption (5) the slow fluxes induce the GGS , where with and . Extending these dissipation potentials to the full flux space by setting and , we can decompose the large-deviation cost function as and accordingly, the macroscopic evolution (16) as:

Due to (22) and the decoupling condition (20), the fast fluxes are orthogonal to the driving force, i.e.,

The quadruple satisfying this condition falls within the class of what is recently coined pre-(Generalised) GENERIC (pGGEN) [9,10].

It was shown in those works that non-interaction condition (25) is a necessary and sufficient condition for the existence of an underlying Generalised GENERIC (GGEN) structure . (One often needs to introduce an auxiliary energy (e.g., a heat bath) to force conservation of energy, which enlarges the degrees of freedom in the system. To keep notation accessible we ignore this issue). This means that the fast flux term is Hamiltonian for some Poisson structure L satisfying the Jacobi identity (see Section 4), is some Hamiltonian energy, and the following two non-interaction conditions are satisfied: for all , and . These two conditions guarantee that along solutions the free energy is non-increasing and the Hamiltonian energy is conserved:

One main message of this section is that the flux large-deviation cost ℒ induces a unique pGGEN system . Although this implies the existence of a GGEN system induced by ℒ, one can not uniquely decide on the basis of ℒ what the “correct” Hamiltonian structure should be. Additional physical or mathematical arguments needed to uniquely fix the Hamiltonian structure are beyond the scope of this paper; for possible constructions of Poisson structures L and energies , we refer the reader to (Section 4, [10]).

It should be noted that the pGGEN system is rather special in that (25) holds for any , since the drift is “divergence-free”. We will see that such systems play a special role in connecting systems in flux and state space.

Another main message is that the GGS part (22) on flux space is strongly related to the GGS (13) on state space, just like both cost functions are related by (21). These observations will be the basis of the general theory of the next section.

4. General Theory

Following the examples of Section 2 and Section 3, we now develop a more abstract framework to study the relation between energy-driven structures in flux and in state space.

4.1. Geometry and Notation

Throughout this section we assume to be given:

A differentiable manifold (“flux space”), where tangents are denoted by and cotangents by . Note that we distinguish between tangents and cotangents; we write the dual pairing between them as or simply ;

A differentiable manifold (“state space”), where tangents are denoted by and cotangents by ;

A surjective differentiable operator with bounded linear differential and adjoint differential . This defines an abstract continuity equation , or in differentiated form , see Figure 1. Contrary to Section 2 and Section 3, the second continuity equation may now also depend on w. In practice, the continuity mapping depends on the initial state , which we assume to be fixed.

Figure 1

We consider a “flux manifold” and a “state manifold” . The “continuity map” maps points to points on ; its differential maps tangents to tangents .

An L-function on flux space (see below for the definition of L-functions). This function could be a dynamic large-deviation cost function corresponding to random fluxes in some microscopic system, but throughout this section it could also be a more general expression.

An L-function on state space, related to the flux space L-function via

This relation is again inspired by large-deviation theory, where such relation holds due to the contraction principle (Theorem 4.2.1, [23]).

Corresponding to the L-functions are their convex duals with respect to their second variable, i.e., and with

We express assumptions in terms of these duals, since in practice they are often more explicitly given than their corresponding L-functions.

Recall that in the previous section we saw that depends on w through only. This condition becomes slightly more complicated in non-flat spaces, for a number of reasons. Firstly, the flux j can not be kept fixed while changing w, unless one has a path-independent notion of parallel transport, i.e., the space is flat. Secondly, even if would not depend on w, this dependence could re-enter through the continuity equation in the infimum . Therefore, the condition that we need is that for all . It is easily seen that this condition is equivalent to the following flux invariance:

All manifolds and functionals are assumed to be sufficiently differentiable wherever needed. For a (differentiable) functional (and similarly on flux space) we write for the derivative, in the sense that on a curve .

4.2. Definitions

We now define the notions of L-functions, dissipation potentials, GGS, pGGEN, Poisson operators and GGEN on flux space; the same concepts on state space are defined analogously. Naturally all notions are compatible with the exposition from Section 2 and Section 3.

We say

Since ℒ is non-negative, the second condition simply means that ℒ should have a unique minimiser; this minimiser corresponds to an evolution equation of the type . Indeed, since is a minimiser, it will generally satisfy the implicit equation . Due to the convexity, ℒ is also the convex dual of ℋ.

Central to GGS, pGGEN and GGEN is the notion of dissipation potentials:

A function dissipation potential whenever for all

is convex;

;

If these conditions hold, then the same conditions hold for the (pre-)dual dissipation potential

We also say that dissipation potential pair whenever

A generalised gradient system (GGS) is a triple

As explained in Section 3.3, by extending a GGS with an orthogonal drift one arrives at

([10]). A Generalised pre-GENERIC system (pGGEN) is a quadruple

We say that an L-function induces a pGGEN

Finally, if the drift has the form of an Hamiltonian system that behaves more or less independently of the GGS part we arrive at a Generalised GENERIC system. In order to define this we first define:

A linear operator for all

Jacobi’s identity implies skew symmetry, i.e., for all , ; in particular one has . Finally,

(Section 2.5, [8]). A generalised GENERIC system (GGEN) is a quintuple

We say that an L-function induces a GGEN

If a GGS, pGGEN or GGEN is given on a manifold

4.3. From L-Functions to GGS, pGGEN and GGEN

We now recall some of the main results from [4,10], that give necessary and sufficient conditions for an L-function to induce a GGS or pGGEN. A similar result for GGEN does not exist, since an L-function does not uniquely determine a Poisson operator and Hamiltonian energy. However, from a pGGEN one can always construct a GGEN (in a non-unique way); for that result we refer the reader to (Section 4, [10]).

Again, the following results are described, but not restricted to flux space.

(Proposition 2.1 & Theorem 2.1, [4]). Let

ℒ induces a GGS for some dissipation potential Ψ,

for some dissipation potential ,

(34)

(35)

In that case

From condition (35) we see that , if it exists, is uniquely given up to constants. Note in particular that conditions (34) and (35) do not involve .

(Theorem 3.6, [10]). Let

ℒ induces a pGGEN for some dissipation potential Ψ,

(37)

,

(38)

In that case

From condition (37) we see that ℋ must consist of a convex part and a linear part , so that the drift b is a priori and uniquely fixed by ℋ. Therefore is again uniquely fixed by condition (38). This is different from the GENERIC setting; uniquely defines and vice versa, but the whole quintuple may not be unique. However, one can still state a GENERIC analogue of Theorems 1 and 2 as follows:

Let

ℒ induces a GGEN

and          (41)

In that case

By Lemma 1 (see below), we can apply Theorem 1 to the shifted L-function and back. This yields the three equivalences, apart from the other non-interaction condition (31). From the explicit formula (43) one finds that the missing non-interaction condition is equivalent to (40) and to (42). ☐

The previous proof made use of the following lemma:

Let

Because of the non-interaction condition (32), the shift transforms relation (33) into (30), and analogously for the other direction. ☐

4.4. Relation between Structures in Flux and State Space

We now consider L-functions ℒ and on flux and state space, and study how their induced structures are related.

Assume that an L-function for some (27) induces a GGS

Since and , we can rewrite and because is a dissipation potential pair, clearly

This is equivalent to , which by Theorem 1 implies that induces a GGS for some . ☐

In the above proposition, also depends on w through only, which is a very physical assumption. It does imply however, that the equilibria of the flux gradient system can only be unique up to the kernel of ; this kernel can be interpreted as a generalisation of divergence-free vector fields.

A natural question is now whether we can turn the statement of Proposition 2 around. Indeed, if the invariance condition (28) holds and we restrict to pGGEN with “divergence-free drifts”, then the statement becomes an equivalence, and we have an explicit relation between the flux and state dissipation potentials. This is a stronger version of the statement in (Proposition 4.7, [15]), where we related GGS to so-called “force structures”.

Assume that an L-function

ℒ induces a pGGEN

If these statement hold, then the dissipation potentials

Assume that ℒ induces a pGGEN with and . Since by assumption does not depend on , by (36) the expression is also invariant under this choice. Therefore we can define by (44); it is easily checked that its convex dual is given by (44). We can write: and hence induces the GGS , which is unique by Theorem (1).

For the other direction, assume that induces a GGS . Define . Then by the invariance condition and by (34):

Now define by (39). In particular since and, by the definition of L-functions, . Then (37) holds and hence by Theorem 2 the flux L-function ℒ induces the pGGEN . ☐

Due to the non-uniqueness of induced GGEN systems, there is no similar “if and only if” statement for the GENERIC setting. Nevertheless, in one direction, the GGEN analogue of Theorem 3 is:

Assume that an L-function for some

Then the L-function

If in addition, ℋ satisfies the invariance principle

We again apply Lemma 1 to transform the problem into a problem of GGSs. Indeed, the L-function induces the GGS . Hence by Proposition 2, a GGS (for some ) is induced by the L-function

If we can now validate that is a Poisson structure, and that the non-interaction conditions are satisfied for , then Lemma 1 concludes the proof.

For the Poisson structure, note that, for any smooth , the Lie bracket remains unaltered:

The non-interaction condition (32) is clearly satisfied as for any we have

To check the other non-interaction condition (31) we use the equivalent formulation (42). Indeed, for any and ,

Finally, if ℋ satisfies the invariance property (28), then Proposition 2 yields relations (44). ☐

Condition (46) is in a sense a natural one as the following result shows:

Assume that L-functions ℒ and

For any and , we may write , and so by (41),  ☐

5. Diffusion

In this section we apply the ideas of the previous section to a model for diffusion. The flux structure related to diffusion is interesting in its own right, and as far as the author is aware, previously unknown. In the next section we show how this model can be coupled with the results of Section 3 to obtain flux and state GGSs/pGGENs for reaction-diffusion systems.

Typical microscopic models of diffusion consist of Brownian particles, or discretised versions thereof, like random walkers or an exclusion process. Since empirical fluxes are a bit easier to define on a lattice, we focus on independent random walkers (With the scaling that we use, the system of independent random walkers is “exponentially equivalent” to a system of Brownian motions, meaning they share the same hydrodynamic limit and large deviations). The model, its many-particle limit and large deviations are similar to e.g., [24,25].

5.1. Diffusing Particle System

This microscopic particle system actually has two scaling parameters: the number of particles, which we denote by V for consistency with the rest of this paper, and the lattice spacing . The speed with which as is irrelevant. For fixed V, let be independent random walkers on the lattice with jump rate . Define the random concentration and (integrated, net) flux by:

As usual denotes a spatial area, possibly a small box surrounding one lattice site, and and are measures.

Now measures the number of particles present in (lattice points in) an area , while measures the net number of particles that have jumped through all midpoints in , in direction , for , see Figure 2. Note that both and are defined as measures on the lattice with shrinking distance between lattice points; this measure-valued formulation is needed to pass to a continuum limit later on.

Figure 2

When a particle jumps from in the direction with unit vector , this event is recorded in the integrated flux in the midpoint . As such the quantity measures the net number of particles that have passed through the (upper-right) boundary of an imaginary box around midpoint x, drawn in dotted lines.

The concentrations and fluxes are related by the V-dependent continuity equation:

Using that , the integrated flux is a Markov process in with generator

Here, the factor comes from the time scaling , together with the fact that we have independent particles to choose from.

5.2. Limit and Large Deviations

For a test function , we set:

The continuity Equation (47) converges to the limit continuity equation (with the usual divergence operator):

With this notation, as the generator (48) converges to (if ):

As in Section 2.3, the limit generator depends on derivatives of the test function only, and so the process converges (pathwise in probability) to the deterministic path satisfying Fick’s Law:

Naturally, combining this equation with the continuity Equation (49) yields the diffusion equation for the empirical measure:

Similarly, we derive the dynamic large deviations by studying the non-linear generator: which follows from expanding the exponentials (with order- exponents) up to second order. Then the following large-deviation principle on flux space holds: with

Note that is independent of , and by (51), the invariance condition (28) is satisfied. Hence by the contraction principle (Theorem 4.2.1, [23]), one obtains the large-deviation principle corresponding to the states (empirical measures): with using the notation .

5.3. Induced GGSs in Flux and State Space

We can now apply Theorem 1 to extract a GGS from the L-function ℒ. We first choose the ‘naive’ flat manifold of non-negative vector measures (equipped with the flat total variation metric). It is easily checked that condition (34) holds for the free energy given by: where we identify , and we implicitly set whenever the measure is not absolutely continuous. (This expression can again be seen as a relative entropy, cf. (22), but now with respect to the Lebesgue measure, where the measure of the whole space—in this case infinity—is omitted. See also (Proposition 3.2, [4]) for a general result in locally finite measure spaces). The dissipation potentials are obtained from (36) and (29), which yields

Theorem 1 states that ℒ induces the GGS on flux space.

For the state space, it is well-known that the state L-function (52) induces the entropy-Wasserstein gradient flow of the entropy functional [4,26,27,28]. By the theory developed in Section 4, we can now see how this gradient structure is related to the flux gradient structure. Indeed, the flux free energy depends on state only, i.e., , where , and so by Proposition 2 the state L-function induces a GGS driven by . Moreover, since the invariance condition (28) holds, the dissipation potentials are related by Equations (44) (recall the norms introduced above in (52)):

5.4. A New Geometry

The form of the dissipation potential and suggests that it is more natural to use different manifolds in the spirit of Remark 2. For the state space this points to the space of probability measures of finite second moment space, equipped with the Monge–Kantorovich–Wasserstein metric [29]: with tangent and cotangent space and . For this setting, the inverse metric tensor is known by the Benamou–Brenier formula (Theorem 8.1, [29]) to be , so that the GGS is indeed the entropy-Wasserstein gradient flow [30]:

Motivated by this observation we can take for the flux manifold the space of signed vector measures of finite first moment . This choice guarantees that the corresponding states have finite second moment (once ):

Moreover, we can now use the dissipation potential to construct a natural metric on : where the infimum runs over paths of fluxes for which remains non-negative. The corresponding tangent and cotangent spaces (in the interior of the domain) are simply and and the inverse metric tensor is . This yields an interesting geometry in flux space, which, as far as the author is aware, is still unknown in the literature.

6. A Simple Reaction-Diffusion Model

We now combine the models from Section 3 and Section 5 to study reaction-diffusion models in flux and state space. The stochastic particle system will now consist of ‘reacting random walkers’. It is known that, if the reaction networks include reactions of different orders (unimolecular, bimolecular, etc.) and we only allow particles to react if the required number of particles are present within the same site/compartment, then the model may not converge to the expected reaction-diffusion Equation [31]. The reason behind this is that for a multimolecular reaction, it becomes very unlikely that the required amount of reactants are all within one site/compartment; different order reactions would require different scalings. This is beyond the scope of the current paper. However we can already illustrate the combination of reaction and transport fluxes for a simple system of unimolecular equations of the type:

In this section we consider GGSs only, hence we shall always consider net rather than one-way fluxes.

6.1. Reacting and Diffusing Particle System

Since we consider unimolecular reactions only, we can take independent reacting random walkers on the scaled lattice , where each reaction occurs locally at each lattice site with rate or respectively, so that

For the transport mechanism, we assume that the two species hop to neighbouring lattice sites with rates and respectively.

As before we consider the random concentrations, as well as the integrated (net) fluxes, where we now distinguish between transport fluxes and reaction fluxes. If is the position and is the species of the i-th particle, then

The concentrations and fluxes are again related by a continuity equation: where the discrete divergence is as in (47), and is the matrix consisting of one state change vector corresponding to a forward reaction.

The pair is then a Markov process with generator where

6.2. Limit and Large Deviations

By the same procedure as in Section 2.3 and Section 5.2, one finds that as and , the continuity operator (56) converges to (assuming ): and the process converges (pathwise in probability) to the solution of the system:

Indeed, putting these together yields the reaction-diffusion equation for the limit concentrations:

To find the corresponding large deviations, we combine (17) and (50) to calculate the non-linear generator:

Let us again abbreviate and . The limiting non-linear generator can now be split into

Since each mechanism corresponds to a separate flux, the corresponding L-function also splits into two parts: using the usual the relative entropy between two measures, i.e.,: if , else .

As before, the calculation above formally shows that the flux large-deviation principle holds (see for example [32] for a similar but rigorous result):

The L-function (59) splits into two parts because the only interaction between the two mechanisms occurs through the state . By contrast, the corresponding state space large-deviation are much more complicated. Observe that the continuity Equation (57) is an affine function of , and so is independent of , and the invariance condition (28) holds. As explained in the beginning of Section 4, this means that one can apply a straightforward contraction principle on the tangents to yield the large deviation cost function for the states/concentrations:

This infimum reintroduces a strong interrelation between the two driving mechanisms. Indeed, for a given tangent , the fluxes in this infimum correspond to an optimal splitting between the two mechanisms, which can be seen as an inf-convolution. Similar interactions also arise when considering multiple reaction pairs, see (Section 3.4, [16]).

6.3. GGSs in Flux and State Space

We now apply Theorem 1 to the reaction-diffusion setting. The symmetry condition (34) holds for the function (58) if we choose the free energy functional

Naturally, this functional can be seen as a combination of (22) and (53), where just like (53), it has the form of a relative entropy with respect to a locally finite invariant measure, namely . We find the corresponding dissipation potentials from (43), which is again a combination of the non-quadratic potentials (23), (24) and the quadratic potential (54): with . Let the flux space be given by , where the first two spaces, corresponding to the transport fluxes, are equipped with the metric (55) introduced in the previous section. By Theorem 1 the flux cost function ℒ induces the GGS . We stress that the dissipation potential splits into two potentials for the transport and reaction mechanisms respectively, but the free energy is one and the same for both mechanisms.

Since a GGS is a special case of a pGGEN, by Theorem 3, the state cost function also induces a GGS , in this case in the space . The same result yields and, using :

We stress that, analogous to the L-function (60), the dissipation potential on state space no longer splits into two parts.

7. Discussion

7.1. General Theory

We studied gradient and (pre-) GENERIC structures induced by flux large deviations, and the relationship between structures induced by state large deviations. At a first glance, the physical interpretation of the resulting GGS/GGEN structures in flux space is not immediate. However, in practice many induced flux GGS/GGEN structures have a free energy and dissipation potential that only depends on the integrated flux through the state of the system. Hence the main difference with state space GGSs/GGENs is that the fluxes rather than velocities are being driven, which seems a very physical assumption.

The general theory that we developed in Section 4 presumes two given L-functions, in flux and state space, where the second is related to the first through an infimum; this is the typical setting for large deviation-based cost functions. A condition that is central to this theory is that this infimum can be taken over the second (time-derivative) argument only. When the two L-functions are indeed large-deviation costs, then this condition means that the jump rate in flux space depends on the state only, and not on the integrated flux, which is a very natural assumption as well.

As described in the introduction, one motivation behind this research was to study whether large deviations can induce a GGS/GGEN structure in the flux space when it fails to induce such structure in the state space. This has to be answered negatively. It turns out that if the flux large deviations induce a GGS/GGEN, then so do the state large deviations. The same principle also holds in the other direction, but the flow in flux space could have an additional ‘divergence’-free Hamiltonian term that is not observed when considering states only.

7.2. New Structures in Flux Space

Nevertheless, we uncovered a number of previously unknown GGSs/GGENs (as far as the author is aware). The interesting feature of their corresponding dissipation potentials in flux space, is that, due to Equations (44), they have much simpler expresssions than their counterparts in state space.

The first new flux structure that we uncovered is the pGGEN for multiscale chemical reaction networks, where the dissipation potential and free energy is given by (24) and (22). In order for this structure to be a pGGEN, we needed to assume that the fast dynamics do not influence the state, i.e., condition (20). From (20) we see that this condition is sufficient but not necessary, so there is room for possible generalisations.

The second new GGS that we identified is the flux counterpart of the Entropy-Wasserstein gradient flow a described in Section 5.3 and Section 5.4. Much of the (well-developped) Wasserstein calculus is based on transforming to the flux space, performing manipulations and limits there, and then transforming back to the state space, see for example [33,34]. Such arguments may be simplified considerably when working in the flux space directly. It would be interesting to see which geometric properties the new manifold has, e.g., in terms of Ricci curvature bounds and geodesic convexity of the free energy [35].

The third new GGS structure that we derived is a combination of the other two, and models reaction-diffusion equations via their reaction and transport fluxes. Here, we restricted to very simple independent unimolecular reactions, which is certainly generalisable. However, more general reactions would require introducing even more notation as well as more different scaling transitions to circumvent the aforementioned convergence problems. This is beyond the scope of the current paper.

7.3. Possible Lines of Future Research

Apart from the open questions mentioned above, the transition from pGGEN to GGEN is worth studying. Throughout the paper we mostly worked with pGGEN rather than GGEN. The reason is that a given L-function can uniquely induce a pGGEN structure, with many different GGEN structures corresponding to it [9,10]. It is still an open question whether and how a meaningful Hamiltonian energy and Poisson structure L in GGEN can be derived from a given drift term b in pGGEN such that (other than by physical arguments). As found in [10], the Hamiltonian part is often related to dynamics that are deterministic even on the microscopic scale, which also coincides with [7]. However, the setting of multiscale reaction fluxes from Section 3 shows that the Hamiltonian part could come from stochastic dynamics that happen on a faster time scale. Therefore, a possible line of future research would be to consider large deviations on the faster time scale (so-called moderate deviations), and study whether the corresponding variational expression can somehow be used to derived a Hamiltonian structure .

Another line of future research could be to study how structures in flux space can be exploited numerically. As we saw in the examples and Section 6 in particular, if a microscopic system consists of multiple driving mechanisms, then the corresponding Markov generator as well as the non-linear generator is a sum over these mechanisms. By considering separate fluxes for each of these mechanisms, the large-deviation L-function and its induced dissipation potential also splits into different terms for each mechanism. By contrast, when working on the state space this would lead to an inf-convolution (compare for example (61) with (62)). The decomposition that occurs in the space of fluxes can be beneficial for numerical purposes, e.g., using operator splitting techniques [36].