Frontiers reaches 6.4 on Journal Impact Factors

Original Research ARTICLE

Front. Phys., 22 January 2014 | https://doi.org/10.3389/fphy.2013.00036

On the relation between the Langmuir and thermodynamic flux equations

  • 1Department of Chemistry, Norwegian University of Science and Technology, Trondheim, Norway
  • 2Process and Energy Laboratory, Delft University of Technology, Leeghwaterstraat, Netherlands
  • 3School of Applied Sciences, Royal Melbourne Institute of Technology, Melbourne, VIC, Australia

Despite its approximative nature, the Langmuir theory has shown to be a very successful approach to describe experimental adsorption isotherms. Langmuir kinetics is based on systems of non-interacting particles that can transfer from the gas phase to the adsorbed phase with a transition flux that depends both on the gas pressure and surface coverage. Recent molecular simulation results suggest, however, that some systems can have isotherms that are apparently Langmuirian while the kinetics are not. This remarkably result seems to question the interpretation of innumerous adsorption experiments. The observed anomalous kinetics were described by thermodynamic rate equations giving exactly the same isotherms. Unidirectional rates, as obtained from mesoscopic non-equilibrium theory, correct for the non-ideality of matter using activities instead of concentrations and seem to suggest that fluxes from phase A to another phase B only depends on the properties of phase A alone. In this article we show, however, that the theories and simulations are actually consistent when the following two points are taken into account. The first point is methodological and related to how one should count crossing events considering the presence of possible correlations. The second point is theoretical and related to the microscopic link between the Langmuir and thermodynamic rate theory. Specifically, we show how to define diffusion and activity coefficients at the border of the gas/solid interface. If both points are taken into account, there is neither a contradiction between both theories, nor with the molecular simulation results.

1. Introduction

The understanding of adsorption and transport phenomena at surfaces is fundamental for optimization of industrial processes. However, despite of more than one century of theoretical research in this field, there are still vigorous debates and conflicting models in this field (see e.g., [1, 2]). Roughly speaking, there are two groups of theories; either based on kinetic or non-equilibrium thermodynamic foundations. Kinetic theories are derived by considering microscopic expressions for the adsorption and desorption fluxes which should be equal under equilibrium conditions (the law of mass action). Henceforth, adsorption isotherm equations can be derived by equating both expressions for the adsorption and desorption fluxes. Famous examples are Langmuir [3] and Brunauer-Emmett-Teller (BET) [4] theory for single layer and multi-layer adsorption, respectively. However, the term “kinetic theory” is somewhat misleading here since Langmuir and BET [5] isotherms can also be derived without considering its kinetics, but by solving the partition functions of the corresponding lattice models using simple Boltzmann statistical physics, often referred to as “statistical thermodynamics.”

The thermodynamic rate theories [69] that we refer to are mainly based on the relation between the non-equilibrium net flux and the gradient of the chemical potential. The use of the chemical potential allows one to go beyond the ideal gas description on which kinetic theories are often based. The net rate of adsorption becomes a function of activities rather than concentrations. However, there are also some hidden approximations in this approach, such as the assumption of local equilibrium. For instance, thermodynamic rate theories assume that the velocities of particles that desorb are rapidly redistributed locally to a Maxwell-Boltzmann distribution while their distance from surface is still very small. Langmuir and BET do not have the same spatial resolution as thermodynamic rate theory since they do not try to describe the time-evolution of the gas-phase densities, but rather treat the gas-phase as a single state. As a result of this, Langmuir and BET can work in some cases where thermodynamic rate theory would fail. For instance, Langmuir or BET can be accurate even if desorbing particles of the real physical system maintain a positive leaving velocity until bouncing to the walls of the container.

Which equation has to be taken is a non-trivial problem [10]. Presently, this choice is mainly based on empirical grounds; the equation, that gives the best fit for the experimental data, is chosen and assumed to give the right physical mechanism of the adsorption process. Seen in this light, the recent paper of Simon et al. [11] is rather worrisome. Molecular dynamics simulations showed that even though adsorption isotherms were accurately modeled using Langmuir, the unidirectional fluxes (adsorption and desorption rates) showed a completely different behavior. The fluxes were then modeled using thermodynamic rate equations and the activity coefficient for the adsorbed phase was constructed in order to agree both with the Langmuir types of isotherms and the observed dynamical behavior. Although based on a fundamental physical theory, the resulting behavior is somewhat difficult to understand by physical intuition. In the resulting set of equations, the adsorption flux is proportional to the gas pressure, but no longer to the fraction of vacant sites as it is in the Langmuir model. In addition, the desorption flux is proportional to the surface coverage just like in Langmuir's expression, but also inversely proportional to the free space. In this article we reanalyze these data and the theoretical approach by making a microscopic link between Langmuir and thermodynamic rate theories. An important part of the analysis is the way in which the separate adsorption and desorption fluxes are calculated. In general, this calculation is based on the counting of how many molecules have crossed a certain border, the so-called crossing events. These events can be correlated, i.e., caused by an immediate crossing back and forth by the same molecule, or uncorrelated i.e., caused by the statistical probability. We will show that the Langmuir model, the thermodynamic rate equations, and the simulation data are in fact consistent, if we disregard correlations in our counting of crossing events and, in addition, use the proper mathematical expression for the diffusion and activity coefficient at the interface between two phases.

2. Langmuir vs. Thermodynamic Rate Equations

We consider the general case of adsorption from a gas to a surface. On the molecular level, the transition for a molecule m can be expressed as a unimolecular rate equation as

m(gas)m(ads.)(1)

The Langmuir model assumes that the non-adsorbed molecules behave as an ideal gas, while the adsorbent consists of energetically equivalent adsorption sites that can be occupied by only one molecule at the time. Hence, the molecules do not interact with each other in neither phase except for an infinite hard repulsion once two molecules try to occupy the same vacant site. Under these assumptions the forward rate, or adsorption flux Ja, is proportional to the gas pressure P and the fraction of vacant sites on the surface (1 − θ), where θ is the surface coverage or the ratio between occupied sites and the total number of adsorption sites. The backward rate, or desorption flux Jd, is simply proportional to the surface coverage.

Ja=kaP(1θ),Jd=kdθ(Langmuir)(2)

Setting in equilibrium Ja = Jd and bringing θ to the left side gives the Langmuir isotherm

θ=(kakd)P(kakd)P+1(3)

It is important to realize that there is, mathematically speaking, an infinite number of alternative sets for the equations of Ja and Jd (Equation 2) that give exactly the same equilibrium statistics, Equation (3). In fact, for any invertible function G the following flux equations Ja = G(kaP (1−θ)), Jd = G(kdθ) will yield the exact same Langmuir isotherm in equilibrium. Henceforth, if experiments show Langmuirian isotherms this does not necessarily imply that the underlying kinetics are Langmuirian as well. This is exactly the conclusion of hydrogen adsorption work of Simon et al. [11].

