Construction of finite rate surface chemistry models from pulsed hyperthermal beam experimental data

Developing accurate surface chemistry models is important for a number of applications such as chemical processes for manufacturing, TPS (thermal protection systems) design, as well as material processing for semiconductors, and technologies for medical sciences, corrosion protection, and lubrication. Fundamental knowledge about the reaction mechanisms and their rates is critical in the modeling and development of various heterogeneous catalytic processes and other aforementioned technologies.

The construction of finite rate chemistry models capable of predicting surface reactions occurring over a wide range of temperatures can be challenging because of the multiple processes and pathways available for adsorption, desorption, and chemical reactions. The values of the rate constants are sensitive to the surface conditions, temperature, and pressure, as a result of the large disparities in the activation energies and power law dependencies.^1–3 Surface chemistry models are often developed with the use of macroscopic experimental data like total product fluxes, heat flux measurements, material recession, and radiative signatures.^3–20 These quantities are used to infer the concentration of the various products near the surface, which are then utilized for obtaining the corresponding reaction rates. Often one or more reaction mechanisms are assumed for each product in order to fit the rate constants to an Arrhenius form. However, these experimentally measured macroscopic quantities are characterized by highly coupled processes on and near the surface, which are almost impossible to isolate.^3,4 This limits the extrapolation of these finite-rate surface chemistry models to conditions different from those of the experiments. Thus, experimental data elucidating the molecular level details of the gas-surface interactions is of paramount importance in constructing a general, physically accurate surface chemistry model.

Molecular beam experiments have been used extensively to study the kinetics and dynamics of gas-surface interactions.^21–36 These experiments are usually performed in high vacuum environments, thus isolating the surface mechanisms from the gas-phase kinetics. The beam is directed at the surface of interest, where the incident atoms/molecules can undergo chemical reactions or simply scatter from the surface non-reactively. The scattered products are typically detected using a mass spectrometer as a function of angle and time. Analyzing the reactive and non-reactively scattered products can help to elucidate the gas-surface collisions at the atomic level. Particularly, the using a super-/hyperthermal beam provides an excellent way of distinguishing the reactive and impulsively scattered products. These beams can introduce the gas-phase reactants to the surface in a continuous manner or in terms of short bursts (pulses). Modulation of the beam into short pulses is ideal for introducing small controlled amounts of the gas onto the surface and further allows the study of the relative rates of different surface processes. The time resolved measurements, commonly known as time-of-flight (TOF) distributions, in addition to the composition of the surface scattered species, can be used to examine the time scales of the surface reactions, which may range from tenths of microseconds up to seconds.^37–78 This information provides a more detailed understanding of the possible mechanisms through which the detected products are formed.

There are three main steps involved in using pulsed hyperthermal beam data to construct finite rate surface chemistry models:⁷⁹ (i) identification of the reaction mechanisms present in the system from the TOF and angular distributions, (ii) computing the total fluxes of each surface reaction product from their respective distributions, and (iii) derivation of the rate constants for each identified mechanism such that the final product fluxes computed from the model are consistent with those obtained from the experimental data. Although the TOF distributions can capture the time scales of the different surface processes, these time scales often overlap and the final TOF measurements represent a cumulation of the products formed via various reaction mechanisms. Accurate decomposition of the experimental TOF distribution is necessary in order to identify the surface mechanisms which produce the observed scattered products. Further, the computation of the total reaction product fluxes from the TOF data is not straightforward for the slowly desorbing products. Finally, due to the introduction of the gas-phase reactants in pulses, the surface concentrations of the different adsorbed species vary significantly with the time evolution of the pulse, thus making it difficult to quantify the rate constants of the reactions.

This paper describes an approach recently developed for the construction of a finite rate carbon oxidation model from pulsed hyperthermal beam scattering experiments.⁸⁰ We elaborate on each of the steps outlined above and use the carbon oxidation model as an example to illustrate the applicability of this general approach. This approach can be applied to analyze surfaces of different types (metallic/non-metallic), initial states (empty/pre-adsorbed); and also experimental setups with multiple gas-phase reactants within the beam. However, it requires the use of pulsed hyperthermal/supersonic beam.

This paper is organized as follows: Section II provides an outline of the pulsed hyperthermal beam experiments. This section also describes direct simulation Monte Carlo (DSMC), which is the gas-kinetic stochastic simulation technique used for simulating the full hyperthermal beam experiments, including the TOF and angular distributions. DSMC is then used to analyze the different components of the experimental distributions that correspond to the various surface mechanisms (Section III). The procedure to compute the total fluxes of each surface reaction product from the experimental distributions (both TOF and angular data) is outlined in Section IV. Section V describes the methodologies to derive the rate constants of the identified reaction mechanisms. First, the steady-state beam approximation is studied and the errors associated with this approximation are analyzed and presented in Section V A. A general calibration methodology for the transient rate calculation is proposed and investigated in Section V B. The procedures involved in the TOF decomposition, computing the total product fluxes, and transient rate calculation are applied to analyze experimental data in Section VI. Finally, the conclusions are provided in Section VII.

The analysis outlined in the following sections is developed for pulsed beam experimental data. In these experiments, a pulsed hyperthermal beam is focused on the material surface of interest, which is placed within a high vacuum chamber. The gas-phase reactants strike the surface with energy characteristic of the beam at a prescribed incidence angle. Upon impact, the particles can adsorb, undergo direct impact surface reactions, or scatter non-reactively. The scattered products are detected and measured using a mass spectrometer. The number density or flux information of the surface scattered products is collected as a function of time (TOF) at various angles in-plane and out-of-plane of the beam. These TOF distributions which are obtained at different scattering angles are then integrated to compute the angular flux distributions. The surface can be empty or pre-adsorbed with one or more reactants, while the incident beam may also contain multiple chemical species.

For the range of surface temperatures considered here, the incident particle energy from the beam is much greater than the fully accommodated energy from the surface. Thus, there is marked energy difference between the products that are impulsively scattered and those formed via thermal mechanisms (based on the surface temperature), which can clearly be distinguished in a TOF distribution. The TOF distributions provide a wealth of information regarding the time scales associated with reaction and desorption of reactively scattered products. This is especially important for developing chemistry models that is capable of predicting the surface reactions over a wide range of conditions, which may be considerably different from the experimental conditions. Experimental data can be acquired at various surface temperatures, incident energies and angles, in order to obtain a comprehensive picture of the reaction system of interest. Details regarding such pulsed molecular beam experiments may be found in the literature.^39,81,82

