Microscopic modeling of gassurface scattering. I. A combined molecular dynamicsrate equation approach
Abstract
A combination of first principle molecular dynamics (MD) simulations with a rate equation model (MDRE approach) is presented to study the trapping and the scattering of rare gas atoms from metal surfaces. The temporal evolution of the atom fractions that are either adsorbed or scattered into the continuum is investigated in detail. We demonstrate that for this description one has to consider trapped, quasitrapped and scattering states, and present an energetic definition of these states. The rate equations contain the transition probabilities between the states. We demonstrate how these rate equations can be derived from kinetic theory. Moreover, we present a rigorous way to determine the transition probabilities from a microscopic analysis of the particle trajectories generated by MD simulations. Once the system reaches quasiequilibrium, the rates converge to stationary values, and the subsequent thermal adsorption/desorption dynamics is completely described by the rate equations without the need to perform further timeconsuming MD simulations.
As a proof of concept of our approach, MD simulations for argon atoms interacting with a platinum (111) surface are presented. A detailed deterministic trajectory analysis is performed, and the transition rates are constructed. The dependence of the rates on the incidence conditions and the lattice temperature is analyzed. Based on this example, we analyze the time scale of the gassurface system to approach the quasistationary state. The MDrate equation model has great relevance for the plasmasurface modeling as it makes an extension of accurate simulations to long, experimentally relevant time scales possible. Its application to the computation of atomic sticking probabilities is given in the second part (paper II).
Keywords: plasmasurface modeling, lowtemperature plasma, gassurface interaction, adsorption and scattering of neutral particles, molecular dynamics, rate equations, transition rates, quasiequilibrium, residence time, argon atom, platinum surface
1 Introduction
Lowtemperature plasma physics has seen remarkable progress over the last decade. This concerns both fundamental science studies and technological applications ranging from etching of solid surfaces to plasma chemistry and plasma medicine. In each of these applications, the contact between particles of the plasma and a solid plays a crucial role. This contact is very complex and includes a large variety of fundamental physical processes such as secondary electron emission, sputtering, neutralization and stopping of ions and adsorption and scattering of neutral particles as well as chemical reactions. In the majority of previous studies in lowtemperature plasma physics, these processes have been treated phenomenologically and a better understanding, which would lead to predictive capability, has been recognized to be a major bottleneck for further progress in the field [1, 2].
The goal of the present work is to present a microscopic theorybased approach to a specific problem of plasmasurface interaction: the scattering, adsorption and sticking of rare gas atoms from a plasma at a metal surface. This is expected to be particularly important at low pressures where the neutral fraction of heavy particles is by far dominating over the ionized component. Our approach presents a combination of microscopic modeling within a Langevin molecular dynamics (MD) approach with an analytical model, formulated as rate equations for relevant surface states. This coupled approach has the advantage of opening the way towards ab initio based longtime simulations, as we will explain in detail below.
We start with a brief overview on previous theoretical works treating the interaction of gas atoms with solid surfaces and summarize their strengths and limitations. The energy transfer and the scattering of rarefied gases from different surfaces have been the subject of multiple studies. Much effort has been devoted to determine an accurate scattering model for gas atoms at a solid surface and a realistic description of the involved collision processes. Already Maxwell [3] proposed a simple model in his studies of gas–surface interactions, where the scattered gas atoms are separated into two fractions: one that is reflected specularly and exchanges no energy, and a second one that is accommodated completely and, eventually, desorbs with an equilibrium distribution. Other studies have been focused on the evaluation of the thermal accommodation coefficients [4].
The angular distributions for different gas atoms (helium, neon, argon, krypton, xenon, and deuterium) scattered by a singlecrystal tungsten (110) surface has been intensively analyzed by Weinberg and Merrill [5]. A direct measurement of the velocity distribution of argon atoms scattered from polycrystalline surfaces has been in the focus of the work of Janda et al [6]. These measurements allowed to reveal a dependence of the average kinetic energy of scattered argon atoms on the average incident kinetic energy and surface temperature.
A theoretical treatment of the interaction of Ar atoms with a selfassembled monolayer on Ag(111) in terms of the stochastic scattering theory, including direct scattering, trapping, and desorption, has been reported by Fan and Manson [7]. Gibson et al [8] presented a detailed study of Ar scattered from an ordered 1decanethiol–Au(111) monolayer.
All these detailed experimental and theoretical studies have proven that the use of atomic probes as scattering projectiles can be a useful experimental tool for studies of structure and dynamical properties of surfaces. The scattered intensities can be measured as a function of the final translation energy, scattered and incident angles. Furthermore, such analyses provide important information on surface corrugation and temperature effects.
Presently, there exist three main theoretical research directions in this area. The first one is based on kinetic theory where important contributions are due to Kreuzer and Teshima [9] and Brenig [10, 11]. More recently Bronold et al. studied the sticking of electrons at a dielectric surface [12] and the neutralization of ions at a gold surface [13]. Overall, the resulting kinetic equations provide a powerful semianalytical tool to compute (mainly) stationary properties of gas atoms or charged particles scattering from a surface within suitable manybody approximations.
The second direction is based on ab initio quantum simulations, most importantly BornOppenheimer density functional theory (DFT) or timedependent DFT (TDDFT). The main advantage of the latter is that the dynamics of the electrons and the internal excitation of adsorbed atoms and molecules can be incorporated on a full quantum level. Among recent applications of TDDFT, we mention the analysis of the chemical reaction dynamics of hydrogen on a silicon surface [14] and the energy loss (stopping power) of ions on graphene and boron nitride sheets [15]. An alternative ab initio approach is based on nonequilibrium Green functions (NEGF) and is advantageous when electronic correlations in the surface are of importance, e.g. [16, 17]. Recently NEGF simulations of the stopping power of hydrogen and helium ions on correlated twodimensional materials were presented [18].
The third approach is semiclassical MD simulations that are extensively used in surface science. Here, the quality depends on the accuracy of effective pair potentials (or force fields) whereas details of the electron dynamics are not resolved. As examples of recent applications, we mention the simulation of metal cluster growth and diffusion [19, 20] and of the dynamics of bimetallic clusters [21].
The temporal evolution of the atoms trapped near the surface is of particular interest for the understanding of the adsorption and scattering of neutral atoms. Their equilibration kinetics is of fundamental importance for the understanding of the dependence of the sticking probabilities on the energy, incident angle and lattice temperature. However, these time dependences are typically out of the range of both kinetic theories and ab initio quantum simulations. While the former usually concentrate on stationary properties, TDDFT and NEGF simulations are computationally extremely expensive and are limited to very short times of the order of a picosecond and small spatial scales. Therefore, MD simulations represent the only approach capable of resolving the dynamics on sufficiently long time and spatial scales being of experimental relevance. Although the typically required time step in these simulations is well below one femtosecond, recent progress in computer power made it possible to investigate the relevant surface effects on a microscopic level reaching simulation times of several hundreds of picoseconds.
In the present work (paper I and paper II) we perform extensive MD simulations to study the scattering and sticking of argon atoms at a platinum (111) surface in a timeresolved fashion. In order to analyze the sticking and thermalization behavior, we introduce a novel approach: the trajectories of gas atoms that are near the surface are subdivided into three classes: trapped (T), quasitrapped (Q) and continuum (C) states. We demonstrate that these states are the relevant observables to analyze the sticking problem with excellent statistics, high accuracy and temporal resolution. Furthermore, we demonstrate that the fractions of atoms in trapped, quasitrapped and continuum states obey a simple system of coupled rate equations. Its solution allows us to significantly reduce the computational cost that is otherwise spent on the temporal resolution of individual particle trajectories. Moreover, we demonstrate how the transition rates between the three states can be accurately extracted from our MD simulations transforming the rate equations into an, in principle, exact description. The accuracy of the latter is only limited by the accuracy of the used pair potentials. Finally, our analysis of the temporal evolution reveals that after a characteristic equilibration time these transition rates become stationary. This means that the rate equations become sufficient to study the dynamics of the systems for longer times without further need of MD simulations. This provides the potential to significantly extend the temporal and spatial scales of the simulations without compromising the accuracy.
The main goal of the present paper is to introduce this new combined molecular dynamicsrate equation (MDRE) approach in detail and to test it thoroughly on a representative example: the scattering of argon atoms at a platinum (111) surface. Here extensive experimental and theoretical data are available for comparison, which allow us to critically assess the validity and limitations of our model. A detailed analysis of the argon sticking probability, its dependence on temperature, incident energy and angle is presented in a separate paper (“paper II”, cf. Ref. [22]).
The paper is organized as follows. In section 2 we introduce the microscopic model of the gassurface interaction. The setup of the MD simulations is explained in section 3. The tracking of the particle trajectories provides a direct access to time dependencies of the energy loss distribution functions and the average momentum (kinetic energy). The timeresolved evolution of the trapped, quasitrapped, and continuum states and its dependence on the lattice temperature and incident angle is discussed in detail in section 4. In sections 5 and 6 we introduce the rate equation model, which allows us to reproduce the MD results of section 4, and, moreover, has a potential to extrapolate these data to much longer times being not accessible by usual MD simulations [see subsection 6.4]. The conclusions are given in section 7.
2 Effective gassurface interaction
In a common approach, the global interaction potential of the th gasatom (a) with the surface (s) is decomposed into a sum of pair potentials according to
(1a)  
(1b) 
where is the distance between the th gasatom and the th atom of the surface and denotes the number of surface atoms. The distance to the surface, , in Eq. (1a) is treated as a parameter
(1b) 
while the lateral atom position is kept fixed to analyze the dependence of the gassurface interaction on the adsorption site (see below). However, this lateral position can be varied for the calculation of the potential energy surface (PES).
Such type of global interaction potential and its decomposition into pairwise terms for the ArPt(111) system has been recently reconstructed using the periodic DFT approach [23]. Similar analyses have been performed for the interaction between an argon atom and gold surfaces [24]. In general, the computation of such interaction potentials is still a challenging task, but it follows a standard scheme. Ab initio calculations are performed for relatively small surface clusters, where the reconstructed pair potential (1b) can be crosschecked to compare their performance with stateoftheart van der Waalscorrected periodic DFT approaches. As a result, accurate pair potentials suitable for molecular dynamics simulations are derived. The parametrization values used for the present ArPt(111) system are listed in table 1.
hcp  fcc  atop  bridge  

