Extending first principle plasma-surface simulations to experimentally relevant scales

Metadynamics is a technique to enhance the computation of a multidimensional free energy surface (FES) of a many-body system. In most cases the FES is far too complicated to be directly computed. Laio and Parrinello in 2002 introduced a method that allows for an efficient computation of the FES by means of molecular dynamics [27] which has become a widely used tool in computational (bio-)physics, chemistry and material science. An overview is given e.g. in [39].

The first key idea (and fundamental assumption) of the method is that the free energy F of a system with a set of coordinates x, the potential V(x) and the inverse temperature β = 1/(k_BT) can be expressed as a function of just a few collective variables ${\bf{S}}=({S}_{1},\,\ldots ,\,{S}_{d})$ according to

$$\begin{eqnarray}&&F=-\displaystyle \frac{1}{\beta }\mathrm{ln}\left(\int \exp [-\beta V({\bf{x}})]\delta [{\bf{S}}-{\bf{S}}({\bf{x}})]\,{\rm{d}}{\bf{x}}\right).\end{eqnarray} \tag{ 1 }$$

If the collective variables provide an adequate representation of the whole configuration space, the FES can be efficiently explored by performing MD simulation where a second key idea is applied: a history-dependent small 'bias potential' ΔV is added, i.e. $V\to V+{\rm{\Delta }}V$. This potential successively enforces the system to leave every occurring minimum in the FES, thus, avoiding bottlenecks due to rare events. In a very simple fashion, this bias potential can be constructed as a sum of weighted Gaussian functions,

$$\begin{eqnarray}&&{\rm{\Delta }}V({\bf{S}},t)=w\displaystyle \sum _{t^{\prime} \in { \mathcal T }}\exp \left(-\displaystyle \sum _{i=1}^{d}\displaystyle \frac{{[{S}_{i}(t)-{S}_{i}(t^{\prime} )]}^{2}}{2{({\rm{\Delta }}{s}_{i})}^{2}}\right),\end{eqnarray} \tag{ 2 }$$

where the set ${ \mathcal T }$ comprises all times t' before the time t, at which the sum of Gaussian functions has been extended by one term. The interval between successive creations of new Gaussian functions as well as their weights w and widths Δs_i should be chosen such that a compromise between computation time and accuracy is achieved. After a sufficiently long simulation time the potential landscape V + ΔV levels out and becomes flat. This can be recognized in the simulation by the 'diffusive' behavior of the considered collective variables [39]. Then the inverse of the final bias potential ΔV provides an accurate estimator of the free energy F.

In contrast to the methods discussed below, the use of metadynamics alone does not yield correct state-to-state dynamics. Nevertheless, it can be utilized for that purpose if it is combined with other methods, such as the CVHD described in section 4.2. Metadynamics is a versatile method and it allows for a relatively easy implementation. However, the choice of the collective variables can be very difficult.

Using hyperdynamics [40], the state-to-state dynamics of an infrequent event system is accelerated by adding a space-dependent bias potential ΔV(r) to the potential energy surface V(r). Thereby, the energy barriers between different states are reduced so that transitions occur more often. For the applicability of the method, it is required that both the unbiased and the biased system dynamics obey the so-called transition state theory (TST) [41, 42]. Furthermore, the bias potential ΔV(r) must be zero at all dividing surfaces, and it must be chosen such that the correct relative probabilities of the transitions are maintained.

The construction of an appropriate bias potential can be an elaborate task in many cases. In the original publication [40], the diffusion of an Ag₁₀ cluster on an Ag(111) surface was investigated by constructing ΔV as a function of the lowest eigenvalue of the Hessian matrix. At this, boost factors of roughly 8 × 103 were achieved.

Another approach introduced by Fichthorn et al [43], the so-called bond-boost method, is to let ΔV be a function of the nearest-neighbor bond lengths in a solid. In a study of the diffusion of Cu atoms on a Cu(001) surface in the temperature range between 230 and 600 K, boost factors of up to 10⁶ could be achieved [44].

Even higher boost factors of up to 10⁹ were obtained for the same system using the CVHD method introduced by Bal and Neyts [29]. The idea of this method is to use the concept of metadynamics for an incremental build-up of a bias potential depending on just one collective variable, i.e., a variable that describes the relevant processes in the system. The CVHD method can be applied whenever the requirements of hyperdynamics are fulfilled and an appropriate collective variable can be found. Even though the latter may be difficult in some cases, the CVHD method has the potential to be applied to many kinds of different systems. For example, it has already been used to study the folding of a polymer chain model [29] and the pyrolysis and oxidation of n-dodecane [45].

The idea of TAD [28] is to make transitions occur more often by performing the simulation at an elevated temperature T_high instead of the temperature of interest T_low. Because this procedure alone would induce wrong ratios of escape probabilities, an additional mechanism is applied to correct for this. The method can only be applied if the system obeys the harmonic TST (HTST). Thus, it is more restrictive than hyperdynamics and parallel replica dynamics (section 4.4), for which the harmonic approximation is not necessary.

TAD is carried out by performing 'basin constrained molecular dynamics' (BCMD) for each system state. Whenever a transition occurs at T_high, the escape path and the corresponding escape time are stored, but the involved particles are reflected back to their initial energy basin and the dynamics is continued. For each observed transition at time t_high, the transition time is extrapolated to the lower temperature according to

$$\begin{eqnarray}&&{t}_{\mathrm{low}}={t}_{\mathrm{high}}\exp \left\{{E}_{{\rm{a}}}\left(\displaystyle \frac{1}{{k}_{{\rm{B}}}{T}_{\mathrm{low}}}-\displaystyle \frac{1}{{k}_{{\rm{B}}}{T}_{\mathrm{high}}}\right)\right\},\end{eqnarray} \tag{ 3 }$$

where E_a is the energy barrier of the transition. The BCMD routine is stopped when the simulation time reaches

$$\begin{eqnarray}&&{\tilde{t}}_{\mathrm{high}}=\displaystyle \frac{\mathrm{ln}(1/\delta )}{{\nu }_{\min }}{\left(\displaystyle \frac{{\nu }_{\min }{t}_{\mathrm{low},\min }}{\mathrm{ln}(1/\delta )}\right)}^{{T}_{\mathrm{low}}/{T}_{\mathrm{high}}},\end{eqnarray} \tag{ 4 }$$

where t_low,min is the minimum of the extrapolated transition times at T_low, ν_min is a guess for a lower bound of the pre-exponential factors occurring in the formulas for all possible transitions, and δ is a pre-defined limit for the probability to observe every transition after ${\tilde{t}}_{\mathrm{high}}$ that would replace the transition at the current minimum t_low,min. As a result of the procedure, the simulation time is advanced by t_low,min, and the process corresponding to this time is executed.