Direct simulation Monte Carlo (DSMC)⁸³ is a probabilistic, particle-based computational kinetic method for solving the Boltzmann equation and is widely used for solving rarefied/non-continuum flows.^84,85 Within DSMC, the real atoms and molecules in a gas are represented by weighted simulator particles, which makes such a molecular level description computationally tractable. The identities, velocities (energies) and positions of these simulator particles are updated at every time-step of the simulation. This allows for the direct reconstruction of the energy (or time) and angular distributions of the scattered products, and also their composition, making DSMC an ideal tool to aid in elucidating the surface interactions characterized by these molecular beam experiments.

The surface reactions in DSMC are modeled using the detailed surface chemistry framework developed by Swaminathan-Gopalan et al.⁸⁶ A comprehensive list of surface reaction mechanisms can be simulated within this framework including adsorption (associative/dissociative, direct/indirect), desorption, collision-induced (CI), Eley-Rideal (ER), various types of Langmuir-Hinshelwood (LH), etc. This framework stochastically models the simultaneous reactions on a set of active sites present on the surface. The reaction probabilities and frequencies are computed using the surface concentrations at the current time-step and user-specified rate constants, sticking coefficients, and other surface and reaction parameters. Multiple surface site sets and multiple adsorbed species can also be incorporated within this framework. The maximum number of adsorbed atoms/molecules is restricted by the surface site density of the material.

The DSMC simulations are configured to emulate the experimental setup. For the case of a pulsed beam, the particles are emitted from within a circular region (source) above the surface at regular intervals based on the frequency of the pulse. For a continuous beam, the particles are released uniformly throughout the duration of the experiment. The velocity (or energy) distribution of the simulated beam particles are drawn directly from the experimental distributions. The gas-phase particles from the beam strike the surface and are processed according to the specified surface scattering/reactions included within the model. The information regarding the scattered reactants and the reaction products are collected as they reach a particular distance from the point of impact, corresponding to the position of the experimental detector. These details are also gathered at various in-plane and out-of-plane angles. In this way, the TOF and angular distributions from the simulations can be directly compared with the experiments.

Although the TOF distributions can aid in distinguishing processes with disparate time scales, identifying the reaction mechanisms from the TOF distributions is not straightforward. Often, there is significant overlap of these timescales resulting in a superposition of several reaction mechanisms within the TOF distribution. In order to identify and isolate the effects of different surface reaction mechanisms, a detailed understanding of the form of TOF distributions corresponding to these different mechanisms is essential. To this effect, we perform numerical simulations of the molecular beam scattering from a smooth surface under a range of surface and beam conditions of the candidate surface processes using DSMC. This parametric study is instrumental in identifying the reaction mechanisms at the surface corresponding to the observed experimental data.

The particles incident upon the surface may not interact reactively with the surface, leading to inelastic (or elastic) scattering. The probability of inelastic or impulsive scattering (IS) can be significant even for atoms/molecules with large adsorption enthalpies owing to the high velocities in a hyperthermal beam. The interaction time with the surface for the impulsively scattered atoms is so short that they do not achieve thermal equilibrium, and hence these atoms have translational energies that are significantly higher than the average energy of a Maxwell-Boltzmann (MB) distribution associated with the surface temperature.^25,26,39,87

The IS component velocity distribution can be expressed using a displaced Gaussian functional form:³⁷ 

The parameter u₀ is the mean and the α is the variance of the distribution. These two parameters are treated as free parameters within the model and can be tuned to match the IS component of the experimental TOF distribution.

The IS atoms can be easily identified within the TOF distribution by a sharp peak at times much shorter than the MB distribution peak as shown in Fig. 1(a). The in-plane angular distribution of the IS atoms is usually peaked at angles larger than the specular angle (away from the normal).^39,88,25,87 The particles tend to have greater accommodation along the surface normal compared to the tangential direction. This causes the angular distribution to shift from the specular angle in the direction away from the surface normal. The distribution also tends to be highly peaked with rapid fall off from the peak angle. A representative IS angular distribution is shown in Fig. 1(b).

When a particle from the beam collides with the surface, it can either adsorb, undergo direct impact reactions or scatter impulsively from the surface. The probability of adsorption is determined by two factors: the surface coverage (θ) and the sticking coefficient (S₀). The adsorption probability is linearly proportional to the sticking coefficient as shown in Fig. 2(a).

The dependence of the adsorption probability on the surface coverage varies based on the dynamic pathway through which the adsorption reaction occurs: direct and indirect adsorption. The direct adsorption process is one in which the particle collides with the surface and immediately bonds with the site at the point of impact. In the indirect adsorption pathway, the particle first adsorbs as a precursor (intermediate). This precursor is weakly bound (possibly through physisorption) and can move freely over the surface, eventually chemisorbing at an appropriate adsorption site.

The direct adsorption model was first proposed by Langmuir,⁹⁰ where the atom/molecule becomes completely immobilized upon striking the surface (immobile adsorption). In this case, the probability of adsorption is determined by two parameters. First, the particle must strike an empty site (not occupied by another adsorbed species); second, the particle must form a bond with the site resulting in adsorption. The probability of striking an empty site is determined based on the current surface coverage, and the probability of adsorbing on an empty site is determined by the sticking coefficient (S₀):

where θ is the total surface coverage and S₀ is the sticking coefficient at zero surface coverage. There can be more than one particle adsorbing on the surface, represented by the parameter α, which takes on a value of unity for associative adsorption, and α = 2 or more for dissociative adsorption. The adsorption probability shows a linear dependence on the surface coverage. This is shown in Fig. 2(b) for a single adsorbate (α = 1).