[Å]  3.3576  3.3576  3.5058  3.3788 
[meV]  76.8951  76.8876  68.8091  75.7618 
[Å]  3.2412  3.2412  3.3894  3.2624 
[meV]  78.1302  78.1246  69.8872  76.9785 
The obtained global ArPt surface interaction potential is presented in figure 1 as a function of the height above the surface plane . For the unrelaxed lattice the uppermost layer is placed at . The lateral position of the Ar atom has been varied in the plane to analyze the corrugation of the potential energy surface. Several sites according to the lattice symmetry (“atop”, “bridge”, “fcc” and “hcp”) have been chosen, as also specified in figure 1 for the fcc(111) lattice. The corresponding potential energy minima and distance to the surface are given in table 1. Their difference with respect to the energy of the fcc site is additionally plotted in figure 1. It characterizes the corrugation of the PES.
In figure 1 (top) we compare the results obtained for the “ideal” (unrelaxed) lattice with the optimized (“relaxed”) lattice. The latter was reconstructed at the surface temperature by minimization of the total energy of the () Pt slab. It accounts for a shift of a few upper layers to more negative values of the coordinate, where the bottom three layers are kept fixed at the same position as in the unrelaxed lattice. This fact can be seen in figure 1 by a shift of the potential energy minima (see the curves labeled “rel”) to smaller heights relative to the plane . The estimated binding energy of meV is in good agreement with the total adsorption energy of about meV reported in the experimental work of HeadGordon et al. [25]. This value is also used as the reference for the optimization of the empirical potentials aimed to reproduce the experimental results on the vertical Pt(111)Ar harmonic vibrational frequency of meV and on the trapping, desorption and scattering data [26, 27, 28]. For an overview of the most commonly used empirical potentials we refer to the recent work of Léonard et al. [23].
In order to perform a diffusive motion on the frozen surface, i.e., lattice atoms are fixed in their equilibrium position, the kinetic energy of the adsorbate atom should exceed the energy barrier at the bridge site. However, no such restriction applies at finite temperatures as the adsorbate atoms can continuously exchange energy with the lattice atoms. In the following, three different lattice temperatures are simulated: , 190, and K, corresponding to 6.90, 16.38, and 25.7 meV, respectively. Hence, the thermal energy supplied from the lattice is significant to overcome the energy barrier for the intersite hops in all cases.
3 MD simulation of the lattice at finite temperature
We perform deterministic molecular dynamics simulations to study the ArPt(111) system. A platinum crystal slab of atoms (consisting of 6 atomic layers) is used. The sample is divided into three parts. Three lower layers form a static crystal. The atoms in this part are frozen at their equilibrium positions and serve as a basement for the dynamical upper layers. The three uppermost layers are an active zone of the crystal. They interact dynamically with the incoming gas atoms. Two bottom layers in the active zone consist of the atoms which are used as a boundary thermostat to realistically treat the removal of energy from the active zone. In this way the excess kinetic energy can dissipate from the active zone. This removal of excess kinetic energy is simulated by restoring the temperature of the heat bath atoms using Langevin dynamics [29]. The Langevin term keeps the kinetic energy of the lattice atoms at the value specified by the lattice temperature . Finally, the Langevin term was switched off for the uppermost layer, and only dynamical correlations due to the gaslattice (ArPt) and lattice atoms (PtPt) binary interactions are retained.
The system of atoms is described by the potential energy
(1c) 
where is the number of gas atoms. The interactions between the gas atoms are neglected, i.e., , assuming sufficiently low gas density. In the directions parallel to the surface ( plane) periodic boundary conditions are applied.
Effective atomatom pair potentials are used for both ArPt [see Eq. (1b)] and PtPt interactions. The interactions between the Pt atoms are modeled using the modified embedded atom potential, which quite accurately reproduces the spectrum of the transverse and longitudinal phonons [30].
At the beginning of the simulation, i.e., before introducing the Ar atom, the crystal is allowed to equilibrate to the surface temperature . This procedure takes about ps. The equilibration is monitored by the instantaneous lattice kinetic energy. After the equilibration phase ( ps) the distribution functions of the tangential () and normal () components of the kinetic energy for the atoms in the three dynamical layers have been evaluated. They quite accurately reproduce the expected (onedimensional and twodimensional) Boltzmann distribution at the temperature
(1d)  
(1e) 
This procedure confirms the correct implementation of the heat bath via Langevin molecular dynamics.
After the lattice has approached steady state, an Ar atom is introduced at a height Å above the surface outside the cutoff radius of the potential . The trajectories are obtained by integrating the classical equations of motion using a fourthorder propagator algorithm [29], to accurately account for the effect of the random force (the Langevin term). The integration time step is fixed at fs. Trapping probabilities, energy exchange calculations and energy distribution functions are evaluated based on samples of trajectories. During an “elementary” event, the impinging gas atom interacts with all atoms within the cutoff radius Å. A simulated trajectory is stopped when (i) the gas atom leaves the surface after undergoing one or more collisions with the surface (called “bounces”) and attains a distance above the surface greater than and (ii) the gas atom experiences more than bounces.
The value of the number of bounces was chosen empirically. By tracking the temporal evolution of trajectories with , we observe in most cases the thermalization of the adsorbate atoms to the lattice temperature and convergence of the energy distribution functions to their quasistationary form on the corresponding time scale. In particular, we analyzed the accommodation of the parallel and perpendicular momentum components. Our MD simulations at low and high temperatures reproduce the general trend of a slower accommodation of the parallel momentum component, as it was first pointed out by Hurst et al [31] and confirmed in many other analyses [32, 34, 35, 23]. This thermalization analysis for the system Ar on Pt (111) is presented in detail in paper II[22].
4 Timeresolved trapped, quasitrapped and scattered fractions
In many physical systems a strong chemical bonding to the surface is the dominant trapping mechanism. In contrast, in the present system the trapping process occurs due to the physisorption potential well described by a van der Waalstype potential, cf. Eq. (1b). Depending on the lattice temperature, the trapped particles desorb after a finite residence time [see subsection 6.4]. The central question is whether this time is sufficient for the adsorbate atoms to equilibrate with the surface so that they eventually desorb with a (quasi)equilibrium energy distribution function. We underline that this is not an academic question but one of practical importance. Indeed, the theoretical description is expected to simplify significantly in cases where thermalization is observed.
This problem was first put forward by Maxwell in his studies of gassurface interactions [3] and was further taken up by Knudsen [36]. The basic concept is based on the thermal accommodation and the efficiency of energy exchange at the gassurface interface. The latter crucially depends on both the incident gas parameters and surface characteristics. In the special case of complete accommodation the desorbed particles leave the surface with a distribution function specified by the Knudsen flux [7]. This type of assumption is frequently used to explain experimental data for the sticking probabilities of the adsorbate and for the distribution functions of energy, momentum and flux vector of the thermally desorbed atoms.
One of our goals is to test this assumption by microscopic modeling of the gassurface interactions and analyze the equilibration kinetics. Using realistic calculations based on the microscopic model introduced in sections 2 and 3, we aim at reproducing the available experimental data for the sticking probabilities [33] and at providing a general framework for the analysis of the temporal evolution of the atom states localized near the surface. While we treat the scattering classically, quantummechanical effects are included by means effective interaction potentials, which are constructed from groundstate DFT calculations [23]. Inelastic effects observed in the scattering events at large energies and high surface temperature are fully taken into account by the Langevin MD scheme [29].
4.1 Classification of particle trajectories: trapped, quasitrapped, and scattering states
The scattering from the surface is modeled by using monoenergetic gas atoms with fixed values of incident kinetic energy and angle . A single collision with the surface introduces a transition from an initial momentum, , to a new momentum state . The final states “f” are distinguished by the surface normal, , and parallel, , components, i.e. . Correspondingly, the kinetic energy of a particle with momentum k is split into two orthogonal contributions,
(1fa)  
(1fb)  
(1fc) 
In addition, every particle moves in the potential landscape of the surface atoms that is characterized by the local surface binding (physisorption) potential, , giving rise to the total energy
(1fg) 
In the following, the dependence of on the local atom position is not explicitly specified. In addition, all energies and their corresponding distribution functions are evaluated at a minimum distance from the surface, where the total energy of a gas atom is approximately conserved. The temporal evolution of the total energy is shown in figure 2. Here and in the following, we use the value as a time unit, where is the Bohr radius, denotes the atomic mass of the argon atoms, and is the Hartree energy. It becomes clear that plateau regions of are well separated from the energy "jumps" corresponding to the inelastic collision processes. During this inelastic process, there is a strong energy exchange with the surface atoms and, therefore, the evaluation of Eq. (1fg) becomes meaningless.
After interaction with the surface a fraction of particles is scattered back, whereas another fraction is (temporally) trapped. A classification of the states is straightforward by analyzing the particle energy.
 I. Scattering states:

the condition to leave the surface due to the momentum exchange is . The backscattered (unbound) particles are referred to as “continuum” (“C”) states.
 II. Bound states:

particles with remain localized near the surface. The localization depends solely on the normal component, while the parallel component can be arbitrary. Therefore, depending on the sign of the total energy (1fg), such states can be further subdivided into
 a. Trapped (“T”) states:

these are particles with .
 b. Quasitrapped (“Q”) states:

these are particles with .
These three categories of particles are characterized by the particle numbers , and , respectively, with the number of atoms . Even for a fixed value , the three contributions can vary with time and with the incidence conditions and the surface parameters. In other words, an analysis of the dynamics of , and should provide detailed information on the gassurface interaction and on the specific system.
It is generally expected that the quasitrapped states (Q states) dominate at large incident angles (with respect to the surface normal) at first. The parallel momentum does not change much during a single reflection event. A particle accelerates towards the surface and gains a large normal kinetic energy of meV by passing the depth of the physisorption well (cf. table 1). As a result, it collides with the surface close to the surface normal when the parallel momentum remains practically unchanged.
If a particle remains localized, its final momentum after every subsequent reflection can be projected sufficiently close to the surface plane, and, hence, the parallel momentum can be strongly perturbed by scattering at the atoms in the upper surface layer. Therefore, it is important to analyze the temporal evolution of the T and Q states and their equilibration mechanism.
As an example, the classification of surface states for three dynamical trajectories is presented in figure 2. In addition to the temporal evolution of the total energy , the distance to the surface is displayed as a function of time. It allows us to uniquely identify reflection events. Notice that a relatively fast energy exchange with the surface takes place at each reflection. Between the reflections the total energy is nearly conserved and correlated with the increase of the height . Once a particle moves in the opposite direction to the surface, it can be reflected back at the turning point specified by the kinetic and potential energy at the right boundary of the physisorption well. As our MD simulations show, here the total energy is conserved and, hence, the reflection can be treated as an elastic process.
4.2 Transitions between T, Q and C states
Furthermore, the transition between the states along a particle trajectory is indicated by the line in figure 2. The three values are used to identify the three states , respectively.
The continuum states () are observed before the first collision () and for the final states with , where both energies are evaluated at the atom position . The corresponding trajectories have no turning point and leave the surface region. In between, such trajectories can experience multiple bounces and transitions between the trapped () and quasitrapped () states depending on the sign of the total energy .
By analyzing a statistical ensemble of the trajectories of , a first “physically” relevant observation can be made. Before a particle is desorbed, it typically gets excited to a quasitrapped state. The direct excitation probability from a trapped to a continuum state is significantly reduced at low lattice temperatures (e.g. K), but it steadily increases with , as will be shown in detail in section 6.
For large statistical ensembles of the trajectory states we can explicitly evaluate the three fractions with , , and , where is the number of trajectories with at a given instant of time. They satisfy the normalization condition .
Figure 3 shows the temporal evolution of these three fractions of Ar atoms influencing a Pt(111) surface at K for two conditions of incidence as well as the corresponding average number of bounces and their variance, , on the same time scale. The value taken at defines the initial sticking probability. Detailed results on the present system can be found in paper II [22]. This quantity remains unchanged until a second bounce takes place. On average, this second bound occurs around .
At low lattice temperature and incident atom energy, the initial sticking fraction remains quite large: for . The continuum states show a saturation already for corresponding to ps. Within this period the trajectories experience on average bounces with the surface.
In contrast, the “pure” trapped states show a convergent behavior at much later times, mainly due to an exponentiallike decay of the quasitrapped states. Within the time interval there is a fast drop of from to for due to a fast conversion into trapped states. This can be seen from the increase in , which is accompanied by a practically constant value of confirming that the particles remain localized near the surface.
The case shows a similar trend, where only the absolute values are different. The fraction of trapped states is large: already after the first reflection. This is to be expected, because less energy is carried by the parallel component at smaller angles, i.e., . Hence, once the trapping condition is satisfied by the normal component, the total energy is negative in most cases (). For the simulations predict that the initial fraction of quasitrapped states can reach only and completely decays into trapped states within .
We conclude from figure 3 that the trapping becomes more efficient at smaller angles. At the lattice temperature K the fraction vanishes after (33 ps), while the fraction remains quasistationary with no noticeable thermal desorption observed on the simulated time scale (40 ps).
5 Derivation of rate equations from kinetic theory
In this section we give a brief derivation of the timedependent integral equations for the population of surface states. The transition probabilities can be expressed via microscopic quantities. Finally, we demonstrate how the kinetic equations can be reduced to a simplified description in terms of a rate equation model and energydistributionaveraged transition coefficients. As we show in section 6, such a model allows for an efficient description of the timeresolved transitions between the trapped and continuum states. It also makes it possible to analyze the time, energy and temperature dependence of the sticking probability and the desorption rates.
Following the derivation presented by Brenig [10], we start with a general definition of the quantum mechanical transition probability between two states of the scattering atom. The transition of an atom from an initial (“i”) state at time to a final (“f”) state at time is given in terms of the transition matrix element of the quantum mechanical time evolution operator . Here, the two indices denote the states of an adsorbate atom, , and of the substrate, . Due to the unitarity of the evolution operator the transition probability obeys the general properties
(1fha)  
(1fhb) 
Starting from the initial state, the transition probability to a final state can be defined as
(1fhi) 
which depends linearly on and contains a delta function to guarantee the energy conservation. While this probability has the familiar form of Fermi’s golden rule, the underlying assumptions should be recalled and critically assessed. First, it is assumed that their are no correlations between the adsorbate atom and the surface. This allows to write the total kinetic energy in the initial and final states as a sum of adsorbate kinetic energy and surface energy. Second, the linear dependence of the transition probability on is a result of the perturbation theory and is valid for small time differences . Third, the fact that the probability depends only on (and is independent of ) implicitly assumes that the system is stationary. Finally, the appearance of the delta function is a consequence of the Markov approximation which assumes that the correlation time has passed, and the energy spectrum has become stationary, see e.g. [16, 37]. All these assumptions are justified in case of a macroscopically stationary surface which is only weakly perturbed by the scattering process of the adsorbate atom.
In particular, we can safely assume for collisions of atoms with thermal and subthermal energy that the initial and final states of the substrate remain close to thermal equilibrium. The recoil energy, initially transferred to the surface layer, is rapidly dissipated to the bulk due to fast atomic vibrations taking place with the Debye frequency and a strong coupling between the substrate atoms when compared to the coupling with the adsorbate atoms. As a result, a saturation to the thermal equilibrium within the substrate is expected on time scales much shorter than the adsorbate thermal accommodation time.
Hence, as a further simplification, we can introduce the transition probabilities averaged over the initial states of the substrate (specified by the Boltzmann factor )
(1fhj) 
This allows us to obtain the temporal evolution of the system of gas atoms alone, which is obtained with the PauliAnsatz [10]
(1fhk) 
which assumes low gas atom density so that all scattering events can be treated as independent of each other. Using the definition of the kinetic coefficients according to , we end up with the system of coupled kinetic equations between the initial and final states
(1fhl) 
where the limit has been taken. This result is well known and directly expresses the connection between microscopic calculations [38, 39] and kinetic theory [40, 41].
As a next step, we explicitly specify the quantum numbers of the gas atoms moving near the surface. Since quantum diffraction effects can be neglected due to the large mass of the gas atoms, a quasiclassical treatment can be used. Thus, it is sufficient to specify the states by the initial and final momenta , which are additionally split into a tangential and parallel component relative to the surface. Consequently, we introduce the following notations
(1fhm)  
(1fhn) 
in the kinetic equations. Note that the averaged transition rate (1fhn) depends solely on the incoming and outgoing momentum.
As a next step, we assume that the transition rate has only a weak directional dependence on the angle between the vectors and . Its main dependence results from two scalars, namely the perpendicular and parallel kinetic energy components of the adsorbate atom, i.e., . This assumption is important for the derivation of the rate equations presented below. It can be justified by the theoretical treatment of a classical collision from vibrating surfaces [42, 43, 44].
In particular, in Refs. [45, 46] an explicit expression for the zerothorder reflection coefficient has been derived for an atomic projectile colliding with a surface consisting of discrete scattering centers (of mass ) whose initial momenta are given by an equilibrium distribution at the lattice temperature . This expression is given by
(1fho)  
and shows a dependence on several key parameters. Here, is the component of the incident momentum, is the form factor of the scattering center, which depends on the interaction potential, and is the kinetic energy after (before) the collision. The only parameter, which contains a directional dependence, is the recoil energy expressed as
Due to the inelastic and Brownianlike character of the scattering processes, the contribution of the first two terms should dominate. The directional dependence (third term) is expected to be weak. Thus, we can use . Finally, it is assumed that the scattering amplitude is a constant, with a value derived for hard sphere scattering.
Now we can proceed and explicitly define the population of surface states and the interstate transition probabilities by the dependence on the kinetic energy components
(1fhp)  
(1fhq) 
In the following we omit the subscripts “i” and “f” and indicate the initial and final states, instead, by using different energy symbols:
(1fhr)  
(1fhs) 
Following the detailed discussion given in section 4, all possible surface and scattering states can be classified into three categories using energy criteria. Now, the trapped, quasitrapped and continuous particle fractions introduced in section 4 can be explicitly defined via integration of the timedependent energy distribution functions, which define the population of different surface states, according to
(1fhta)  
(1fhtb)  
(1fhtc) 
Here, we introduced the shorthand notation with T, Q, C for the distinction of different states and the integration limits. This energy integral always comprises a double integration over the tangential and normal kinetic energy components. The upper integration limit for the tangential component specifies that the trapped and quasitrapped states stay localized (bound) near the surface due to the condition (or ). Here, is the potential energy in the physisorption well where the total particle energy is nearly conserved. The second integral, i.e., that over the parallel component, introduces a distinction between bound and continuous states.
Using the definitions (1fhr) and (1fhs) for the notation of the initial and final state energies, we can rewrite the kinetic equation (1fhl) in the form
This equation specifies the temporal evolution of the three types of states introduced in Eqs. (1fhta)(1fhtc) by taking the time derivative. This yields the result
(1fhtu) 
with = Q, T, C. The first term defines the incoming flux from all possible states and the second term represent the outgoing flux from the state . The inner integral over can be splitted into the contribution of different surface states, . This allows to introduce the energyresolved interstate transition rates, which carry the energydependence of initial state ( = Q, T, C) and are integrated over the energy of the final state ( = Q, T, C) according to
for . The diagonal terms with are mutually cancelled in Eq. (1fhtu) and, therefore, they can be excluded from the consideration. We rewrite Eq. (1fhtu) to the form
(1fhtv)  
and end up with a set of rate equations
(1fhtw) 
for the final states with = Q, T, C, where
(1fhtx) 
are the energy distributionaveraged transition coefficients.
This set of coupled equations can be further simplified, once the system reaches a quasiequilibrium state. This regime can be identified from the convergence of the momentum and energy distribution functions to a quasistationary form, . The distribution converges to the shape specified by the quasiequilibration time , and stays unchanged up to some timedependent scaling factor . In this regime the transition rates (1fhtx) also become timeindependent with the constant values
(1fhty) 
defined by the quasiequilibrium energy distribution of the system in the different surface states (e.g. = Q, T, C.). This result will be used in the next section, where we demonstrate how these specific values can be derived from MD simulations.
6 Rate equation model
Based on Eq. (1fhtw), the temporal evolution of the three types of surface states can be analyzed on the quantitative level by the set of rate equations
(1fhtza)  
(1fhtzb)  
(1fhtzc) 
where the backwards transitions from the continuous to the bound states are assumed to be negligible. The fractions of atoms in the three states change with time due to different decay channels. For instance, the first two terms in Eq. (1fhtza) take into account the decay of the Q states into T and C states with the transition rate and , respectively. This means we use the notation for transitions of the type . According to (1fhtx), the transition rates in the set of equations (1fhtza)(1fhtzc) are generally timedependent, i.e., . They crucially depend on the energy distribution function of the initial and final state in general. Hence, the energy distribution function is nonstationary as well until the system reaches quasiequilibrium.
6.1 Reconstruction of the transition rates from the MD simulations
Now we demonstrate how the rates can be extracted from the MD simulation data. In figure 4 we plot the change in the population of three states (Q, T and C) due to the net transition fluxes
(1fhtzaaa)  
(1fhtzaab)  
(1fhtzaac) 
Two incident conditions are compared, which correspond to the energy meV at the angle and meV at , respectively, for the lattice temperature of 80 K. Each of the decay channels can be uniquely identified by analyzing the initial and final state along the trajectory for every bounce event shown in figure 2 for several examples. The transition rate can be extracted by performing the numerical differentiation
(1fhtzaaab) 
for the curves shown in figure 4.
The following procedure has been used. First, the MD data have been smoothed by a Gaussian kernel
(1fhtzaaac) 
where (), and is the total number of points on the curve. The new values are evaluated using the weighted contribution on neighboring points, and, hence, the statistical fluctuations at each point are suppressed. The differentiation of with respect to time leads to
(1fhtzaaad) 
The smoothness of the derivative can be controlled by the parameter and adjusted to give better agreement with the MD data. Typically, the value was found to be a reasonable choice.
Alternatively, the transition rates can be expressed via the integral form
(1fhtzaaae)  
Here, we used the assumption that the system state has reached a quasiequilibrium state for . Then, the transition rate has only a weak timedependence, i.e., . For the example shown in figure 3, the equilibration time is about . Finally, the quasiequilibrium transition rate can be determined from Eq. (1fhtzaaae) as the ratio of the integrated population of states
(1fhtzaaaf) 
By a proper choice of the equilibration time , the estimated rate should exhibit only a weak time dependence and represent the asymptotic limit of the more general timedependent rate in Eq. (1fhtzaaab).
The data are provided by the MD simulations. To determine the transition rates and the population of states more accurately, we have used the statistical averages over several thousand trajectories. We performed the analysis similar to figure 2 for each set of initial parameters () and determined the fluxes between different states.
6.2 Test of the approach
The comparison of the transition rates given by Eqs. (1fhtzaaab) and (1fhtzaaaf) is demonstrated in figure 4. We choose the time of equilibration for the rate and for . The results of Eq. (1fhtzaaab) are represented by the dots, and those from Eq. (1fhtzaaaf) by the solid lines. Both should match at . The choice of in each case needs some adjustments, such that should be nearly constant for . The quality of the quasiequilibrium approximation can be checked on the right panel. The rates , reconstructed from Eq. (1fhtzaaab), accurately fit the MD data (presented by the dotted curves) for all simulation times and the changes in state populations (, and ).
Due to a statistical noise present in the MD data, the extracted rates exhibit artificial oscillations, see e.g. . At longer times they can be successfully removed using the equilibrium estimator Eq. (1fhtzaaaf), which is much less influenced by the statistical noise. The corresponding equilibrium rate is shown by the horizontal dashed curves on the left panels. In particular, the rate for the incident angle () can reproduce the change in the population of quasitrapped states, , for (). The left panel in figure 4 shows that the corresponding nonequilibrium rate actually converges slowly to on a time scale, which depends on the incident angle.
A similar analysis can be performed for and . Here, we observe that the quasiequilirium assumption can be introduced at a significantly earlier moment. Using we can nearly exactly reproduce the MD data (see the right panel) for the conversion of the quasitrapped states to the trapped state, i.e., , and with some deviations observed at the desorption of the quasitrapped states to the continuum, i.e., . The rates and stay practically constant over the entire simulation period, starting at . However, at the rates drop to smaller values. The later behaviour is artificial and should be explained by the smoothing procedure applied in Eq. (1fhtzaaac). It perturbs the slope of the MD curve around the time of second reflection from the surface ().
In summary, we have demonstrated the efficiency of the rate equation model and provided the proof of convergence to quasiequilibrium after some transient time .
The next important point concerns the analysis of the impact of the lattice temperature . In figures 5 and 6 we compare results for the two temperatures K and K. The incident angle is and the initial gas energy is meV. Compared to the case with meV shown in figure 3, we observe an increase by a factor of in at in figure 5 for K, i.e., a much higher fraction of particles is reflected after the first bounce. At the same time, the initial population of the trapped states is reduced by factor 6 (from to ), while the population of Q states is increased by (from to ). During the Q states practically vanish due to the decay channel Q T. Furthermore, the trapped state demonstrate stability against the thermal desorption at this lattice temperature of 80 K. However, the temporal evolution of the state populations is very similar to the case with the incident energy meV (figure 3).
The situation changes at the higher lattice temperature K shown in figure 5. The initial reflection coefficient is similar to the lowtemperature case ( K), but then it rapidly increases with two characteristic rates. As shown in figure 6, a higher rate is found for due to a fast decay of the Q states into the C and T states: and . Because is larger than , the conversion to the trapped states is the dominant process when the fraction (t) is large.
This trend changes at (or after 23 bounces with the surface). As it is becomes clear from figure 5 the fraction is reduced by a factor of 2, while the fraction reaches a local maximum. For the trapped states dominate and are steadily converted to the continuum states. Note that the decay of the T states takes place much faster than the decay of the Q states.
The fraction first saturates around and then slowly decays to at . This behavior can be explained on basis of the rate equations. As shown in figure 6 for K, the rates of the mutual conversion T Q between the states T and Q differ by a factor of 4. However, they contribute in the rate equations (1fhtza) and (1fhtzb) being multiplied by the population factors and, hence, all incoming and outgoing net fluxes can be compensated for a given state (e.g. ) if the detailed balance condition
(1fhtzaaag) 
approximately holds. Indeed, this condition can be satisfied for the case shown in figure 5 when the ratio of both populations increases to for . The proof why the relation (1fhtzaaag) should hold in general will be given in section 6.4.
One can use the rates shown in figure 6 to analyze how fast the system reaches quasiequilibrium. Using the data for and , we estimate for K and for K. Hence, a higher lattice temperature favors a faster adsorbate equilibration. A more quantitative discussion of the convergence to quasiequilibrium is presented in paper II [22].
The corresponding analysis of the temperature effects for the larger incident angle is shown in figures 7 and 8 for the two temperatures and K. Due to the larger incident angle, the initial population of quasitrapped states is above and slightly decreases with increasing temperature. The population of trapped states is around for both the surface temperatures. While the values specified at the time of the first bounce () are similar, their temporal evolution is quite different due to an enhanced thermal desorption at the larger temperature K. This can be clearly resolved from the temporal evolution of the continuum states, given by the curves . The trapped fraction seems to saturate during the timeinterval , and it subsequently decreases for due to thermal desorption going faster at K. This can also be identified by a faster increase of the continuum fraction . During the period the fraction is reduced by a factor of from to , while it is reduced only by at K.
In contrast, the states show only a weak dependence on . They decay with a similar slope for both lattice temperatures. As in figure 5, we observe a fast conversion via the two channels Q T and Q C at first. The first channel dominates due to a high rate (cf. figure. 8). Its value decreases only slightly with increasing temperature, whereas the asymptotic value of the two other rates and at quasiequilibrium drops by a factor of 2 by lowering the temperature from to K. This finding seems quite reasonable. The deexcitation transition takes place, when an excited state (a quasitrapped trajectory) releases an energy to go into a lower energy state (a trapped trajectory).
In contrast, the two transitions T Q and Q C require that a finite portion of energy should be supplied from the lattice. In this case, the excitation probability should scale with the population of phonon modes and depend on multiphonon excitations. Here, a strong temperature dependence is expected. Indeed, the data for and presented in figures 6 and 8 confirm this expectation and demonstrate a clear temperature dependence. However, if the rates are compared at the same lattice temperature, but different incident angles (see and in figure 4), the difference in the rates does not exceed .
We can conclude that the rate equations provide a very useful tool to analyze the nonequilibrium kinetics during the first few picoseconds. The temporal evolution of the population of different states can be successfully described by the net fluxes in terms of a set of statistically averaged parameters – the transition probabilities . These transition rates can be accurately extracted from the MD data by analyzing the temporal behavior of the particle trajectories. Typically, we observe that the saturation of the transition rates at their equilibrium values can be reached within , i.e., ps, for incident energies below meV. The equilibration takes longer for larger incident angles/energies and lower lattice temperatures.
6.3 Temperature dependence of the transition rates
Results for the transition rates at quasiequilibrium for various incidence conditions, i.e., different incident energies , incident angles and surface temperatures , are summarized in figure 9.
First, we discuss the deexcitation transition from Q to T states (QT). The corresponding rate exhibits only a weak dependence on temperature and incident angle/energy. It stays practically unchanged for incident energies meV and when the lattice temperature is varied in the range . Hence, the binary atomatom collisions play a major role here, while surface temperature effects are secondary.
The value increases when the incident angle is closer to the surface normal. It is about larger for than for . However, this comparison is performed at different incident energies to guarantee that the initial sticking probability is the same for both the incident angles. If the rates are compared at similar incident energies, the observed difference is reduced. This becomes obvious e.g. when comparing the cases meV at versus meV at or meV at versus meV at in figure 9. A more detailed analysis of the ()dependence would require a larger set of data.
Now, we consider the transition rates to a higher energy state: T Q, T C and Q C. They reveal a strong, partially linear dependence on the lattice temperature. In addition, the rate shows a clear dependence on the initial energy , whereas the rate is almost independent of . Such behavior originates from a different portion of energy transferred from the lattice and required for each type of excitation. Much less energy is required for the T Q transition, when only the parallel component of kinetic energy needs to be changed to make the total energy positive. The most probable contribution is expected from the trapped states in the highenergy tail of the distribution function. Here, some correlations with in the incident beam should be present. In contrast, during the T C excitation the normal component needs to be changed significantly by an amount comparable with the depth of the physisorption well being about 80 meV for Ar on Pt(111) to bring a particle to the continuum. Here, the correlations with the incident energies below should be small. Note also that the transition probability is a factor of smaller than that of and on the average. Therefore, the most probable excitation to the continuum is a two stage process: T Q followed by Q C.
The transition rates presented in figure 9 have a direct practical application. In combination with the rate equations, they can be used to extrapolate the temporal evolution of the state populations to longer time scales being not accessible by usual MD simulations. In particular, this is important for the analysis of thermal desorption at low lattice temperatures, such as K, when a significant depletion of the trapped states can be observed only on the time scales exceeding those that are used in the present simulations, i.e., for (40 ps). Some applications of this idea will be discussed more in detail below.
6.4 Analytical solution of the rate equations
In the present section we present the analytical solution of the rate equations introduced in section 6.1 and demonstrate its efficiency to predict the temporal evolution for longer times. In the end, we derive the estimator of the average residence time of the adsorbate atoms trapped on the surface prior to their thermal desorption.
In case the transition rates are timeindependent, the rate equation model reduces to a system of homogeneous linear differential equations of first order and can be solved analytically. Therefore, we rewrite Eqs. (1fhtza)(1fhtzb) in matrix notation according to
(1fhtzaaah) 
Here is a column vector with the elements , , and is a matrix of the transition rates with the elements
(1fhtzaaai) 
Notice that the fraction of continuum states is not considered here as it follows directly from particle number conservation.
Once the matrix is diagonalized, the eigenvalues and the eigenvectors with , define the complete solution which can be written in the form