The boost factor that can be achieved by means of TAD depends critically on the ratio of T_high to T_low. As one cannot choose arbitrarily high values of T_high without breaking the requirements of HTST, the method is particularly effective for systems at low-temperature. For example, a boost factor of 10⁷ was achieved in a simulation of the growth of a Cu(100) surface at T_low = 77 K [38]. A recent example for a simulation at a higher temperature of ${T}_{\mathrm{low}}=500\,{\rm{K}}$ can be found in [46], where the sputter deposition of Mg–Al–O films was studied.

While standard parallelization techniques are usually applied to extend the accessible system sizes, parallel replica dynamics (ParRep) allows one to use parallel computing to extend the time scales, too [47]. Among the three methods hyperdynamics, TAD and ParRep, ParRep is the most accurate one, and a higher boost can be trivially achieved by increasing the number of processors N_p [36, 37]. ParRep can be applied to any infrequent event system with first-order kinetics, i.e., with exponentially distributed first-escape times of all occurring processes,

$$\begin{eqnarray}&&f(t)=\lambda \exp (-\lambda t).\end{eqnarray} \tag{ 5 }$$

The ParRep procedure starts by replicating and dephasing the system on each available processor. Each copy of the system is propagated independently and in parallel on each processor until a transition is detected on one of the processors. Then, the system clock is advanced by the sum of the N_p individual simulation times, and the global system state is set to the state reached after the observed transition. Subsequently, a short serial run is performed to allow for the occurrence of correlated events. After that, the whole procedure is repeated.

As the number of processors is one limiting factor of the achievable speed-up, the ParRep method is often less effective than hyperdynamics and TAD. Nevertheless, it is simple to implement, and it can be combined with other acceleration methods. Therefore, ParRep has become a valuable tool in the field of computational materials science [48]. For example, it has been applied to simulate the diffusion of H₂ in crystalline C₆₀ [49], the diffusion of lithium in amorphous polyethylene oxide [50], and the crack-tip behavior in metals [51].

In some cases, it is reasonable to assume that the motion of atoms on a surface or in a medium is approximately Brownian. This type of motion can be generated by solely performing Langevin dynamics for the particles of interest, while the other atoms and molecules of the background medium do not have to be explicitly included in the simulations. Here we consider, as an example, the deposition of gold atoms onto a polymer surface. The MD simulations tracked each individual atom, its diffusion on the surface, the emergence and growth of clusters and, eventually the coalescence of the latter. A typical example is presented in figure 4 and shows the cluster configuration at an early moment (top) and a later time point (bottom).

Zoom In Zoom Out Reset image size

The influence of the plasma environment is mostly due to the impact of energetic ions. This leads to the formation of surface defects that trap incoming atoms and prevent their diffusion. The figure compares the cases of a weak plasma effect (right column, the fraction of atoms trapped equals γ = 0.001) and a stronger effect (left column, γ = 0.05.) In the former case, a small number of large clusters is being formed, due to cluster coalescence, whereas in the latter case the film is much more homogeneous, containing a much larger number of smaller clusters [52].

We underline that the lower snapshots in the figure refer to a film thickness of about one nanometer which requires a deposition time of about two minutes. This is impossible to achieve with first principle MD simulations. The key to achieve this extreme simulation duration and to compare to experiments was SPA. The main ideas of the approach are explained in the following.

In the MD simulations the isotropic Langevin equation of motion for all gold particles with the mass m and spatial coordinates ${\bf{r}}=({{\bf{r}}}_{1},{{\bf{r}}}_{2},\,\ldots )$ are solved:

$$\begin{eqnarray}&&m\ddot{{\bf{r}}}=-{\rm{\nabla }}U({\bf{r}})-\displaystyle \frac{m}{{t}_{\mathrm{damp}}}\dot{{\bf{r}}}+\sqrt{\displaystyle \frac{2\,{{mk}}_{{\rm{B}}}T}{{t}_{\mathrm{damp}}}}{\bf{R}},\end{eqnarray} \tag{ 6 }$$

where the potential U describes the interaction between gold particles. For this potential ab initio force field data are being used (the MD simulations used the LAMMPS package). Further, t_damp has the role of a damping parameter, and R is a delta-correlated Gaussian random process. This random force and the viscous damping simulate the effect of the polymer on the heavy gold particles. If only neighborless atoms are considered, i.e., ∇U ≡ 0 , the combination of the last two terms on the right side of equation (6) induces a diffusive motion with the diffusion coefficient

$$\begin{eqnarray}&&D=\displaystyle \frac{1}{m}{k}_{{\rm{B}}}{{Tt}}_{\mathrm{damp}}.\end{eqnarray} \tag{ 7 }$$

Thus, it is clear that the utilization of the Langevin dynamics allows one to control the speed of the surface diffusion and bulk by choosing a specific combination of the temperature T and the damping parameter t_damp. An anisotropic diffusive motion can be generated if one generalizes equation (6) by separately defining damping parameters ${t}_{\mathrm{damp}}^{x}$, ${t}_{\mathrm{damp}}^{y}$ and ${t}_{\mathrm{damp}}^{z}$ for each of the three spatial directions. Beyond that, it is possible to add a spatial dependence to the diffusion coefficient if one lets the damping parameters depend on the position of the particle.

Based on the above considerations, Abraham et al [30] developed a procedure to simulate the growth of nanogranular gold structures on a thin polymer film. Instead of simulating the polymer with explicit particle models, their method relies on performing Langevin dynamics for the gold atoms with the simulation box being partitioned into three parts representing the upper part of the polymer bulk (I), the surface of the polymer (II) and the region above the surface (III). By choosing appropriate ratios of the damping parameters, one can make sure that the atoms spend most of the time in the surface layer (II), where they perform a random walk which is restricted to a small range of possible z-coordinates. The use of Langevin dynamics is restricted to regions (I) and (II); in region (III), the dynamics is purely microscopic. This allows one to add particles to the system by creating particles at the top of the simulation box and assigning them a negative initial velocity. Therefore, it is possible to perform the simulation with values of the deposition rates J_sim and diffusion coefficients D_sim that are much higher than the values in typical experiments.