In the indirect adsorption model, the particle gets trapped and forms a loosely bound intermediate or precursor that is short lived. This intermediate can move freely on the surface like a two-dimensional gas. If a suitable site is encountered, the precursor is chemisorbed, otherwise the particle desorbs from the surface. A simple model for indirect adsorption was presented by Kisliuk,⁸⁹ which describes the variation of the adsorption probability as a function of surface coverage θ according to:

Here, K_eq is the equilibrium constant of the adsorption-desorption process of the intermediate:⁹¹ 

The variation of adsorption probability with surface coverage θ for different values of K_eq are shown in Fig. 2(c). For an ideal precursor (K = ∞), the adsorption probability is equal to S₀ for all values of surface coverage except 1. This model reduces to the Langmuir model when K = 0. However, K is expected to have a temperature dependence. Although this model is not derived from first principles, it is well grounded in physics and reproduces many experimental observations.^89,92–95 In general, the indirect or precursor-mediated adsorption is much more likely than direct adsorption. This is due to the lower entropy loss encountered during the process of indirect adsorption compared to direct adsorption.^90,93,94,96

The Eley-Rideal mechanism is a non-thermal direct impact mechanism. This results from the fact that the products formed through the Eley-Ridel mechanism are formed immediately after the gas phase reactant strikes the surface, and the interaction time is not long enough for thermal equilibration with the surface.^25,26,97,98 The products have translational energies that are significantly higher than the average MB distribution, but whose energies are much lower than those of the impulsively scattered (IS) atoms, as shown in Fig. 3(a). Similar to IS atoms, the velocity distribution of the ER products are modeled as a Gaussian distribution (Eq. (1)), where the mean u₀ and variance α are treated as free parameters to fit the observed TOF distribution.³⁷ The in-plane angular distribution of the ER products do not follow a particular distribution and are reported in literature to vary based on the reaction. The only consistent feature is that a peak in angular distributions occurs between the normal and the specular angle^37,99 as shown in Fig. 3(b).

The hot-atom (HA) mechanism is another direct impact mechanism that is possible in a hyperthermal beam. Similar to ER mechanism, HA is a non-thermal mechanism where the products formed on the surface retain a portion of the high beam energy before reacting with adsorbed species.^100–102 These adsorbed species can be several atomic distances away from the point of impact. The resulting products immediately desorb into the gas-phase and do not become accommodated with the surface owing to the slow energy dissipation. Thus, the TOF of the products are more energetic than the average MB distribution, but have lower energies compared to the IS atoms (See Fig. 4(a) and Fig. 1(a) green curves). The angular distributions of the HA products are similar to the ER products with a peak occurring between the normal and the specular angle^103,104 (Fig. 4(b)).

A CI mechanism is similar to an ER mechanism as both involve gas-phase and surface adsorbed reactants. However, the two reactants within the CI mechanism do not chemically react with each other. The energy from the collision of the gas-phase reactant induces desorption of adsorbed atoms/molecules from the surface, and hence these reactions usually involve high-speed particles (super-/hyperthermal velocities).^105–107,25 The incident particle may or may not adsorb onto the surface. If the incident particle does not adsorb, its scattering will be similar to an IS particle. The TOF or velocity distribution of the desorbing particles usually follows a MB distribution at time t=0 (Fig. 5(a)), while the angular distributions is described by the cosine law^108,109 (Fig. 5(b)). However, if a desorption energy barrier or additional energy transfer mechanisms are present, they can alter the shape of both angular and TOF distribution of the desorbing particle (Section III F 2).

Desorption is a thermal mechanism and the products exit the surface with a MB distribution.^25,39,93,95 However, since the reactants are all adsorbed species, this mechanism is not directly dependent on the beam. Hence, the products are not formed immediately after the beam hits the surface, if the reaction rate is relatively low (corresponding to a relatively low surface temperature) and do not follow a MB distribution at time t=0. The desorption of the products into the gas phase occurs over an extended period of time, leading to a cumulative MB distribution over that time period. If the surface were continuously replenished with adsorbed reactants, the TOF signal would be constant with time. However, for a pulsed beam, the TOF signal of these reactions will be based on a convolution of MB distribution with an exponential or power decay depending on the order of the reaction as shown in Fig. 6(a).³⁹ As the rate constant increases (corresponding to increasing temperature), the final TOF distribution shifts closer to the MB distribution at time t=0. When the reaction rate is relatively high (corresponding to a relatively high surface temperature), the products are formed almost immediately after the beam collides with the surface and they follow a MB distribution at time t=0 (Fig. 6(c)). Regardless of the rate of the reaction, the angular distribution of the products follow a cosine distribution (shown in Fig. 6(d)) provided that the species has no desorption barrier.

Due to the nature of interaction between the products and the surface, the desorption may occur over an energy barrier, characterized by a desorption activation energy. For the surface desorption process, this energy barrier exists only in the normal direction.^25,110–112 Thus, only the products having enough energy to overcome the barrier will be able to desorb from the surface. This will change the velocity distribution of the observed products along the normal direction, thus leading to the distortion of the speed/energy distribution and thereby the TOF distribution.^25,113 Fig. 7(a) shows a representative TOF distribution of a thermal process with and without an energy barrier. In a system with an energy barrier, the TOF distribution of the desorbing product is shifted left from a MB distribution. The angular distribution, which is determined by the proportion of the normal to the tangential velocity is also altered due to the energy barrier. The angular distribution is much more peaked towards the normal and follows a cosine power distribution (cosⁿθ).^25,111

Besides the desorption barrier, there are several other energy transfer mechanisms that can alter the observed distributions of the desorbing product. In the local hot-spot mechanism,^105,114,115 the local surface environment surrounding the desorbing molecule has energies greater than the bulk surface temperature. These local hot-spots are formed due to the energy deposited by the incoming gas-phase particles. This is highly likely in a hyperthermal beam, especially for the CI mechanism. In this case, although the products are in thermal equilibrium with local surroundings, the final desorbing energies of the products might not be characterized only by the equilibrium bulk surface temperature, but may also depend on the incoming velocities of the particles. Another energy transfer pathway is due to dissociation or adsorption bond energy corresponding to the incident particle, causing the product to desorb with superthermal energies.^108,109 Although these mechanisms cause the TOF distribution to be different from a MB distribution at time t=0, the “extra” energy does not show a preference to the normal/tangential directions or a dependence on the azimuthal angle, and hence the angular distribution will be a cosine distribution (Fig. 8).

