Rates of adsorption and desorption: Entropic contributions and errors due to mean-field approximations

Adsorption and desorption are processes of interest to surface science and catalysis.^1,2 In this article, we present simulations whose purpose is to examine how the rates of adsorption and desorption depend on the number of particles N_a adsorbed per unit area (the coverage). We are particularly interested in the manner in which the activation entropy affects the adsorption and desorption rate coefficients and also in the validity of the transition state theory (TST) and of several mean-field approximations.

It is common practice to calculate, by using density functional theory (DFT), energy barriers for desorption and estimate the desorption rate by choosing a reasonable pre-exponential. This is useful for qualitative studies, especially when comparing similar systems (for example, the same molecule desorbing or adsorbing on different solid surfaces). To calculate the rate constants more accurately, one needs to compute the activation entropy. The entropy of the transition state is that of an ideal two-dimensional gas and it can be calculated analytically. The entropy of the adsorbates is most often determined by using simple models: either a two-dimensional lattice gas (2DLG) or a two-dimensional ideal gas (2DIG).^3–6 Campbell and co-workers have shown that the measured entropy of adsorption differs significantly from that predicted by these models^7,8 and provided a linear relation between the entropy of the adsorbates and that of a gas. A linear relation to estimate adsorbed-phase entropy has also been suggested by Dumesic and co-workers.^9–11 

A simple analysis (see the supplementary material) shows that the desorption rate constants calculated with these two models can differ from each other by five orders-of-magnitude for the diatomic molecule considered in this work. Since the two calculations use the same activation energy, this difference originates from the way the models estimate the activation entropy. One expects that the lattice model gives an adequate approximation when the barrier to diffusion is much higher than k_BT. The two-dimensional ideal gas model might work when the temperature is high and the barrier to diffusion is lower than k_BT. There is no simple model for intermediate situations. In addition, neither model covers the case when the adsorbed molecules have a tendency to aggregate and form “islands” whose size and existence depend on temperature and coverage.

This situation raises a number of questions. First, is the lattice-gas model correct at low temperatures or the two-dimensional gas model at high temperatures? Second, what can one say about the magnitude of the activation entropy in the case of intermediate temperatures? How does it depend on temperature and coverage? In this article, we answer these questions, by calculating exactly, for a generic model system, the adsorption and desorption rate constant.

The model we use is not meant to simulate a specific system. It studies the adsorption and desorption rate constants for a generic diatomic molecule, interacting with the surface of a generic solid. We use a diatomic molecule, rather than an atom, because the energy transfer from kinetic energy to rotational energy affects the probability that the molecule is trapped at the surface. The parameters in the interaction energies were chosen so that the model is representative of reality. The calculation of the rate constants includes the reaction of all participants: the adsorbed molecules, the incoming (or departing) molecule, and all the atoms of the solid. The infinite nature of the surface is captured by using periodic boundary conditions in two directions. The supercell (the cell repeated periodically) is large enough to prevent artificial size dependence.

We collect here some elementary facts about adsorption and desorption in order to define the terminology and prevent confusion. If the adsorbate is not in equilibrium with a gas, the number of particles per unit area, $Na$⁠, will change with a rate given by

Here, p is the pressure of the gas in contact with the surface, k_d is the desorption rate constant, and k_a is the adsorption rate constant. Equation (1) is empirical and defines the constants k_d and k_a. The equation assumes, correctly, that adsorption and desorption are statistically independent events. The adsorption rate depends on the number of available sites. We include this dependence in the adsorption rate coefficient. Therefore, both k_d and k_a are functions of the coverage $Na$ (the rate constants are not constant).

The following thought experiment clarifies the operational meaning of k_d. Imagine that a surface is held at constant temperature T and it is in contact with a gas having pressure p. The surface is equilibrated with the gas and has the coverage (at equilibrium) $Nae(p,T)$⁠. We can control the equilibrium surface coverage while keeping the temperature T constant, by changing p. After equilibration, we pump out the gas so that the gas pressure is instantly equal to zero and is maintained this way until all atoms have desorbed. Throughout this, the temperature of the surface is held constant. We monitor the time dependence of the coverage $Na(t)$⁠. Since the gas pressure is zero,

The statement that the adsorption rate and desorption rate are uncorrelated means that k_d defined by Eq. (2) is identical to the one present in Eq. (1).

The following “thought” experiment defines the adsorption rate constant. Start with a clean surface in perfect vacuum. Increase the gas pressure suddenly and monitor the time evolution of the coverage $Na(t)$⁠. The adsorption rate constant is defined by

The rate constant k_a defined by Eq. (3) is the same as the one present in Eq. (1).

The rate of desorption is the same whether a gas is present or not. This is not true at very high gas pressures or for desorption from a solid in contact with a liquid.

Desorption takes place with the adsorbate in thermal equilibrium with the surface, and the probability of different states in the system is given by a canonical ensemble. The only quantity that is out of equilibrium, when the gas is pumped out, is surface coverage.

We explain next the procedure used for exact calculations of the adsorption and desorption rate constants and of other quantities that characterize the dynamics of adsorption or desorption processes.

