# On Oscillations and Noise in Multicomponent Adsorption: The Nature of Multiple Stationary States

Authors: Olga M. Jakšić; Zoran Jakšić; Milena B. Rašljić; Ljiljana Z. Kolar-Anić
Journal: Advances in Mathematical Physics
 (2019)

---

## Abstract

Starting from the fact that monocomponent adsorption, whether modeled by Lagergren or nonlinear Riccati equation, does not sustain oscillations, we speculate about the nature of multiple steady state states in multicomponent adsorption with second-order kinetics and about the possibility that multicomponent adsorption might exhibit oscillating behavior, in order to provide a tool for better discerning possible oscillations from inevitable fluctuations in experimental results or a tool for a better control of adsorption process far from equilibrium. We perform an analysis of stability of binary adsorption with second-order kinetics in multiple ways. We address perturbations around the steady state analytically, first in a classical way, then by introducing Langevin forces and analyzing the reaction flux and cross-correlations, then by applying the stochastic chemical master equation approach, and finally, numerically, by using stochastic simulation algorithms. Our results show that stationary states in this model are stable nodes. Hence, experimental results with purported oscillations in response should be addressed from the point of view of fluctuations and noise analysis.


---

## Body

## 1. Introduction

Adsorption processes are surface phenomena that are well understood, explored, and employed in various scientific and industrial fields [1–4]. Some mathematical models used for the interpretation of adsorption kinetics are phenomenological, obtained through analyzing experimental data, while some are derived analytically. Either way, mathematically, adsorption can be modeled either as a process with pseudo-first-order kinetics (described by linear first-order differential equations, first proposed by Lagergren [5]) or as a process with second-order kinetics (described by nonlinear Riccati equations [6].) However, every mathematical model is a simplification of what is observed in real life experiments, and real life experiments are sometimes hard to interpret due to fluctuations and noise.Natural systems are genuinely stochastic in nature. There are inherent intrinsic noise sources, extrinsic fluctuations originating from the surrounding equipment, and signal artifacts that also affect observed results. If deviations from ideal response in practical applications are reproducible, with some sort of purported regularity, we need a proper tool for interpretation of their origin with respect to difference between the oscillations which might emerge if proper conditions are met and the fluctuations which are inevitable. Oscillations and noise may relate in many different ways. If present in real life systems, they coexist. Sometimes small fluctuations and noise may invoke escalating oscillations in highly unstable systems, and sometimes oscillations are small and latently hidden in the noise, but sometimes noise may mimic oscillations and it is hard to deduce the origin of temporal deviations from the steady state, especially when the number of the samples is not very high. Stability is an important issue in many applications, as well as the appropriate modeling of the process under investigation.Stability and chemical oscillations were analyzed in a classical way by Prigogine and coworkers in the early sixties. In [7] Nicolis and Portnow addressed the general criteria about the type of chemical reactions exhibiting sustained oscillations, available even before. Although these methods are still employed in practical situations [8], new approaches emerged with time. Very complex systems can be analyzed with greater ease numerically rather than analytically, and algorithms for simulations of systems, modeled as either deterministic or stochastic ones, enabled a more informative insight into specific reactions observed. Gillespie and colleagues proposed algorithms [9–16] that are frequently applied and also implemented into a standalone software solution for multiple platforms (StochKit2). These algorithms for simulation of stochastic reaction kinetics are based on the master equation for the probabilities that the system has a given composition as a function of time, i.e., chemical master equation (CME), the mandatory method in literature on stochastic analysis [17]. However, despite the completeness of the insight into the system one could gain by knowing any given composition at any instant, it is seldom the case that the CME can be solved analytically. Then, it is often the case that Langevin approach is applied, where a known deterministic model of the system is augmented with stochastic terms representing intrinsic noise sources and random fluctuations, with the main issue being the modeling of variances of the Gaussian zero mean added terms called the Langevin forces.Apart from the focus on the reactants and products in reactions, there is also an approach where the reaction flux, the chemical reaction rate of the flow of reactants and products through the whole reaction network, is in the focus of the analysis, analogously to the analysis of generalized flux and its conjugated generalized force in the framework of nonlinear thermodynamics. The flux balance analysis approach enables the prediction of the growth of an organism in population dynamics or prediction of the rate of production of biotechnologically or chemically important metabolites and chemicals. The stochastic treatment of the flux is even more insightful because it reveals statistics about species creation and depletion [18]. Works of Bianca and Lemarchand on the reaction fluxes of the steady state and the oscillating biological and chemical systems showed how the relation between the reaction flux and the cross-correlation functions of the concentration fluctuations may be utilized in the analysis of biological and chemical systems [19–21]. This approach has been so far the most appropriate for the analysis in this paper since it allows for addressing both oscillations and fluctuations at the same time.Our motivation to contribute to this field emerged from the contradictory interpretations of the experimental results, where fluctuations can be misunderstood for oscillations and vice versa, and from the fact that adsorption is important beyond laboratories, and it is also an industrial process used in chemical, petrochemical, oil-refining, pharmaceutical, and related industries. In those complex systems, due to nonlinearities, small oscillations may escalate enough to endanger the whole operation of the chemical plant so additional investments in improvements of the automatized control of the system are necessary and economically justified [22]. Despite the fact that the phenomenon of oscillating adsorption process has been reported in praxis recently [23], we present here our doubts on the possibility of sustained oscillations in monolayer adsorption process itself, in spite of intrinsic nonlinearities in the model of the process [24].We perform a theoretical analysis of the nature of multiple steady states in multicomponent monolayer adsorption with second-order kinetics in order to provide arguments for proper interpretation of seemingly repeatable deviations from expected values in experimental results. We perform the analysis in multiple ways. We addressed perturbations around the steady state analytically, first in a classical way, then by introducing Langevin forces and analyzing reaction flux and cross-correlations, then by the use of stochastic chemical master equation approach, and finally, numerically, by the use of stochastic simulation algorithms.


## 2. Materials and Methods