In Simon et al. [11] the adsorption of hydrogen was analyzed using Langmuir and thermodynamic rate theory. The isotherms showed Langmuirian behavior, but the adsorption flux (and therefore also the desorption flux) continued to increase linearly with pressure. The simulation data were then modeled using the thermodynamic rate equation in which the adsorption and desorption fluxes are considered to be proportional to the activities of the corresponding phases. The gas phase was still considered to be ideal, so the gas activity was proportional to the pressure agasP which satisfies the observed flux behavior. Note that in this case the adsorption flux is independent of the fraction of the vacant sites on the surface. In contrast, the adsorbed phase was considered non-ideal and the activity of the adsorbed phase can be written as aads. ∝ γθ where γ is the activity coefficient that corrects for the non-ideal behavior of the adsorbed phase. The desorption flux is Jd = kdaads. Equating Ja = JdkaP = kdγθ shows that if γ ∝ 1/(1 − θ) the exact same Langmuir isotherm, Equation (3), is re-obtained. The thermodynamic adsorption and desorption fluxes in that case equal

Ja=kaP,   Jd=kdθ1θ     (4)

                                                       (Simon et al. [11])


This set of equations suggests a completely different mechanism than the Langmuir model despite having exactly the same isotherm. The adsorption kinetics were interpreted as the consequence of a mobile adsorbed phase in which molecules are not trapped on specific adsorption sites, but can even move away from the surface by several angstroms. Therefore, this kinetics resembles the adsorption of a gas into an ideal liquid. In such a case the unidirectional flux from the gas to the absorbed state can still increase linearly with pressure even beyond the saturation point due to gas molecules not requiring vacant sites, but adsorbing anywhere at the liquid surface and getting emerged into the liquid. Simultaneously, the liquid has to release another solute molecule to the gas phase since it can not maintain a number of solute molecules that is larger than its saturation limit. The desorption rate was explained by the non-ideal factor γ that enhances the desorption due to the increase of repulsion between adsorbed molecules when the density at the surface increases.

The non-equilibrium thermodynamic approach has another starting point than kinetic theories which is the phenomenological flux-force relations for determining the expression for the net flux through the system [12]. One of the flux-force relations expresses the net flux in a system as a function of the gradient of a chemical potential [13, 14]

J(z)=l(z)Tdμ(z)dz(5)

Here z is the reaction coordinate or, in case of adsorption, the Cartesian coordinate that is orthogonal to the surface of the adsorbent. T is the temperature and μ(z) is the chemical potential along the coordinate z. In equilibrium, the chemical potential becomes a constant and the net flux vanishes. The coefficient l in this flux-force relation is proportional to the diffusion coefficient D and concentration C, l = DC(z)/kB where kB is the Boltzmann constant, so that

J(z)=βDC(z)dμ(z)dz,(6)

with β =1/kBT. This equation is often referred to as the Kramer's or Teorell's diffusion equation. In the case of an ideal gas μ(z) = kBT ln C(z) + const. so that dμ(z)/dz = (kBT/C)dC(z)/dz and the above equation transforms into Fick's law: J = −DdC/dz. The expression 6 that contains the gradient in chemical potential rather than the gradient in concentration is an important improvement compared to Fick's law in order to describe transport in realistic systems. The equation can further be generalized by taking the diffusion coefficient to be position dependent: D = D(z). In section 5 we give a derivation of Equation (6) based on microscopic principles which also sheds light on how D(z) can be defined in the case of non-uniform density or a position dependent external potential. In section 5 this will be used to derive an expression for this diffusion and activity coefficient at the solid gas interface, which is crucial for a correct interpretation of surface adsorption in terms of thermodynamic rate theory.

Now, the equation for the unidirectional fluxes follows from the application of the above equation, Equation (6), to describe the diffusion over a chemical barrier. In fact, it can be shown, assuming quasi-stationary state behavior, that the net rate J is identical to the difference of two flux-terms that are proportional to the activities of the states at either side of the barrier.

J=kfaRkraP,with kf=Kkr(7)

Here, kf and kr the forward and reverse rate constants and aR and aP the activities of the reactant and product state, respectively. K is the generalized equilibrium constant. For a reaction with n chemical species (reactants and products) this constant is given by

K=kfkb=exp(i=1nβνiμi°)(8)

where μ°i is the standard chemical potential of species i and νi the stoichiometric coefficient. Hence, a unimolecular reaction like Equation (1) results into a generalized equilibrium constant that is an exponential of the difference between the reference chemical potentials μ°R and μ°P of the reactant and product state, respectively

K=exp(β[μR°μP°])(9)

We come back to Equations (6–9) in section 4 where we derive them from microscopic principles. Since the rate or net flux of Equation (7) consists of two terms with opposite sign, it is logical to assume that these correspond to the forward and backward flux. In other words, the thermodynamic adsorption and desorption fluxes are

Ja=kfagas,Jd=kraads.(10)

Henceforth, the thermodynamic rate theory seems to assume that the adsorption rate only depends on the properties of the gas, while the desorption rate only depends on the properties of the adsorbent. However, there are a few points that we need to consider here. The use of chemical potentials and activities can account for the non-ideality of the phases and is, therefore, an improvement on simple idealistic systems like Langmuir. Yet, also the non-equilibrium thermodynamic approach bears some hidden assumptions and approximations whose validity needs to be evaluated carefully. In addition, also the separation of the net flux Equation (7) into adsorption and desorption rate is not a trivial step. While the net flux is always well defined, this is not always true for the individual forward- and backward fluxes. Specifically, if we include or exclude correlated crossings, this will not change the difference, the net flux, but can give huge changes to the forward and backward fluxes. Finally, we will also have a closer look at the microscopic expressions for kf and kr in the forthcoming section 4 to analyze the conditions under which they can really be viewed as plain constants that do not change during the process of a reaction.

3. Chemical Potentials, Activities and Activity Coefficients for the Langmuir Model

In order to compare the Langmuir and thermodynamic rate expressions we will treat an ideal Langmuir model using the concepts of equilibrium statistical thermodynamics. A first step in this approach is to calculate the chemical potentials of the two phases in this model and the corresponding activities and activity coefficients. In specific, the activity coefficient of the adsorbed phase is of interest since the phenomenological analysis of the hydrogen adsorption study [11] suggested that it is proportional to 1/(1 − θ). Here, we will do this analysis analytically for an ideal adsorbent. For this purpose, we will use the microscopic expression of the chemical potential which for a system at constant temperature T, volume V, and number of particles N is the following

μ=kBTln[Q(N1,V,T)Q(N,V,T)](11)

Here, Q is the partition function:

Q(N,V,T)=1Λ3NN!dRNexp(βW(RN))(12)

where W is the total potential energy function of the system and Λ the thermal de Broglie wave length and RN the complete set of coordinates of all N molecules in the system. Substitution of this expression results in

μ=kBTln[NΛ3dRN1exp(βW(RN1))dRNexp(βW(RN))](13)

The definition of activity is often given by the implicit relation:

μ=μ°+kBTlna(14)

where μ° is the chemical potential of a reference state. This implies that the activity can be written as

a=vNdRN1eβW(RN1)dRNeβW(RN)(15)

