Non-Equilibrium Thermodynamics and Stochastic Dynamics of a Bistable Catalytic Surface Reaction

Catalytic surface reaction networks exhibit nonlinear dissipative phenomena, such as bistability. Macroscopic rate law descriptions predict that the reaction system resides on one of the two steady-state branches of the bistable region for an indefinite period of time. However, the smaller the catalytic surface, the greater the influence of coverage fluctuations, given that their amplitude normally scales as the square root of the system size. Thus, one can observe fluctuation-induced transitions between the steady-states. In this work, a model for the bistable catalytic CO oxidation on small surfaces is studied. After a brief introduction of the average stochastic modelling framework and its corresponding deterministic limit, we discuss the non-equilibrium conditions necessary for bistability. The entropy production rate, an important thermodynamic quantity measuring dissipation in a system, is compared across the two approaches. We conclude that, in our catalytic model, the most favorable non-equilibrium steady state is not necessary the state with the maximum or minimum entropy production rate.


Introduction
The study of dissipative systems has been an active research topic during many years, and a plethora of interesting studies have been published [1]. Nowadays, the notion of dissipative structures appears in a wide range of fields from physics to the economics and social sciences [2]. These dissipative structures, which are characterised by the existence of non-equilibrium steady states, are created by irreversible processes that dissipate energy and generate entropy [1,2]. However, a question of current interest is whether the entropy production rate reaches an extremum (maximum/minimum) or not at the non-equilibrium steady state (NESS) [3][4][5][6][7][8][9][10][11][12][13][14]. In this work, we explore this problem from the viewpoint of bistable catalytic reaction networks.
The idea that universal extremal principles determine many of the phenomena that occur in nature has been around for many years. The point is that, independent of the form of the specific system, among all possible trajectories that may link an initial state from the final one, the trajectory that will be actually followed, under the laws governing the system, extremizes a certain quantity. One of the most famous extremal principle is the principle of least action for conservative systems, which is used to derive many physical theories, like for example Maxwell's equations, Newton's laws, and quantum mechanics [15]. The extend to which this type of principles can be expected in systems operating far from the state of thermodynamic equilibrium (dissipative systems) has attracted interest during many decades. However, no general principles of this kind have been rigorously proven, and many conflicting results exist [3][4][5][6][7][8][9][10][11][12][13][14]. A well known extremal principle is the so-called Prigogine's minimum entropy production principle (MinEPP) which states that, through appropriate constraints and close O 2 (gas) + 2 * 2O(ads), CO 2 (gas) + 2 * CO(ads) + O(ads).
with * and ads denoting a vacant site on the surface and adsorbed atoms or molecules, respectively [30]. For thermodynamic consistency, each forward reaction is accompanied by the corresponding reverse reaction. On surfaces with high adsorbate motility, the catalytic system is well mixed and the phenomenon of bistability is one of its most prominent non-equilibrium features. Thus, two different stable non-equilibrium steady states coexist for the same control parameter values (we refer to Table 1 for a description of the control parameters of the model). The less active NESS where, a high CO coverage inhibits the dissociative O 2 adsorption, and the active NESS, where the surface is predominantly oxygen covered. The active (less active) NESS is characterised by a high (low) CO 2 production rate. [22,23].
In the so-called mean-field approach, where spatial correlations are ignored, the bistable catalytic system is mathematically described with ODEs for the coverage of chemical species involved. This approach predicts that, depending on the initial conditions, the system will reside on one of the two non-equilibrium steady states for an indefinite period of time [23]. However, this prediction breaks down when the surface area is of meso or nanoscale dimensions. By decreasing the area of the catalytic surface, random coverage fluctuations become important and transitions, between the two non-equilibrium steady states, occur [25][26][27][28]. Thus, in order to take into account these random aspects, a stochastic mean-field description is needed. Such a stochastic approach and its macroscopic deterministic limit are discussed in the following sections.