In order to analyze multicomponent adsorption using various approaches, we form the stoichiometric equation. In chemisorption, there is a chemical reaction between an adsorbed particle and an adsorption center, whereas it is not the case with physisorption, where no chemical bonds are made. However, in both cases adsorption is considered here as the process that is well described by the same mathematical formalism. The stoichiometric equation valid for general multicomponent adsorption ofr different adsorbate species in gas/liquid phase (Ag,i) competing for free adsorption sites (Af) on the same surface, reversibly forming adsorbed particles (Aa,i), without dissociation or mutual interaction, is(1)Ag,i+Af⇄kd,ika,iAa,ii=1,…,rwhere ka,i and kd,i are adsorption and desorption rate constants, respectively. We assume a homogeneous surface with equivalent adsorption sites with no dissociative adsorption. Let us denote the total number of different species in the overall system by N01,N02,…,N0r. In any given moment, N1,N2,…,Nr of them are adsorbed whereas (N01-N1),…,N0r-Nr are in gaseous phase. There are no interactions between species in the gas phase. Let NM be the overall number of adsorption sites on the surface. There is a conservation relation that says that at any instant the sum of number of overall adsorbed molecules of all species on the surface and the number of free adsorption sites must equal NM. Considering all that, the set of deterministic differential equations that correspond to the stoichiometric equation (1) is as follows.(2)dN1dt=ka1Ng1Nf-kd1N1=ka1N01-N1NM-N1-N2-⋯-Nr-kd1N1dN2dt=ka2Ng2Nf-kd2N2=ka2N02-N2NM-N1-N2-⋯-Nr-kd2N2⋮dNrdt=karNgrNf-kdrNr=karN0r-NrNM-N1-N2-⋯-Nr-kdrNrThese relations are also valid when the system reaches steady state. In the steady state the number of adsorbed molecules of each mixture component on the surface is constant and the adsorption rates equal the desorption rates. As a matter of fact, adsorption and desorption of molecules keep on occurring, making topological changes of adsorbate molecules on the surface but keeping the mean percentual surface coverage and mean number of adsorbed particles constant. These perturbations to a steady state will be addressed later. Let us now denote the numbers of adsorbate molecules on the surface for every component in the mixture in the steady state with N1s,N2s,…,Nrs. In the steady state the reaction flux of all relations in (2) equals zero. Equation (2) becomes a set of r interconnected quadratic equations where each of the steady state numbers can be written in terms of the first.(3)ka1N01-N1skaiN0i-Nis=kd1N1skdiNisi=2,…,rAfter replacing all relations from (3) into the first equation from (2), keeping the left hand side equal to zero, the steady state value N1s appears to be a solution of a polynomial of the order of 1+r.(4a)ka1NM-N1sN01-N1s-kd1N1ska1N01-N1s-∑j=2rkaj/kdjN0jN1ska1/kd1N01+N1skaj/kdj-ka1/kd1=Pr+1N1sQrN1s=0(4b)ka1N01-Na1kaiN0i-Nai=kd1Na1kdiNaii=2,…,rEach solution to the polynomial Pr+1N1s, referring to the steady state number of adsorbate species denoted by one, generates the steady state numbers of all other adsorbate species, so it appears that the adsorption process of mixtures with r adsorbate species may exhibitr+1 steady states in that particular system.This result holds for an arbitrary mixture, even for the monocomponent adsorption, which is modeled by second-order Riccati differential equation. It is shown in [24] that only one of those two potential steady states obtained by mathematical calculus is practically possible; only one solution of the two is the steady state, with the other one being a mathematical artifact as it refers to a value greater than NM, the maximal number of adsorption centers on the surface, and also the maximal number of adsorbed molecules at any instant during the process of monolayer adsorption. The article [24] also addresses the criteria for modeling adsorption with approximate pseudo-first-order kinetic model while the article [28] addresses the stochastic analysis of adsorption of an arbitrary mixture for which the use of the approximate first-order kinetic modeling is appropriate. According to the classical criteria on stability [7], the linear closed systems are stable and do not exhibit sustained oscillations; the results from [28] are in agreement and also show that the mean value of the instantaneous number of adsorbed molecules of every species equals its deterministic value and that all the higher moments can be expressed in terms of the first moment (the mean value). Since the analysis of multicomponent adsorption stability for linear systems is resolved, the purpose of this work now is to investigate the nature of the steady states in adsorption of an arbitrary multicomponent mixture modeled as a process with the second-order kinetics only.Inspired by the fact that any mixture may be approximated with fictive parameters referring to a single component gas, like using UDG option (user defined gas) or “one gas analyzer” option while operating binary gas analyzers [29, 30], we start the analysis of the possibility that the adsorption process of an arbitrary mixture may exhibit multiple steady states, with speculations referring to adsorption process of binary mixtures.Stability is addressed firstly in a classical way, by analyzing the Jacobian of the Taylor expansion of the equation set (2), and then the procedure based on the Langevin approach is applied. Langevin forces and the relation between reaction flux and cross-correlation functions are modeled according to the procedures presented in [19–21]. Then stochastic analysis, based on the chemical master equation, is applied. Final conclusion on the stability of adsorption of an arbitrary mixture is based on the fact that the analysis is general, with the only constraints being that the rate constants are positive real numbers and all numbers of molecules are integers. In the case kinetics and surface occupancy by any mixture can be modeled as if it were a single gas, just like when operating a binary gas analyzer, conclusions valid for the binary mixture can be generalized to more complex systems.


## 3. Results and Discussion

### 3.1. Classical Approach

The set of interdependent deterministic matrix Riccati differential equations valid for binary adsorption is as follows.(5a)dN1dt=fN1,N2=ka1N01-N1NM-N1-N2-kd1N1(5b)dN2dt=gN1,N2=ka2N02-N2NM-N1-N2-kd2N2Here, just like in general solution (4a) and (4b), the steady states, N1s and Nr2, are interdependent. There are three possible sets of steady state numbers of adsorbate molecules in adsorption of a binary mixture. They can be calculated from the following equations.(6)N1s3+b2N1s2+b1N1s+b0=0b0=ka1kd2NMN012kd1ka2-ka1kd2b1=NMN01-N01ka1kd2N01+kd1ka2N02+ka1kd2NM+kd1kd2kd1ka2-ka1kd2b2=ka1kd2N01+kd1ka2N02kd1ka2-ka1kd2-NM-N01-kd1ka1(7)N2s=kd1ka2N02N1ska1kd2N01+kd1ka2-ka1kd2N1sWe seek the solutions in the form of positive real numbers smaller than or equal to the maximal possible number of the adsorbed particles, NM, keeping all other solutions neglected as a mathematical result without physical meaning. We start the analysis of the stability of steady states considering small perturbation of the system from a steady state. The kinetic equations set (5a) and (5b) written in the form of Taylor expansion around the steady state is as follows.(8a)dN1s+δN1dt=dδN1dt=fN1s,N2s+∂f∂N1sδN1+∂f∂N2sδN2+⋯(8b)dN2s+δN2dt=dδN2dt=gN1s,N2s+∂g∂N1sδN1+∂g∂N2sδN2+⋯The functions f and g equal zero in a steady state and the higher order terms may be omitted if perturbation from the steady state is small, so equation (8a) and (8b), written in the matrix form, is as follows.(9a)dδNdt=J×δNJ=∂f∂N1∂f∂N2∂g∂N1∂g∂N2s=ABCDN is the column matrix of the elementsN1 andN2, and the expressions that correspond to the elements of the Jacobian J are as follows.(9b)A=∂f∂N1s=2ka1N1s+ka1N2s-ka1N01+NM+kd1B=∂f∂N2s=-ka1N01-N1sC=∂g∂N1s=-ka2N02-N2sD=∂g∂N2s=2ka2N2s+ka2N1s-ka2N02+NM+kd2The solution to (9a) and (9b) implies an exponential regression to/from steady state or equilibrium.(10a)δN1t=C11eλ1t+C12eλ2t(10b)δN2t=C21eλ1t+C22eλ2tHere, Cij are constants related to the initial perturbation from the steady state, while the time constants λ1,2 written in terms of the Jacobian J are as follows.(11)λ1,2=12trJ±trJ2-4det⁡JBearing in mind that the Jacobian refers to the steady state, we analyze now the nature of the time constants. In order to determine the signs of tr(J), det⁡(J) and of the discriminant in expression for lambda, we investigate the signs of A, B, C, and D first.If we employ relations valid for steady states, we transform the expressions forA and D into the following.(12)A=ka1N1+ka1N01N2-NMN1D=ka2N2+ka2N02N1-NMN2In monolayer adsorption the instantaneous number of adsorbed particles of any species can never be greater than the initial number of particles, and the instantaneous number of all occupied adsorption sites on the surface can never exceed NM, so(13)A≤ka1N011+N2-NMN1=ka1N01N1N1+N2-NM≤0(14)D≤ka2N021+N1-NMN2=ka2N02N2N1+N2-NM≤0and, hence, the Jacobian trace is always negative, B and C elements are always negative, and the discriminant in the expression for lambda is always positive.Slightly more cumbersome algebra shows that the Jacobian determinant is always positive, which means that the discriminant root in (11) is always less than the Jacobian trace. Both of the time constants λ1,2 (the eigenvalues of the Jacobian matrix) are always real negative numbers. That means that the steady state (N1s,N2s) is astable nodeand oscillations cannot be sustained. All departures from the only practically feasible steady state end up with exponential regression back towards that one and the only steady state without the chance to reach the state of a Hopf bifurcation point under any circumstances.By analyzing the stability of binary adsorption with this procedure, we came to conclusion about the nature of the steady states without calculating or knowing their exact value. We analyzed the deviations from the steady state and the relations valid for the steady state, but the value of the steady state was not needed for the application of criteria in this procedure. This is important because the analysis of the cubic polynomial is not needed.So far we did not address the stochastic nature of the system, its intrinsic fluctuations, and causes or frequencies of the deviations from the steady state.