We can distinguish an adsorbed molecule from a gaseous one by using a dividing surface that is parallel to the solid surface and located at a distance z_d from it. A molecule is adsorbed if the position of its center of mass is between the solid surface and the dividing surface. We chose z_d so that a molecule whose center of mass is located on the dividing surface does not interact with the solid or with the molecules adsorbed on it. This choice of the dividing surface ensures that transition state theory (TST) gives the exact value for the rate constant k_d because a desorbing molecule passing through the dividing surface will never recross it. One must, however, keep in mind that this nonrecrossing condition is not satisfied if the pressure of the gas is very high since a desorbing molecule can collide with a gas phase molecule and readsorb. We do not examine this situation here.

The calculation of the rate constant at this level of detail, without making any approximations, allows us to understand the aspects of adsorption and desorption rates that are not accessible to experiments.

Since we calculate k_d at various values of T, we can make an Arrhenius plot. If the plot is linear, it will provide a numerical value for the activation energy. We define this to be the true activation energy and compare it to the values obtained by making various approximations.

It is conventional to assume that the pre-exponential provides the activation entropy and that the activation energy is the energy of the transition state minus the energy of the adsorbates. The exact rate constant calculations we performed allow us to test and understand these assumptions.

During our calculation, we evaluate by Monte Carlo (MC) the energy E^† of $Na$ adsorbed molecules and one molecule on the dividing surface, as well as the energy E_R of $Na+1$ adsorbed molecules. The difference ΔE^† ≡ E^† − E_R defines an activation energy ΔE^†. Is this energy equal to the energy obtained from the Arrhenius plot?

While ΔE^† is an average quantity, each desorbing molecule interacts with the $Na$ adsorbed molecules having specific instantaneous positions on the surface. The energy needed for desorbing a specific molecule, when the rest of the adsorbates have specific positions on the surface, is different from the average energy ΔE^†. This means that the desorption energy fluctuates. To get an indication of the magnitude of these fluctuations, we calculate, by Monte Carlo, the standard deviation of the desorption energy.

The implementation of TST requires the computation of the ratio of two partition functions: one for $Na$ adsorbed molecules and one molecule on the divided surface and the other for $Na+1$ adsorbed molecules. There are several ways in which this ratio can be calculated. We use the multiple-window umbrella-sampling¹² combined with weighted-histogram-analysis method¹³ (WHAM) because this provides the two partition functions and also the potential of mean force (which, in turn, provides a definition of the activation free energy).

We use the calculated activation energy and the activation free energy to calculate an activation entropy. Calculating how this entropy depends on coverage and temperature will allow us to evaluate the adequacy of the 2-dimensional lattice gas and 2-dimensional ideal gas models.

The adsorption rate constant, however, has to be calculated by the methods that take into account the recrossing of the dividing surface. It is possible that TST is a good approximation but we do not know this a priori. We calculate k_a exactly by using the method explained below by describing the steps of the simulation.

1. Start with $Na$ molecules adsorbed on the surface. Use a sufficient number of Monte Carlo moves to generate an equilibrium configuration of these adsorbates and of the atoms in the solid. Memorize the positions of all atoms in the system obtained in the last Monte Carlo move.

2. Place the center of mass of one diatomic molecule at a random location on the dividing surface. To distinguish this diatomic molecule from the ones already adsorbed on the surface, we call it the incoming diatomic.

3. Solve Newton’s equation for all atoms in the system. Use the initial positions for the solid and adsorbate atoms generated by the Monte Carlo simulation. The initial velocities for all atoms are generated by sampling Maxwell distribution. Make sure that the initial velocity of the center-of-mass of the incoming diatomic (starting on the dividing surface) is such that the molecule moves toward the solid surface.

4. The flux of the molecules approaching the surface is $vz$(0)ρ, where ρ is the density of the gas and $vz$(0) is the z-component of the velocity of the center of mass of the incoming diatomic. We use a coordinate system with the z-axis perpendicular to the solid and pointing toward it so that $vz$(0) > 0. Assuming an ideal gas gives ρ = p/k_BT. The incident flux is, therefore,

Multiply this by the function

which differs from zero only if the molecule is adsorbed. Record, for each trajectory, the quantity

The index α labels the result obtained by running one trajectory. The time τ should be long enough to allow one to decide whether the center of mass recrosses the dividing surface (moving toward the gas side) or is stuck to the surface with no chance of getting away on a time scale much shorter than the inverse of the desorption rate constant. Here, $X[zα(τ)]$ is the sticking probability of the αth trajectory.

1. Go back to step 1 to generate a new equilibrium configuration and repeat this procedure N_MD times.

The quantity

is the adsorption rate for the coverage, $Na$⁠, and the temperature, T. The adsorption rate coefficient is

The rate constant defined by Eq. (8) is the same as the rate constant defined by the phenomenological equation (1). This formula is a correlation function between the incident flux and a function that is equal to 1 when z(t) < z_d and equal to 0 otherwise; it is called a flux-position correlation function. The rate constant is given by the asymptotic values of this function when t = τ is larger than the collision time but smaller than the time 1/k_a.

With no extra cost, one can calculate

to have a quantitative measure of the dispersion of the quantity J_α, which is a measure of the effect of fluctuations on the rate. This quantity is of interest if we want to know the chance that any mean-field method calculation will give large errors.