The Langmuir-Hinshelwood (LH) mechanism is a pure-surface mechanism (PS), where the reaction occurs strictly between adsorbed particles and/or bulk atoms which constitute the solid material. The LH products become accommodated with the temperature of the surface, resulting in desorption based on a MB distribution. As discussed in the previous section, the desorption of these species are not directly dependent on the beam. The TOF distribution of the desorbing products may follow a MB distribution at time t=0 or a cumulative distribution corresponding to the desorption rate.

The LH mechanism can be sub-divided into three parts or steps: (i) adsorption, (ii) formation, and (iii) desorption. First, all the reactants must adsorb onto the surface. In the second step, formation, the adsorbed species diffuse on the surface and react with other adsorbed species and surface atoms leading to the formation of the reaction product, which remains adsorbed. In the final desorption step, the adsorbed product desorbs into the gas-phase. Each of these steps may contain many elementary sub-steps. The three steps and the final reaction for a simple CO formation on carbon surfaces with adsorbed oxygen are shown below:

In this system, O(s), CO(s) and CO(g) are the reactant, surface intermediate, and product respectively. The quantities t_f and t_d are the characteristic times for formation and desorption, respectively. If τ is the time scale of interest (which may be the experimentally observable time scale or time scale related to other processes in the system, etc.), an order of magnitude comparison between the reaction times (t_f and t_d) and the time scale of interest τ, results in the following four types of LH mechanisms:

1. t_f ≪ τ, t_d ≪ τ: Prompt LH mechanism

2. t_f ∼ τ, t_d ≪ τ: LH limited by formation

3. t_f ≪ τ, t_d ∼ τ: LH limited by desorption

4. t_f ∼ τ, t_d ∼ τ: LH limited by both formation and desorption

The LH-1 mechanism corresponds to a situation where both processes of formation and desorption are extremely fast compared to the time scale of interest. This means that the surface intermediate is formed immediately after the reactant strikes the surface. Desorption may be a part of the formation step or an independent process that is also rapid. The variation in the concentration of all the species are shown in Fig. 9(a). The concentration of the reactants falls rapidly to zero in a single time step (dt = 10⁻⁶). All the reactants are instantly converted to products, and the concentration of the surface intermediate remains zero. The gas-phase products exit the surface promptly at all temperatures and the TOF distribution follows a Maxwell-Boltzmann (MB) distribution at time t=0 corresponding to the surface temperature as shown in Fig. 9(b). For a LH-1 (prompt LH mechanism), both formation and desorption processes are rapid, leading to the collapse of Eqs. (5)(b) and (c) to a single step.

LH mechanism of type 2 occurs when the time of desorption is much smaller than the time scale of interest τ, but the formation time is of the same order as τ. The formation reaction is slow and is the rate determining step (RDS) in this type of mechanism. The reactants remain on the surface for a relatively long time before forming the surface intermediate. The surface intermediate, however, desorbs promptly as a gas-phase product since the desorption process is very rapid. Fig. 10(a) shows the variation in the concentration of all the species as a function of time. The concentration of the reactant decays either exponentially or based on a power law expression, depending on the order of the reaction. The product shows a corresponding increase with time, while the concentration of the intermediate remains zero. Although the gas-phase products exit the surface immediately after they are formed, the formation step is not instantaneous and follows a decay rate equation. Thus, the observed TOF distribution does not follow a MB distribution, but a convolution of MB distribution and the corresponding decay rate equation as shown in Fig. 10(b).³⁹ 

A LH reaction system follows the characteristics of a type 2 mechanism only for a certain temperature range. If the rate constant of each step follows an Arrhenius form, then the time of desorption and formation reduces with temperature. Thus, for a fixed time scale of interest τ, a reaction system might transition from type 2 to type 1 after a certain temperature threshold. For a LH-2 mechanism (LH limited by formation), only the desorption step is quick, resulting in the collapse of Eq. (5)(c) into Eq. (5)(b) and giving rise to a two-step LH mechanism:

The LH type 3 mechanism corresponds to a situation where the formation process is rapid, and the desorption time is on the same order of the time scale of interest τ. Since the formation time is very small, the surface intermediate is formed almost instantly after the reactants strike the surface. But, the desorption process, however, is slow and is the RDS in this type of LH mechanism. The intermediates that are formed stay adsorbed on the surface for a relatively long time before desorbing. The variation in the concentration of all the species are shown in Fig. 11(a). The concentration of the reactant falls rapidly to zero in a single time step to form the surface intermediate. The concentration of the intermediate decays (either as an exponential or power law based on the order of the reaction) with time, and the product concentration shows a corresponding increase. Since the desorption occurs slowly, the TOF distribution of the gas-phase product follows a convolution of MB distribution at time t=0 and the corresponding decay rate equation as shown in Fig. 11(b). As mentioned previously, a LH reaction system might transition from type 3 to type 1 after a certain temperature threshold. In the case of a LH-3 mechanism (LH limited by desorption), the formation step is quick, leading to the collapse of Eq. (5)(b) into Eq. (5)(a) and giving rise to a two-step LH mechanism again:

The fourth and the final type of LH mechanism occurs when both the desorption and formation times are on the order of the time scale of interest τ. Either of these reactions could be the RDS. Since there is no rapid process involved, the variation of the concentration of the reactant, intermediate and product is more complex. Fig. 12(a) shows the concentration of the species for a typical LH-4 mechanism. Although the exact shape of the curves will vary widely depending on the precise values of the individual rates of formation and desorption, the same general trends are observed. The reactant concentration decreases steadily based on the rate of formation. The concentration of the intermediate shows an initial increase, where the rate of formation is greater than the rate of desorption, owing to the greater concentration of the reactant compared to the intermediate. As time progresses, the concentration of the intermediate reaches a peak and then starts to decrease. This results from the declining formation rate (lower concentration of reactant) and increasing desorption rate (higher concentration of intermediate). The product concentration shows a monotonic increase with time, however its rate varies widely. Initially the rate is small, which increases with time to reach a peak and then starts to fall again, following the trend of the intermediate concentration. Owing to the time-varying rate of the gas-phase product desorption, the TOF distribution has a complex form as shown in Fig. 12(b). Similar to the previous cases, a LH reaction system might transition from type 4 to either type 2 or 3, and finally to type 1 after certain temperature thresholds. For a LH-4 (LH limited by both formation and desorption), none of the processes are rapid, hence no simplifications can be made and no steps can be collapsed:

It is observed that the TOF distribution for a LH-2 and LH-3 mechanism have the same form. Hence, the TOF distribution alone is not enough to distinguish between type 2 and type 3 LH mechanisms if there is only one LH reaction. However, the existence of a parallel reaction system with two or more LH reactions enables us to gain further insight into the exact type of the LH mechanism utilizing only the TOF distribution.

Consider a simple system of two first-order parallel reactions with the same reactant R, and having different surface intermediates (I₁ and I₂) and products (P₁ and P₂) with corresponding rate constants as shown in Fig. 13. If the system of reactions are of type 2 (LH limited by formation), i.e., formation rates (k_f1, k_f2) are slow and desorption rates (k_d1, k_d2) are high, then the concentration of each species would follow the trends shown in Fig. 14(a). The concentration of both the intermediates are zero, and the decay rate of the reactant is determined by the sum of k_f1 and k_f2. The selectivity of products formed is determined by the ratio of k_f1 and k_f2.

If the LH reaction system is of type 2, the values of k_f1 and k_f2 can be determined by the rate of decay of the reactant and the relative concentration of the products using Eq. (10).

The TOF distribution of both products (normalized independently) are shown in Fig. 14(b). The distributions follow a convolution of MB distribution at time t=0 and exponential decay (first-order). Since the product formation rate is determined by the same reactant concentration, the TOF distribution of the products are aligned exactly.

If the reaction system follows a LH-3 mechanism (LH limited by desorption), i.e., formation rates (k_f1, k_f2) are high and desorption rates (k_d1, k_d2) are low, then the variation of the concentration with time of each species would be very different from the previous case. This is shown in Fig. 15(a). The reactant immediately forms the surface intermediates and its concentration instantaneously falls to zero. The concentration of both intermediates slowly decay over time (exponential in this case), while the products show a corresponding increase. The TOF distribution of products (normalized independently) are presented in Fig. 15(b). Similar to the LH-2 mechanism, the TOF distributions follow the convolution of MB distribution and an exponential decay (first-order). However, the TOF distribution of the products from the type 3 mechanism is not the same. This results due to the fact that the formation rate of the products are based on the concentration of their respective intermediates, rather than the concentration of the reactant. This key difference between the decay rates of the products can be used to distinguish the products formed via LH-2 or LH-3 mechanisms.

Similar to the LH-2 mechanism, the final relative concentration of the products is determined by the ratio of k_f1 and k_f2 (Eq. (10)(b)). If the system of LH reactions are of type 3, the values of k_f1 and k_f2 can be determined by the ratio of the concentration of the products to the total influx of reactants. The decay rate of the two products can be used to determine the value of desorption rates k_d1 and k_d2, respectively.

It is important to keep in mind that the desorption rates do not reflect the total concentration of the products, but rather the rate at which these products are observed in the gas-phase. Hence, a reaction with higher decay rate does not necessarily imply that the total concentration of that particular product is higher. The final concentration of the products are determined only by the formation rates. This can be seen more clearly from Fig. 16. The rate of desorption of reaction 1 is ten times greater than that of reaction 2. Thus, a higher amount of product 1 is initially observed from the TOF distribution (Fig. 16(b)). However, the distribution dies off quickly in comparison with reaction 2. Hence, the total final concentration (obtained from integrating the TOF distribution) of product 2 is higher than that of product 1. In fact, the final concentration of product 2 will be ten times greater than of product 1, as can be inferred from the ratio of the corresponding formation rates (Eq. (10)(b)).

The total TOF distribution obtained from the experiments will often have significant overlap of the product distribution from various mechanisms. However, the information from the previous parts of this section can be used to identify and isolate the different mechanisms. The different components comprising the total TOF distribution are identified in Fig. 17. Component I is IS made up of inelastically or impulsively scattered atoms, which have the greatest energy and the shortest flight time. The non-thermal (NT) mechanisms (ER and HA) form products that constitute component II. The mechanisms which produce thermal products immediately after the beam strikes the surface are mainly concentrated in component III. LH-1 and CI mechanisms fall under this category. Finally, the slowly desorbing products resulting from desorption and LH mechanisms of type 2, 3 and 4 are present in component IV.

Although the exact procedure for decomposing the TOF data will vary from one system to another, the procedure outlined here has recently been employed in Swaminathan-Gopalan et al,⁸⁰ and is briefly described with the intent of providing a guideline. In the first step, we take advantage of the fact that all the distributions except for the slow component decays rapidly to zero and are not present at long times. Thus, the profile of the slow portion is determined promptly by fitting to the TOF distribution at long times. In the next step, we identify the thermal (TD) component of the TOF distribution. The desorption of the thermal products is based on a MB distribution at time t=0 and are preferentially scattered towards the surface normal following a cosine law. Both the velocity and angular distribution of the TD products might be altered due to the existence of an energy barrier for desorption. The polar angular distributions can be used to calculate the energy barrier by fitting to a cosine power law. The energy barrier value can then be used to compute the modified TD velocity distribution.⁸⁰ 

Once the slow and thermal components have been captured, the component of the TOF distribution containing the impulsively scattered particles is characterized. A displaced Gaussian profile (Eq. (1)) can be used to describe the distribution of velocity for the IS component. There are two variable parameters: mean u₀ and variance α. Due to the hyperthermal energies of the particles, the IS particle collision with the surface belongs in the structural regime.⁸⁸ In this regime, the soft-sphere scattering model^87,116,39 can be utilized to compute the mean energy of the inelastically scattered particles. The exact procedure to obtain the free parameters (variance and soft-sphere model parameters) from the experimental data has been outlined in Ref. 80.