### 3.2. Langevin Forces and Correlation Functions

In this approach we focus on the part of the system we are interested in as if it were purely deterministic, and we consider the rest as a statistical bath [17]. We do not consider the Taylor expansion of functions in kinetic equations now; instead we rewrite the original kinetic equations replacing the instantaneous number of adsorbed molecules N with the sum of its steady state value Ns and the instantaneous deviation from the steady state δN. As before, the indices refer to a particular species in the binary mixture.(15a)dN1s+δN1dt=dδN1dt=ka1N01-N1s+δN1NM-N1s+δN1-N2s+δN2-kd1N1s+δN1(15b)dN2s+δN2dt=dδN2dt=ka2N02-N2s+δN2NM-N1s+δN1-N2s+δN2-kd2N2s+δN2We assume that the departures from the steady state are small so it is convenient to neglect all the products resulting from their multiplication. We also neglect all terms related to the reaction flux, which does not comprise infinitesimally small variables but vanishes at the steady state. However, we add the Langevin forces, random terms to interpret randomness in departures and as a compensation for all of the neglected terms. The added terms are zero mean Gaussian noise terms ξN1 and ξN2 whose variances will be addressed later in this section. (16a)dδN1dt=κ1δN1+κ2δN2+ξN1t(16b)dδN2dt=κ3δN1+κ4δN2+ξN2tBy doing so we linearized the equation set (15a) and (15b). The constants in (16a) and (16b) are defined by the following expressions.(17a)κ1=-ka1N01-N1s+NM-N1s-N2s-kd1κ2=-ka1N01-N1s(17b)κ3=-ka2N02-N2sκ4=-ka2N02-N2s+NM-N1s-N2s-kd2In order to investigate the stability of the steady state, we rewrite the equation (16a) and (16b) in the matrix form and compute the eigenvalues of the matrix M.(18)M=κ1κ2κ3κ4Now eigenvalues are as follows.(19)λ±=12trM±trM2-4detMBearing in mind that expressions in (18) refer to the sorption rate constants, the overall numbers of molecules of reactant speciesN01 andN02, and the overall number of adsorption centers on the surface NM that are all real positive numbers, we see that the trace of the matrix M is always negative as it never reaches zero, which means that a Hopf bifurcation can never occur. That is in accord with the stability analysis performed without Langevin linearization.Langevin forces have zero mean but nonzero variance. In order to determine the variances of the Langevin forces, we focus on the stoichiometric equations and the reaction flux in equilibrium. In case the reactions are reversible, as they are in the model of adsorption, forward and backward fluxes are treated separately. In the case of binary adsorption we have the following.(20)Ag1+Af→ka1Aa1Aa1→kd1Ag1+AfAg2+Af→ka2Aa2Aa2→kd2Ag2+AfAlthough movements of molecules of different species are independent in the same manner as they are in any Brownian motion with Maxwell-Boltzmann distribution, on the surface they fight for the same area so we have the additional conservation equation that says that at any moment the sum of the number of adsorbed molecules Ag1 and Aa2 and the number of free adsorption sites Af must equal the overall number of adsorption sites on the surface, NM. Additionally, for each of both species the initial number of molecules in the system,N0, must equal the sum of the number of its adsorbed molecules Aa and the number of the molecules in the gas phase Ag.For every stoichiometric equation and every variable we are interested in, we define one stochastic term; let us denote it byζ multiplied by the square root of adsorption or desorption rate in the steady state. Then we form the Langevin force for every stochastic variable, taking into account all stoichiometric equations in which it appears.Langevin forces in (16a) and (16b) are given by the following.(21a)ξN1t=ka1N01-N1sNM-N1s-N2sζ1t-kd1N1sζ2t(21b)ξN2t=ka2N02-N2sNM-N1s-N2sζ3t-kd2N2sζ4tAll stochastic terms ζ are zero mean Gaussian white noise terms with unit variance, mutually uncorrelated, with their cross-correlations being zero. Hence, the variances of the Langevin forces are(22a)ξN1tξN1t′=F11δt-t′(22b)ξN1tξN2t′=F12δt-t′(22c)ξN2tξN2t′=F22δt-t′where(23a)F11=ka1N01-N1sNM-N1s-N2s+kd1N1s(23b)F12=0(23c)F22=ka2N02-N2sNM-N1s-N2s+kd2N2s.In order to gain a better insight into the adsorption process and how it behaves at very short times and high frequencies, we focus now on the correlation functions. We use the procedure in [19–21] and determine the correlation functions after the change of the basis by the use of the following matrix. (24)P=P11=κ2κ22+λ+-κ12P12=κ2κ22+λ--κ12P21=λ+-κ1κ22+λ+-κ12P22=λ--κ1κ22+λ--κ12The matrix P relates δN1 and δN2 with new coordinates.(25)δN1tδN2t=PXtYtAfter finding the correlation functions of the new coordinates and using the change of the basis, in the long time limit when all transients are gone and only the dependence on the time difference τ=t–t′ remains, then the difference of cross-correlations of the deviations to the steady state is given by the following.(26)Iτ=δN10δN2τ-δN1τδN20=-P21P22F11-P11P22+P12P21F12+P11P12F22P11P22-P12P21eλ+τ-eλ-τλ++λ-After multiple rearrangements, by using all expressions valid for the matrix elements, the variances of the Langevin forces, and time constants in (26) and by using all relations valid for the steady state, the expression for the difference of cross-correlation functions is obtained as follows.(27)Iτ=Φkd1N1s-kd2N2sNM-N1s-N2sκ1-κ42+4κ2κ3eλ+τ-eλ-τκ1+κ4The reaction flux, related to the original equation set (5a) and (5b), before linearization, vanishes at equilibrium and equals(28)Φ=ka1N01-N1sNM-N1s-N2s-kd1N1s=ka2N02-N2sNM-N1s-N2s-kd2N2s.Although the Langevin approach is widely used in practical applications, it does not give the full insight into all statistical parameters of the process as a fully stochastic approach by the use of the chemical master equation would do.


### 3.3. Stochastic Approach: Chemical Master Equation