This simulation also provides the information needed for evaluating the TST approximation for k_a. All we have to do is exclude from the sum in Eq. (7) those trajectories in which the incident molecule bounces off the solid surface and recrosses the dividing surface. This is not the most efficient way of performing TST calculations but when we calculate the exact rate constant k_a, the TST is obtained as a side product, by taking into account only trajectories that do not recross the dividing surface. Obviously, this simulation also allows the definition of a sticking coefficient and of its dependence on the energy of the incoming molecules.

This calculation of k_a assumes that the equilibration rate of the adsorbate with the solid is fast enough so that every incident diatomic molecule comes in contact with the equilibrated surface. This means that the equilibration rate is faster than the rate of deposition of the molecules to the surface. It is also assumed that the heat of adsorption is removed so fast, by the thermostat, that it exceeds the rate of heating by adsorption so that the temperature is constant. Both conditions are satisfied in experiments in which the incident flux is low.

This simulation will answer a number of interesting questions. Is the rate constant calculated with transition state theory correct? It is common to describe experiments in terms of a sticking coefficient, which is the fraction of the incident flux that sticks to the surface. This is essentially the average of the quantity, $X[zα(τ)]$⁠. Our simulations allow us to determine how the sticking coefficient depends on the kinetic energy and orientation of the incident molecules, surface temperature, and surface coverage. It is assumed that the sticking coefficient is proportional to the fraction of the empty surface but this is not obvious and needs to be tested.

We have mentioned here only the calculation of the adsorption rate coefficient. However, other quantities that describe the dynamics of adsorption, as well as their fluctuations can be recorded during the simulation. These will be mentioned later.

The remainder of this article is organized as follows: in Sec. IV, we summarize the methods used in this work. In Sec. V, we present the results of rate calculations and entropic contribution to desorption rates. Finally, in Sec. VI, we offer concluding remarks.

We consider the desorption and the adsorption of a generic heteronuclear diatomic molecule, AB, which is adsorbed, without dissociating, on a generic solid surface. The solid is simulated by using a slab that contains 6 atomic layers. We fix the positions of the atoms in the bottom layer of the slab, which is located in the plane z = 0, as shown in Fig. 1. When the molecule hits the slab, it creates a “phonon wave packet” which moves away from the impact site, is reflected by the bottom layer, and travels back to the solid-vacuum interface. Our slab is sufficiently thick so that the time of this round trip is shorter than the time it takes an incoming molecule to be reflected or to settle as an adsorbate. Moreover, the impact energy is distributed over many atoms of the solid, not just on the atoms with which the newly adsorbed molecule is in contact. In other words, the phonon packet produced by the impact does not affect the fate of the impacting molecule since we equilibrate the system prior to sending a molecule in the simulation at the same temperature in each step.

We used Morse potentials¹⁴ 

to model the interaction between the solid atoms and bonded interactions in a diatomic molecule. The Morse potential parameters are given in Table I. For the solid, these parameters yield atomization energies of 3.3 eV, which falls in the range of atomization energies of metals.¹⁵ For a diatomic molecule, using the parameters given in Table I, we get a vibrational frequency of 2157 cm⁻¹ which is close to the fundamental frequency of C-O (2143 cm⁻¹). We note here that although Morse potential parameters are inspired from solid Cu and CO molecules,^16,17 we make no attempt to simulate the behavior of CO molecules on a Cu surface. We use these parameters to give our generic model a behavior similar to that of real systems.

The interaction of the atoms in the solid with the atoms of the diatomic, as well as the interaction between the atoms belonging to two different diatomics, are given by the Lennard-Jones (LJ) potential¹⁸ 

with the parameters given in Table I. The lowest adsorption energy is obtained when the molecule binds to the hollow site. The barrier to diffusion is at a bridging site and its magnitude is 0.2 eV. This was determined by using the climbing image nudged elastic method^19,20 with the implementation previously reported.²¹ The A atom in the diatomic AB interacts more strongly with the atoms of the solid. The Lennard-Jones parameters for two dissimilar, nonbonded atoms of two diatomic molecules were obtained using the Lorentz-Berthelot^22,23 mixing rules. Furthermore, we included dipole-dipole interactions between the diatomic molecules using the function^24,25

where μ_i is the dipole moment of the ith molecule and $pi→^$ is the unit vector along the dipole direction of the ith molecule. In this work, we used μ = 0.1 D which is close to the dipole moment of CO in gas (which is 0.122 D). A cutoff radius of 12 Å was used for all Lennard-Jones and for the dipole-dipole interactions. These potential functions yield a desorption energy of 0.9 eV for one molecule at 0 K.

We begin by reviewing a few concepts from statistical mechanics needed for describing the calculations performed here. The transition state theory provides an exact expression for the desorption rate constant k_d, if the dividing surface is properly chosen. This dividing surface is parallel to the solid surface and its distance z_d from the solid surface must be such that a molecule whose center of mass moves away from the solid surface and crosses the dividing surface will not return to the solid. If such a choice is made, then

Here, Q^† is the partition function of a system of N molecules, one of which has the center of mass on the dividing surface while the remaining N − 1 are adsorbed. Therefore, Q^† is given by

Γ represents the Cartesian coordinates of all atoms in the system. z is the distance of the center of mass of the diatomic from the solid surface.

It is also useful to define the potential of mean force $w$(ξ) through

where ξ is the distance from the solid surface and $w$(ξ) is the free energy of the system when the center of mass of the diatomic is constrained to stay at a distance ξ from the surface. With this notation, Q^† = exp[−$w$(z_d)/k_BT].