In [30], it was argued that the simulations yield an adequate description of a real experimental deposition process if the ratio J_sim/D_sim is equal to the ratio J_exp/D_exp of the corresponding quantities in the experiment. The idea behind that is that—at least at the early stage of the deposition process—the growth should be essentially determined by the average distance an atom travels on the surface between successive depositions of atoms. Hence, the absolute time of the process is assumed to be irrelevant. The results presented in [30] were obtained with a time step of 1 fs, and a damping parameter for the diffusion in x- and y-directions of 1 ps. The temperature and the deposition rate were set to match the conditions of the experimental results in [53] for the sputter deposition of gold on polystyrene. Using these parameters, the direct MD simulation time for the growth of a thin gold film is roughly 10⁹ times shorter than the corresponding time in the experiment. Or in other words, the duration of the MD simulations could be extended by nine orders of magnitude. To verify the validity of such a dramatic shift of the time scales, comprehensive tests of the method were performed, see also [34, 52] for a discussion. In particular, as one accelerates only selected processes, i.e., the deposition of atoms and the diffusion of atoms on the surface, one has to make sure that the neglect of other processes, e.g., the relaxation of a cluster structure, does not lead to artifacts in the simulation results. In [30], the method was tested by comparing several quantities describing the evolution of the gold film morphology with the results of time-resolved in situ grazing incidence x-ray scattering experiments of Schwartzkopf et al [53]. It turned out that many of the experimentally observed features could be reproduced for film thicknesses up to 3 nm. This thickness corresponds to an impressive effective simulation time of 367 s which is directly suited for comparison with measurements. As an example of the compared quantities, figure 5 shows the number density of metal clusters on the polymer surface. The comparison between experimental data and simulation results shows very good agreement, at least up to a time of about 350 s. For longer times, the simulation start to deviate from the measurements indicating that the procedure is no longer applicable.

Zoom In Zoom Out Reset image size

The present approach of selective acceleration of dominant processes can be generalized to other systems as well. A recent application concerned the deposition and growth of bi-metallic clusters on a polymer surface [31] where the acceleration allowed one to study the very slow process of demixing of the two metals. Applications of this approach to various plasma processes should be possible as well. One effect that has already been studied is the creation of defects by ion impact. The main effect is trapping of clusters [54] case, at the defect locations which reduces their diffusion and limits cluster coalescence, see figure 4 above. In addition to the deposition of neutral atoms, the method also allows one to describe the impact of ions and the growth of charged clusters.

In concluding this section we mention that similar problems of rare events appear not only in MD simulations but also in statistical approaches such as Monte Carlo simulations. Some strategies are discussed in [55], where further references are given as well.

We now discuss the approaches listed in the right column of figure 3. The main idea of the coarse graining approaches is to perform an analysis of the different length and time scales that exist in a many-particle system—such as the plasma-surface system—driven out of equilibrium (one example could be a dense system of ions that is excited by an electric field pulse. Another example could be an ensemble of neutrals from the plasma that impact a solid surface and equilibrate there due to collisions with the lattice atoms). Depending on the type of excitation and on the properties of the system, the relaxation towards equilibrium typically proceeds in a number of steps. Even though this is a highly complex process in general, it is often possible to identify a sequence of relaxation processes or even a hierarchy.

A situation typical of gases and plasmas is sketched in table 1. Here, four relaxation stages are distinguished which are separated by the following characteristic time scales: the hydrodynamic time scale t_hyd, as well as the kinetic time scales t_rel (relaxation time) and τ_cor (correlation time). The latter is directly related to an equilibration of pair, g_ab (triple, g_abc, and higher) correlations, t_rel denotes the time necessary to establish a Maxwell (equilibrium) distribution f^EQ, and t_hyd is associated with the decay of density, velocity or temperature fluctuations or inhomogeneities or similar large-scale excitations. At each of these stages the system is adequately described by a specific set of quantities, and their dynamics is governed by specific equations: hydrodynamic equations (Stage III), kinetic equations (Stage II) and generalized non-Markovian kinetic equations or ab initio simulations (Stage I), respectively.

Table 1.  Characteristic scales and relaxation processes in correlated many-particle systems (schematic). Typical examples are the relaxation of electrons in a plasma following local ionization or excitation by a short electric field pulse or the thermalization of atoms from the plasma on a solid surface. Beginning at the initial time t₀, the evolution goes (from bottom to top) through several time stages and extends from small to larger length scales. This can be viewed as successive coarse graining, see figure 3. Accordingly, the relevant observables and concepts for a statistical description change. For explanations and details, see section 5. Adapted from [56].

Hierarchy of scales and relaxation processes in many-body systems
Time and length scales	Stage, Effects	Quantities	Theory
IV	Equilibrium	${n}_{a}^{\mathrm{EQ}}$, TEQ, pEQ	Equilibrium theory
		$p=p({n}_{1},{n}_{2}...,T)$	Equation of state
t > thyd		${n}_{a}={n}_{a}({n}_{1},{n}_{2}...,T)$	Mass action law
	Correlated equilibrium, or	p = pideal + pcor, etc
l > lhyd	Stationary nonequilibrium state	${n}_{a}({ \mathcal U })$, $T({ \mathcal U })$, $p({ \mathcal U })$	Quasi-equilibrium theory
	in an external field ${ \mathcal U }({\bf{R}})$
III	Hydrodynamic stage	na(Rt), ua(Rt), Ta(Rt)	Hydrodynamic equations
$t\in [{t}_{\mathrm{rel}},{t}_{\mathrm{hyd}}]$	Local equilibrium	${f}_{a}={f}_{a}^{\mathrm{EQ}}(n({\bf{R}}t),{\bf{u}}({\bf{R}}t),T({\bf{R}}t))$	Gas-dynamic equations
$l\in [{l}_{\mathrm{mfp}},{l}_{\mathrm{hyd}}]$			Reaction-diffusion eqs.
	Correlation corrections	$\ {n}_{a}({\bf{R}}t)={n}_{a}^{\mathrm{ideal}}({\bf{R}}t)+{n}_{a}^{\mathrm{cor}}({\bf{R}}t)\,$	Rate equations
		etc	Master equation
II	Kinetic stage	fa(pR, t)	Kinetic theory/
		fa(t0)	Relaxation time
t ∈ [tcor, trel]	Functional hypothesis		approximation
	Equilibrium correlations	${g}_{{ab}}={g}_{{ab}}^{\mathrm{EQ}}(\{f(t)\})\,$	Markov limit (M) +
l ∈ [lcor, lmfp]	Kinetic energy conservation	$=\,{g}_{{ab},0}^{{\rm{M}}}+{g}_{{ab},1}^{{\rm{M}}}+...$	Correlation corrections
I	Initial stage	gab(paRapbRb, t)	Generalized
	Initial correlations	${g}_{{ab}}({t}_{0})$	kinetic equations