Instead of investigating deterministic equations in the vicinity of a steady state or stochastic Langevin differential equations with added random white noise, here we analyze adsorption as a stochastic process and work with probabilities that instantaneous numbers of adsorbate molecules on the surface may have a specified value. By doing so, we apply the conventional procedure of stochastic analysis and assume that the time span for our system observations,Δt, is short enough, so that only one transition between the neighboring states may occur. That means that the system can reach the state withN1 adsorbed molecules of the first species andN2 adsorbed molecules of the second species in five ways: one molecule of the species denoted by 1 has been adsorbed, one molecule of the species denoted by 2 has been adsorbed, one molecule of the species denoted by 1 has been desorbed, one molecule of the species denoted by 2 has been desorbed, and no transition occurred. The probability that the system will be in a state withN1 adsorbed molecules of the first gas species andN2 adsorbed molecules of the second gas species in the moment t+Δt is the sum of these five probabilities.(29)PN1,N2t+Δt=PN1-1,N2tpN1-1→N1+PN1,N2-1tpN2-1→N2+PN1+1,N2tpN1+1→N1+PN1,N2+1tpN2+1→N2+PN1,N2tpN1→N1,N2→N2The probability that no transition occurs is as follows.(30)pN1→N1,N2→N2=1-pN1→N1-1-pN1→N1+1-pN2→N2-1-pN2→N2+1The transition probabilities among neighboring states are proportional to the instantaneous sorption rates, so(31)pN1-1→N1=ka1N01-N1-1M-N1-1-N2ΔtpN2-1→N2=ka2N01-N1-1M-N1-N2-1ΔtpN1+1→N1=kd1N1+1ΔtpN2+1→N2=kd2N2+1ΔtpN1→N1+1=ka1N01-N1M-N1-N2ΔtpN1→N1-1=kd1N1ΔtpN2→N2+1=ka2N02-N2M-N1-N2ΔtpN2→N2-1=kd2N2Δt.In the short time limit, after rearranging terms in (29), we obtain the differential equation for the probabilities, i.e., the chemical master equation.(32)dPN1,N2tdt=PN1-1,N2tka1N01-N1-1M-N1-1-N2+PN1,N2-1tka2N01-N1-1M-N1-N2-1+PN1+1,N2tkd1N1+1+PN1,N2+1tkd2N2+1-kd1N1+ka1N01-N1M-N1-N2+kd2N2+ka2N02-N2M-N1-N2In order to find the solution to the chemical master equation governing the adsorption process of binary gas mixture, we need a bivariant probability generating function. (33)Fs1,s2,t=∑N2=0∞∑N1=0∞PN1,N2ts1N1s2N2s1,s2<1Its first derivative with respect to s1, calculated fors1 = 1 ands2 = 1, is then used for the calculations of the first moment of the number of adsorbed molecules for the first gas species.(34)∂Fs1,s2,t∂s1=∑N2=0∞∑N1=0∞N1PN1,N2ts1N1-1s2N2s1=s2=1=N1¯Its first derivative with respect tos2, calculated fors1 = 1 ands2 = 1, equals the first moment of the number of adsorbed molecules for the second gas species.(35)∂Fs1,s2,t∂s2=∑N2=0∞∑N1=0∞N2PN1,N2ts1N1s2N2-1s1=s2=1=N2¯Its second derivative with respect tos1 ands2 is as follows.(36)∂2Fs1,s2,t∂s1∂s2=∑N2=0∞∑N1=0∞N1N2PN1,N2ts1N1-1s2N2-1s1=s2=1=N1N2¯≠N1¯·N2¯Sums go to infinity but their terms are zero beyond physical limits. In order to find the solution to CME (32) we first multiply it by ∑N1=0∞s1N1 and then by ∑N2=0∞s2N2 and then rearrange its members grouping together the terms that form the probability generating function (PGF) or the derivatives of PGF. By doing so, we transform CME, differential equation on probabilities (32), into a partial differential equation on PGF,F(s1,s2,t). Once solved, the partial differential equation onF(s1,s2,t), would reveal all statistical data on the process with all higher moments.Since it is a cumbersome and possibly unfeasible task, with no great importance for the applications where knowing of the first moments only is crucial, it is often much faster to seek for the solution by simulations or the analysis of the moments only. Here, we will address it using simulation algorithms based on the chemical master equation.


### 3.4. Stochastic Approach: Simulation Algorithms

Our routines for visualization of stochastic transitions in multicomponent adsorption of molecules are based on stochastic simulation algorithms (SSA) developed by Gillespie and colleagues [9–16]. The routines are part of our custom, in-house designed, permanently growing software solution. Features of its current version, ADmoND, written for MathWorks MATLAB release 2013, will be briefly outlined in next section and the whole package is available at Mendeley Data repository [31].The propensity functions, used for modeling the probabilities of transitions between states and for the modeling of the time sequence of random moments for every transition, are the same as in the chemical master equation (32). Here, the mean time between transitions of the molecules of each species in the mixture is a random variable with the exponential distribution whose parameters are related to propensity function, i.e., sorption rate of the relevant transition (adsorption or desorption of the first or second component in the mixture). The exponential distribution is designed using in-built function for uniform distribution in MATLAB and inverse CDF sampling technique method [32]. After generating the sequences of time and the instantaneous numbers of adsorbed molecules for each component in the mixture, it is easy to obtain the time evolutions of the number of adsorbed molecules for every component in the mixture, the instantaneous changes in the phase space, and the correlation functions.Here we give results obtained by the simulation algorithm valid for binary adsorption of the mixture of Sarin and Sulphur Mustard on graphene surface. Table1 lists parameters of the mixture. The sorption rate constants are calculated by the procedure from [24], and the desorption energy and the surface density of adsorbed particles, needed for the calculations of the sorption rates, are taken from the online Harvard Dataverse repository [25] and from [26, 27]. Assumed volume of the reaction chamber is 1 dm3, temperature is 282 K, and the number of adsorption centers on the surface, NM, is 38.17.1014 (the surface density of graphene atoms being 38.17.1018 and surface area being 1 cm2). Figure 1(a) shows phase space assuming that, at the initial moment, the surface was clean with no adsorbed amount. Figure 1(b) shows the corresponding cross-correlation function. Due to very high desorption rate constants (Table 1) and high number of adsorption sites on the surface, 1 million transitions were not enough to show the reaching of equilibrium. The corresponding correlation function is symmetric.Table 1
Parameters used in simulations of binary adsorption on graphene in Figure1.
Binary mixture component pressure [Pa] k a [molecules/s] k d [1/s] References Sarin 46 2.89.10−14 1.7.109 [24–26] Sulphur Mustard 360 1.17.10−17 1.57.108 [24, 25, 27]Figure 1
Phase space (a) and the corresponding cross-correlation function (b) of numbers of adsorbed molecules of Sarin and Sulphur Mustard on graphene assuming clean surface in the beginning.
(a) (b)Figure2 shows simulations of adsorption of fictive binary mixture assuming the number of adsorption centers on the surface, NM, is 38.17.105, adsorption rate constant and desorption rate constant of the first component 2.89 and 1.7, respectively, and adsorption rate constant and desorption rate constant of the second component 1.17 and 1.57, respectively, assuming the initial condition is the clean surface with no adsorbed molecules. Other parameters are the same as for Figure 1. On the left, up, and down, on time evolutions of numbers of adsorbed molecules, we see either that reaching the equilibrium was fast and that there was just one transient or that we need more than 10 million transitions in simulations in order to see vanishing of the other transient. In order to see exact equilibrium levels, we would have to check with the analytical solutions. On the right, we see that every simulation led to the same node and that the cross-correlation functions were symmetrical.Figure 2
Time evolutions of the number of adsorbed molecules (left), phase space (up right), and the corresponding cross-correlation function (down right), fictive mixture, clean surface at the beginning.Figure3 shows results for the same fictive mixture as Figure 2, with the only difference being the starting moment of observation, which is now random. Phase space, shown on the upper diagram on the right side, shows that all different starting points are leading to the same node. The diagram on the lower right part in Figure 3 shows that the cross-correlation function is not ideally symmetric in this case. Time evolutions in Figure 3 (top and bottom diagram on the left) demonstrate now two changes in the process kinetics. Due to our analysis we now know that there are no two consecutive steady states, and change in kinetics is the consequence of the vanishing of multiple transients (for binary mixture the exact solution would give at least two transients). Exact solution to multicomponent adsorption modeled with linear Lagergren model, published in article [28], states that in linear systems the number of transients equals the number of components in the mixture. The exact solution to nonlinear systems would reveal that there are more transients than there are components in the mixture but they may not be always visible in the observed time and amplitude span. The quality of results also depends on the quality and features of used software.Figure 3
Time evolutions of the number of adsorbed molecules (left), phase space (up right), and the corresponding cross-correlation function (down right), fictive mixture, random initial moment.