where the dependency on μ° and Λ have been adsorbed into the constant v* = exp(−βμ*) Λ3 that has the dimension of volume. For the ideal gas we can simply set W = 0 so that Q(N, V, T) = VN/N! Λ3N. Hence, μ = kBT ln [CΛ3], a = v*C, and γ = a/C = v*, where C = N/V is the concentration.

In the adsorbed layer the particles are assumed not to interact with each other, but only with the adsorbent in which each adsorbed molecule gives a contribution of −|ua| to the total potential energy. The adsorbed molecules can occupy vacant sites only and they have no interaction with each other except that they cannot occupy the same vacancy. Let us assume that the volume of the vacancy pocket equals Vvac and M is the number of vacant sites. Then the integration is replaced by counting the number of states that are not doubly occupied. There are M vacancies for the first molecule to choose from, but once this one is fixed, the second molecule has only M − 1 sites left, etc. which results in:

dRN1eβW(RN1)dRNeβW(RN)      =eβ(N1)(|ua|)eβN(|ua|)VvacN1VvacN                    ×M(M1)(M2)(MN+2)M(M1)(M2)(MN+1)                      =eβ|ua|Vvac(MN+1)eβ|ua|Vvac1(MN)(16)

We can substitute the last equation in Equations (13) and (15) and use the following relations θ = N/M = C/Csite with Csite = M/V the concentration of vacancies on the adsorbent. With this, we can list the chemical potentials, activities and activity coefficients for the two phases in the ideal Langmuir model:

μgas=kBTln[CΛ3]=kBTln[βPΛ3],agas=v*C=v*βP,γgas=v*μads.=|ua|+kBTln[θ(1θ)Λ3Vvac],aads.=v*eβ|ua|Vvacθ(1θ),γads.=v*eβ|ua|VvacCsite(1θ)(17)

Here, we also used the ideal gas law PV = NkBTC = βP. In equilibrium, μgas = μads. or equivalently agas = aads. which gives back again the Langmuir isotherm

θ=(βVvaceβ|ua|)P(βVvaceβ|ua|)P+1(18)

which is identical to Equation (3) for ka/kd = βVvaceβ|ua|. This shows again that the Langmuir theory is maybe not well categorized as being a kinetic theory since it can be derived from statistical principles as well. In equilibrium, the ideal gas in contact with the ideal adsorbent gives the Langmuir isotherm regardless of the underlying kinetic mechanism. The only restriction is that the dynamics should obey microscopic reversibility and, therefore, conserve the equilibrium distribution. It is also interesting to realize that the γads. ∝ 1/(1 − θ) is exactly the same behavior that Simon et al. [11] found empirically. This behavior has, therefore, little to do with non-ideality or repulsive forces between the adsorbed molecules. The non-overlap condition is sufficient.

Substitution of the activities of Equation (17) into Equation (10) results into the thermodynamic adsorption and desorption fluxes, Equations (4). This creates a paradox; thermodynamic rate theory seems to imply that Langmuir kinetics (Equation 2) is incorrect even if the underlying model is the ideal Langmuir model. This, of course, can't be true and the only explanation is that at least one of the assumptions of thermodynamic rate theory is false for this case. In order to detect this flaw, we will derive Equations (10) from microscopic principles and scrutinize each assumption in the derivation with respect to the ideal Langmuir model.

4. Microscopic Derivation of Thermodynamic Rate Equation

In this section we will derive Equations (10) from microscopic principles. This equation leads to the thermodynamic adsorption fluxes, Equation (4), after substitution of Equations (17). The analysis, carried out in this section, will be used in section 5 in order to show that that ka, kd in Equation (4) are in fact not always constants.

The derivation of Equations (10) is divided in two stages. First, we will derive the general diffusion equation 6. The importance of this step is that we get a microscopic definition for the diffusion coefficient. In particular, the diffusion coefficient at the boundary of two phases is a nontrivial issue. The second step is to use that expression to derive Equations (10) and give explicit expressions for the constants in term of system parameters.

4.1. Diffusion Equation

Equation (6) contains the diffusion coefficient which, for a homogeneous system, is defined as

D=limt(r(t)r(0))22 dt(19)

where d the dimension of the displacement and r are the corresponding coordinates of the molecules for which the self-diffusion is measured. The brackets imply an average over all molecules and initial conditions. As in the case of adsorption, D refers to the one-dimensional diffusion along the coordinate z that is orthogonal to the surface of the adsorbent. We can relate this diffusion to that of a discrete one-dimensional random walk in which at discrete time intervals of Δt there is an equal probability to move Δz to the right or to the left. The diffusion coefficient for such 1D random walk equals

D=12Δz2Δt(20)

To describe diffusion in a real system such as a liquid or gas we consider slabs with a width of Δz and measure the frequency of particle exchanges between the slabs. In this model Δz can be small, though the limit lim Δz → 0 cannot be taken. Δz should be large enough so that transition between the slabs can be consider uncorrelated. Though, for a finite Δz we can match the corresponding Δt that is the typical time interval in which we can expect a transition from one slab to the slab at the left or to the slab at the right for a given molecule. Suppose that we simulate the realistic system using molecular dynamics with a timestep δt « Δt. Then Δt equals δt times the number of atomistic steps that are required on average to expect a transition from one slab to another for a given particle. Moreover, the average number of atomistic steps, needed to observe a transition for a specific particle, must be inversely proportional to the chance that a single molecular dynamics time step creates such a transition. In other words

Δt=δt/{2P(III)}(21)

where P(III) is the chance that a specific particle that is in slab I will move to slab II within δt. Henceforth, the chance to move either to the left or to the right must be the double if we assume that these chances are equal, P(III) = P(III). This explains the factor 2 in Equation (21). Substitution of Equation (21) into Equation (20) gives

D=Δz2δtP(III)(22)

We should note, however, that the diffusion coefficient D is only well defined if the slabs have similar conditions, for instance have the same chemical potential. This condition implies that P(III) = P(III). Though, it is certainly possible to give meaningful definitions for a position dependent D(z) for inhomogeneous systems as we will do in the following.

Also the fluxes from slab I to slab II and the reverse can be expressed by the transition probabilities P(III) and P(III). These fluxes measure the number of crossings between the two slabs per unit of time per surface area. Henceforth, these must be proportional to P(III) or P(III), the chance that a single particle makes the transition within δt, times the number of molecules in the respective slabs or

yes

where yes = Vz is the surface area of the slabs. For the net rate we can write

J=JIIIJIII(24)

So, in order to express both diffusion and fluxes in terms of microscopic properties, we need to develop microscopic relations to express P(III) and P(III).

In general, both the number density and the external potential in each slab might differ. Though, within each slab we assume local equilibrium which implies a fast diffusion in the x and y directions and a relatively fast decay of the velocity autocorrelation function along the z-direction. Let NI and NII be the number of particles in each slab and let WI and WII be the potentials that are felt by these particles in each slab that includes both intra-particle interactions and the effect of a possible external potential. At the microscopic level we can express any possible one-particle III transition as a transition between two microstates i and j that belong to two different classes, X ≡ {NI, NII − 1} and Y ≡ {NI − 1, NII} respectively. The class or macrostate {n, m} refers to the collection of points in phasespace, describing the two slabs as collective system, having n and m particles in each slab. Now, by introducing t(ij) as the dynamical transition probability to make a physical transition from microstate i to j using a δt timestep, we can write that