The partition function Q_R is defined by

where $X(z−zd)=1$ if z < z_d and $X(z−zd)=0$ if z > z_d. This is the partition function when all, N, molecules are adsorbed.

The probability that the center of mass of a molecule is located at a distance from the surface between ξ and ξ + dξ is

In order to compute P(ξ), we partition the space between the solid surface and the dividing surface into bins defined by planes parallel to the surface. The width, Δ_b, of each bin is much smaller than z_d and is chosen so that we can assume that the free energy is constant within each bin. We use the Monte Carlo method to create a histogram that gives the number n_i of particles whose center of mass is in bin i. If the number of Monte Carlo moves N_MC is sufficiently large, then the probability P(ξ_i) that a molecule in bin i is

where n_i is the total number of configurations for which one molecule is in the ith bin. A straight implementation of the Monte Carlo method does not frequently sample enough positions corresponding to high energy. To overcome this sampling problem, we use a multiple-window umbrella sampling method^12,13 with the bias potential

Here, $Nw$ is the number of windows and ξ_j = ξ₁ + jΔW. The effect of the bias potential is to increase sampling around ξ_j. ΔW and k_j are chosen to achieve a significant overlap between the distributions of adjacent windows. The unbiased distributions are obtained by combining biased distributions from each “umbrella window” by using the maximum likelihood of the multinomial distribution function.^13,26 This method of obtaining free energy is known as the weighted histogram analysis method (WHAM).¹³ If n_ij is the total number of observations in the ith bin and jth window, the unbiased distribution, $Pi$⁠, is

Here, $Nw$ is the total number of windows, N_b is the total number of bins, $Nj=∑i=1Nbnij$⁠, and

Equations (20) and (21) are solved self-consistently, starting with an initial guess of the μ_j. Finally, we obtain the ratio of partition functions by using

where z_d is the position of the dividing surface and the value of $X(ξ−zd)$ is 1 when ξ < z_d and 0 otherwise.

The activation Helmholtz free energy ΔA^† is computed from

where the factor Δ_b is needed to make the argument of the logarithm dimensionless.

The activation energy ΔE^† = E^† − E_R is obtained by two separate Monte Carlo runs: one calculates the mean energy E^†, by restricting the center of mass of one molecule to be located on the dividing surface, and the other calculates E_R, when all molecules are adsorbed. We obtain the activation entropy change ΔS^† by using the thermodynamic relation

If the number of molecules included in the simulation is N, the calculations of averages pertaining to the transition state are performed with one molecule on the dividing surface and N − 1 molecules adsorbed on the surface. All calculations of the properties of A are performed with N adsorbed molecules. Independent simulations with different values of N correspond to different surface coverages.

Desorption was studied at 200 K, 500 K, and 700 K for the coverages θ of 1/50, 7/50, 11/50, and 14/50. We considered the (100) surface of a face-centered close-packed structure with a unit cell parameter of 3.615 Å. Calculations were performed on a supercell containing a 5 × 5 surface unit. The slab was 6 atomic layers thick and we fixed the atoms in the bottom layer of the slab. The slab has a total of 300 atoms with 50 hollow-sites on the surface for adsorption. We applied periodic boundary conditions in all directions according to the minimum image convention.^27,28 A vacuum layer of 40 Å was used to decouple interaction of the slab with its images in the z-direction. All atoms in the slab, except those in the bottom layer, were allowed to move, both in the Monte Carlo and in the molecular dynamics simulation.

For the Monte Carlo part of the computations, we used a fast Mersenne-Twister algorithm based on 128-bit operations for generating random numbers because it produces statistically good equidistribution over long periods.²⁹ Two types of MC moves were attempted: (1) the translation of atoms and (2) the translation of the center of mass of diatomic molecules with fixed orientation. The atomic translations were attempted with 90% probability and molecular translations with 10% probability. In principle, any combination of attempted probabilities could be used, but each will have its own efficiency. We obtained these probabilities by performing MC simulations on a test system for which maximum efficiency was obtained (data not shown here). In these moves, an atom (or molecule) is selected at random and a random displacement between −λ/2 and λ/2 is added to the three Cartesian coordinates (or center-of-mass coordinates). The value for λ was tuned to have an acceptance probability between 0.2 and 0.4. For atomic moves, we used λ = 0.2–0.4 Å, and for the center molecular moves, we used λ = 6.0–7.0 Å. These large values for λ are needed for sampling jumps from one binding site to another to generate new configurations and get the correct entropy of mixing. At 700 K, only atomic moves were used. The moves were accepted or rejected based on the Metropolis algorithm.³⁰ 

The system was equilibrated at the target temperature for 20 000 sweeps before performing production runs. Each sweep is equal to N_T moves, where N_T is the total number of atoms in the system.

For the molecular dynamics part of the computations, we used the velocity-Verlet algorithm³¹ to integrate the Newton’s equations of motion, with a time step of 1.0 fs. Maxwell-Boltzmann velocity distribution was obtained using the Box-Muller method.^32,33 Molecular dynamics trajectories for computing adsorption flux were computed for a total time of 15 ps, which is typically the time scale for collision. Within this time scale, the average sticking coefficient converged to a constant value for all coverages and temperatures. For low coverages and low temperatures, it was sufficient to run the simulations for 10 ps.