t ∈ [t0, τcor]	Correlation buildup		Correlation time approx [74].
	Total energy conservation
l < lcor	Higher correlations	gabc, ${g}_{{abcd}},\,...$	First principle simulations

Even though each of these equations is an approximation to the full many-body equations, these equations are accurate within their respective stages and time scales. Thereby, the accuracy of these models is typically not limited by the equations themselves (e.g. resulting from the omission of higher-order terms) but by the parameters entering these models. For example, these parameters are the transport or hydrodynamic coefficients (determined by the distribution function) in hydrodynamic equations. In the case of kinetic equations, these parameters are cross sections or collision integrals. Of course, approximation schemes are used in practice for these parameters. But here we consider a different (although hypothetical) situation: if these input quantities would be known exactly, the underlying model equations would be exact, within their respective range of applicability⁴ . Of course, this requires the application of rigorous coarse graining procedures for the derivation of these equations, some examples and properties of which we will discuss in sections 5.3–5.5.

The existence of formally exact coarse grained equations is at the heart of this work and is explored in the remainder of this paper. In particular, our main idea is to obtain these exact input quantities for models from first principle MD simulations and thereby realize our goal to significantly extend the duration of first principle simulations. Since the coarse grained equations (Stages II–IV) emerge dynamically in the course of the equilibration, this novel method is called DFDM. We demonstrate it for simple examples in section 6.

An observation of interest for the following discussion is that the distribution function f_a(R, p, t) (a could denote ions or neutrals, in the examples above) is still time- and space-dependent, even after its equilibration. However, this dependence is only implicit and arises exclusively from the (slower) time evolution and weaker space dependence of the macroscopic fields entering the function

$$\begin{eqnarray}&&{f}_{a}({\bf{R}},{\bf{p}},t)={f}_{a}^{\mathrm{EQ}}[n({\bf{R}},t),{\bf{u}}({\bf{R}},t),T({\bf{R}},t)],\end{eqnarray} \tag{ 8 }$$

i.e., the density n, the mean velocity u and the mean energy (temperature T). By contrast, the momentum dependence is fixed by the Maxwellian form. This applies to the hydrodynamic Stage (III), see table 1, whereas the distribution carries, in addition, an explicit time dependence, on the kinetic stage (II). This concept is based on the functional hypothesis of Bogolyubov [57]. Close to the border of the two stages, i.e. for t slightly below ${t}_{\mathrm{rel}}$, the deviations of f from the distribution function f^EQ at equilibrium are small, and the time derivative of f (the collision integral, see equation (14)) can be approximated by the linear relation

$$\begin{eqnarray}&&{I}_{f}(t)\approx -\displaystyle \frac{1}{{t}_{\mathrm{rel}}}[f(t)-{f}^{\mathrm{EQ}}].\end{eqnarray} \tag{ 9 }$$

This is nothing but the familiar relaxation time approximation (or Bhatnagar–Gross–Krook (BGK) collision integral [58]), where the relaxation time is usually taken from equilibrium models of the system which have limited accuracy. In cases where ab initio simulation data are available, as we assume here, the MD simulation results can be used to extract an ab initio relaxation time which transforms equation (9) into an exact relation. This also requires to use the MD data for the equilibrium distribution that may be modified in the presence of correlations between the particles. This improved description also allows for a more accurate modeling of the hydrodynamic stage. We mention for the sake of completeness, that a similar crossover can also be studied between stages I and II. In that case, the pair correlation function reaches its equilibrium form ${g}^{\mathrm{EQ}}$ around the correlation time t ∼ τ_cor and retains only an implicit time dependence via the distribution functions g(t) = g^EQ[f(t)] (see table 1) thereafter. Similar to the relaxation time approximation, here one can use a correlation time approximation to describe the final approach to the equilibrium correlations [56, 59].

As a simple example, we consider the second order differential equation

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}\hat{A}}{{\rm{d}}{t}^{2}}+\gamma \displaystyle \frac{{\rm{d}}\hat{A}}{{\rm{d}}t}+\hat{k}(t)\hat{A}=0,\end{eqnarray} \tag{ 10 }$$

which represents e.g. a Newtonian equation of motion. Here, γ is a dissipation coefficient, and $\hat{k}(t)$ denotes a quickly fluctuating force constant or frequency. This could be a dust particle in a complex plasma that experiences collisions with plasma neutrals, ions or electrons. Another example could be a large molecule on a solid surface that undergoes collisions with the lattice atoms. A frequent situation is that $\hat{k}(t)\,=k+\delta \hat{k}$, where $k\equiv \langle \hat{k}\rangle$, and the brackets denote time averaging over the period of the fast oscillations or rapid dynamics of the light particles. Correspondingly, the dynamics of the variable $\hat{A}$ can be split into a slowly changing term and a rapidly oscillating contribution according to $\hat{A}(t)=A+\delta \hat{A}$, where $A=\langle \hat{A}\rangle$ obeys the equation of motion

$$\begin{eqnarray}&&\displaystyle \frac{{{\rm{d}}}^{2}A}{{\rm{d}}{t}^{2}}+\gamma \displaystyle \frac{{\rm{d}}A}{{\rm{d}}t}+k(t)A=-\langle \delta \hat{k}(t)\delta \hat{A}(t)\rangle ={I}_{A}.\end{eqnarray} \tag{ 11 }$$

This equation describes the 'coarse grained dynamics' where the fast 'random fluctuations' seem to be eliminated, and indeed, the left-hand side of this equation coincides with the original equation (10). However, the fast processes leave a trace on the slow dynamics via the new term on the right side, which is the correlation function of two fluctuating quantities. If the right-hand side is known, equation (11) will still be exact. In fact, it is not difficult to derive the equation of motion for $\langle \delta \hat{k}(t)\delta \hat{A}(t)\rangle$ by first finding the equation for $\delta \hat{A}(t)$, see e.g. [60, 61]. However, it is easy to see that this equation is not closed as well, and couples to a new quantity, $\langle \delta \hat{k}(t)\delta \hat{k}(t)\delta \hat{A}(t)\rangle$, which obeys its own equation of motion. Thus, an infinite hierarchy of equations emerges which is the direct consequence of the averaging procedure.

The common solution is to decouple this hierarchy by invoking additional assumptions on the fast dynamics. A common approximation is to assume ${\tau }_{{\rm{cor}}}^{\delta \hat{k}}/{t}_{A}\to 0$ which means that the correlation time of the rapid process (i.e. of the light particles) is vanishingly small compared to the characteristic time scale of A (the heavy particle). In that case, the term on the right side of equation (11) becomes a delta correlated random process, and one recovers a Langevin-type equation with a Gaussian white noise term.