P(III)=iXP(i|iX)jYt(ij)P(III)=jYP(j|jY)iXt(ji)(25)

Using the Boltzmann statistics that says that the equilibrium probability of state i is proportional e−βUi with Ui the total energy of microstate i and assuming local equilibrium we can express the conditional probability P(i|iX) as the equilibrium probability Peq(i) ∝ e−βUi divided by the probability P(X) ∝ ∑i′ ∈ X e −βUi of macrostate X. Therefore P(i|iX) = e−βUi/ ∑i′ ∈ X e−βUi′, which leads to

P(III)=iXeβUijYt(ij)iXeβUiP(III)=jYeβUjiXt(ji)jYeβUj(26)

Naturally t(ij) will depend on the details of the dynamics of the system, but it should obey microscopic reversibility

   Peq(i)t(ij)=Peq(j)t(ji)eβUit(ij)=eβUjt(ji)(27)

and therefore we can write

iXjYeβUit(ij)=iXjYeβUjt(ji)      (28)

For further use we will define

τNII,UIINI,UI=iXjYeβUit(ij)      =iXjYeβUjt(ji)(29)

which obeys the symmetry relation τNI,UINII,UII. = τNII,UII NI,UI.

Here we should reflect upon the microscopic reversibility equation Equation (27) regarding the issue of whether we should consider States i and j to be defined in phase space or configuration space. Both treatments are valid, but for our purpose configuration space is more convenient. Hence Ui in Equations (25) and (27) should be considered as the potential energy. This implies that t(ij) also includes the probability that configuration state i has the right velocities, both in magnitude and direction, to make the transition ij possible, in addition to possible stochastic parameters. The relation 27 is therefore valid for a wide range of dynamics ranging from deterministic NVE dynamics, thermostated dynamics like Andersen or Nosé-Hoover, Langevin, and Brownian dynamics. However, it is also a valid remark whether we should consider all transitions. Δz has a finite value, but δt can, in principle, be taken to the limit δt → 0. Therefore, only particles close to the edges can make the transition within a single timestep. Since it takes many timesteps for this particle to diffuse to the middle of the new slab, it might well recross the boundary several times before this happens. These types of correlations can be avoided by only counting the so-called effective crossings [15, 16] (See Figure 1].

FIGURE 1
www.frontiersin.org

Figure 1. Illustration of how the counting of diffusive correlations can be avoided using the effective positive flux formalism [15]. The approach is based on the principle that for each trajectory that starts at the middle (point a) of one slab and ends at the middle of neighboring slab (point f) only one crossing should be considered as “effective.” For instance, the last time that it crosses the interface between the slabs (crossing ij). This implies that the crossings at points b, c, d and e are not considered. For the reverse trajectories, that move from right to left, we should then consider the first crossings of these trajectories. The microscopic reversibility condition, Equation (27), is then still valid since the statistical probability of trajectory (abcdeijf) equals the probability of the reversed path (fjiedcba) in which time order and velocities are reversed [15].

Now, we will assume that the total potential energy Ui and Uj of states i and j, that describe both slabs I and II as a whole, can be written as the sum of energies of the individual slabs. In other words, we assume that, for the calculation of the total energy, we can neglect the interaction between particles that are in different slabs. This implies that

iXeβUi=ZI(NI)ZII(NII1)jYeβUj=ZI(NI1)ZII(NII)(30)

where we introduced the partition function integrals ZK(N)

ZK(N)=dRNeβWK(RN)(31)

Hence,

  P(III)=τNII,UIINI,UIZI(NI)ZII(NII1),P(III)=τNII,UIINI,UIZI(NI1)ZII(NII)(32)

We can now substitute these expressions into the earlier derived expressions for the flux and diffusion coefficient:

yes

Though, for a non-uniform system the diffusion coefficient is, in principle, not exactly defined. In fact, the transition probability to jump to the right might differ from to the transition probability to jump the left. We will show that, if we define the diffusion coefficient in the middle of two slabs as the geometric mean of these two possibilities

DI,II=Δz2δtP(III)P(III)(non-uniform) =Δz2δtτNII,UIINI,UIZI(NI)ZII(NII1)ZI(NI1)ZII(NII),(34)

we can reobtain the diffusion equation, Equation (6). For this purpose, we also need to define the activity coefficient at the interface. Using the partition function integral 15, we get simple expressions for the activities, Equation (15), and activity coefficients

aK=v*NKZK(NK1)ZK(NK),γK=aK/C=v*VZK(NK1)ZK(NK)(35)

Now, similar to Equation (34), we will also use the geometric mean to define the activity coefficient of the intermediate region between two slabs

γI,II=γIγII=v*VZI(NI1)ZII(NII1)ZI(NI)ZII(NII)(36)

With these definitions at hand, we can derive a neat expression for the net rate

yes

Hence, we reobtained the general diffusion equation Equation (6) from microscopic principles. In the last line we used βμ = ln a + const. ⇒ Cβdμ/dz = (C/a)da/dz = (1/γ)da/dz. In retrospect, this also justifies the use of the geometric mean instead of a normal average to express the diffusion coefficient and activity coefficient at the interface of two slabs. We will use these definitions, Equations. (34) and (36), in the following two sections to analyze the possible density dependence of ka and kd of Equation (4).

4.2. Activity Based Flux Expression

Based on the diffusion equation 6 we can derive the activity based flux expression 10. The derivation given here is slightly different from the one given in Kjelstrup et al [13]. Instead of using expression 6 we use the one-but-last expression of 37 that is equivalent

J(z,t)=D(z,t)γ(z,t)da(z,t)dz(38)

Using the steady state approximation, we can assume that after a short initialization the flux no longer depends on position.

for t>tinit and zbarrier region:  J(zt) = J(t)(39)

Now, we can divide Equation (38) by D/γ and integrate both expressions at either side of the equals sign from z = zR to z = zP with zR and zP being the reactant state and product state values for the reaction coordinate.

zRzPJ(t)γ(z,t)D(z,t)dz=zRzPda(z,t)dzdzJ(t)=[zRzP(γ(z,t)D(z,t))dz]1(a(zR,t)a(zP,t))(40)

The largest contribution to the integral is generally around the transition state. Since the transition state is generally a state of low density, in many cases the ideal gas law applies to this region which gives an activity coefficient γ(zTS) ∝ eβUTS with UTS the potential barrier that each molecule acquires when it is at the transition state [See also Equation (A17) of the appendix]. In addition, the fact that the transition state has a low density during the complete course of the reaction, also little variations in the diffusion coefficient are expected: D(zTS,t) ≈ DTS. Under these circumstances the integral becomes independent of time

zRzP(γ(z,t)D(z,t))dzA(eβUTSDTS)(41)

where A is a frequency factor related to the curvature of the potential barrier at the transition state. We will come back to Equation (41) later on. As we will show, it is exactly this approximation that fails in the case of surface adsorption and, therefore, explains the contradiction between Equation (2) and 4.