For the free energy part of the calculations, we considered 8 different initial guesses picked randomly from a large sample of equilibrated configurations. Each configuration was run for a total of 20 000 sweeps. We collected data at every sweep, thus making the total sample size equal to 160 000. There was no appreciable change in the free-energy surface on increasing the sample size by another 20 000. The umbrella windows were placed at an interval of 0.2–0.3 Å and the k_j for umbrella Gaussian potential was chosen to be 800 kJ/(mol/Ų). These values were obtained by performing benchmarking calculations to obtain a significant overlap between adjacent windows. The self-consistency loop in Eq. (20) was terminated when $∑j=1Nw1−μjnewμjold2<10−8$⁠. Numerical integration to compute the ratio of partition functions was performed using composite Simpson’s rule.³² The width of the bin (Δ_b) for computing probability distributions was 0.01 Å.

The dividing surface was placed at z = 21 Å from the bottom layer of the slab, which is fixed at the plane z = 0. This choice corresponds to nearly 11 Å away from the surface. We say “nearly” because this distance will depend on temperature. However, because the bottom layer is fixed, the distance of the molecule from the bottom layer will not change with the temperature. Adsorption rates were computed using 50 trajectories (with velocities obtained from the Maxwell-Boltzmann velocity distribution at the target temperature) fired from each of the 400 random configurations drawn from a sample of configurations on the dividing surface generated from MC simulations; in all, this corresponds to a total of 20 000 trajectories.

The potential of mean force was calculated from Eq. (15) by the procedure explained in Sec. IV C. Figure 2 shows a plot of the potential of mean force as a function of the distance ξ between the center of mass of the diatomic and the solid surface.

As shown in Fig. 2, the restricted free-energy barrier for desorption decreases when the coverage or the temperature is increased.

An interesting feature of these free-energy profiles is the minimum at ξ_d ∼ 16 Å for coverages greater than 0.02. The minima are deeper when the coverage or the temperature is increased. A similar minimum has been previously observed by Doren and Tully,^34,35 who suggested that it corresponds to a “precursor state.” The potential of mean force includes contributions from entropy, and therefore, this minimum might not be present in the potential energy surface.

It is common to use density functional theory to calculate the energy barrier for desorption. To study the coverage dependence of the barrier, one assumes a certain distribution of the adsorbates on the surface. Because of computer-power limitations, one cannot examine all possible configurations of the adsorbed molecules. We use a simple model for calculating the potential energy, so we are not subject to this limitation. To explore the fluctuation in the barrier height caused by fluctuations in the positions of the adsorbates, we calculated the energy barriers for 14 different adsorbate configurations generated by Monte Carlo simulations.

Figure 3 shows that there are substantial variations in the energy barriers, the difference between the largest and the smallest being 20 kJ/mol. If the pre-exponential factor is assumed to be the same for both processes, the desorption rates from mean-field theory will differ by three orders-of-magnitude between the two barriers at 500 K.

The numerical values of Q^†/Q_R are given in Table II. The units are Å⁻¹ because when Q^† is calculated, we do not integrate over the distance between the center of mass and the surface, while that integral is included in Q_R. In the calculation of Q^†, the center of mass of one molecule is fixed on the dividing surface and the remaining (N − 1) molecules are adsorbed. The ratio Q^†/Q_R increases with coverage, which we expect because the interaction between adsorbates is repulsive. However, this interpretation is qualitative: Q^†/Q_R depends on the entropy change ΔS^†, which is also coverage dependent. Finally, we point out that to calculate the change in Helmholtz free energy ΔA^† from Q^†/Q_R, we need to use Eq. (23) in which Q^† is multiplied by the bin width Δ_b. Doing so makes the quantity Δ_bQ^†/Q_R dimensionless, which is necessary for calculating ΔA^† = −k_BT ln(Δ_bQ^†/Q_R).

ΔE^† was calculated by performing two independent Monte Carlo simulations. In one, we calculated the mean energy E^† for a system in which one molecule had the center of mass on the dividing surface and the other N − 1 molecules were adsorbed. In the other, we calculated the mean energy E_R when all N molecules were adsorbed. ΔE^† is defined by ΔE^† = E^† − E_R. Because the adsorbate-adsorbate interactions are repulsive, ΔE decreases as the coverage increases.

It is interesting to compare ΔE^† to the desorption energies given in Fig. 3 for a variety of configurations. At a coverage of 0.28 and a temperature of 500 K, the desorption energies in Fig. 3 range between 73 kJ and 90 kJ. Under the same conditions, ΔE^† = 74.7 kJ; this value is sensible because the higher desorption energies in Fig. 3 are weighted by a Boltzmann factor and contribute less to the average ΔE^†.

In order to define the activation energy E_act, we make Arrhenius plots of the logarithm of the desorption rate constant vs the inverse temperature. The plots are linear at all coverages. The results are shown in Table III.

It is interesting to compare the activation energy to the mean desorption energy ΔE^†. Because the rate depends exponentially on energy, we also tabulated the values of exp[(ΔE^† − E_act)/k_BT] in Table III. At a temperature of 500 K, these values are close to 1 at all coverages. This result means that under these conditions, one can approximate E_act by ΔE^†. This is a good news because ΔE^† is easier to calculate. However, the nearness of ΔE^† to E_act is accidental. Using ΔE^† as a proxy for E_act leads to large errors in the rate constant at 200 K and 700 K.