The above model is representative for a large variety of problems containing multiple time scales. A similar situation appears in the case of multiple length scales, where a spatial averaging leads to a coarse grained description.

A further example for a coarse grained description is many-body dynamics in phase space, which is described by a generalized distribution function $\hat{N}({\bf{r}},{\bf{p}},t)$—the microscopic phase space density introduced by Klimontovich [60]. This distribution function obeys a mean field type equation

$$\begin{eqnarray}&&\left\{\displaystyle \frac{\partial }{\partial t}+{\bf{v}}\cdot {{\rm{\nabla }}}_{{\bf{r}}}+\hat{{\bf{F}}}\cdot {{\rm{\nabla }}}_{{\bf{p}}}\right\}\hat{N}({\bf{r}},{\bf{p}},t)=0,\end{eqnarray} \tag{ 12 }$$

where $\hat{{\bf{F}}}$ denotes the total force on the particles which includes all external and induced forces. In general, the force also contains rapid and slow contributions. Thus, we can write again $\hat{{\bf{F}}}={\bf{F}}+\delta \hat{{\bf{F}}}$, and the above coarse graining procedure can be repeated. In fact, this procedure leads to the well-known Boltzmann-type kinetic equations for the one-particle distribution function $f({\bf{r}},{\bf{p}},t)\equiv \langle \hat{N}({\bf{r}},{\bf{p}},t)\rangle$. It has the same form as equation (12) and reads

$$\begin{eqnarray}&&\left\{\displaystyle \frac{\partial }{\partial t}+{\bf{v}}\cdot {{\rm{\nabla }}}_{{\bf{r}}}+{\bf{F}}\cdot {{\rm{\nabla }}}_{{\bf{p}}}\right\}f({\bf{r}},{\bf{p}},t)={I}_{f}({\bf{r}},{\bf{p}},t).\end{eqnarray} \tag{ 13 }$$

However, it contains the collision integral I_f on the right-hand side in addition, which is directly determined by the correlation function

$$\begin{eqnarray}&&{I}_{f}({\bf{r}},{\bf{p}},t)=-\langle \delta \hat{{\bf{F}}}\cdot {{\rm{\nabla }}}_{{\bf{p}}}\delta \hat{N}\rangle \end{eqnarray} \tag{ 14 }$$

of the fluctuations. As before, the collision term is an unknown quantity, and one can derive an equation of motion for it. Again it includes a higher-order correlation function, and the whole system turns into an infinite hierarchy of equations. The standard solution is to use physical approximations for the choice of the collision integral I_f using e.g. the Boltzmann collision integral or improvements such as the Balescu–Lenard integral [56].

Now, let us discuss how the idea of the present paper can be applied here. In cases where we have first principle simulation data at our disposal, there exists a straightforward way how the hierarchies of equations discussed above can be decoupled. When performing MD simulations of a classical system, every observable can be computed. For example, it is possible to evaluate the right-hand side of equation (11) during a MD simulation. However, the result for $\delta \hat{k}(t)\delta \hat{A}(t)$ is random, depending on the choice of the initial conditions for the particle trajectories. The simple solution consists in running many independent MD simulations. Thereby, the initial conditions have to be chosen with a proper probability distribution given by the ensemble over which the averaging in the collision term I_A is performed. Then, the first principle MD result for I_A follows by averaging over M realizations according to

$$\begin{eqnarray}&&{I}_{A}(t)=\mathop{\mathrm{lim}}\limits_{M\to \infty }\displaystyle \frac{1}{M}\displaystyle \sum _{l=1}^{M}{(\delta \hat{k}(t)\delta \hat{A}(t))}^{(l)}.\end{eqnarray} \tag{ 15 }$$

This method is in principle exact for a classical system. The efficiency crucially depends on the cost of the MD simulation and on the number M of trajectories needed to obtain a converged result. At the same time, the fluctuation around the result (15) gives an estimate of the statistical uncertainty. A similar approach can be applied to the kinetic equation (13) and the collision integral I_f, equation (14).

Of course, the question arises what is gained by this approach compared to performing only MD simulations without resorting to the many-body equation (11) or the kinetic equation (13) at all. The point is that accurate MD simulations are typically substantially more costly than the latter approaches. Therefore, an advantageous compromise between accuracy and computational effort consists in performing MD simulations, for short time scales, and in continuing the simulation by solving a kinetic equation, for longer times. For the latter, the use of the MD input for the collision integral, as explained by equation (15), could yield a significant increase in accuracy and, it could eventually offer a way to extend first principle simulations to longer times. We return to this issue in section 6.

Notice that similar approaches exist also for quantum systems. An example is the 'Stochastic mean field' approach, see [62] and references therein. Instead of MD simulations, here time-dependent Hartree–Fock simulations are performed over which the averaging is carried out. The result of this procedure turned out to be very encouraging when compared to exact simulation results, at least for short and intermediate time scales [62].

Another application of coarse graining concepts is frequently used for particles in contact with a reservoir or 'bath'. Again this can be heavy dust particles in contact with lighter plasma particles or a macromolecule in a plasma or on a surface. The complete state of the system of N particles and N_B bath particles is described by the phase space distribution function F(R, P; R_B, P_B), where we use the compact notation ${\bf{R}}={{\bf{r}}}_{1},{{\bf{r}}}_{2}...{{\bf{r}}}_{N}$ and ${\bf{P}}={{\bf{p}}}_{1},{{\bf{p}}}_{2}...{{\bf{p}}}_{N}$ for the particles and similar notations for the bath particles and restrict ourselves to classical particles. It is easy to formulate the equations of motion for the whole system, but the solution for the N + N_B particles is extremely costly. More importantly, one is typically not interested in the details of the dynamics of the bath particles, where usually ${N}_{B}\gg N$.

Therefore, the standard procedure is to switch to a 'reduced' description, which resolves only the degrees of freedom of the system, by integrating over the bath parameters according to

$$\begin{eqnarray*}&&f({\bf{R}},{\bf{P}})\equiv \int {\rm{d}}{{\bf{R}}}_{B}{\rm{d}}{{\bf{P}}}_{B}F({\bf{R}},{\bf{P}};{{\bf{R}}}_{B},{{\bf{P}}}_{B}).\end{eqnarray*}$$

If the bath is in thermodynamic equilibrium at temperature T_B, this coarse graining procedure transforms the dynamics of the system of particles from the microcanonical ensemble into the canonical or grand-canonical ensemble. The corresponding equations of motion can then be solved by Langevin MD simulation or by using a Nose–Hoover thermostat [63].