### 3.5. Algorithms and ADmoND Software for Adsorption Analysis

The software package ADmoND we used is aimed at the analysis of monolayer adsorption and desorption processes in novel devices. It has been developed for the software platform MathWorks MATLAB, 2013a release. ADmoND can be used as a MATLAB toolbox or MATLAB application with graphical user interface (GUI) that is adapted MATLAB GUI for figure handling, extended with additional tabs (System, Evolutions, and Equilibrium). All calculations are based upon analytically derived expressions, implemented in mathematical forms optimized to avoid or minimize numerical error propagation. Theoretical background is based on well-known laws or new original solutions published in peer-reviewed journals. The full software package is openly available and can be downloaded from Mendeley Data repository [31].The 'System' tab is aimed at handling parameters saved outside MATLAB and interactive setting of system parameters (temperature, adsorbent area, reaction chamber volume, and composition of an adsorbate mixture). Figure4, generated with ADmoND, illustrates adsorption of binary gas mixture on a solid surface.Figure 4
Competitive adsorption between particles of two species for equivalent adsorption sites on homogeneous surface.The 'Evolutions' tab of ADmoND is dedicated to calculations related to adsorption kinetics in time domain. Figure5 illustrates possible time evolutions of monocomponent adsorption (green lines) obtained by stochastic simulation algorithm (SSA, inspired by [9], adapted for ADmoND case studies), their mean (magenta line), and the ideal deterministic solution (red dashed line).Figure 5
Time evolution of monocomponent adsorption obtained by stochastic simulation algorithm (SSA, green lines), mean of SSA ensemble (magenta line), and deterministic solution (red dashed line).The 'Equilibrium' tab is dedicated to calculations related to adsorption dynamics in equilibrium. As we saw, for every mixture ofr components, there arer+1 possible sets of r numerical values,r+1 possible steady states. Each possible steady state is a vector ofr numbers, where each number corresponds to the number of adsorbed molecules for every component in the mixture (Nis). ADmoND findsr+1 roots of the polynomial of ther+1 order (candidates for the steady state of component 1) then computes all other possible steady states for all other components in the mixture. Afterwards, it performs inspection if exactly one set of solutions for the number of adsorbed particles of mixture components in steady stateNs is realistic and satisfies three conditions. The first condition is that the steady state values for the number of adsorbed molecules of particular components must be positive numbers, lower thanNM. The second condition is that the sum of all steady state values for the number of adsorbed molecules of particular components on the surface may never exceed NM. The third condition is that the steady state values for the number of adsorbed molecules of each component in the mixture must be positive number, lower than the overall number of its molecules in the system, N0i.The theoretical proof that there is just one steady state, the stable node, ensures that there is just one set of solutions and thus enables the optimization of the code (the first set that is found is the only set; there is no need for inspection of all other sets).


### 3.6. Discussion of Results

The stability analysis, described above, helps in a better understanding of adsorption dynamics in equilibrium and a more accurate interpretation of results. In research span of almost a decade, in simulations of multicomponent adsorption, we did not encounter a case of practically feasible multiple steady states obtained by MATLAB; only one set of solutions satisfied required conditions (steady state adsorbate molecules on the surface are positive numbers with overall sum that is smaller than or equal to the maximal possible number of adsorbed particles,NM). After thorough investigation, some complex solutions that appeared once proved to be the consequence of the numerical error propagation, which resulted in our additional efforts towards code improvements.All of the methods we applied here for stability analysis in this paper result with the same conclusion: There is no combination of system parameters that ensures the appearance of sustained oscillations. On the contrary, the steady state is the stable node. Thus, if there is a seemingly oscillatory response in the experiment, the cause of oscillations is not the adsorption process (modeled with the nonlinear set of Riccati equations). Our conclusion is that the possible explanations are as follows:(i)
There are no oscillations; noise is mimicking them.(ii)
The reactions are more complex and are not well modeled by the set of Riccati equations.(iii)
The oscillations emerged due to some other cause in the control loop of complex systems used in chemical, petrochemical, oil-refining, pharmaceutical, and related plants.An example that supports this reasoning is presented in the paper [23], where self-oscillating adsorption–desorption of silver nanoparticles on the gel formed by incorporating graphene oxide into a poly (methacrylic acid)–polyethylene glycol copolymer film is reported. Here, adsorbent is a pH-responsive hydrogel which acts as a pH oscillator in a closed reaction and the oscillations in reversible adsorption–desorption of Ag nanoparticles on the patterned hydrogel surface are just response to pH oscillations.


## 3.1. Classical Approach

The set of interdependent deterministic matrix Riccati differential equations valid for binary adsorption is as follows.(5a)dN1dt=fN1,N2=ka1N01-N1NM-N1-N2-kd1N1(5b)dN2dt=gN1,N2=ka2N02-N2NM-N1-N2-kd2N2Here, just like in general solution (4a) and (4b), the steady states, N1s and Nr2, are interdependent. There are three possible sets of steady state numbers of adsorbate molecules in adsorption of a binary mixture. They can be calculated from the following equations.(6)N1s3+b2N1s2+b1N1s+b0=0b0=ka1kd2NMN012kd1ka2-ka1kd2b1=NMN01-N01ka1kd2N01+kd1ka2N02+ka1kd2NM+kd1kd2kd1ka2-ka1kd2b2=ka1kd2N01+kd1ka2N02kd1ka2-ka1kd2-NM-N01-kd1ka1(7)N2s=kd1ka2N02N1ska1kd2N01+kd1ka2-ka1kd2N1sWe seek the solutions in the form of positive real numbers smaller than or equal to the maximal possible number of the adsorbed particles, NM, keeping all other solutions neglected as a mathematical result without physical meaning. We start the analysis of the stability of steady states considering small perturbation of the system from a steady state. The kinetic equations set (5a) and (5b) written in the form of Taylor expansion around the steady state is as follows.(8a)dN1s+δN1dt=dδN1dt=fN1s,N2s+∂f∂N1sδN1+∂f∂N2sδN2+⋯(8b)dN2s+δN2dt=dδN2dt=gN1s,N2s+∂g∂N1sδN1+∂g∂N2sδN2+⋯The functions f and g equal zero in a steady state and the higher order terms may be omitted if perturbation from the steady state is small, so equation (8a) and (8b), written in the matrix form, is as follows.(9a)dδNdt=J×δNJ=∂f∂N1∂f∂N2∂g∂N1∂g∂N2s=ABCDN is the column matrix of the elementsN1 andN2, and the expressions that correspond to the elements of the Jacobian J are as follows.(9b)A=∂f∂N1s=2ka1N1s+ka1N2s-ka1N01+NM+kd1B=∂f∂N2s=-ka1N01-N1sC=∂g∂N1s=-ka2N02-N2sD=∂g∂N2s=2ka2N2s+ka2N1s-ka2N02+NM+kd2The solution to (9a) and (9b) implies an exponential regression to/from steady state or equilibrium.(10a)δN1t=C11eλ1t+C12eλ2t(10b)δN2t=C21eλ1t+C22eλ2tHere, Cij are constants related to the initial perturbation from the steady state, while the time constants λ1,2 written in terms of the Jacobian J are as follows.(11)λ1,2=12trJ±trJ2-4det⁡JBearing in mind that the Jacobian refers to the steady state, we analyze now the nature of the time constants. In order to determine the signs of tr(J), det⁡(J) and of the discriminant in expression for lambda, we investigate the signs of A, B, C, and D first.If we employ relations valid for steady states, we transform the expressions forA and D into the following.(12)A=ka1N1+ka1N01N2-NMN1D=ka2N2+ka2N02N1-NMN2In monolayer adsorption the instantaneous number of adsorbed particles of any species can never be greater than the initial number of particles, and the instantaneous number of all occupied adsorption sites on the surface can never exceed NM, so(13)A≤ka1N011+N2-NMN1=ka1N01N1N1+N2-NM≤0(14)D≤ka2N021+N1-NMN2=ka2N02N2N1+N2-NM≤0and, hence, the Jacobian trace is always negative, B and C elements are always negative, and the discriminant in the expression for lambda is always positive.Slightly more cumbersome algebra shows that the Jacobian determinant is always positive, which means that the discriminant root in (11) is always less than the Jacobian trace. Both of the time constants λ1,2 (the eigenvalues of the Jacobian matrix) are always real negative numbers. That means that the steady state (N1s,N2s) is astable nodeand oscillations cannot be sustained. All departures from the only practically feasible steady state end up with exponential regression back towards that one and the only steady state without the chance to reach the state of a Hopf bifurcation point under any circumstances.By analyzing the stability of binary adsorption with this procedure, we came to conclusion about the nature of the steady states without calculating or knowing their exact value. We analyzed the deviations from the steady state and the relations valid for the steady state, but the value of the steady state was not needed for the application of criteria in this procedure. This is important because the analysis of the cubic polynomial is not needed.So far we did not address the stochastic nature of the system, its intrinsic fluctuations, and causes or frequencies of the deviations from the steady state.