We conclude that the equation

is approximate. Because ΔS^† and ΔH^† depend on temperature, a fit of k_d to P exp[−E_act/RT] gives P ≠ (k_BT/h) exp[ΔS^†/R] and E_act ≠ ΔE^†.

Moreover, the barrier in the potential of mean force also gives a poor approximation to the activation energy. For example, in Fig. 2(b), the desorption free energy at 500 K, for a coverage θ = 0.02, is equal to 60 kJ/mol. The activation energy is 73.7 kJ/mol. This discrepancy means that using TST with the free energy surface is not recommended. In part, the reason for the differences in these energies is the large fluctuations in the desorption energies caused by the fluctuations in the adsorbate distribution on the surface (see Fig. 3).

It is customary to use density functional theory (DFT) to obtain an estimate of the desorption energy. In order to make contact with rate theory, one must estimate the activation entropy. This estimation is most often done by using either a lattice gas model or a two-dimensional ideal gas model for the adsorbed molecules. One hopes that the two-dimensional gas model is adequate at high temperatures and the lattice gas model is adequate at low temperatures. There is no reasonable candidate to model the intermediate temperatures.

In this work, the activation entropy was calculated from ΔS^† = (ΔE^† − ΔA^†)/T, with ΔE^† and ΔA^† obtained from our simulations. The results are given in Table III.

The change in ΔS^† with temperature is remarkably large. In a simplified model for the rate constant, it is assumed that in the empirical Arrhenius formula k_d = A_p exp[−E_act/k_BT], one should identify E_a with ΔE^† and take A_p = (k_BT/h) exp[ΔS^†/k_B]. It is further assumed that ΔS^† and ΔE^† are independent of temperature because A_p and E_a in the empirical formula are temperature-independent. Unfortunately, the situation is not so simple. We find that k_d calculated in our simulation does satisfy a Arrhenius formula: the plot of ln k_d vs 1/T is linear (Fig. 6). From this plot, we extract an activation energy E_a and we find that it is different from ΔE^† (see Table III). Moreover, the pre-exponential A_p is temperature-independent and differs from (k_BT/h) exp[ΔS^†/k_B].

The reason for the mismatch is that ΔH^† and ΔS^† depend on temperature. In the case of ΔS^†, the dependence is fairly strong. For a coverage of θ = 0.02, ΔS^† changes from 44.6 J/mol K at 200 K, to −7.0 J/mol K at 500 K, and to −17 J/mol K at 700 K.

It is common to obtain an approximation to the desorption rate by using DFT to calculate desorption energy for one adsorbate concentration and use it as the activation energy. The entropy is calculated by simple models: either an ideal lattice gas (2DLG) or a two-dimensional ideal gas (2DIG). Since there are many such models, we explain in detail, in the  Appendix, the versions used here.

We now compare our computed results for ΔS^† obtained from exact rate constant theory to those obtained from assuming simple models of adsorption. The 2DIG does not depend on coverage because of the assumption that adsorbate-adsorbate interaction is negligible. 2DLG depends on coverage through the configurational entropy. As shown in Table III, we find that the 2DIG ideal gas model is very close in predicting the behavior of adsorbed molecules at high temperatures. We do not expect a perfect match since ν_⊥ (see the  Appendix) was obtained from experiments and, as mentioned in Sec. IV, we did not try to build a “perfect” model of CO adsorption on Cu. Nonetheless, even this level of matching is encouraging. We find our results to depend slightly on the coverage while the 2DIG ideal gas model is independent of coverage. As mentioned above, we find ΔS^† to increase by ∼10 J/(K/mol) when θ increases from 0.02 to 0.28, which results in a rate-coefficient increase by ∼30%.

At T = 500 K, the calculated ΔS^† matches well with that obtained from 2DIG; however, at a high coverage (θ = 0.28), the calculated ΔS^† lies in between 2DIG and 2DLG. Because the calculated ΔS^† is ∼40 J/(K/mol) away from that for either of the two models, the rate coefficients obtained from the two models will differ from the exact rate by nearly two orders of magnitude. Similarly, at T = 200 K, ΔS^† lies between the values given by the limiting models and picking one of them could lead to huge errors in the rate constant.

Finally, we compare our computed ΔS^† to simplified relations provided by the Campbell^7,8 and Dumesic^9–11 groups. As shown in Table III, the relation provided by Campbell and co-workers is closer to our exact values at T = 200 K, whereas the predictions from Dumesic’s group are not close to the exact values for any temperatures. It is not a surprise that the relation provided by Campbell and co-workers is closer to the values obtained at lower temperatures as it is obtained by fitting to experiments which were performed at lower temperatures. We do not expect an exact match because experiments have been performed on different systems. Moreover, solid surfaces in real systems are different from the one considered here. We cannot exclude the possibility of island formation in experiments, which has not been considered in this work. It may be fortuitous that we get a close match to that obtained by Campbell and co-workers; however, as they showed, the relation is valid for a large number of systems and may also be applicable to our system.

The outcome of gas-surface collision (which is nondissociative) depends on several factors: the surface corrugation, the impact point on the surface, the energy exchange between phonons and colliding molecule, and the conversion, during the collision, of kinetic energy perpendicular to the surface into rotational energy or kinetic energy parallel to the surface. Loss of kinetic energy in the direction perpendicular to the surface leads to trapping of the molecule at the surface.