Here, the main assumption is of course the equilibrium of the bath. This neglects the influence of the dynamics of the system particles on the bath particles which might be questionable, in particular at strong excitation conditions. In that case, an alternative strategy consists in performing short-time MD simulations including the bath dynamics by using different realizations of the bath over which an averaging is performed, in analogy to equation (15). This accurate procedure is computationally expensive and cannot be extended to long times. Thus, strategies can be developed, where one switches to the equilibrium description of the bath at times exceeding the thermalization time ${t}_{B}^{\mathrm{rel}}$ of the bath. A similar procedure is used in the simulation of disordered systems where an average over different realizations of the disorder is performed.

Such a strategy should be applicable to certain plasma-surface simulations such as scattering of plasma particles from a surface or diffusion of an adsorbate atom or molecule on a surface. In that case it can be justified to treat the surface as a 'bath' at sufficiently long time scales. In contrast, the initial time period of the interaction of a plasma particle with the surface requires a full dynamic treatment of the adsorbate and the surface atoms. A similar idea is realized in section 6.3.

Finally, we discuss an important approach to treat correlated many-particle systems which is based on the concept of reduced distribution functions. A typical example is the dynamics of electrons following a strong excitation and their subsequent thermalization, due to electron–electron collisions. The dynamics of a classical N-particle system can be treated by first principle MD simulations or, equivalently, by the N-particle distribution function ${F}_{N}({x}_{1},{x}_{2},\,\mathrm{...}\,{x}_{N})$ where x_i = (r_i, _pi) and F_N obeys the Liouville equation

$$\begin{eqnarray}{f}_{s}({x}_{1}...{x}_{s}) & = & \displaystyle \frac{N!}{(N-s)!}\\ & & \times \int {{\rm{d}}{x}}_{s+1}...{{\rm{d}}{x}}_{N}\,{F}_{N}({x}_{1}...{x}_{N}).\end{eqnarray} \tag{ 16 }$$

The equation of motion for f_s follows directly from the original equation for F_N by integration over the remaining variables, as in the definition (16) of f_s. The resulting equation for f_s is not closed but involves contributions from f_s+1 giving rise to a hierarchy of equations—the BBGKY hierarchy, which is discussed in detail e.g. in [56].

If the system is non-interacting, the s-particle distribution is a product of s single-particle factors, e.g. ${f}_{2}({x}_{1},{x}_{2})\,={f}_{1}({x}_{1}){f}_{1}({x}_{2})$. Here, we do not consider quantum exchange effects. In case of correlations between the particles, this relation is generalized to

$$\begin{eqnarray}&&{f}_{2}({x}_{1},{x}_{2})={f}_{1}({x}_{1}){f}_{1}({x}_{2})+{g}_{2}({x}_{1},{x}_{2}),\end{eqnarray} \tag{ 17 }$$

where g₂ is the pair correlation function. The BBGKY hierarchy is decoupled by invoking physically motivated approximations e.g. for the pair correlation function

$$\begin{eqnarray}&&{g}_{2}({x}_{1},{x}_{2},t)\to {g}_{2}^{{\rm{app}}}([{f}_{1}];{x}_{1},{x}_{2},t),\end{eqnarray} \tag{ 18 }$$

which is a given functional of the one-particle distribution function. As a result one obtains a closed equation for f₁—the kinetic equation—where the pair correlation function determines the collision integral ${I}_{f}={I}_{f}[{g}_{2}^{{\rm{app}}}]$, which coincides with the previous result given by equation (14).

In the spirit of the present paper, we underline that the resulting equation for f₁ is a dramatic simplification compared to the description of the full dynamics of all N particles, in particular when N is macroscopically large. This is the result of a very efficient coarse graining procedure. However, the quality of the result depends on the accuracy of the approximation, ${g}_{2}^{{\rm{app}}}([{f}_{1}];{x}_{1},{x}_{2},t)$. As discussed above, reliable approximations exist for limiting cases, e.g. when the problem contains small parameters (such as for weak interaction) and for long times $t\geqslant {\tau }_{{\rm{cor}}}$. However, for the initial time period, $0\leqslant t\lt {\tau }_{{\rm{cor}}}$, the standard (Markovian) results for ${g}_{2}^{{\rm{app}}}$ are known to fail, e.g. [56, 59]. Here, the concept of the present paper can be utilized again: perform MD simulations at short time scales, use their result to reconstruct the exact functions g₂([f₁], t), and extend this result to longer times using the kinetic equation for f₁.

Such a procedure has not yet been realized so far for kinetic equations. For this reason it is interesting to look at other examples where this concept has been tested.

Franke and Pehlke performed extensive DFT simulations of the diffusion of a 1,4-butaneditiol molecule on a gold surface [32]. They found the local adsorption energy minima of the molecule and then applied the nudged elastic band approach [64] to compute the transition rates between them. This allowed them to record the entire 'network' of atomic scale diffusion paths of the molecule on the surface. There is a large number of processes that are analyzed in [32]. Here, we are not concerned with the details of the associated diffusion hops, [32]—but focus on the coarse graining idea.

In order to connect the ab initio short-time diffusion simulations to the long time behavior Franke and Pehlke considered a master equation [32]

$$\begin{eqnarray}\displaystyle \frac{{{\rm{d}}{p}}_{i}(t)}{{\rm{d}}{t}} & = & \displaystyle \sum _{j\ne i}\{{{\rm{\Gamma }}}_{j\to i}\,{p}_{j}(t)-{{\rm{\Gamma }}}_{i\to j}\,{p}_{i}(t)\},\\ & & 0\leqslant {p}_{i}(t)\leqslant 1,\quad \displaystyle \sum _{i}{p}_{i}(t)=1,\end{eqnarray} \tag{ 19 }$$

where i is a multi-index numbering the configurations of the molecule which have a probability _pi(t), and Γ_a→b are the transition rates (probability per unit time) from state a to b. The first term on the right side of equation (19) describes processes which increase the probability to realize state i ('gain'), whereas the second term describes the analogous loss processes. For the computation of the transition rates the authors used standard TST [65],

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{a\to b}={\nu }_{a\to b}^{0}\,{{\rm{e}}}^{-{\rm{\Delta }}{E}_{a\to b}/{k}_{{\rm{B}}}T}.\end{eqnarray} \tag{ 20 }$$

Here, ${\rm{\Delta }}{E}_{a\to b}$ is the energy barrier for the transition between states a and b which is computed using DFT, and ν⁰ is the attempt frequency which is of the order of 10¹² s⁻¹ [32].