Finally, we will have to take care for the fact that activities of different phases are generally defined using different reference states. Hence, instead of v* = exp(−βμ°)Λ3 we should use the appropriate values vR = exp(−βμ°R3 and vP = exp(−βμ°P3 that are used for the reactant and product state which implies

aR(t)=a(zR,t)v*exp(βμR°Λ3),aP(t)=a(zP,t)v*exp(βμP°Λ3)(42)

which finally gives the result of Equation (7): J = kf aRkr aP, with constants

kf=(v*DTSeβμR°A)eβUTS,kr=(v*DTSeβμP°A)eβUTS,(43)

and the ratio kf/kb corresponds to Equation (9). Splitting the forward and backward part and setting aR = aads. and aP = agas results in Equation (10).

In the appendix we show how this works for simple three-level system. An important lesson from that analysis is that the theoretical forward rate kfaR and backward rates kraP are indeed the uncorrelated fluxes. Henceforth, despite that the mathematical derivation involves quite a lot operations with the sole aim to get expressions for the net rate rather than the separate forward or back rates, the theoretical expressions for Jf and Jr should, in principle, reflect the same fluxes as the ones considered in kinetic models. However, the only microscopic term that is really of dynamical nature are the t(ij) transition probabilities and everything that depends on these. Therefore, all dynamical properties in the thermodynamic rate equation are somehow adsorbed in the diffusion coefficient D. As we will see, the apparent contradiction between kinetic and thermodynamic rate theory is due to the fact that kf and kb are in fact not constants but depending on θ because the the approximation given in Equation (41) is invalid for adsorption processes.

5. Diffusion and Activity Coefficients at the Adsorbent-Gas Interface

As stated in the previous section, the key assumption that determines the validity of the thermodynamic rate equations (Equation 10) is the approximation of Equation (41). In order to address the validity of Equation (41) for surface adsorption, we need to determine the activity and diffusion coefficient in the whole system, i.e., also on the boundary of the gas-solid interface. For this, we will use the definitions given by Equations (34) and (36).

First, we need to derive an expression for τNads., Uads. Ngas,Ugas that is given by Equation (29)

τNgas,UgasNads.,Uads.=iXjYeβUit(ij)            =iXjYeβUjt(ji)(44)

where X and Y are descriptions of a state consisting of two layers; an adsorbent layer in contact with a layer of gas, X having {Nads., Ngas − 1} and Y having {Nads. − 1, Ngas} molecules in either phase. Now, it is most convenient to use the first expression that includes t(ij) and assume that for each iX there is a fixed number of states jY with non-zero transition probability t(ij) with j being a state that is almost identical to state i except for the z coordinate of the target particle. Then, by stating that ∑j t(ij) is a constant equal to k′δt we get

τNgas,UgasNads.,Uads.=iXeβUikδt             =Zads.(Nads.)Zgas(Ngas1)kδt(45)

and for the diffusion coefficient

Dgas,ads.  =Δz2kZads.(Nads.)Zgas(Ngas1)Zads.(Nads.1)Zgas(Ngas)          =Δz2keβ|ua|Vvac(MvacNads.)1V          =Δz2keβ|ua|CsiteVvac(1θ)(46)

where we used again Equation (16), θ = Nads./Mvac, and Csite = Mvac/V. As a result, we have obtained the interesting conclusion that the diffusion constant depends on the coverage as ∝ 1θ.

The activity coefficient at the interface follows directly from γgas,ads.=γgasγads. and Equations (17)

γgas,ads.=v*eβ|ua|VvacCsite(1θ)(47)

and hence

Dgas,ads.γgas,ads.=Δz2kv*eβ|ua|VvacCsite(1θ)(48)

Let us now consider Equation (40) for the adsorption process

J/(agasaads.)=[zgaszads.(γ(z,t)D(z,t))dz ]1[(Dgas,ads.γgas,ads.)1Δz+(Dgasγgas)1|zgaszads.Δz| ]1(49)

This relation is in fact not completely valid since we cannot assume steady state behavior (position independent flux) in the gas phase z ∈ [zads. + Δz : zgas] where the potential is flat. Still it is useful to consider the above equation qualitatively. It shows that for low coverage the system might well behave like the thermodynamic rate equation which implies that the fluxes are proportional to activities. The fluxes mainly depend on the diffusion coefficient in the gas phase. However, when the coverage increases the (Dgas,ads.gas,ads.)−1 always becomes dominant since it grows as 1/(1 − θ).

Since the kinetic models consider simply the fluxes between one phase to the other, we will just consider the fluxes between two layers, one layer of gas and one layer describing the adsorbed state. In other words |zgaszads.| = Δz and therefore

J=(Dgas,ads.γgas,ads.Δz)(agasaads.)(50)

This implies that for surface adsorption kf = kr = (Dgas,ads.gas,ads. Δz) in Equation (10) and these are not constant but depend on the coverage θ. The generalized equilibrium constant K is unity; since the hydrogens do not dissociate or undergo another type of chemical reaction, we use the same reference chemical potential for product and reactant in Equation (9). Substitution of the activities of Equations (17) yields the “constants” of Equation (4) as

ka=Dgas,ads.γgas,ads.Δzv*βkd =Dgas,ads.γgas,ads.Δzv*exp(β|ua|)Vvac(51)

or after substitution of Equation (48)

ka=Δzkeβ|ua|VvacCsite(1θ)βkd=ΔzkCsite(1θ)(52)

This shows that since Dgas,ads.1θ and γgas,ads. ∝ 1/1θ the 1/(1 − θ) dependence of aads. is exactly cancelled while the missing (1 − θ) dependence of the adsorption flux is exactly compensated. In fact, Ja = kaP and Jd = kd θ/(1 − θ) is identical to the Langmuir equations Equations (2), Ja = kaP (1 − θ), Jd = kd θ with

ka=ΔzkCsiteβeβ|ua|Vvackd=ΔzkCsite(53)

Moreover, ka/kd = βVvaceβ|ua| which is in agreement with Equation (18). The notion that the Langmuir equations do not automatically follow from the standard non-equilibrium thermodynamics approach has been made before by Pagonabarraga and Rubí in Pagonabarraga and Rubi [14]. They state that non-linear relations between fluxes and forces are indispensable to derive the Langmuir equation from thermodynamic rate theory. We are the first to derive these relations from microscopic principles and show that this non-linearity is related to the diffusion and activity coefficients at the solid gas interface for which we can give exact expressions.

6. Correlated Fluxes

We derived that the thermodynamic rate equations and kinetic rate equations are consistent if the right expressions for diffusion and activity coefficients are used. This has the implication that the thermodynamic rate constants are in fact not constant for the cases in which the diffusion and activity coefficients depend on the densities. And in fact, this is exactly the case for a system that behaves according to the Langmuir model. However, we still need to explain the results of Simon et al. [11] who found that kinetic data of hydrogen adsorption simulations agreed with thermodynamic rate equations in which fluxes are simply proportional to activities, i.e., with constants being truly constant.

We will reason that measured fluxes contain the reactive flux (which is approximately Langmuirian) but also correlated fluxes. These are basically the particles that approach the adsorbent sufficiently close at a spot that is already occupied by another molecule. Therefore, this particle is repelled back almost instantly into the gas. If we just count all events in which molecules cross a certain threshold distance away from the adsorbent, this process would be considered as adsorption followed by an immediate desorption. Although this is a valid way to measure the net flux (where these processes cancel), we should disregard these correlated events in our theoretical analysis since these do not contribute to the rate of an reaction nor to the diffusion.

Now, we will decompose the fluxes for adsorption and desorption that are just based on the distance and not corrected for correlations Jdist.a and Jdist.d, into two components: the Langmuir-flux and the correlated-flux

Jadist.=JaL+Jac,Jddist.=JdL+Jdc  (54)

Since each correlated adsorption is directly followed by a desorption, we can write Jca = Jcd = Jc. We can assume that these bounce-back events occur with a frequency that is proportional to the fraction of occupied sites and the pressure.

Jc=k^aPθ(55)

Now, if we take k^a = ka we reobtain the expression for Ja in Equation (4)

Jadist.=JaL+Jc=kaP(1θ)+kaPθ=kaP(56)

and for the desorption we can write

Jddist.=JdL+Jc=kdθ+kaPθ,(57)

which still looks rather different from Equation (4). However, since Langmuir statistics is conserved, we can simply invert Equation (3) to get pressure as function of θ:

P=θkakd(1θ)(58)

Substitution of this expression into Equation (57) gives

Jddist.=kdθ+kaθ2kakd(1θ)=kdθ1θ(59)

which exactly agrees with the simulation data of Simon et al. [11] and thermodynamic rate theory, Equation (4), with the (false) assumption that kf and kr are constants.

Often, there is a large separation in timescales which makes it easy to distinguish between correlated or uncorrelated events. However, there are also cases where this is less obvious. In the next section, we will apply a correlated-flux analysis on the simulation data of Simon et al. [11] based on the residence time at the surface.

7. Hydrogen on Graphite

Finally, we come back to the results of Simon et al. [11] and show that the main disagreement between Langmuir kinetics and the fluxes measured in the simulations is due to the counting of correlated data. We repeated the simulations of Simon et al. [11] using the same parameters. 1 ns NVE simulations consisting of 100, 200, 300, 350, 400, 450, 500, 550, 600, 650, and 700 hydrogen molecules were simulated at an average kinetic energy corresponding to temperature of T = 180 K. The isotherm is given in Figure 2 and shows that it is reasonably well described by the Langmuir model where the surface density equals Cads.(P) = Cmaxθ(P) = CmaxαP/(1 + αP) with Cmax being the maximum surface density and α = ka/kd. The fitted Langmuir parameters are Cmax = 0.034 mmol/m2 and α = 0.013. The maximum surface density is somewhat larger than the one reported by Simon et al. [11] (≈ 0.015 at 180 K). This might be the effect of small differences in the simulation set-up. Simon et al. used a smaller simulation box which might lead to larger finite size effects. More significant is probably the flexible graphite surface that was used in Simon et al. [11] whereas we adopted a fixed graphite surface. On the other hand, the pressure range in both studies is still very much inside the linear regime of the adsorption isotherm which make it difficult to fit accurately Cmax (see bottom-right corner of Figure 2]. However, the most important conclusion is that the Langmuir model can be fitted to our data, just like the results of Simon et al. [11].

FIGURE 2
www.frontiersin.org

Figure 2. Isotherm (red) with corresponding Langmuir fit (green). The inset in the top-left corner shows 1/Cads.(P) vs 1/P which should be linear for Langmuirian types of isotherms: (1/Cads.(P)) = (1/Cmaxα) (1/P) + (1/Cmax). Langmuir parameters were obtained from the slope and intersection of the linear fit. The bottom-right corner shows an extrapolation of the fitted curve at a larger pressure domain.

In Figure 3 (red line) we show the flux using the same counting strategy of Simon et al. [11]. This implies that a hydrogen is considered adsorbed when it is within 3 Å of the surface. The flux is, hence, measured by counting the molecules crossing the z = 3 Å line. The desorption and adsorption fluxes are equal since we consider a system in equilibrium in which the net flux is zero, but the individual desorption/adsorption flux increases linearly as function of pressure and faster than linearly as function of coverage θ. As argued in previous section, this behavior might be due to correlations that correspond to events that do not settle but are repelled back to the gas phase. In order to quantify the fraction of correlated crossings we examined the distribution of residence times on the surface after an adsorption event. In the case of ideal two-state kinetics, this distribution should be simple exponential P(t) ∝ exp(−kt) with k being the rate constant. In the presence of correlations, we expect the distribution to be a superposition of the ideal distribution and a distribution that is peaked at the small t. Therefore, we assume that we can model this overall time distribution as

P(t)=q1ek1t+q2ek2t with k1k2(60)
0P(t)dt=1q1=k1k2(k2q2),(61)

where we used the fact that the distribution is normalized. From an accurate estimation of P(t) it is relatively easy to fit k1, k2 and q2 since in a log(P(t)) vs. t two straight lines should be visible; a first line at small t with slope −k1 and a second with slope −k2. The determination of P(t) is technically a little problematic since it requires to set up bins of finite width describing small time-intervals and creating a histogram by counting the number of residence times for each bin. However, taking a too small bin-size results into spikes since many bins will contain no points at all. Too large bins will result in a graph having a huge first block that contain all correlated crossings and is therefore not so informative. One could think of setting up bins of variable length, but our solution is more convenient. Instead of computing P(t) we compute the integral

T(t)10tP(t)dtq1k1ek1t+q2k2ek2t(62)

The construction of T(t) is much less sensitive to the bin-size similar to the construction of a crossing probability [17]. In Figure 4 we constructed the histogram for T(t) using a bin size of 0.2 ps and the plot looks smooth upto at least 60 ps. The parameters k2 and q2 are directly obtained by the slope and intersection of a straight line fitted to the linear part (10 ps < t <30 ps) in the log plot. After subtracting this exponential function q2 exp(−k2t)/k2 from the histogram T(t), the remaining part can be fitted again at the small timescale to obtain k1. In Figure 4 we show the graph T(t) obtained from simulations together with the fit. The insets of Figure 4 show the non-integrated time distribution P(t) both in normal-scale and in log-scale together with the fitted double exponential that was obtained from the integrated curve. The two timescales corresponding to correlated and uncorrelated adsorbing events are clearly visible, though they are not completely non-overlapping. Therefore, taking a fixed time, say 5 ps, and assuming that all crossings that remain a shorter time on the surface should be considered as a correlated crossing, won't work very well. Even a significant part of the uncorrelated distribution q2 exp(−k2t) is on this short timescale. In our approach we get a better separation between correlated and uncorrelated crossings when the timescale k−11 and k−12 are still within one order of magnitude. The fraction of correlations can also be obtained as (k2q2)/k2 which increases more or less linearly with increasing coverage just as expected (see Figure 5]. We can repeat the same analysis on the residence time in the gas phase. Surprisingly, that analysis shows an even stronger separation and a fraction of correlations that is around 80%. The reason for this is that the distance z = 3 Å is still very much inside the potential of attraction. Hence, many particles crossing this line from the adsorbed phase have not enough energy to really escape and quickly fall back to the surface. The fraction of these short time recrossings is therefore more or less constant, independent of coverage and pressure. In fact, it is following an opposite trend with a small decrease in the number of correlations. In our analysis we do not correct for these recrossings. Hence, crossing the z = 3 Å from the lower distance is assumed to be a successful transition to the gas phase. This way of counting makes a consistent story with the Langmuir kinetic model. The calculation of the population correlation function [18, 19] would be a way to account for the multiple correlation events, but that is not the aim of this article. What is important is that k2 corresponds to the frequency that particles coming from the gas phase become long-lived surface particles. If multiplied with the average number of surface particles, this gives the uncorrelated flux that must be the same for desorption and adsorption. Figure 3 also shows the uncorrelated flux that goes to a constant for increasing pressure which is in agreement with Langmuir kinetics. When considered as function of coverage, the flux is increasing a bit less than linear due to a k2 that is somewhat smaller at the higher pressures. This is also not illogical since the gas is not an ideal gas. Particles from it will collide with surface atoms that try to escape from surface but are then kicked back without crossing the z = 3 Å line. Hence, our results suggest that the biggest approximation of the Langmuir model for this system lies in the treatment of the gas phase as an ideal gas rather than the neglect of repulsions between adsorbed molecules.

FIGURE 3
www.frontiersin.org

Figure 3. Absolute flux and the flux corrected for correlations.

FIGURE 4
www.frontiersin.org

Figure 4. Residence time analysis for the hydrogen adsorption system with 700 hydrogen molecules. This system corresponds to a gas pressure of 89 bar. The insets shows the probability distribution of the residence times that the hydrogens spend on the surface. The main panel shows the integrated curve, Equation (62), that was used to fit the double exponential which are, for this pressure, k1 = 1.37 and k2 = 0.106 ps−1.

FIGURE 5
www.frontiersin.org

Figure 5. Fraction of correlations against coverage.

8. Conclusions

We have shown that apparent discrepancies between kinetic models and thermodynamic rate models for the description of surface sorption is related to the microscopic definition for the diffusion and activity coefficient at the interface between the gas phase and adsorbed phase. We show that the diffusion coefficient D and the activity coefficient γ depend on the coverage θ of the adsorbent: D1θ and γ ∝ 1/1θ. If the proper expressions are applied in the thermodynamic rate equations for the adsorption and desorption fluxes, these are no longer linearly proportional to the activities. The “constants” that are the proportionality coefficients between the sorption fluxes and corresponding activities depend on D and γ and are, therefore, also dependent on the coverage θ. The resulting equations are identical to the ones derived from kinetic principles. We also analyzed molecular simulation data that showed Langmuirian isotherms, but with adsorption and desorption fluxes that were proportional to the activities of each phase, in agreement with thermodynamic rate theory with constant D and γ. We have interpreted these results by making a distinction between correlated fluxes and uncorrelated fluxes. A molecule approaching the surface might discover that the site on the surface is already occupied after which it is bounced back into the gas phase. This can be counted as an extra adsorption event followed by an extra desorption event. By eliminating the excess of events which have a too short residence time on the surface, one can correct for these correlated recrossings. The final corrected data are more consistent with Langmuir kinetics. Elimination of correlated recrossings is also essential for mapping the microscopic full-atom system on more mesoscopic or macroscopic descriptions in which the dynamics are fully described by diffusion coefficients and rate constants.

The results presented in this article are useful since they connect two classes of theories that so far have been used by different groups independently, sometimes leading to contradicting results. Our approach shows how thermodynamic rate theory can be used when the gradient of the chemical potential is difficult to evaluate due to discontinuities or when the diffusion coefficient is physically not clearly defined. By doing so, the thermodynamic rate theory becomes identical to Langmuir kinetics if the system behaves as an ideal Langmuir model. In addition, it allows to derive more general flux and rate equations since the use of chemical potentials can correct for the non-ideality of matter. For instance, we showed that the adsorption fluxes are not simply proportional to the activities (Jaagas, Jdaads.) but proportional to a more complex expression: Ja ∝ (Dgas,ads.gas,ads.)agas, Jd ∝ (Dgas,ads.gas,ads.)aads. where Dgas,ads. and γgas,ads. can in principle depend on both surface coverage and gas pressure. This makes it possible to use thermodynamic rate theory for the modeling of transport and chemical reactions in complex systems which are not well described by the Langmuir kinetics.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the Research Council of Norway and the Faculty of Natural Sciences and Technology (NTNU) for support via the institution-based strategic projects (ISP): From molecules to process application.

References

1. Zhdanov VP. Adsorption-desorption kinetics and chemical potential of adsorbed and gas-phase particles. J Chem Phys. (2001) 114:4746–8. doi: 10.1063/1.1349178

CrossRef Full Text

2. Ward CA. Liquid-vapour phase change rates and interfacial entropy production. J Non Equilib Thermodyn. (2002) 27:289–303. doi: 10.1515/JNETDY.2002.017

CrossRef Full Text

3. Langmuir I. The adsorption of gases on plane surfaces of glass, mica and platinum. J Am Chem Soc. (1918) 40:1361–403. doi: 10.1021/ja02242a004

CrossRef Full Text

4. Brunauer S, Emmett PH, Teller E. Adsorption of gases in multimolecular layers. J Am Chem Soc. (1938) 60:309–19. doi: 10.1021/ja01269a023

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

5. Keil T. A statistical derivation of the bet equation. J Chem Phys. (1954) 22:1617–8. doi: 10.1063/1.1740480

CrossRef Full Text

6. Ward CA, Findlay RD, Rizk M. Statistical rate theory of interfacial transport.1. theoretical development. J Chem Phys. (1982) 76:5599–605. doi: 10.1063/1.442865

CrossRef Full Text

7. Kreuzer HJ. Thermal-desorption kinetics. Langmuir (1992) 8:774–81. doi: 10.1021/la00039a009

CrossRef Full Text

8. Payne SH, Kreuzer HJ. Analysis of thermal-desorption data. Surf Sci. 1989 222:404–29. doi: 10.1016/0039-6028(89)90369-5

CrossRef Full Text

9. Kreuzer HJ, Payne SH. Nonequilibrium thermodynamics of a 2-phase adsorbate. Surf Sci. (1988) 198:235–62. doi: 10.1016/0039-6028(88)90482-7

CrossRef Full Text

10. Ho YS, Ng JCY, McKay G. Kinetics of pollutant sorption by biosorbents: review. Sep Purif Methods. (2000) 29:189–232. doi: 10.1081/SPM-100100009

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

11. Simon JM, Haas OE, Kjelstrup S. Adsorption and desorption of H2 on graphite by molecular dynamics simulations. J Phys Chem. (2010) 114:10212–20. doi: 10.1021/jp1011022

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

12. de Groot SR, Mazur P. Non-Equilibrium Thermodynamics. New York: Dover (1984).

13. Kjelstrup S, Bedeaux D, Johannessen E. Non-equilibrium Thermodynamics for Engineers. Science and culture series (Singapore): Physics. World Scientific; (2010). Available online at: http://books.google.no/books?id=U3VaEONIDv0C.

14. Pagonabarraga I, Rubí JM. Derivation of the Langmuir adsorption equation from non-equilibrium thermodynamics. Physica A (1992) 188:553–67. doi: 10.1016/0378-4371(92)90331-J

CrossRef Full Text

15. van Erp TS, Bolhuis PG. Elaborating transition interface sampling methods. J Comput Phys. (2005) 205:157–81. doi: 10.1016/j.jcp.2004.11.003

CrossRef Full Text

16. van Erp TS. Dynamical Rare event simulation techniques for equilibrium and non-equilibrium systems. Adv Chem Phys. (2012) 151:27. doi: 10.1002/9781118309513.ch2

CrossRef Full Text

17. van Erp TS, Moroni D, Bolhuis PG. A Novel path sampling method for the sampling of rate constants. J Chem Phys. (2003) 118:7762–74. doi: 10.1063/1.1562614

CrossRef Full Text

18. Chandler D. Introduction to Modern Statistical Mechanics. New York: Oxford University Press (1987).

19. Frenkel D, Smit B. Understanding Molecular Simulation, 2nd Edn. San Diego, CA: Academic Press (2002).

20. Moroni D, Bolhuis PG, van Erp TS. Rate constants for diffusive processes by partial path sampling. J Chem Phys. (2004) 120:4055–65. doi: 10.1063/1.1644537

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Appendix

To illustrate the thermodynamic rate equation we consider a simple three-level system that is an abstract model for a reaction barrier with an intermediate energy state (See Figure A1). The overall rate constants k13 and k31 are given by [20]

k13=k12k23k21+k23=k12p2+k31=k31k21k21+k23=k32p2(A1)

where p+2 = k23/(k21 + k23) and p2 = k21/(k21 + k23) are the probabilities that, if the system is at the barrier, it will go the right or left side, respectively. The reaction constants are therefore determined by the frequency that the system jumps on the barrier times the probability to fall off at the other side rather than falling back. The fluxes J1 → 3 and J3 → 1 are, henceforth, simply

yes

and the net flux equals

yes
FIGURE A1
www.frontiersin.org

Figure A1. Three-level system: transitions between two low energy states (state 1 and 3) have to pass through an intermediate energy state (state 2) that is higher in energy. kij denotes the rate constant from state i to j. The barrier between state 1 and 2 is the transition state and therefore k21 < k23.

Now, let us use thermodynamic rate theory to derive this expression. From Equation (15) we can write down the partition function integral for each each energy state:

Zi(N)=(VeβEi)N(A4)

and derive the activities and activity coefficients of each level by Equation (35)

ai=NieβEi,γi=VeβEi(A5)

where we have gauged the reference chemical potential μ° = kBT ln(V3) so that v*/V = 1. The detailed balance condition implies that

eβEiki,i+1=eβEi+1ki+1,i(A6)

which is independent to the number of particles Ni, Ni + 1 at each level since we assume non-interacting particles. Moreover, the chance to move to right or left with a time-interval δt is simply given by t(ii ± 1) = ki, i± 1 δt and, therefore, similar to Equation (45)

τi,i+1=Zi(Ni)Zi+1(Ni+11)ki,i+1δt(A7)

Hence, the diffusion coefficients are (Equation 34)

   Di,i+1=Δz2δtZi(Ni)Zi+1(Ni+11)Zi(Ni1)Zi+1(Ni+1)ki,i+1δt      =Δz2ki,i+1ki+1,i(A8)

The activity coefficient at the boundaries are

γi,i+1=Vki+1,iki,i+1eβEi(A9)

So that

yes

Now we can approximate the derivative da/dz

dadz|i,i+1(ai+1ai)Δz=1Δz(eβEi+1Ni+1eβEiNi)      =eβEiΔz(ki+1,iNi+1ki,i+1Ni)ki,i+1)(A11)

Therefore, the net rate between level 1 and 2 is:

yes

This result is not surprising since it is straightforwardly obtained from a kinetic viewpoint. Though, it shows that thermodynamic rate theory is al least equivalent provided that one uses the right definition for the diffusion and activity coefficients. The net rate between level 2 and 3 follows directly by changing the labels in the equation above

yes

Now, the steady state approximation implies setting dN2/dt = 0 in the following equation

dN2dt=k12N1+k32N3k21N2k23N20(A14)

which results in

N2(t)=k12N1(t)+k32N3(t)k21+k23(A15)

Substitution of this expression into Equations (A12) and (A13) shows that the net rates at the barrier are indeed the same for steady state conditions

yes

and from the integral expression we get

yes

where we used eβE1/k12 = eβE2/k21. The sum k−121 + k−123 is mostly determined by the smallest of the two rate constant. For instance if k21 « k32 the thermodynamic rate equations is simply J13 ≈ exp(−βE2)k21/yes. In addition, the highest points of the potential barrier, the transition state, is likely to be found between states 1 and 2. This all agrees with Equation (41). Using again the microscopic reversibility relation, we can express activities as

a1=N1eβE2k21k12,a3=N3eβE2k23k32,(A18)

which gives

yes

The overall net flux J13 has the same form as J12 or J23 and the kinetic theory based equation, Equation (A3). Another way to express the forward flux is J1 → 3 = k12p2+N1/yes = k12 (1 − p2) N1/yes = J1 → 2k12p2N1/yes which is the total flux from state 1 to state 2 minus the correlated flux (crossings that recross the transition state before entering state 3). Hence, it is not that correlated recrossings cancel in the net flux when taking the difference J1 → 3J3 → 1, but the individual forward and backward fluxes, as obtained from thermodynamic rate theory, are already corrected for fast recrossings.

Keywords: Langmuir kinetics, mesoscopic non-equilibrium thermodynamics, interface diffusion, correlated crossing events

Citation: van Erp TS, Trinh T, Kjelstrup S and Glavatskiy KS (2014) On the relation between the Langmuir and thermodynamic flux equations. Front. Physics 1:36. doi: 10.3389/fphy.2013.00036

Received: 03 November 2013; Accepted: 23 December 2013;
Published online: 22 January 2014.

Edited by:

Sotiris Sotiropoulos, Aristotle University of Thessaloniki, Greece

Reviewed by:

Fernando Bresme, Imperial college London, UK
Miguel Rubi, Universitat de Barcelona, Spain
Dick Bedeaux, Norwegian University of Science and Technology, Norway
Margaritis Kostoglou, Aristotle University of Thessaloniki, Greece

Copyright © 2014 van Erp, Trinh, Kjelstrup and Glavatskiy. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Signe Kjelstrup, Department of Chemistry, Norwegian University of Science and Technology, Høgskoleringen 5, 7941 Trondheim, Norway e-mail: signe.kjelstrup@ntnu.no