The efficiency of trapping, as a function of coverage and temperature, is described by the sticking coefficient, given in Table II. We also plot the average sticking coefficient as a function of time of molecular dynamics trajectory for various coverages and temperatures in Fig. 4. Trapping can also be seen in Fig. 5, which displays 500 trajectories, chosen at random, for several surface coverages and temperatures.

At the lowest coverage (θ = 0.02), trapping is very efficient. Even at 700 K, very few trajectories hit the surface and bounce back into the vacuum. Trapping takes place because, during the collision, some of the kinetic energy in the direction perpendicular to the surface is transferred to rotation, to kinetic energy parallel to the surface, or to phonons. As a result, the molecule does not have enough kinetic energy to escape the attraction of the surface. Of course, trapping is also helped by the fact that the desorption energy is much larger than k_BT.

As the coverage increases, a smaller fraction of trajectories is trapped. This happens, primarily, because the attraction between molecules is much weaker than the attraction between a molecule and the surface. Not all incoming molecules that collide with an adsorbed molecule are repulsed into the vacuum. Some acquire velocity in a direction parallel to the surface, do not have enough energy to escape, travel sideways, and may find an empty surface site to stick to.

It is somewhat surprising that the temperature does not have a strong effect on the sticking coefficient. Additionally, we find that the sticking coefficient is weakly correlated with velocity, i.e., mean field rate is equivalent to the exact rate obtained from flux-correlation theory. This is because rotational anisotropy is very large and the barriers are much larger than k_BT.

We show a sample of 500 trajectories for each temperature and coverage studied in this work in Fig. S1 of the supplementary material. Several incoming trajectories are displayed in detail for illustration through movies. Figure S2 of the supplementary material shows a trajectory, for θ = 0.18 and T = 500 K, in which the molecule is trapped as soon as it arrives at the surface. The incoming molecule lands on a surface site that is not covered by an adsorbate and sticks there. Its center-of-mass performs a small-amplitude oscillation in the direction perpendicular to the surface and its rotational motion is suppressed.

Figure S3 of the supplementary material shows a trajectory for θ = 0.28 and T = 500 K. The molecule hits the surface and manages to return to the vacuum after several attempts. Alternatively, Fig. S4 of the supplementary material shows a movie of a molecule colliding with the adsorbed molecule, loosing energy, after several failed attempts find a vacant site and gets trapped on the surface.

We have performed exact calculations of the adsorption and desorption rate constants for a model that contains all essential features of a real system. If the dividing surface is taken sufficiently far from the surface, no molecule whose center of mass reaches it will ever return to the surface: the transition state theory gives the exact rate coefficient. The adsorption rate is more interesting because not all incoming molecules will stick to the surface. Besides the usual thermodynamic partition function, the theory contains a dynamic factor that reflects the sticking probability. We found that the sticking probability is very sensitive to coverage: most incoming molecules stick when the coverage is low. At high coverages, many incident molecules interact with adsorbed molecules. The interaction between molecules is much weaker than the molecule-surface interaction, and as a result, the adsorption probability is diminished. This lessening is not strictly proportional to the fraction of surface sites not covered with adsorbates because some incoming molecules collide with an adsorbed molecule, lose energy, and move along the surface until they find an empty site to which they bind.

The entropic contribution to the desorption rate increases with the temperature and decreases with the coverage. At high temperatures (T = 700 K), the exact entropic contribution is close to that provided by the two-dimensional gas model. However, at lower temperatures, both the two-dimensional gas model and the two-dimensional lattice gas model give errors in the rate constant of two orders-of magnitude.

The mean-field model, which applies the transition state theory to the potential of mean force, also gives errors of two orders-of-magnitude. The reason is that each desorbing molecule sees a different potential energy surface because of different adsorbate configurations. These fluctuations in the energy barrier for desorption can be quite large and have a big effect on the desorption rate.

See supplementary material for (a) the ratio of 2D ideal gas and 2D lattice gas desorption rates for a diatomic molecule, (b) sample trajectories from the dividing surface for different coverages and temperatures, and (c) selected embedded movies of molecular dynamics trajectories from the dividing surface.

We acknowledge support from the High Performance Computing center at IITK and the Center for Scientific Computing from the California NanoSystems Institute, UCSB Materials Research Laboratory: an NSF MRSEC (Grant No. DMR-1720256) and NSF Grant No. CNS-1725797. V.A. acknowledges the funding support from IITK (initiation grant) and SERB (Ramanujan Fellowship No. SB/S2/RJN-013/2017). H.M. acknowledges the financial support from the Department of Energy, Office of Science, Office of Basic Energy Sciences, DE-FG03-89ER14048. We also thank Dr. Celia Wrathall for extensive proofreading and typesetting the manuscript.

Here, we discuss simple models—2D lattice gas (2DLG) and 2D ideal gas (2DIG)—which are used in the literature to calculate the entropies of adsorbed molecules. In the case of desorption, the entropy of a molecule in the transition state is given exactly by a two-dimensional ideal gas. Only the entropy of the adsorbates is treated by using a simplified model.

The 2DLG model assumes that the adsorbed molecules behave like a two-dimensional ideal lattice gas, i.e., they remain localized on an adsorption site where they vibrate. The 2DIG model assumes that the adsorbates translate parallel to the surface and are also free to rotate in two dimensions; the model ignores the corrugation of the surface.