The authors of this reference consider two stages of the evolution: the initial stage, corresponding to stage I in table 1, and the asymptotic hydrodynamic state, corresponding to stage III. In stage I the dynamics depend strongly on the initial configuration of the molecule, and the diffusion retains a memory of the initial state being anisotropic. In contrast, one expects spatially isotropic motion of the molecule in stage III, where averaging over many initial configurations is assumed, since all memory of the initial state has been lost. Correspondingly, it is expected that the standard diffusion equation

$$\begin{eqnarray}&&\displaystyle \frac{\partial n({\bf{r}},t)}{\partial t}=D{\rm{\Delta }}n({\bf{r}},t),\end{eqnarray} \tag{ 21 }$$

holds with the well-known time-dependent solution of the initial value problem for the initial condition n(r, 0) = Nδ(r − r₀),

$$\begin{eqnarray}&&n({\bf{r}},t)=\displaystyle \frac{N}{4\pi {Dt}}\,{{\rm{e}}}^{-{({\bf{r}}-{{\bf{r}}}_{0})}^{2}/(4{Dt})}.\end{eqnarray} \tag{ 22 }$$

The solution (22) was recovered in [32] by mapping of the _pi(t) on the associated spatial coordinates of the center of mass of the molecule, giving rise to the space-dependent probability density P(r, t). P(r, t) = n(r, t)/N is proportional to the particle density n(r, t), where N is the total number of molecules (number of initial configurations). In the long time limit, an exponential temperature dependence of the diffusion coefficient results, i.e., an Arrhenius law

$$\begin{eqnarray}&&D(T)={a}_{0}^{2}\,\nu \,{{\rm{e}}}^{-{\rm{\Delta }}E/({k}_{{\rm{B}}}T)},\end{eqnarray} \tag{ 23 }$$

which made it possible to recover the effective attempt frequency ν and effective diffusion energy barrier ΔE. In (23), a₀ is the surface lattice constant.

To summarize, ab initio results for the elementary diffusion motions of a molecule on a surface have been obtained in this example. While it provides a complete microscopic picture of surface diffusion, this information is far too detailed for many purposes, in particular, for comparisons with measurements. In order to characterize the mobility of the molecules, the main interest concerns the long time behavior, which is governed by much simpler physics described by the classical diffusion equation (21). This reduced dynamics emerges dynamically during the course of the evolution of the system due to self-averaging effects. The main advantage of this approach is that the involved diffusion coefficient can be obtained exactly from ab initio DFT data instead of using standard approximate results from transport models.

At the same time the present combination of DFT and a master equation approach is in principle able to provide additional information, beyond that presented in [32]. First of all, it would be possible to establish the equilibration time scale t_rel, when the system reaches the isotropic diffusion regime (stage III in table 1). Furthermore, it should be possible to investigate the transient behavior being relevant at shorter time (stage II) as well, and this might give rise to a modified diffusion equation.

Finally, we note that similar master equation based approaches have been applied to the computation of chemical reaction rates, see [66, 67] and references therein. In the context of plasma-surface interaction, diffusion and chemical reaction rates in the presence of a plasma are of high interest. Therefore, the present approach might be useful to derive improved surface diffusion and reaction models that take into the influence of a plasma.

At this point, a first assessment of the concept to couple first principles data with an analytical model is given. The central question concerns, of course, the validity limits of the master equation. Let us summarize the corresponding requirements.

• (i)  
  First, it is known that TST is the basis for the transition rates (20), and it assumes that the surface is in thermal equilibrium.
• (ii)  
  All energy barriers have to be large compared to the thermal energy.
• (iii)  
  The transition probabilities in the master equation (19) depend linearly on the current occupation probabilities, i.e., the rates Γ_a→b are independent of all _pi(t).
• (iv)  
  The master equation is Markovian, i.e., the transition probabilities depend only on the current state of the system (no memory).
• (v)  
  The transition rates Γ_a→b in the master equation are time-independent.

Let us ask the question now, how important these conditions are and which of them can eventually be relaxed, if we would have first principle simulation data at our disposal. First of all, conditions (i) and (ii) can be easily dropped and the rates can be replaced by numerical simulation results, ${{\rm{\Gamma }}}_{a\to b}\longrightarrow {{\rm{\Gamma }}}_{a\to b}^{\mathrm{sim}}$. Second, the restrictions (iii)–(v) on the master equation can be dropped as well in favor of numerical results ${{\rm{\Gamma }}}_{a\to b}^{\mathrm{sim}}(\{{p}_{i}\},t)$, which are updated during the simulation. Obviously, these rates depend generally on the probabilities as well.

In fact, such a generalized version of the master equation which takes into account memory of the previous system states, is nothing but a generalized non-Markovian kinetic equation [56], which can be derived rigorously from the fundamental equations of many-body physics, such as the BBGKY-hierarchy, see section 5.5. The necessary requirements to such a generalized master equation are:

• A.  
  The microstates of the many-body system are mapped onto a complete set of events 'i' with well defined probabilities, i.e. ${\sum }_{i}{p}_{i}(t)=1$, and $0\leqslant {p}_{i}(t)\leqslant 1,$ for all i and all times.
• B.  
  There has to be a rigorous and numerically stable procedure how to identify these states and to assure ergodicity.
• C.  
  There has to be a consistent, stable and sufficiently accurate procedure how to determine the corresponding probabilities and transition rates 'on the fly' during a simulation.

A first example how to realize such a procedure is given in the following section 6.3.

Filinov et al [33, 68] presented a first application of the procedure outlined above to the adsorption dynamics and sticking probability of argon atoms on a platinum surface. They performed semi-classical MD simulations of the atom dynamics using ab initio pair potentials. These simulations yield the complete information on the particle trajectories x_i(t) = {r_i(t), _pi(t)} with $i=1...N$. These trajectories depend on the initial conditions such as the incident energy E_i(0) and angle θ_i(0), and position r_i(0). All observables of interest can in principle be computed from these trajectories (microstates) without resorting to additional approximations such as in TST (20).

Using the procedure described below, Filinov et al computed the sticking coefficient of argon atoms, as a function of time and obtained very good agreement with experiments [68], as is illustrated in figure 6. Interestingly, the results provide energy and angle resolved data for the probability that an impacting atom will be adsorbed at the surface (or be reflected). This is valuable input information for microscopic plasma simulations.

Zoom In Zoom Out Reset image size