Theoretical Framework
In this section, we present the chemical master equation (CME) of the catalytic reaction model, valid in the limit of very fast diffusion and Markovian dynamics [31,32]. We also introduce the average stochastic expression for the entropy production rate. From the stochastic description, we derive a set of ODEs for the CO and oxygen coverages on the surface. The associated macroscopic entropy production rate is also presented. In the next section, we explore whether the entropy production is maximised, minimised, or achieves no extremum in the non-equilibrium steady states.
The random jumps and waiting times are normally obtained using the so-called Gillespie algorithm [24,33]. When the system is in state Z, the next reaction σ to occur is a random variable with probability given by while the waiting time τ is also a new random variable which is exponentially distributed with the probability density. where Please note that the successive jumps are statistically independent. In Equations (6) and (8), the sum over the reactions runs over both the forward and reverse reactions σ = ±1, ±2, ±3. The vector Z follows as stochastic trajectory, and therefore, it is characterised by the probability P(Z; t) of having a number Z of particles on the surface at time t. The temporal evolution of this probability distribution is given by the following CME d dt where the sum runs over both the forward and reverse reactions (for this notation see [34,35]). The terms W σ stand for the transition rates for the forward and reverse reactions (see Table 1 and Appendix). The validity of the Markovian approximation is based on the idea that there exists a clear separation of time scales between the dynamics of the stochastic observables and that of the faster (microscopic) degrees of freedom that are not included in the stochastic description. Because of this time scale separation, the microscopic degrees of freedom follow an equilibrium distribution all along the stochastic trajectories. Table 1. Processes, population changes, and transition rates for our well-mixed CME treatment of the dynamics of Z = {N CO , N O } for a surface with N L available sites. Parameters k ads co and k des co represent the rate constants for CO(gas) adsorption and CO(ads) desorption, respectively. k ads o 2 and k des o are the rate constants for O 2 (gas) dissociative adsorption and O(ads) associative desorption. The parameter k co 2 is the rate constant for CO 2 (gas) dissociative adsorption, and k r is the rate constant for CO(ads) + O(ads) reaction. ζ is the coordination number or the number of nearest neighbours of a site. In a 2D regular lattice ζ can take one of the following values: 3 for honeycomb-type lattice, 4 for a square lattice, and 6 for a hexagonal lattice [31,32]. Please note that the number of free sites on the surface is In Appendix A, we demonstrate the consistency of these transition rates.

Process
Population Change Transition Rate The probability distribution relaxes toward a stationary state value, P(Z; t) = P st (Z), characterised by a NESS or the state of thermodynamic equilibrium. If the stationary state becomes the equilibrium one (P st (Z) = P eq (Z)), the detailed balance (for all the reactions σ ∈ {±1, ±2, ±3}) applies. In equilibrium, the solution of Equation (9) is a multinomial distribution (see Appendix A). For a NESS, detailed balance does not hold and an analytical expression for P st (Z; t) is not available. However, in the following, the steady state probability distribution will be obtained after averaging over individual trajectories generated by the Gillespie stochastic simulation algorithm. In the next section, we introduce the corresponding stochastic entropy production rate.