At last the portion between the peaks of IS and TD, termed as non-thermal (NT), is characterized. The reason for analyzing this NT portion towards the end is because of the several pathways through which this component could be produced: (i) direct non-thermal reactions like ER or HA, (ii) local hot-spots on the surface leading to superthermal energies of TD products,¹⁰⁵ (iii) additional energy transfer due to bond breaking,¹⁰⁸ (iv) multiple bounces of the IS atoms on the surface resulting in high energy loss.¹¹⁷ The profile of the velocity and angular distribution of this TOF portion may be used to identify the exact pathway. The displaced Gaussian functional form (Eq. (1)) can be used to represent the velocity distribution similar to the IS part. Again, the procedure to determine the fitting parameters from the experimental data can be found in Ref. 80.

Once the TOF distributions have been decomposed, the total product flux from each of the mechanisms must be computed in order to calculate the rate constants of the reactions. First, the fluxes must be computed from the TOF distributions at each angle, then these fluxes must be integrated over the angular distributions to obtain the final product fluxes.

The TOF distribution detected using a mass spectrometer can either record the number density^118,87,39 or the flux¹¹⁹ of the products as a function of time. If the flux distribution is measured, it can directly be integrated over the length of time between two pulses to obtain the integrated flux of the reaction products at that particular angle. However, if the number density is measured, then these must be converted to a flux distribution before integration. In order to obtain the flux or intensity distribution (I(t)) from the number density distribution (N(t)), the following formula is usually employed:^39,68,72,120

However, this formula only applies for the products that scatter or desorb from the surface immediately after the beam strikes. If the products desorb slowly over a period of time, this density-to-flux conversion is not applicable since the time t is not only indicative of the time needed for the particle to reach the detector from the surface. The time t also includes the diffusion, formation and residence times of the slowly desorbing products. In order to obtain the flux distribution of the slowly desorbing products, first the desorption rate constant is determined by fitting the convolution of the MB number density distribution at time t=0 and the decay rate equation to the experimental distribution. The flux distribution is then given by the convolution of the MB flux distribution at time t=0 and the decay rate equation (with the rate constant determined in the previous step) according to:

Within the experimental setup, the detection of products desorbing very slowly over extended periods of time is challenging owing to the difficulty in distinguishing their low signal from the background noise.

This becomes especially important at low temperatures when the reaction rate constants are lower. In addition, the sampling window in the experiments during which the TOF data is collected may not cover the entire time between the beam pulses,^39,68,72,120 thus further increasing the fraction of products that goes undetected. Hence, care must be taken so that the sum of all the desorbing product fluxes account for the entire incident flux from the beam. Another important consideration, especially when multiple products are involved, is the mass sensitivity of the detector. Typically, the detector efficiency is higher for larger masses.¹¹⁸ Thus, this mass sensitivity factor of the detector must be taken into account for obtaining accurate fluxes from experimental data.

Integration of the TOF distribution provides the fluxes at a particular reflected angle. This must be repeated for each angle to compute the angular distribution. Typically, these measurements are carried out over a limited number of angles in the plane of the beam (and the surface normal).^{108,37,39,68,72,120} Both the in-plane and out-of-plane distributions is required for each mechanism to obtain the total fluxes of the products from the TOF integrated fluxes at different angles. The polar (in-plane) angular distribution of the TD products is characterized by a cosine law, while the azimuthal (out-of-plane) distribution is usually uniform. However, this can be altered due to the presence of energy barriers and additional heat transfer mechanisms. The angular distributions of other direct impact reaction mechanisms and inelastically scattered species must be explicitly determined either in the current or previous experiments or theoretical studies. While integrating the distributions over the angles, it is important to use analytical or functional forms and perform the integration over the complete range of angles: [0, π/2] for in-plane, and (−π, π] for out-of-plane angle.

In this section, we describe methodologies to compute the rate constants for all the surface reactions identified in the TOF and angular distributions. First, the rate constant of desorption reactions can be directly determined from fitting the TOF distributions of the slowly desorbing products. As described previously, the form of the slowly desorbing product TOF distribution follows a convolution of MB distribution at time t=0 and the decay rate equation. The best fit to the slow component of the experimental TOF distribution with the analytical form gives the desorption rate constants. It is important to note that the desorption reactions need not necessarily be first order reactions.^93,95 Further, their order might vary with coverage of the different species on the surface.

Once the rate constant of desorption reactions have been computed, the final step is to obtain the Gas-Surface (GS) reaction rate constants. The GS reactions are comprised of all reactions which contain gas-phase and surface (or bulk) reactants. This includes direct impact reactions (ER, HA) as well as some thermal mechanisms such as LH-1, LH-3 formation, and adsorption reaction.

The GS reaction rate constants for adsorption-mediated mechanisms are a function of only the surface properties such as coverage and temperature. The rate constants of direct impact and adsorption (sticking coefficient) are functions of surface properties and the incoming gas-phase particle properties such as energy and angle. The variation of the sticking coefficient with the surface coverage and temperature is reasonably captured by the Kisliuk model.^89,92 The parameters within the Kisliuk’s model can be obtained from experimental data or from detailed computational techniques like Molecular Dynamics (MD).

The value of the sticking coefficient may also have a strong dependence on the incident energy of the atoms/molecules. This variation of sticking coefficient is expected to be highly dependent on the system of interest,^110,90,93 and therefore the sticking coefficient must be obtained for each particular case by using experimental data. In the molecular beam experiments, the sticking coefficient can be estimated using the fraction of the impulsively scattered particles.

In order to obtain these rate constants, usually the steady-state approximation is used to estimate the surface conditions to simplify the analysis.^{1,34,35,121,122,36} With this approximation, the analytical solution for the rate constant values can be computed directly from the steady state equations for the product fluxes, surface coverage and site conservation equation. However, this approach neglects the transient characteristic of the pulsed beam, in which the surface conditions still vary widely over the course of a pulse. For instance, based on the temporal location of each particle within the pulse, the surface coverage it encounters will vary, thereby leading to different reaction rates. After the system has reached steady state (where “steady-state” refers to a long exposure of the surface to the pulsed beam prior to the acquisition of the experimental data) the incoming particles at the start of each pulse will experience the same surface coverage at time t=0. As the pulse progresses, the surface coverage increases till the end of the pulse. Then, the desorption of the particles from the surface leads to steady decrease in surface coverage till the next pulse. By the start of the next pulse, the surface coverage recovers to the value at t=0 and this cycle continues. This is in complete contrast to the steady-state approximation, where it is assumed that each particle experiences identical surface conditions, thus neglecting the variations in the surface within the duration of the pulse.