## 3.2. Langevin Forces and Correlation Functions

In this approach we focus on the part of the system we are interested in as if it were purely deterministic, and we consider the rest as a statistical bath [17]. We do not consider the Taylor expansion of functions in kinetic equations now; instead we rewrite the original kinetic equations replacing the instantaneous number of adsorbed molecules N with the sum of its steady state value Ns and the instantaneous deviation from the steady state δN. As before, the indices refer to a particular species in the binary mixture.(15a)dN1s+δN1dt=dδN1dt=ka1N01-N1s+δN1NM-N1s+δN1-N2s+δN2-kd1N1s+δN1(15b)dN2s+δN2dt=dδN2dt=ka2N02-N2s+δN2NM-N1s+δN1-N2s+δN2-kd2N2s+δN2We assume that the departures from the steady state are small so it is convenient to neglect all the products resulting from their multiplication. We also neglect all terms related to the reaction flux, which does not comprise infinitesimally small variables but vanishes at the steady state. However, we add the Langevin forces, random terms to interpret randomness in departures and as a compensation for all of the neglected terms. The added terms are zero mean Gaussian noise terms ξN1 and ξN2 whose variances will be addressed later in this section. (16a)dδN1dt=κ1δN1+κ2δN2+ξN1t(16b)dδN2dt=κ3δN1+κ4δN2+ξN2tBy doing so we linearized the equation set (15a) and (15b). The constants in (16a) and (16b) are defined by the following expressions.(17a)κ1=-ka1N01-N1s+NM-N1s-N2s-kd1κ2=-ka1N01-N1s(17b)κ3=-ka2N02-N2sκ4=-ka2N02-N2s+NM-N1s-N2s-kd2In order to investigate the stability of the steady state, we rewrite the equation (16a) and (16b) in the matrix form and compute the eigenvalues of the matrix M.(18)M=κ1κ2κ3κ4Now eigenvalues are as follows.(19)λ±=12trM±trM2-4detMBearing in mind that expressions in (18) refer to the sorption rate constants, the overall numbers of molecules of reactant speciesN01 andN02, and the overall number of adsorption centers on the surface NM that are all real positive numbers, we see that the trace of the matrix M is always negative as it never reaches zero, which means that a Hopf bifurcation can never occur. That is in accord with the stability analysis performed without Langevin linearization.Langevin forces have zero mean but nonzero variance. In order to determine the variances of the Langevin forces, we focus on the stoichiometric equations and the reaction flux in equilibrium. In case the reactions are reversible, as they are in the model of adsorption, forward and backward fluxes are treated separately. In the case of binary adsorption we have the following.(20)Ag1+Af→ka1Aa1Aa1→kd1Ag1+AfAg2+Af→ka2Aa2Aa2→kd2Ag2+AfAlthough movements of molecules of different species are independent in the same manner as they are in any Brownian motion with Maxwell-Boltzmann distribution, on the surface they fight for the same area so we have the additional conservation equation that says that at any moment the sum of the number of adsorbed molecules Ag1 and Aa2 and the number of free adsorption sites Af must equal the overall number of adsorption sites on the surface, NM. Additionally, for each of both species the initial number of molecules in the system,N0, must equal the sum of the number of its adsorbed molecules Aa and the number of the molecules in the gas phase Ag.For every stoichiometric equation and every variable we are interested in, we define one stochastic term; let us denote it byζ multiplied by the square root of adsorption or desorption rate in the steady state. Then we form the Langevin force for every stochastic variable, taking into account all stoichiometric equations in which it appears.Langevin forces in (16a) and (16b) are given by the following.(21a)ξN1t=ka1N01-N1sNM-N1s-N2sζ1t-kd1N1sζ2t(21b)ξN2t=ka2N02-N2sNM-N1s-N2sζ3t-kd2N2sζ4tAll stochastic terms ζ are zero mean Gaussian white noise terms with unit variance, mutually uncorrelated, with their cross-correlations being zero. Hence, the variances of the Langevin forces are(22a)ξN1tξN1t′=F11δt-t′(22b)ξN1tξN2t′=F12δt-t′(22c)ξN2tξN2t′=F22δt-t′where(23a)F11=ka1N01-N1sNM-N1s-N2s+kd1N1s(23b)F12=0(23c)F22=ka2N02-N2sNM-N1s-N2s+kd2N2s.In order to gain a better insight into the adsorption process and how it behaves at very short times and high frequencies, we focus now on the correlation functions. We use the procedure in [19–21] and determine the correlation functions after the change of the basis by the use of the following matrix. (24)P=P11=κ2κ22+λ+-κ12P12=κ2κ22+λ--κ12P21=λ+-κ1κ22+λ+-κ12P22=λ--κ1κ22+λ--κ12The matrix P relates δN1 and δN2 with new coordinates.(25)δN1tδN2t=PXtYtAfter finding the correlation functions of the new coordinates and using the change of the basis, in the long time limit when all transients are gone and only the dependence on the time difference τ=t–t′ remains, then the difference of cross-correlations of the deviations to the steady state is given by the following.(26)Iτ=δN10δN2τ-δN1τδN20=-P21P22F11-P11P22+P12P21F12+P11P12F22P11P22-P12P21eλ+τ-eλ-τλ++λ-After multiple rearrangements, by using all expressions valid for the matrix elements, the variances of the Langevin forces, and time constants in (26) and by using all relations valid for the steady state, the expression for the difference of cross-correlation functions is obtained as follows.(27)Iτ=Φkd1N1s-kd2N2sNM-N1s-N2sκ1-κ42+4κ2κ3eλ+τ-eλ-τκ1+κ4The reaction flux, related to the original equation set (5a) and (5b), before linearization, vanishes at equilibrium and equals(28)Φ=ka1N01-N1sNM-N1s-N2s-kd1N1s=ka2N02-N2sNM-N1s-N2s-kd2N2s.Although the Langevin approach is widely used in practical applications, it does not give the full insight into all statistical parameters of the process as a fully stochastic approach by the use of the chemical master equation would do.