Stochastic Entropy Production Rate
The most important thermodynamic quantity measuring dissipation is the so-called entropy production rate. The entropy of the Markov process described by the probability distribution P(Z; t) is given by in units where the Boltzmann constant is equal to one, k B = 1 [34,[36][37][38]. Equation (11) is simply a contribution to entropy due to the probability distribution of the number of CO molecules and oxygen atoms on the surface. It is assumed to be valid in equilibrium and as well as in non-equilibrium situations [34,[36][37][38]. The sum over the vector Z means a double sum over N CO and N O . The time evolution of Equation (11) can be written as [34,[36][37][38] where we define the net reaction rates or mesoscopic thermodynamic fluxes from state Z − ω σ to Z as As above, in Equation (12), the sum over the reactions runs over both the forward and reversed reactions σ = ±1, ±2, ±3. Equation (12) represents the total variation of entropy, which is also obtained from where dS i dt is the entropy production rate due to internal processes and dS e dt is the net entropy flow rate, which can be positive or negative [1]. This net entropy flow rate is calculated from the exchange fluxes of entropy into and out of the system. The entropy production should meet two conditions: it should be non-negative and vanish in equilibrium. An expression that satisfies these requirements is provided by Schnakenberg [39], which relates the entropy production rate to the transition rates of the master equation as where are the so-called affinities or mesoscopic thermodynamic forces associated with the reactions σ = ±1, ±2, ±3. Because the inequality (R σ − R −σ ) ln R σ /R −σ ≥ 0, the entropy production rate is always positive in agreement with the second law of thermodynamics [34,[36][37][38]. Please note that The entropy flow is obtained by replacing Equation (12) and Equation (15) into Equation (14).
The equations presented above define the ensemble stochastic thermodynamics of our catalytic system for equilibrium and non-equilibrium conditions [40]. In a NESS, dS/dt = 0, implying which means that the net entropy exchanged with the environment must be negative. Please note that the entropy flow is easier to calculate and can be used in place of the stochastic entropy production rate, when dealing with non-equilibrium steady states. As expected, at thermodynamic equilibrium, the entropy flow and entropy production are zero. After having presented the stochastic framework to be implemented in this work, we proceed in the following section to discuss the deterministic limit of it.

Deterministic Mean-Field Description
In the macroscopic limit, the dynamics of the CO coverage (u) and the oxygen coverage (ν) are described by the following two coupled first-order nonlinear ODEs which are obtained from the so-called chemical Fokker-Planck equation (CFPE) approximation of the CME, in the limit of N L → ∞. This CFPE is obtained after truncating the Kramers-Moyal expansion of the CME (Equation (9)) [27,41,42]. The forward and reverse reaction rates or fluxes are given in Table 2. At steady state, where du/dt = dν/dt = 0, there can be one, two or three non-negative steady state solutions of the system of equations. A solution corresponds to certain coverage of CO and oxygen on the surface. The macroscopic thermodynamic equilibrium condition, where each forward reaction rate is equated to its corresponding reverse reaction (w σ = w −σ , for all the reactions σ ∈ {±1, ±2, ±3}), is characterised by one single solution (the equilibrium one). This macroscopic thermodynamic equilibrium condition leads to [7] K ads co K des co K ads as expected. However, far from equilibrium up to three non-negative steady state solutions can exist.
Here, we are interested in the bistable case where two non-equilibrium steady states are separated by an unstable one. However, first, let us introduce the associated macroscopic entropy production rate of the catalytic system. Table 2. Reaction rates or fluxes pertaining to Equations (19) and (20). Please note that the macroscopic rate constants are K ads co = k ads co , K des co = k des co , K ads o 2 = ζk ads o 2 /2, K des o = ζk des o /2, K co 2 = ζk co 2 /2, and K r = ζk r .

Process
Reaction Rates or Fluxes

Macroscopic Entropy Production Rate
In the macroscopic limit (i.e., N L → ∞), the entropy production rate to maintain the NESS is obtained from Equations (17) and (18). It is given by the following sum in units in which k B = 1 [34]. We refer to [1,40] for a derivation of this expression. However, note that Equation (22) should be evaluated at the corresponding steady state value. The sum over the reactions runs over both the forward and reverse reactions σ = ±1, ±2, ±3. Each term is the product of the net reaction rate or macroscopic thermodynamic flux (w +σ − w −σ ) and the macroscopic thermodynamic force or affinity (ln w +σ w −σ ) [34]. At thermodynamic equilibrium, there is no net entropy production rate because the thermodynamic fluxes and forces vanish. Because of the inequality (w +σ − w −σ ) ln w +σ /w −σ ≥ 0, the macroscopic entropy production rate is always positive in agreement with the second law of thermodynamics. It is also important to mention that, in the limit N L → ∞, EP = N −1 L dS i /dt [43]. Please note that this last statement follows from Equation (18).

Results and Discussion
Let us start with the deterministic analysis of the catalytic system. Then, we continue with the corresponding stochastic framework. It is known from experiments that, in the range where bistability exists, the oxygen coverage is small [23]. Therefore, for all cases, k ads o 2 = 0.2, k des o = 20, k co 2 = 0.5, k r = 50, and ζ = 4 (square lattice). The system of ODEs is solved by using the so-called Euler method and the probability distribution for the stochastic description is constructed from stochastic trajectories generated by Gillespie's algorithm together with the transition rates presented in Table 1.

Deterministic Analysis
In the macroscopic regime, where the deterministic approach is valid, the condition of thermodynamic equilibrium leads to the following equilibrium relations and k ads co = u eq 1 − u eq − ν eq k des co , as expected.
This allows us to determine the equilibrium or thermodynamic branch in the parameter space (k ads co , k des co ). Note also that, from Equation (21), one gets In Figure 1, we plot in red dashed line the equilibrium line given by Equation (25), with u eq and ν eq obtained from Equations (23) and (24) (or Equation (26)). For parameter values along this line, the system always relaxes to the equilibrium steady state, (u eq ≈ 0.044, ν eq ≈ 0.087). The figure also shows the far-from-equilibrium bistable region obtained when du/dt = dν/dt = 0 (black zone). The boundaries of this region are given by two saddle node (sn) lines meeting each other in a cusp. Inside this region, the system of ODEs admits two stable solutions and one unstable one, for the same control parameters. Each of these solutions is characterised by a certain amount of CO molecules and oxygen atoms on the surface. However, the unstable one is never observed in simulations.
For clarity, Figure 2a and b show separately the steady state behaviour of u and ν as a function of k ads co , for k des co = 0.15. Figure 2a shows a branch along which the CO coverage is small (denoted as u − ) and a branch along which the CO coverage is high (denoted as u + ). Figure 2b shows the corresponding bifurcation diagram for the oxygen coverage, in which we can also observe a lower branch denoted as ν − and an upper branch denoted as ν + . The figures also show the existence of a region inside which the two NESS branches coexist with a middle one (red dots), which is always unstable (denoted as u saddle in the case of CO coverage and ν saddle in the case of oxygen coverage). The boundaries of these bistable regions are characterised by the annihilation of the corresponding stable NESS branches with the saddle or unstable one (a saddle node bifurcation occurs at the sn l and sn h points) [23]. This unstable branch is not observed in simulations. The bistable region constitutes a so-called hysteresis loop. When k ads co increases from zero, the system leads simultaneously along u − and ν + until they disappear simultaneously in a saddle node bifurcation at the sn h point, then switches to u + and ν − , respectively. As k ads co decreases from above, the system remains along u + and ν − until it turns around the low sn l point and jumps to u − and ν + , respectively. Thus, the bistable region consists of the stable branches (u − , ν + ) and (u + , ν − ) coexisting with the unstable one given by (u saddle , ν saddle ).
The stability of the branches under weak perturbations is normally analysed in terms of the Jacobian matrix of Equations (19) and (20) evaluated at each of the steady state points along the branches [44]. The points along the branches (u − , ν + ) and (u + , ν − ) are stable, while the points along the saddle branch (u saddle , ν saddle ) are unstable. However, as Figure 3 shows, the selection between the two stable NESS branches depends on the initial condition of the simulation. As an example, the figure shows that starting from (u o = 0.25, ν o = 0.02) leads to a point along the branch (u − , ν + ), while starting from (u o = 0.5, ν o = 0.004) leds to a point along the branch (u + , ν − ). Outside the bistable region only one of the branches dominates no matter what the initial condition is.
We turn now to the properties of the macroscopic entropy production rate or dissipation. In Figure 4a the macroscopic entropy production rate, EP (Equation (22)), is plotted as a function of k ads co . We can observe the monostable and bistable behaviours associated with the branches presented in Figure 2. The entropy production vanishes at the thermodynamic equilibrium. As Figure 4 show, in the entire bistable interval, the dissipation along the points defined by the branch (u − , ν + ) is always larger than its value along the points associated with the branch (u + , ν − ). Thus, the state eventually selected for a given set of initial conditions is not necessarily the state where the system dissipates the most (which in the bistable region is the branch (u − , ν + ), for all values of k ads co ).

Stochastic Analysis
The results presented above suggest that the macroscopic entropy production rate is not necessarily at its minimum or maximum value when the catalytic system is in a NESS. Depending on the initial conditions, the system can end up in a NESS with high or low entropy production rate.
In this subsection, we explore whether this conclusion is still valid when the dynamics of the catalytic system are under the influence of the coverage fluctuations.  (19) and (20). Inside black region the system exhibits the bistable phenomenon. Its boundaries are given by the upper saddle node (sn h ) line and the lower saddle node (sn l ) line meeting each other in a cusp. Along the dashed-red line (k ads co ≈ 0.05k des co ), the system relaxes to thermodynamic equilibrium. Other parameters are k ads o 2 = 0.2, k des o = 20, k co 2 = 0.5, k r = 50, and ζ = 4 (square lattice). When decreasing the catalytic surface area, the predictions of the CME (Equation (9)) simulated by the Gillespie algorithm deviate from the predictions of the deterministic approach. One of the most prominent features is the phenomenon of fluctuations-induced transitions between the two non-equilibrium steady states of the bistable region [45]. As an example, Figure 5 shows time series pertaining to CO coverage inside the bistable region as obtained from both approaches, the deterministic and stochastic one. Blue dashed lines are the two non-equilibrium steady states predicted by the deterministic approach. As expected, depending on the initial condition, the catalytic system relaxes towards a NESS with a large amount of CO molecules on the surface or a NESS with a small number of CO molecules. In relation to the stochastic predictions, Figure 5a shows that, if the area of the surface is large enough, the stochastic time series follow on average the deterministic trajectories. However, Figure 5b shows that random transitions between the two non-equilibrium steady states may occur when the surface area is small enough.   (22)), corresponding to the steady-state solutions of Equations (19) and (20) as a function of k ads co , for k des co = 0.15. (b) The same EP but only for values of k ads co around the bistable region. Full blue lines are the EP associated with the (u − , ν + ) and (u + , ν − ) stable branches and dotted line to the one of the unstable or saddle branch or (u saddle , ν saddle ). Other parameters are k ads o 2 = 0.2, k des o = 20, k co 2 = 0.5, k r = 50, and ζ = 4 (square lattice). In a stochastic system, a useful quantity to study bistability is the normalised steady state probability distribution, P st (Z; t), of finding a population vector Z = {N CO , N O } [27]. Throughout this work, P(Z; t), is obtained by using the so-called Gillespie algorithm together with the transition rates presented in Table 1. Figure 6 shows the behaviour of P st (Z; t) projected on the plane (N CO , N O ), for three different values of k ads co , and k des co = 0.15. P st (Z; t) exhibits bimodal distribution in all three cases, but one of the peaks is usually much more dominant.  The stationary probability distribution allows us to calculate the steady state average quantities u = ∑ Z uP st (Z), (27) and where u = N CO /N L and u = N O /N L . Figure 7 shows u , and ν versus k ads co for different system sizes. The figure shows that the average oxygen coverage and CO coverage exhibit single-valued rather than the S-shaped multivalued behaviour of the deterministic approach. As in the case of equilibrium first-order phase transition, a transition from the steady state of low CO coverage and high oxygen coverage to the steady state of high CO coverage and low oxygen coverage occurs at a critical point, and it becomes sharper as N L increases. Hence, the state with low CO coverage and high oxygen coverage is selected below the critical point, while above this point the state with high CO coverage and low oxygen coverage is selected. Such kinetic phase transitions were first observed in the 80's for a CO oxidation model that did not, however, obey microscopic reversibility (as it contained irreversible reactions) [46].
In Figure 8 we show the corresponding stochastic entropy production rate (Equation (15)) as a function of k ads co . As expected, the stochastic entropy production rate also vanishes at the thermodynamic equilibrium. Note also that, outside the bistable region, the stochastic entropy production rate follows the trend of the macroscopic entropy production rate, in particular close to equilibrium (for a comparison see Figure 4 and corresponding discussions). It also shows that, in contrast to the macroscopic entropy production rate, this quantity always follows a monostable behaviour along the bistable region (see Figure 4a,b for a comparison with the deterministic entropy production, EP). However, there is a critical point around which a transition from a state with high entropy production (associated with the branch with low CO and high oxygen coverage) to a state with low entropy production (associated with the branch with high CO and low oxygen coverage) occurs. Thus, a state of high (low) entropy production is selected below (above) the critical point.

Overall CO 2 Production Rate
Finally, it is interesting to compare the overall CO production rate as calculated from the deterministic and stochastic approaches. For our deterministic mean-field description, with ζ = 4, the deterministic overall reaction rate is given by (see Figure 9a) As Equation (29) shows, this expression can be zero, positive, or negative. When it is negative the CO 2 (gas) dissociative adsorption dominates, while when it is positive the CO(ads) + O(ads) reaction dominates. It is to equal zero only when the system is at thermodynamic equilibrium (in Figure 9a this occurs at k ads co ≈ 0.0075). When the parameter k ads co is varied, the two steady states define two separate stable branches that we call upper rate (UR) and lower rate (LR). The UR correspond to the points along the (u − , ν + ) branch, while the LR correspond to the points along the (u + , ν − ) branch. When k ads co increases from zero, the reaction rate leads along the UR until it disappears in a saddle node bifurcation at the sn h point, then switches to LR. As k ads co decreases from above, the reaction rate remains along LR until the branch turn around the sn l point and the reaction rate jumps to the UR (a hysteresis loop occurs). From Figures 2 and 4, it is clear that, all along the bistable region, the UR has a higher entropy production rate.
The stochastic overall CO production rate is calculated using the Gillespie's algorithm. We generate a reaction rate time series where each reading of the reaction rate is R(t k ) = (r −3 − r +3 )/(N L ∆t), with ∆t = t k − t k−1 = 0.1. The terms r −3 and r +3 are the numbers of elementary events of reaction CO(ads) + O(ads) and CO 2 (gas) dissociative adsorption, that occur during the time ∆t, respectively. Then, to obtain r CO 2 , we take the average over the number of readings. Instead of the hysteresis behaviour presented in Figure 9a, Figure 9b shows that, around a critical point, r CO 2 exhibits a monotonic transition from the UR to LR, which becomes sharper as the system size N L increases. Below the critical point the UR is selected, while above that point the LR is the one selected. These correspond respectively to the high and low entropy production rates calculated in Figure 8. (b) Figure 9. (a) Deterministic overall CO 2 production rate; (b) Stochastic overall CO 2 production rate at stationary state. In both cases, k ads co was variated and k des co was fixed at 0.15. Other parameters are k ads o 2 = 0.2, k des o = 20, k co 2 = 0.5, k r = 50, and ζ = 4 (square lattice).

Summary and Conclusions
In this work, we explored the non-equilibrium thermodynamics and non-linear dynamics of a model for a catalytic surface reaction with emphasis on the question whether the entropy production rate is maximised, minimised, or does not achieve any extremum at a NESS. We considered a minimalistic model for the bistable catalytic CO oxidation reaction on single crystal surfaces, and explored the deterministic and stochastic aspects of it. From the deterministic approach, we derived the region of the control-parameter space where bistability is observed and we identified the equilibrium branch of the system. This phase diagram allows us to identify, inside the bistable region, a hysteresis loop in which a NESS characterised by a high CO 2 production rate (reactive state) coexists with a NESS with a low CO 2 production rate (less reactive state). We identified the characteristic S-shaped behaviour of bistability. The deterministic approach also predicts that, depending on the initial conditions, the system will reside in one of the two stable non-equilibrium states for an indefinite period of time.
The corresponding stochastic model exhibits, due to coverage fluctuations, random transitions between those non-equillibrium states. These transitions were analysed using the whole probability distribution of the catalytic system. We used this probability distribution to calculate the CO coverage and the oxygen coverage inside and outside the bistable region. We found that, instead of the deterministic S-shaped behaviour observed when k ads co is varied, the stochastic framework predicts monotononic variations of these quantities. This indicates a first-order phase transition at a critical value of k ads co = k which becomes sharper as the size of the system, N L , increases. Hence, for k ads co < k the high reactive state with a low CO coverage is selected, while for k ads co > k the less reactive state with a high CO coverage is selected.
To analyse the non-equilibrium thermodynamic behaviour of the system, we calculated the macroscopic (deterministic) entropy production rate as well as its corresponding stochastic version. The macroscopic entropy production predicts that inside the bistable region, there are three entropy production rates, one for each steady state. Depending on the control parameters and initial conditions, the system can reach a stable state of low or high entropy production rate. Inside this bistable region, the state with a high entropy production rate corresponds to a NESS with a high CO 2 production rate, while the state of low entropy production rate corresponds to a NEES characterised by a low CO 2 production rate. In the stochastic description, there is a unique entropy production rate. Inside the bistable region, the stochastic entropy production rate exhibits a first-order phase transition at a critical point, which becomes sharper as N L increases. For k ads co below the critical point (for k ads co above the critical point), a state with high (low) entropy production is selected. These observations suggest that, in our model, the entropy production rate is neither maximised nor minimised for all non-equilibrium states.
It is important to emphasise that the results presented in this work only correspond to deterministic and average stochastic interpretations of entropy production. A more formal study will require an interpretation in terms of trajectory stochastic thermodynamics, information theory, and the recently established notion of stochastic least-action principle for dissipative systems [3,8,37,47]. This is an important and interesting extension of our work that will be explored in the near future. For example, one can wonder about the relation, as far as entropy production is concerned, between the deterministic trajectory and the set of trajectories around it associated with the presence of coverage fluctuations. In the present stochastic approach, these stochastic trajectories are averaged out. It should be also very convenient to extend this study to other type of dissipative structures observed in catalytic surface reaction systems, like for example oscillatory and excitable dynamics [48,49].
As a final note on the future extensions of this work, we would like to emphasise that, when dealing with non-equilibrium mesoscopic systems, we need to make a proper distinction between the predictions of the stochastic approach and those of the ODEs based on the traditional law of mass action kinetics. For example, it has been recognised that the size of the system, which does not appear in the ODEs, is an important bifurcation parameter [50,51]. In the context of catalytic surface reactions, a work along these lines was presented in [52]. Using an irreversible minimalistic model for the catalytic CO oxidation, a shift of the bistable region was predicted as a function of the size of the catalytic surface.

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

Abbreviations
The following abbreviations are used in this manuscript:

NESS
Non-equilibrium steady state ODEs Ordinary differential equations CO Carbon monoxide CO 2 Carbon dioxide O 2 Oxygen LH Langmuir-Hinshelwood CME Chemical master equations 2D two dimensional EP Macroscopic entropy production rate sn Saddle node UR Upper rate LR Lower rate MinEPP Minimium entropy production principle MaxEPP Maximium entropy production principle

Appendix A. Transition Rates and Detailed Balance
Here, we show that the transition rates, W σ , satisfy the following constraint [49] with where N * = N L − N CO − N O . It is easy to show that N CO eq = N L u eq , N O eq = N L ν eq , and N * eq = N L − N CO eq − N O eq = N L (1 − u eq − ν eq ). For simplicity, we omitted the arrow used in Equation (10). As before, in Z = {N CO , N O }, N CO is the number of adsorbed CO species on the surface and N O is the number of adsorbed oxygen atoms. The stoichiometric vectors for forward reactions are ω +1 = {1, 0}, ω +2 = {0, 2}, and ω +3 = {1, 1}, and the once for the reversed reactions ω −1 = −ω +1 , ω −2 = −ω +2 , and ω −3 = −ω +3 .
satisfy Equation (A6), where K co eq = K des co /K ads co , K des co = k des co , and K ads co = k ads co Appendix A.2. Dissociative O 2 (gas) Adsorption and Associative O(ads) Desorption.