First, a study is conducted to quantify the errors introduced by the steady-state approximation using a baseline surface reaction system with the given mechanisms and rate constants as shown in Table I. These mechanisms and rate constants are taken as the baseline values which serve as the standard for comparison to assess the accuracy of the steady-state assumption vs. pulsed beam approach for determining the model rate constants. DSMC simulations of a pulsed hyperthermal beam of O atoms directed at the surface of interest are performed at various temperatures, and the final compositions of the products that are observed at steady state are shown in Fig. 18. The frequency of the pulse is 2 Hz, and the flux is 5 × 10¹⁸ O atoms/m²/s. The width of the pulse is approximately 20 μs, and the total site density on the surface (Φ) is taken as 6 × 10¹⁸ m⁻². These compositions are then given as input to the steady-state approximation to compute the rate constants, which are then compared with the baseline rate constant values in Table I.

The list of equations included within this study are shown below:

As described in the previous section, the desorption rate constant (k₂) can be directly obtained from the TOF distributions. Hence assuming that the values of k₁ and k₂ are known, the reaction rate constants of k₃, k₄, and k₅ are computed using the above equations with the product fluxes (Fig. 18) as input. The percent error between the computed rate constant values and the baseline values from Table I are plotted in Fig. 19(a). Significant discrepancies are observed between the rate constants. Major differences are observed at higher temperatures, when the surface coverage is very low at steady state. The errors in rate constants k₃ and k₅ are identical owing to the similar dependence of the corresponding products on the surface coverage. These errors also reach very high values at high temperatures as a result of low product fluxes. In order to understand the impact of the differences in the computed rate constants, the steady-state rate constant values are used to perform DSMC simulations of the pulsed hyperthermal O beam. The product fluxes obtained from the rate constants computed using the steady-state approximation are compared with the baseline product fluxes (Fig. 18) and the errors are plotted in Fig. 19(b). The errors in the mole fractions of O₂ and CO₂ reach values greater than 100%, while the error in CO reaches a maximum value of 25%. Due to the coupled nature of the reaction system, the errors in O flux also reach 100% at higher temperatures.

Using the steady-state approximation is analogous to assuming a continuous beam, in which the particles are continuously incident on the surface. In a pulsed beam, however, all the particles strike the surface in a very short period of time (20-100 μs). In the remaining time between pulses, in which the scattered product information is collected, there are no particles from the beam colliding with the surface as shown in Fig. 20(a). Thus for a given overall incident flux, the continuous beam experiences a much smaller instantaneous flux spread uniformly throughout the time between pulses, while the pulsed beam has a much greater flux incident over a short amount of time (duration of the pulse). Fig. 20(b) presents the corresponding variation with time of the surface coverage values under both of these beams. The surface upon which the continuous beam is incident has a surface coverage value that is time invariant. However, within the pulsed beam setup, the surface experiences a sharp increase in the coverage as the pulse progresses. Once the pulse is completed, the surface coverage value drops gradually as the adsorbed particles slowly desorb from the surface. This difference in the time evolution of the surface coverage is the reason for the observed differences between the final product fluxes. Hence, analyzing the experimental data using the steady-state approximation will lead to over/under prediction of finite rate chemistry at the surface.

In the previous section it was established that the steady-state approximation leads to significant errors in the predicted products fluxes and the transient nature of the pulsed beam must be considered. However, this is not straightforward since the rate constant information is not directly available, but must be determined from experimental product flux data. Here we propose a general methodology to precisely compute the rate constants from pulsed molecular beam data.

At a given temperature, the rate constant of a reaction is simply a constant value. The instantaneous rate (or probability in DSMC) is then determined solely by these rate constant values and the current surface conditions. Thus, for each temperature, the problem breaks down to finding the value of these constants for each reaction mechanism for which the final product compositions reproduce the experimentally observed values. In order to construct this calibration problem, the initial values are prescribed for each rate constant, and the complete beam simulations are carried out in DSMC to obtain the final product composition formed via each reaction mechanism. The error is then defined as the difference between the resulting composition and the experimental composition. The goal of the calibration procedure is to minimize this error by calibrating the reaction rate constant values. This process is then repeated at each temperature to obtain the reaction rate constants, which can then be fitted to any form (preferably Arrhenius). A major advantage of such an approach is that no assumption is made about the surface coverage; rather, it is obtained as a result of this calibration procedure. This is of particular importance since the experimental determination of surface coverage is challenging, and it is also expected to widely vary over a range of temperatures. However, this procedure does require as an input the maximum available adsorption sites per unit area of the surface material. This inherently affects the surface coverage and hence the final rates and rate constants of the reactions. Thus, it is important to provide accurate value for this parameter.

Since the objective function to be minimized (i.e., the error between the DSMC and experimental product composition) cannot be explicitly expressed in a functional form, heuristic or meta-heuristic algorithms like the particle swarm optimization algorithm^123,124 is ideally suited for this calibration problem. In addition, this algorithm can search over very large multi-dimensional spaces of candidate solutions, which is required for this calibration process. Furthermore, the particle swarm algorithm is stochastic in nature and is shown to consistently converge to the global minimum.¹²⁵ 

This methodology is applied to the sample system described in the previous section. Initial values are guessed for the three rate constants (k₁, k₂, and k₃). The product composition corresponding to these guess values are computed and compared with the baseline values (Fig. 18). Based on the error between these product fluxes, the guess values are updated. This process continues until the error between the two sets of product fluxes falls below the specified threshold. The rate constant obtained from this methodology is compared to the baseline rate constants (Table I) and the errors are plotted in Fig. 21(a). Excellent agreement is observed between the sets of rate constants and the maximum error is around 0.5%. The corresponding errors in the mole fractions are shown in Fig. 21(b) and as expected the maximum errors are less than 0.3%.