## 3.3. Stochastic Approach: Chemical Master Equation

Instead of investigating deterministic equations in the vicinity of a steady state or stochastic Langevin differential equations with added random white noise, here we analyze adsorption as a stochastic process and work with probabilities that instantaneous numbers of adsorbate molecules on the surface may have a specified value. By doing so, we apply the conventional procedure of stochastic analysis and assume that the time span for our system observations,Δt, is short enough, so that only one transition between the neighboring states may occur. That means that the system can reach the state withN1 adsorbed molecules of the first species andN2 adsorbed molecules of the second species in five ways: one molecule of the species denoted by 1 has been adsorbed, one molecule of the species denoted by 2 has been adsorbed, one molecule of the species denoted by 1 has been desorbed, one molecule of the species denoted by 2 has been desorbed, and no transition occurred. The probability that the system will be in a state withN1 adsorbed molecules of the first gas species andN2 adsorbed molecules of the second gas species in the moment t+Δt is the sum of these five probabilities.(29)PN1,N2t+Δt=PN1-1,N2tpN1-1→N1+PN1,N2-1tpN2-1→N2+PN1+1,N2tpN1+1→N1+PN1,N2+1tpN2+1→N2+PN1,N2tpN1→N1,N2→N2The probability that no transition occurs is as follows.(30)pN1→N1,N2→N2=1-pN1→N1-1-pN1→N1+1-pN2→N2-1-pN2→N2+1The transition probabilities among neighboring states are proportional to the instantaneous sorption rates, so(31)pN1-1→N1=ka1N01-N1-1M-N1-1-N2ΔtpN2-1→N2=ka2N01-N1-1M-N1-N2-1ΔtpN1+1→N1=kd1N1+1ΔtpN2+1→N2=kd2N2+1ΔtpN1→N1+1=ka1N01-N1M-N1-N2ΔtpN1→N1-1=kd1N1ΔtpN2→N2+1=ka2N02-N2M-N1-N2ΔtpN2→N2-1=kd2N2Δt.In the short time limit, after rearranging terms in (29), we obtain the differential equation for the probabilities, i.e., the chemical master equation.(32)dPN1,N2tdt=PN1-1,N2tka1N01-N1-1M-N1-1-N2+PN1,N2-1tka2N01-N1-1M-N1-N2-1+PN1+1,N2tkd1N1+1+PN1,N2+1tkd2N2+1-kd1N1+ka1N01-N1M-N1-N2+kd2N2+ka2N02-N2M-N1-N2In order to find the solution to the chemical master equation governing the adsorption process of binary gas mixture, we need a bivariant probability generating function. (33)Fs1,s2,t=∑N2=0∞∑N1=0∞PN1,N2ts1N1s2N2s1,s2<1Its first derivative with respect to s1, calculated fors1 = 1 ands2 = 1, is then used for the calculations of the first moment of the number of adsorbed molecules for the first gas species.(34)∂Fs1,s2,t∂s1=∑N2=0∞∑N1=0∞N1PN1,N2ts1N1-1s2N2s1=s2=1=N1¯Its first derivative with respect tos2, calculated fors1 = 1 ands2 = 1, equals the first moment of the number of adsorbed molecules for the second gas species.(35)∂Fs1,s2,t∂s2=∑N2=0∞∑N1=0∞N2PN1,N2ts1N1s2N2-1s1=s2=1=N2¯Its second derivative with respect tos1 ands2 is as follows.(36)∂2Fs1,s2,t∂s1∂s2=∑N2=0∞∑N1=0∞N1N2PN1,N2ts1N1-1s2N2-1s1=s2=1=N1N2¯≠N1¯·N2¯Sums go to infinity but their terms are zero beyond physical limits. In order to find the solution to CME (32) we first multiply it by ∑N1=0∞s1N1 and then by ∑N2=0∞s2N2 and then rearrange its members grouping together the terms that form the probability generating function (PGF) or the derivatives of PGF. By doing so, we transform CME, differential equation on probabilities (32), into a partial differential equation on PGF,F(s1,s2,t). Once solved, the partial differential equation onF(s1,s2,t), would reveal all statistical data on the process with all higher moments.Since it is a cumbersome and possibly unfeasible task, with no great importance for the applications where knowing of the first moments only is crucial, it is often much faster to seek for the solution by simulations or the analysis of the moments only. Here, we will address it using simulation algorithms based on the chemical master equation.


## 3.4. Stochastic Approach: Simulation Algorithms

Our routines for visualization of stochastic transitions in multicomponent adsorption of molecules are based on stochastic simulation algorithms (SSA) developed by Gillespie and colleagues [9–16]. The routines are part of our custom, in-house designed, permanently growing software solution. Features of its current version, ADmoND, written for MathWorks MATLAB release 2013, will be briefly outlined in next section and the whole package is available at Mendeley Data repository [31].The propensity functions, used for modeling the probabilities of transitions between states and for the modeling of the time sequence of random moments for every transition, are the same as in the chemical master equation (32). Here, the mean time between transitions of the molecules of each species in the mixture is a random variable with the exponential distribution whose parameters are related to propensity function, i.e., sorption rate of the relevant transition (adsorption or desorption of the first or second component in the mixture). The exponential distribution is designed using in-built function for uniform distribution in MATLAB and inverse CDF sampling technique method [32]. After generating the sequences of time and the instantaneous numbers of adsorbed molecules for each component in the mixture, it is easy to obtain the time evolutions of the number of adsorbed molecules for every component in the mixture, the instantaneous changes in the phase space, and the correlation functions.Here we give results obtained by the simulation algorithm valid for binary adsorption of the mixture of Sarin and Sulphur Mustard on graphene surface. Table1 lists parameters of the mixture. The sorption rate constants are calculated by the procedure from [24], and the desorption energy and the surface density of adsorbed particles, needed for the calculations of the sorption rates, are taken from the online Harvard Dataverse repository [25] and from [26, 27]. Assumed volume of the reaction chamber is 1 dm3, temperature is 282 K, and the number of adsorption centers on the surface, NM, is 38.17.1014 (the surface density of graphene atoms being 38.17.1018 and surface area being 1 cm2). Figure 1(a) shows phase space assuming that, at the initial moment, the surface was clean with no adsorbed amount. Figure 1(b) shows the corresponding cross-correlation function. Due to very high desorption rate constants (Table 1) and high number of adsorption sites on the surface, 1 million transitions were not enough to show the reaching of equilibrium. The corresponding correlation function is symmetric.Table 1
Parameters used in simulations of binary adsorption on graphene in Figure1.
Binary mixture component pressure [Pa] k a [molecules/s] k d [1/s] References Sarin 46 2.89.10−14 1.7.109 [24–26] Sulphur Mustard 360 1.17.10−17 1.57.108 [24, 25, 27]Figure 1
Phase space (a) and the corresponding cross-correlation function (b) of numbers of adsorbed molecules of Sarin and Sulphur Mustard on graphene assuming clean surface in the beginning.
(a) (b)Figure2 shows simulations of adsorption of fictive binary mixture assuming the number of adsorption centers on the surface, NM, is 38.17.105, adsorption rate constant and desorption rate constant of the first component 2.89 and 1.7, respectively, and adsorption rate constant and desorption rate constant of the second component 1.17 and 1.57, respectively, assuming the initial condition is the clean surface with no adsorbed molecules. Other parameters are the same as for Figure 1. On the left, up, and down, on time evolutions of numbers of adsorbed molecules, we see either that reaching the equilibrium was fast and that there was just one transient or that we need more than 10 million transitions in simulations in order to see vanishing of the other transient. In order to see exact equilibrium levels, we would have to check with the analytical solutions. On the right, we see that every simulation led to the same node and that the cross-correlation functions were symmetrical.Figure 2
Time evolutions of the number of adsorbed molecules (left), phase space (up right), and the corresponding cross-correlation function (down right), fictive mixture, clean surface at the beginning.Figure3 shows results for the same fictive mixture as Figure 2, with the only difference being the starting moment of observation, which is now random. Phase space, shown on the upper diagram on the right side, shows that all different starting points are leading to the same node. The diagram on the lower right part in Figure 3 shows that the cross-correlation function is not ideally symmetric in this case. Time evolutions in Figure 3 (top and bottom diagram on the left) demonstrate now two changes in the process kinetics. Due to our analysis we now know that there are no two consecutive steady states, and change in kinetics is the consequence of the vanishing of multiple transients (for binary mixture the exact solution would give at least two transients). Exact solution to multicomponent adsorption modeled with linear Lagergren model, published in article [28], states that in linear systems the number of transients equals the number of components in the mixture. The exact solution to nonlinear systems would reveal that there are more transients than there are components in the mixture but they may not be always visible in the observed time and amplitude span. The quality of results also depends on the quality and features of used software.Figure 3
Time evolutions of the number of adsorbed molecules (left), phase space (up right), and the corresponding cross-correlation function (down right), fictive mixture, random initial moment.


