Advances in Mathematical Physics

Volume 2019, Article ID 7687643, 12 pages

https://doi.org/10.1155/2019/7687643

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

^{1}Centre of Microelectronic Technologies, Institute of Chemistry, Technology and Metallurgy, University of Belgrade, Njegoševa 12, Belgrade, Serbia^{2}Faculty of Physical Chemistry, University of Belgrade, Studentski trg 12-16, 11000 Belgrade, Serbia

Correspondence should be addressed to Olga M. Jakšić; sr.ca.gb.mthi.sysonan@aglo

Received 7 August 2018; Revised 25 November 2018; Accepted 16 December 2018; Published 1 January 2019

Academic Editor: Carlo Bianca

Copyright © 2019 Olga M. Jakšić et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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.

#### 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 of different adsorbate species in gas/liquid phase competing for free adsorption sites on the same surface, reversibly forming adsorbed particles , without dissociation or mutual interaction, iswhere and 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 . In any given moment, of them are adsorbed whereas are in gaseous phase. There are no interactions between species in the gas phase. Let 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 . Considering all that, the set of deterministic differential equations that correspond to the stoichiometric equation (1) is as follows.These 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 . In the steady state the reaction flux of all relations in (2) equals zero. Equation (2) becomes a set of interconnected quadratic equations where each of the steady state numbers can be written in terms of the first.After replacing all relations from (3) into the first equation from (2), keeping the left hand side equal to zero, the steady state value appears to be a solution of a polynomial of the order of .Each solution to the polynomial , 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 adsorbate species may exhibit* r*+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 , 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.Here, just like in general solution (4a) and (4b), the steady states, and , 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.We seek the solutions in the form of positive real numbers smaller than or equal to the maximal possible number of the adsorbed particles, , 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.The functions and 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. is the column matrix of the elements* N*_{1} and* N*_{2}, and the expressions that correspond to the elements of the Jacobian are as follows.The solution to (9a) and (9b) implies an exponential regression to/from steady state or equilibrium.Here, are constants related to the initial perturbation from the steady state, while the time constants written in terms of the Jacobian are as follows.Bearing 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 , and of the discriminant in expression for lambda, we investigate the signs of , , , and first.

If we employ relations valid for steady states, we transform the expressions for and into the following.In 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 , soand, hence, the Jacobian trace is always negative, and 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 (the eigenvalues of the Jacobian matrix) are always real negative numbers. That means that the steady state is a** stable node **and 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 with the sum of its steady state value and the instantaneous deviation from the steady state . As before, the indices refer to a particular species in the binary mixture.We 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 and whose variances will be addressed later in this section. By doing so we linearized the equation set (15a) and (15b). The constants in (16a) and (16b) are defined by the following expressions.In 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 .Now eigenvalues are as follows.Bearing in mind that expressions in (18) refer to the sorption rate constants, the overall numbers of molecules of reactant species* N*_{01} and* N*_{02}, and the overall number of adsorption centers on the surface that are all real positive numbers, we see that the trace of the matrix 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.Although 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 and and the number of free adsorption sites must equal the overall number of adsorption sites on the surface, . Additionally, for each of both species the initial number of molecules in the system,* N*_{0}, must equal the sum of the number of its adsorbed molecules and the number of the molecules in the gas phase .

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.All 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 arewhereIn 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. The matrix relates and with new coordinates.After 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 remains, then the difference of cross-correlations of the deviations to the steady state is given by the following.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.The reaction flux, related to the original equation set (5a) and (5b), before linearization, vanishes at equilibrium and equalsAlthough 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, , is short enough, so that only one transition between the neighboring states may occur. That means that the system can reach the state with* N*_{1} adsorbed molecules of the first species and* N*_{2} 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 with* N*_{1} adsorbed molecules of the first gas species and* N*_{2} adsorbed molecules of the second gas species in the moment is the sum of these five probabilities.The probability that no transition occurs is as follows.The transition probabilities among neighboring states are proportional to the instantaneous sorption rates, soIn the short time limit, after rearranging terms in (29), we obtain the differential equation for the probabilities, i.e., the chemical master equation.In 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. Its first derivative with respect to , calculated for* s*_{1} = 1 and* s*_{2} = 1, is then used for the calculations of the first moment of the number of adsorbed molecules for the first gas species.Its first derivative with respect to* s*_{2}, calculated for* s*_{1} = 1 and* s*_{2} = 1, equals the first moment of the number of adsorbed molecules for the second gas species.Its second derivative with respect to* s*_{1} and* s*_{2} is as follows.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 and then by 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*(*s*_{1},*s*_{2},*t*). Once solved, the partial differential equation on* F*(*s*_{1},*s*_{2},*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. Table 1 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 dm^{3}, temperature is 282 K, and the number of adsorption centers on the surface, , is 38.17^{.}10^{14} (the surface density of graphene atoms being 38.17^{.}10^{18} and surface area being 1 cm^{2}). 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.