In this section, both the TOF decomposition strategy (presented in Section III) and the transient rate calculation methodology (described in Section V B) are applied to the experimental data of Murray et al.,³⁹ for oxidation of vitreous carbon. First, the analysis of TOF and angular distributions over the range of temperatures is used to identify the set of 9 reaction mechanisms shown in Table II. Notice that there are two different types of CO being produced from the surface labeled as CO{a} and CO{b}. CO{a} is observed to be loosely bound and exits the surface quickly, while CO{b} has a much stronger bond the surface and exited at much slower rates.⁸⁰ Such a distinction between the different types of CO having varying time scales is possible only from the use of TOF distributions. The information of the total composition of different products will not be able to distinguish between these different mechanisms. Further, the analysis outlined in Section III G 2 was used to identify the type of LH mechanism present in the system.

Once the different reactions are identified, the TOF and angular distributions are decomposed into the individual components following the procedure described in Section III H. Fig. 22 presents the decomposition of experimental distributions (both TOF and angular) into the different components (IS, NT, TD, and slow) for O at 1875 K.

Next, the product fluxes for each reaction mechanism are calculated as described in Section IV. These product fluxes are then given as inputs to the particle swarm algorithm to obtain the reaction probabilities at each temperature (Section V B), which are fit to a desired form (Table II). The sticking coefficient of O was found to be independent of the surface coverage and temperature. Thus, the oxygen atom behaves like an ideal precursor (K_eq = ∞), where the particles continue to adsorb on the surface with the same probability (sticking coefficient at zero surface coverage) until saturation. With this assumption, the sticking coefficient for adsorption was also treated as a fitting parameter in the calibration procedure, which converged to a value of 0.85. It is noted that the sticking coefficient may have a dependence on the incident energy of the particle. The current value of 0.85 is valid only at the incident energy of 5 eV (experimental beam energy). The sticking coefficient for oxygen atom is expected to increase with decreasing beam energy, approaching unity for sufficiently low beam energies. Fig. 23 shows the percent error between the analyzed experimental TOF data and the DSMC simulation performed using the derived rate constants⁸⁰ at three different temperatures of 800 K, 1000 K, and 1875 K. The oscillations observed in the error values are due to the stochastic nature of DSMC. The maximum error is less than 7.5% for all the temperature cases. This error is a result of two factors: (i) error from the calibration methodology to obtain the rate constant values at each temperature and (ii) error from fitting the rate constant values to an Arrhenius form as a function of temperature. The error from the calibration procedure is around 0.1%, and hence almost all the error observed is due to the Arrhenius form fits. Finally, Fig. 24 provides the percent error between the analyzed experimental flux and the constructed finite-rate model⁸⁰ from DSMC for different types of O, CO and total fluxes. Here, again all the error values are below 7.5%, and occur because of the Arrhenius fits. Excellent agreement observed for the TOF and angular distributions, and total product fluxes, establishes the validity of the proposed methodologies.

We have presented a general approach for constructing finite rate surface chemistry models using pulsed hyperthermal beam experimental data, which is of great interest in a broad range of fields. First, detailed DSMC simulations of various reaction mechanisms for a wide range of conditions are performed to obtain the TOF and angular distributions. This step is crucial in understanding and interpreting the experimental data and isolating the effects of different mechanisms contained in the pulsed beam experimental data. The Langmuir-Hinshelwood mechanisms were categorized into four types based on the comparison of time scales for formation and desorption steps with the system time scale. The characteristics of each type of mechanism were analyzed in detail and a method to identify and distinguish between them using only the experimental TOF data is provided. A new procedure to accurately compute the total integrated flux of the products, specifically the slowly desorbing species from the TOF and angular distribution is described.

The final and most crucial step of this approach is the development of a general methodology to accurately derive the rate constants of the identified surface reaction mechanisms under transient (pulsed beam) conditions. An analysis of the commonly used steady-state approximation revealed significant differences in terms of both the rate constants and total product composition, specifically when the surface coverage at steady state is low. Modulation of the beam into short pulses inherently introduces a transient component which must be accounted for in order to reliably obtain the rates. Within the proposed methodology, the changing surface conditions are explicitly considered while computing the product fluxes. The fitting procedure uses particle swarm optimization to systematically search through the rate constant parameter space in order to match the final product composition with the experimental values. In addition, the surface coverage of the different species is readily obtained as output and are consistent with the experimental product fluxes. This procedure is performed at each temperature, and the computed rate constant values can be expressed as a function of temperature in any desired form.

The applicability of the TOF decomposition procedure, total product flux computation, and transient rate calculation methodology is demonstrated by applying it to actual experimental data of Murray et al.,³⁹ for oxidation of vitreous carbon. Excellent agreement was observed between the constructed model and the experimental data, thus validating the methods proposed in this work. Although one particular case is considered in this work, the methodology described in this paper is general and can be used to construct finite rate models for chemical reactions at any gas-solid interface from pulsed hyperthermal beam experimental data. This approach can be applied to analyze surfaces of different types (metallic/non-metallic), initial states (empty/pre-adsorbed), experimental setups with multiple gas-phase reactants within the beam, and experimental setups with multiple beams.

This work was performed under the Entry System Modeling Project (M. J. Wright Project Manager) at the NASA Game Changing Development (GCD) Program and supported by NASA Grants NNX15AU92A and the University of Illinois. KAS also acknowledges support from an Early Career Faculty award from NASA’s Space Technology Research Grants Program.

B

total available sites for adsorption

I

surface intermediate

K

equilibrium constant

k

reaction rate constant

M(b)

bulk species on the surface

M(g)

chemical species in the gas phase

M(s)

chemical species adsorbed on the surface

P

product

R

reactant

S₀

sticking coefficient

(s)

empty surface site

t

time

Abbreviations

ad

adsorption

DSMC

direct simulation Monte Carlo

des, d

desorption

ER

Eley-Rideal reaction mechanism

f

formation

HA

hot-atom reaction mechanism

LH

Langmuir-Hinshelwood reaction mechanism

MB

Maxwell-Boltzmann distribution

PT

prompt thermal reaction mechanism

TOF

time-of-flight

TPS

thermal protection system

Greek

θ

surface coverage, angle

τ

characteristic time, s