## 3.5. Algorithms and ADmoND Software for Adsorption Analysis

The software package ADmoND we used is aimed at the analysis of monolayer adsorption and desorption processes in novel devices. It has been developed for the software platform MathWorks MATLAB, 2013a release. ADmoND can be used as a MATLAB toolbox or MATLAB application with graphical user interface (GUI) that is adapted MATLAB GUI for figure handling, extended with additional tabs (System, Evolutions, and Equilibrium). All calculations are based upon analytically derived expressions, implemented in mathematical forms optimized to avoid or minimize numerical error propagation. Theoretical background is based on well-known laws or new original solutions published in peer-reviewed journals. The full software package is openly available and can be downloaded from Mendeley Data repository [31].The 'System' tab is aimed at handling parameters saved outside MATLAB and interactive setting of system parameters (temperature, adsorbent area, reaction chamber volume, and composition of an adsorbate mixture). Figure4, generated with ADmoND, illustrates adsorption of binary gas mixture on a solid surface.Figure 4
Competitive adsorption between particles of two species for equivalent adsorption sites on homogeneous surface.The 'Evolutions' tab of ADmoND is dedicated to calculations related to adsorption kinetics in time domain. Figure5 illustrates possible time evolutions of monocomponent adsorption (green lines) obtained by stochastic simulation algorithm (SSA, inspired by [9], adapted for ADmoND case studies), their mean (magenta line), and the ideal deterministic solution (red dashed line).Figure 5
Time evolution of monocomponent adsorption obtained by stochastic simulation algorithm (SSA, green lines), mean of SSA ensemble (magenta line), and deterministic solution (red dashed line).The 'Equilibrium' tab is dedicated to calculations related to adsorption dynamics in equilibrium. As we saw, for every mixture ofr components, there arer+1 possible sets of r numerical values,r+1 possible steady states. Each possible steady state is a vector ofr numbers, where each number corresponds to the number of adsorbed molecules for every component in the mixture (Nis). ADmoND findsr+1 roots of the polynomial of ther+1 order (candidates for the steady state of component 1) then computes all other possible steady states for all other components in the mixture. Afterwards, it performs inspection if exactly one set of solutions for the number of adsorbed particles of mixture components in steady stateNs is realistic and satisfies three conditions. The first condition is that the steady state values for the number of adsorbed molecules of particular components must be positive numbers, lower thanNM. The second condition is that the sum of all steady state values for the number of adsorbed molecules of particular components on the surface may never exceed NM. The third condition is that the steady state values for the number of adsorbed molecules of each component in the mixture must be positive number, lower than the overall number of its molecules in the system, N0i.The theoretical proof that there is just one steady state, the stable node, ensures that there is just one set of solutions and thus enables the optimization of the code (the first set that is found is the only set; there is no need for inspection of all other sets).


## 3.6. Discussion of Results

The stability analysis, described above, helps in a better understanding of adsorption dynamics in equilibrium and a more accurate interpretation of results. In research span of almost a decade, in simulations of multicomponent adsorption, we did not encounter a case of practically feasible multiple steady states obtained by MATLAB; only one set of solutions satisfied required conditions (steady state adsorbate molecules on the surface are positive numbers with overall sum that is smaller than or equal to the maximal possible number of adsorbed particles,NM). After thorough investigation, some complex solutions that appeared once proved to be the consequence of the numerical error propagation, which resulted in our additional efforts towards code improvements.All of the methods we applied here for stability analysis in this paper result with the same conclusion: There is no combination of system parameters that ensures the appearance of sustained oscillations. On the contrary, the steady state is the stable node. Thus, if there is a seemingly oscillatory response in the experiment, the cause of oscillations is not the adsorption process (modeled with the nonlinear set of Riccati equations). Our conclusion is that the possible explanations are as follows:(i)
There are no oscillations; noise is mimicking them.(ii)
The reactions are more complex and are not well modeled by the set of Riccati equations.(iii)
The oscillations emerged due to some other cause in the control loop of complex systems used in chemical, petrochemical, oil-refining, pharmaceutical, and related plants.An example that supports this reasoning is presented in the paper [23], where self-oscillating adsorption–desorption of silver nanoparticles on the gel formed by incorporating graphene oxide into a poly (methacrylic acid)–polyethylene glycol copolymer film is reported. Here, adsorbent is a pH-responsive hydrogel which acts as a pH oscillator in a closed reaction and the oscillations in reversible adsorption–desorption of Ag nanoparticles on the patterned hydrogel surface are just response to pH oscillations.


## 4. Conclusions

The presented results prove that a binary adsorption process, moved from the equilibrium, goes through different transition states while reaching a new equilibrium, without experiencing oscillations. It is valid for all processes governed by the analyzed equations or their mathematical equivalents, i.e., for adsorption of gases but also for biological antibody-antigen kinetics or adsorption of macromolecules observed in commercial surface plasmon resonance-based instruments. If we extend this reasoning to mixtures with more than two components, we see that this result is also valid. It is valid whenever we can divide a mixture into a target gas and a fictive gas with user defined parameters (analogously to “user defined gas” option in commercial binary gas analyzers [30]), if the before mentioned assumptions are valid (monolayer adsorption, noninteracting species, and positive real sorption rates).This result can be used for better understanding of seemingly oscillating experimental values, for better judgment of experimental setup, or for better control of adsorption systems far from equilibrium. Further efforts in research and development of adsorption based devices useful for operation of commercial surface plasmon resonance instruments or chemical plants in petrochemical, oil-refining, paper-producing, pharmaceutical, and related industries might be directed to the field of analysis of noisy response.



---

*Source: 7687643-2019-01-01.xml*