We want to compute entropy to remove one molecule from the surface and place it on the dividing surface, i.e., we want to calculate

where S_R and Q_R are the entropy and canonical partition functions of the adsorbed molecules, respectively, and S^† and Q^† are the entropy and the canonical partition function of the molecule when its center of mass is restricted to stay on the dividing surface.

We also want to compute the transition state theory (TST) rate constant for the two limiting cases. The TST rate constant is given by

Consequently, we assume that the adsorbates are randomly distributed on the surface and they adsorb nondissociatively. We also assume that gas molecules are ideal heteronuclear diatomics (because our system consists of diatomic molecules). Additionally, we use a classic, rigid-rotor, harmonic oscillator approximation to separate rotational and vibrational degrees of freedom. Therefore, we assume that canonical partition function of each molecule can be approximated as

where q_trans, $qvib$⁠, q_rot, q_elec, and q_nucl are the translational, vibrational, rotational, electronic, and nuclear partition functions, respectively. We further assume that the solid surface is homogeneous. N molecules are adsorbed on the surface, and the vibrational stretching frequency of the diatomic molecule is retained on adsorption. The vibrational frequency of the CO molecule changes from 2157 cm⁻¹ in the gas-phase to 2050 cm⁻¹ upon adsorption on Cu.^36,37 We want to compute the change in entropy when one molecule is removed from the surface and placed on the dividing surface.

We first evaluate ΔS^† when the adsorbed molecule is considered a 2DIG. In this case, the ratio of partition is given by

where Q_L is the total partition function of the solid. For the ideal case, we assume that the phonons are not changed by the adsorption of the molecule and therefore their effect will cancel in the ratio of partition functions. ΔE^† is the electronic adsorption energy. The factor N! is a correction that comes from the fact that at sufficiently high temperatures, identical molecules obey “Boltzmann” statistics. Note that for transition state there is one less molecule on the surface and hence this term would become $1(N−1)!$⁠. The factor 1/N appears to preserve the extensivity of $lnQ†$⁠. For the 2D ideal gas model, vibrational, rotational, and 2D translational degrees are retained on adsorption; hence, the partition functions will cancel in the ratio. $qν⊥$ is the vibrational partition function for the motion perpendicular to the surface. In principle, one should use a quantum vibrational partition function. However, as we are performing classical simulations, we take the classical approximation to the vibrational partition function, i.e., $qν⊥=1/βhν⊥$⁠. This is not a bad approximation since the frequencies associated with adsorption are usually small and there is a very small error in using the classical approximation. For example, for a vibrational frequency of 200 cm⁻¹, the error in vibrational partition is less than 10% between quantum and classical harmonic approximations at 200 K. The nuclear partition function contributes through a multiplicative constant which cancels on computing the ratio of partition functions.

Finally, for a 2DIG model, the entropy change (⁠$ΔS2DIG†$⁠) can be written as follows:

and the TST rate constant for the 2DIG is given by

We now evaluate ΔS^† when the adsorbates are considered as 2DLG. We assume that there are M sites available for the molecule to adsorb. After adsorption, the two translational degrees of freedom are converted to two low-frequency vibrational modes (ν_t). For a diatomic molecule, there are two rotational degrees of freedom which we assume are converted to two low-frequency vibrational modes (ν_r). For the rest of degrees of freedom, we use arguments as before. For the 2DLG model, there is an additional term that takes into account the number of ways of arranging N molecules on M adsorption sites. The ratio of partition function is given by

where $M!(M−N)!N!$ is the number of ways for arranging N molecules on M adsorption sites. Note that for the transition state, there is one molecule on the surface and hence this term would become $M!(M−N+1)!(N−1)!$⁠. Again, the factor 1/N is needed to make $lnQ†$ extensive. The value of a is the area of per site (or A/M). Θ_ν is the vibrational temperature and is equal to hν/k_B. Θ_r is the rotational temperature and is equal to h²/8π²Ik_B, where I is the moment of inertia. a is the area/site, and θ is the coverage. We note here that in the last equality, we have assumed M ≫ 1; therefore, 1/M ≈ 0. The ratio of the partition function obtained above has a singularity when θ → 1 and should not be used for θ close to one. We note here that this relates with the fact that Langmuir’s isotherm has complete coverage only at infinite pressures.

Now, for the 2DLG model, the entropy change is given by

The transition state theory rate constant is given by

The dependence of the rate constant on 1/(1 − θ) has previously been shown by many authors.^38–40 We note here that it might seem that the dependence of the rate constant on 1/(1 − θ) would mean that the rate of desorption will be proportional to θ/(1 − θ), which would, in turn, get the adsorption isotherm incorrect. However, as has been pointed out,⁴⁰ one needs to multiply the TST rate with the appropriate recrossing term. The 1/(1 − θ) occurs from the “excluded-area” effect similar to the “excluded-volume” term in van der Waal’s equation. To put this into perspective, let us consider the van der Waal’s equation of state with only the excluded volume correction, i.e., we consider the following equation of state:

where N is number of molecules, R is the universal gas constant, and b is the volume per molecule. The chemical potential of this gas is

This equation is similar to the chemical potential of the 2D ideal lattice gas

where a is the area/site, and thus, Ma is the total area of the solid.