Let us now return to the idea of this paper—how to extend these simulations to longer times. Even though for the present problem of computing the sticking probability MD simulations of $10...40$ ps duration are sufficient and no extensions are required, it is very instructive to analyze the potential of this scheme. First, it is clear that macroscopic properties such as the sticking coefficient do not require the complete microscopic information of all particle trajectories. Therefore, we attempt to map this microscopic information on a finite set of many-body states that are distinguished by the energy of the atoms.

The kinetic energy E_p contains two orthogonal contributions and reads

$$\begin{eqnarray}&&{E}_{{\rm{p}}}={E}_{{\rm{p}}}^{\perp }+{E}_{{\rm{p}}}^{\parallel },\quad {E}_{{\rm{p}}}^{\perp ,\parallel }=\displaystyle \frac{{({p}^{\perp ,\parallel })}^{2}}{2m}.\end{eqnarray} \tag{ 24a }$$

In addition, each particle moves in the potential landscape of all surface atoms V giving rise to the total energy

$$\begin{eqnarray}&&{E}_{t}({\boldsymbol{r}})={E}_{{\rm{p}}}({\boldsymbol{r}})+V({\boldsymbol{r}}).\end{eqnarray} \tag{ 25 }$$

All trajectories of the atoms can be uniquely classified by their energy at every time: (1) There are atoms with positive surface-normal energy, i.e., ${E}_{{\rm{p}}}^{\perp }+V\gt 0$. These particles are desorbed from the surface and their fraction is denoted as N_C (continuum states). (2) The remaining atoms have ${E}_{{\rm{p}}}^{\perp }+V\leqslant 0$ and belong into two fractions. The first fraction has a positive total energy, i.e., E_t > 0. These atoms can freely move across the surface and are denoted by N_Q (fraction of 'quasi-trapped' atoms). (3) The remaining atoms have a negative total energy (E_t < 0). That is they are trapped in local potential minima, and their fraction is denoted by N_T.

An example of such a trajectory is depicted in figure 7(a). There, an atom approaches the surface (being in a continuum state) and by colliding with the surface atoms, rapidly looses a large fraction of its energy, becoming trapped. Afterwards it gains again energy from another collision with the surface (thereby being transferred to a quasi-trapped state) until it is eventually desorbed (returning to a continuum state C).

Zoom In Zoom Out Reset image size

The main advantage of this mapping approach is that a large statistical ensemble of trajectories is available for sufficiently many atoms leading to an excellent accuracy of the results. We underline that this classification of all atoms into just three categories (three 'states') is unique, i.e., ${N}_{Q}(t)+{N}_{T}(t)+{N}_{C}(t)=1$, in agreement with the condition A listed in section 6.2. In fact, these fractions can be understood as occupation probabilities p_i with i = {C, Q, T} of the three distinct macro-states.

The three fractions of atoms are time-dependent and change during the time the atoms spend on the surface. The time dependence is governed by the system of rate equations [33]

$$\begin{eqnarray}&&{\dot{N}}_{Q}=-({T}_{{TQ}}+{T}_{{CQ}}){N}_{Q}+{T}_{{QT}}{N}_{T},\end{eqnarray} \tag{ 26a }$$

$$\begin{eqnarray}&&{\dot{N}}_{T}=-({T}_{{QT}}+{T}_{{CT}}){N}_{T}+{T}_{{TQ}}{N}_{Q},\end{eqnarray} \tag{ 26b }$$

$$\begin{eqnarray}&&{\dot{N}}_{C}=-({\dot{N}}_{Q}+{\dot{N}}_{T})={T}_{{CT}}{N}_{T}+{T}_{{CQ}}{N}_{Q},\end{eqnarray} \tag{ 26c }$$

that is, in fact, just an example of the master equation discussed in section 6.2, and the transition rates T_αβ [α, β = {C, Q, T}] are just the coefficients ${{\rm{\Gamma }}}_{\alpha \to \beta }$ occurring in the latter.

In the terms used in table 1, this system of rate equations corresponds to the methods listed for the hydrodynamic stage (stage III). But, nothing prevents us of course from using these equations also for earlier time scales, i.e., for the kinetic and initial stages (stage II and I, respectively). As discussed in section 6.2, to be applicable at earlier times, we have to permit a dependence of the rates on time and on the individual probabilities in that case, i.e., ${T}_{\alpha \beta }={T}_{\alpha \beta }({N}_{C},{N}_{Q},{N}_{T};t)$. Furthermore, it was shown in [33] that these rates depend sensitively on the energy distribution function of the gas atoms $F({E}^{\perp },{E}^{\parallel };t)$, which themselves evolve with the time duration atoms spend on the surface before they are desorbed.

In [33] it was demonstrated in detail, how the complete information from the MD trajectories and it was demonstrated how the MD data can be used to explicitly reconstruct the rates T_αβ for all times. This means that the conditions B and C of section 6.2 are also realized. Furthermore, the analysis revealed that the gas atom–surface interaction proceeds in two stages. At first, the energy distribution functions thermalize within the relaxation time ${t}_{\mathrm{rel}}\approx (20...40)$ ps, which depends on the energy and angle of incidence of the atom. During this period of time the transition rates also reach their equilibrium form, T_αβ(t) → T^EQ_αβ, and remain practically constant for times t > t_rel. This can be seen in figure 7(b). This first time interval corresponds to the stages I and II in table 1. The subsequent temporal evolution corresponds to stage III and turns out to be accurately described by the rate equations (26a)–(26c) with constant transition rates.

This behavior is verified by a comparison of MD simulation data with the analytical solutions of the rate equations (26a)–(26c). The corresponding results are shown in figure 7(c). While the analytical results differ qualitatively from the MD simulation results, for short times, the analytical fractions of trapped and quasi-trapped atoms practically coincide with the MD data, for times larger than about 5t₀ ≈ 33 ps. This is approximately the time, where the transition rates have saturated (figure 7(b)).

In other words, the set of three occupation probabilities (fractions) ${N}_{C},{N}_{Q}$ and N_T is sufficient to capture the entire sticking and desorption properties of the gas atoms, for times $t\geqslant {t}_{\mathrm{rel}}$. These three collective variables have emerged dynamically during the temporal evolution and the associated coarse graining dynamics. Thus, the system of rate equations is sufficient to describe the mean adsorption/desorption dynamics for a sufficiently large ensemble of atoms at longer times. The solution of these equations is computationally cheap and allows one to propagate the system, in principle, to arbitrarily long times. This could be of relevance for experiments where the surface properties change in time, e.g. due to an AC field or in the course of continuous sputter deposition. Note that the transition rates ${T}_{\alpha \beta }^{\mathrm{EQ}}$ obtained from the MD simulations are not approximations but have first principle quality, in principle.