Concurrent Characterization of Surface Diffusion and Intermixing of Ge on Si: A Classical Molecular Dynamics Study

The surface diffusion and intermixing of Ge ad‐atoms over Si (001) 2 × 1 substrates using classical molecular dynamics (CMD) simulations are characterized here. Several interatomic potentials, parametrizations, and parameter mixing rules are contemplated. A novel simulation scheme is devised to characterize the effective frequency of surface diffusion and intermixing events overcoming the inherent difficulties related to their interdependency in heteroepitaxial systems. The effective energy barriers of these events encompass different atomistic mechanisms weighted by their occurrence probabilities. The overall description of surface diffusion and intermixing based on Stillinger–Weber (SW) potential is in agreement with ab initio calculations and experimental observations, though some atomistic details differ. This study is extended to Si(001) substrates with stressed Ge monolayers grown on top. It is found that Ge ad‐atom dynamics is accelerated with respect to the case of the pure Si substrate and that diffusion across dimer rows is mainly mediated by the atomic exchange of the Ge ad‐atom with a Ge atom on the surface.


Introduction
The progress in Si nanotechnology involves the growth of SiGe layers with high Ge content or pure Ge on Si, to fabricate faster CMOS transistors, [1] to create quantum dots for photonic devices, [2] or to improve optic absorption in photodetectors. [3][4][5] The 4.2% lattice mismatch between the two materials generates compressive stress in the grown layer, which has some benefits (enhanced carrier mobility) but makes the fabrication of perfect flat layers a tremendous challenge. Depending on the deposition conditions and Ge content in the film, the compressive stress DOI: 10.1002/adts.202200848 relaxes by atomic interdiffusion, formation of dislocations, or generation of islands. [6] Experimental and theoretical studies investigated heteroepitaxy of SiGe and Ge films on Si, [7][8][9][10][11][12][13] but some aspects are still unknown, especially those regarding atomic interactions at the early stages of deposition. Their direct experimental characterization is hampered by the small size and time scales involved. To shed some light on this matter, atomistic physical modeling becomes a useful tool.
Epitaxial growth can be described using lattice-gas stochastic atomistic models, in which atoms are randomly deposited in a periodic distribution of adsorption positions, and their evolution is simulated using kinetic Monte Carlo (KMC) [14] or continuum [15] methods. To develop predictive and efficient computational models within these simulation schemes, all possible atomistic events and their associated probabilities (energies) must be defined a priori. This information is extracted from more fundamental atomistic simulation techniques, such as ab initio or classical molecular dynamics (CMD). Ab initio methods have been used to determine adsorption positions of Ge ad-atoms on the Si surface, [8,[16][17][18] and to characterize Ge/Si intermixing [9,19] and surface diffusion [8,17,18,20] events. Due to the heavy computational cost of this method, ab initio results come mainly from static calculations where particular atomic paths are assumed a priori, and energy barriers are extracted using the nudged elastic band method (NEBM). [21] CMD has the advantage of allowing direct simulation of dynamics in larger systems than ab initio calculations, although results depend on empirical interatomic potentials. Using CMD techniques, various local minima for Ge adsorption on the Si surface were determined and some insight into possible diffusion paths was offered, [10,22] but no full dynamic characterization of atomistic surface diffusion and intermixing mechanisms has been performed. The lack of CMD results on this matter could be related to the inherent difficulties for such a characterization in heterosystems, as we shall expose later. There are also no dynamical studies about the behavior of Ge ad-atoms over Ge layers already deposited on Si. The description of atomic mechanisms provided by an empirical potential, although it may not be as accurate as in ab initio, is essential to interpret macroscopic structures and results obtained in extensive CMD simulations where ab initio methods are not affordable. [23] Here we present a comprehensive CMD study aimed to characterize surface diffusion and intermixing of Ge ad-atoms on Si(001) 2 × 1 surfaces (the most technologically relevant). To carry on this study, we selected the most suitable interaction potential to describe Si-Ge mixtures, and we developed a novel CMD simulation scheme to extract the energetics of Ge ad-atom surface diffusion and intermixing processes. We explored pure Si substrates as well as Si substrates with Ge perfect layers deposited on top.

Selection of the Interatomic Potentials for Si-Ge
The CMD simulation technique consists of the numerical resolution of Newton's equations of motion for a system of atoms. [24] Atoms are modeled as point-like particles, and interactions among them are evaluated from empirical potentials, which include parameters fitted to experimental data and/or ab initio calculations. A careful selection of empirical potentials describing Si, Ge, and their mixtures, is essential to have meaningful CMD results. Recently, machine learning interatomic potentials (MLIPs) have emerged as a powerful tool for materials modeling, being able to predict energy and atomic forces with nearly ab initio accuracy and orders of magnitude faster. [25] MLIPs have been developed for pure Si and Ge, [26,27] but not for Si-Ge mixtures. For this reason, our study resorts to traditional and extensively tested potentials for semiconductors, such as Stillinger-Weber (SW) [28] and Tersoff, [29] for which several parametrizations and mixing rules for Si-Ge systems are available in the literature.
For pure Si or Ge, most studies using Tersoff potential employ the parameters originally defined in ref. [29]. In contrast, a wide variety of parametrizations have been proposed for the SW potential. We tested two of them for Si [28,30] and five for Ge. [31][32][33][34][35] For Si-Ge systems, the parameters to describe atomic interactions are determined according to specific mixing rules involving the pure Si and Ge parameters. For the Tersoff potential, mixing rules were already provided in the author's original work. [29] For the SW potential, different approaches were considered. [36][37][38] We calculated some basic properties of Si, Ge, and SiGe to be utilized as figures of merit: the lattice parameter (l o ), cohesive energy (E coh ), melting temperature (T m ), and elastic properties such as the elastic constants, bulk modulus (K), shear modulus (G), and Poisson ratio ( ). l 0 is key to properly describe the lattice mismatch on Si-Ge heterostructures, and E coh is the depth of the potential well at the bonding distance and influences atom dynamics (i.e., how easily atoms can escape from their lattice sites and diffuse). A correct T m assures the agreement between simulated and experimental temperatures, which is necessary to calculate diffusion coefficients and migration energies in dynamic simulations. Elastic properties determine how the material behaves under stress, as it occurs in heteroepitaxial films. For cubic crystals, elastic constants C ij are reduced to three independent values: C 11 = C 22 = C 33 , C 12 = C 13 = C 23 , and C 44 = C 55 = C 66 (C ij = 0 for the remaining elastic constants). [39] In all simulations, we used the open-source MD software package LAMMPS. [40] To calculate l o and E coh (at 0 K), we employed a system consisting of 4096 atoms, arranged in 8 × 8 × 8 unit cells of the diamond lattice with X, Y, and Z cell axes lying along the [100], [010], and [001] directions, respectively. Periodic boundary conditions (PBCs) were applied in all spatial directions. We used Ref. [29] Ref. [28] Ref. [30] l a) From ref. [45]; b) from ref. [46]; c) from ref. [47]; d) from ref. [48]; and e) from ref. [49].
the conjugate-gradient method to relax the system with different lattice parameters, and we estimated l o and E coh from the minimum energy-zero pressure point. [41] To determine T m at zero external pressure, we applied the twophase molecular dynamics method. [42] Our simulation system consisted of 32 000 atoms arranged in 40 × 10 × 10 unit cells of the diamond lattice with X, Y, and Z cell axes lying along the [100], [010], and [001] directions, respectively. PBCs were applied in Y and Z-axes, and the lattice parameter in these two directions was scaled to account for thermal expansion at a temperature close to the expected T m . Free boundary conditions were applied in the X-axis to allow the system to relax at zero external pressure in that direction. Initially, half of the simulation cell in the X-direction was melted by heating it to a temperature much higher than the expected T m . Once melting was completed, the liquid was cooled down to a temperature T * close to the expected value of T m . The other half of the simulation cell was kept in the crystalline phase, and it was also equilibrated at T * . Then, the system evolved freely in the NVE ensemble. If T * < T m , some of the liquid recrystallizes liberating the corresponding latent heat to the system and increasing the temperature within the simulation cell. If T * > T m , some of the crystal melts absorbing the corresponding latent heat leading to a decrease of the temperature in the cell. In this way, the temperature of the simulated system automatically tends to T m .
To evaluate the elastic constants C 11 , C 12 , and C 44 at 300 K, we created deformations in the simulation cell in different directions and calculated the change in the stress tensor. [43] We employed a system consisting of 512 atoms arranged in a diamond lattice of size 4 × 4 × 4 unit cells with X, Y, and Z axes oriented along the main crystallographic directions. PBCs were applied in all spatial directions. Elastic parameters K, G, and were obtained by applying the following expressions  Ref. [29] Ref. [31] Ref. [32] Ref. [33] Ref. [34] Ref. [35] l a) From ref. [50]; b) from ref. [51]; c) from ref. [52]; d) from ref. [48]; and e) from ref. [49]. a well-known fact, so temperature re-scaling is commonly performed in order to compare real with simulated temperatures in CMD. [44] The original parametrization [28] of the SW potential underestimates E coh . The parametrization of Balamane et al. [30] adjusts the value of E coh at the expense of a small 6% deviation in T m . None of the potentials and parametrizations give a set of values for the elastic properties clearly better than the rest and, in all cases, values could be considered acceptable. In view of these results, among the SW parametrizations for Si, we chose that of Balamane et al. [30] for our study. Table 2 summarizes the calculated and experimental values for l o , E coh , T m and elastic properties for Ge. We only evaluated elastic properties with SW parametrizations of refs. [32,34,35] that provide l o , E coh , and T m in good agreement with experimental values. All of them predict fair values for the calculated parameters, but the parametrization of Posselt et al. [34] gives a better value for T m (overestimates the experimental value by 7%, analogously to the case of the SW parametrization for Si of Balamane et al. [30] ) For this reason, we chose the parametrization of Posselt et al. [34] to describe Ge using SW. Tersoff potential gives adequate values for l o and E coh , but T m is again highly overestimated, with a value very close to the one predicted for Si melting, contrary to experiments where there is a difference of 476 K between the two T m values. This fact precludes the possibility of temperature rescaling in simulations of Si and Ge mixtures, which constitutes a serious drawback for the Tersoff potential. Nevertheless, Tersoff gives better values for the elastic properties. Table 3 shows the experimental and calculated values for l o , E coh , T m and elastic properties for a SiGe system. In our simulations, SiGe is a perfect Si 0.5 Ge 0.5 alloy in the diamond crystal structure, where each Si atom is surrounded by four Ge atoms, and vice versa. Tersoff potential provides fair values for l o and E coh , it describes elastic properties better than SW, but T m is again highly overestimated (and in fact, it is not even between the Ge and Si melting temperatures). With respect to the SW potential, we applied the previously cited mixing rules [36][37][38] to the selected parameter sets for Si [30] and Ge [34] and the evaluated properties are similar for all of them. Values obtained for l o and E coh compare fairly well with experiments, while T m is slightly overestimated (between 6% and 9%), as it also happens for SW models of Si and Ge. In view of these results, we chose Gilmer and Grabow mixing rules, [36] since they are more consistent with the standard combination rules in binary systems, [53] and they have been commonly employed in the literature. [37,38,[54][55][56]

Methodology to Evaluate Surface Diffusion and Intermixing Parameters in Heterosystems
Ad-atom surface diffusion and intermixing mechanisms in heterosystems are interdependent, and this complicates their characterization. To monitor surface intermixing events, a large CMD sample with a great amount of ad-atoms would be required to have good statistics; but, before surface intermixing occurs, adatoms diffuse and encounter each other giving rise to the formation of ad-dimers that preclude surface intermixing. To characterize ad-atom surface diffusion, long CMD simulation times would be required to extract a clear slope in the mean square displacement (MSD) curves; however, the ad-atom frequently interdiffuses or interacts with defects that spontaneously form on the substrate surface, and diffusion stops. To overcome these difficulties, we propose a novel CMD simulation scheme (depicted in Figure 1), based on a collective analysis of individual simulations (CAIS). Each simulation with a single Ge ad-atom follows its diffusion in a small cell until an event that alters system dynamics occurs. From the collective treatment of data from all individual simulations run in parallel, we concurrently extract frequencies and activation energies of surface diffusion and intermixing.
In our study, the simulation cell consisted of 1200 Si atoms arranged in a diamond lattice of size 3l T × 10 where l T is the Si equilibrium lattice constant at temperature T and zero applied stress. Cell axes X, Y, and Z were oriented along the [100], [011], and [0-11] directions, respectively, with PBCs in Yand Z-directions and free conditions in X. The whole system was composed of twelve Si monolayers in the X-direction. To simulate the bulk behavior, atoms in the four bottom layers were held fixed, while the rest of the atoms were free to move. Atoms in the middle four layers were used to create a thermal bath (thermostat) to control the system temperature. In the topmost layer, atoms were arranged along Z to form the 2 × 1 dimer rows typical of the Si(001) surface reconstruction. Initially, the atoms in the thermostat layers were given random velocities in a Maxwell-Boltzmann distribution corresponding to the simulated temperature. To equilibrate the substrate, the cell was relaxed in the NVE ensemble for 50 000 steps, with the time step fixed at 1 fs, scaling the atom velocities within the thermostat layers every 2000 steps at the desired temperature.
In every simulation, a single Ge atom was deposited on a random (Y, Z) position at a height of 10 Å above a clean Si(001) 2 × 1 surface, with an initial velocity of −0.1 nm ps −1 in the Xdirection (this condition assures that the impact of the Ge atom does not damage the Si surface). The system evolved in the NVE ensemble at the target temperature, with 1 fs timesteps, and the Ge ad-atom MSD i was computed. Atom positions were averaged for 10 000 time steps to remove thermal vibrations. Each averaged configuration was evaluated to assess whether a surface intermixing event occurred or the Si surface degraded (dimers broke down). The use of averaged instead of relaxed instantaneous configurations filters unsuccessful transition attempts. [58] We considered that surface intermixing occurred when the following two conditions were fulfilled: 1) the distance of the Ge adatom to a Si lattice site is smaller than 0.4 Å; and 2) the shortest distance of any of the Si atoms to Si lattice sites is greater than 2.2 Å. When such an event occurred, we recorded the time elapsed from the deposition of the Ge atom (surface intermixing time t i ) and stopped the simulation. We considered that surface degradation occurred when at least two Si substrate atoms were displaced from perfect lattice sites more than 2 Å, which happened more frequently at higher temperatures. When this condition occurred, we stopped the simulation because surface defects act as trapping centers for the Ge ad-atom, thus affecting the calculation of MSD i and t i . Hundreds of parallel simulations were run to achieve enough statistics.
The total amount of simulated individual surface intermixing events, N 0 , is statistically equivalent to N 0 Ge ad-atoms simultaneously deposited at time t = 0 on a large Si substrate where distances among them are long enough to assure that they do not interact before surface intermixing occurs. To extract the surface intermixing parameters, for every simulated temperature we sorted the surface intermixing times obtained in individual simulations t i in increasing order: t i1 was the shortest surface intermixing time, t i2 the second shortest time, and so on. After t i1 , one of the Ge ad-atoms intermixed, so the total amount of remaining Ge ad-atoms was N(t i1 ) = N 0 − 1. Again, that after t i2 another Ge adatom intermixed, being N(t i2 ) = N 0 − 2 the number of remaining Ge ad-atoms, and so on. As an example, in Figure 2 we represent N(t) in logarithmic scale as a function of time t for the simulated temperature of 1350 K using the SW potential. As predicted by transition-state theory for first-order processes, N(t) follows an exponential decay, characterized by a well-defined reaction rate i (T) at each temperature T according to: This scheme was repeated for several temperatures to extract the surface intermixing activation energy E i from the Arrhenius equation where i0 is the frequency prefactor and k B the Boltzmann constant.
Since diffusion processes are equivalent in all the individual simulations, the total MSD of a single Ge ad-atom that had not been disrupted by surface intermixing or degradation events is equivalent to the cumulative summation of the MSD i during the time corresponding to the summation of all simulation times. As an example, Figure 3 represents the cumulative MSD obtained using the SW potential for a temperature of 1350 K. Even though the cumulative MSD is composed by the MSD i extracted from 1200 parallel CMD runs, it shows a clear slope which remains constant along all the simulations (the validity of this approach is also addressed in Section 4.2). The MSD slope is related to the temperature-dependent diffusion coefficient, D(T), according to the Einstein relation in an n-dimensional Brownian motion We extracted D(T) from the cumulative summation of all MSD i , thus improving statistics in its calculation compared to that obtained from short individual CMD runs. D(T) follows an Arrhenius dependence, where D 0 is the diffusivity prefactor and E d the activation energy for diffusion. To quantitatively compare the frequencies of diffusion and surface intermixing processes, we extracted the diffusion jump frequency d (T) from D(T)  (3) and (6).
where is the jump length for ad-atom diffusion and d 0 = 2n 2 D 0 the jump frequency prefactor. By splitting the MSD into parallel (Z-direction) and perpendicular (Y-direction) components (see Figure 3), we characterized Ge ad-atom surface diffusion frequencies along and across dimer rows, respectively (n = 1).

Dynamics of Ge Ad-Atoms on a Si(001) Substrate
To characterize Ge ad-atom surface diffusion and intermixing on a Si(001) 2 × 1 surface, we carried out CMD simulations at temperatures between 0.60T Si m,SW and 0.90T Si m,SW for the SW potential, and between 0.60T Si m,T and 0.80T Si m,T for the Tersoff potential (T Si m,SW = 1795 K and T Si m,T = 2444 K according to the calculations shown in Section 2). We selected these specific ranges for each potential to maintain a compromise between a fast system dynamics (to accelerate the simulations as much as possible) and a sufficiently stable substrate surface (dimers break down at high temperatures). Figure 4 shows the Arrhenius plot for the frequencies of surface intermixing, i (T), and surface diffusion, d (T), split into parallel (diffusion along dimer rows) and perpendicular (diffusion across dimer rows) components, obtained with SW and Tersoff potentials. Table 4 compares the extracted activation energies E d∕∕ , E d⟂ , and E i with those reported by other authors using ab initio calculations, [8,9,[17][18][19][20] lattice-gas techniques, [22] and experimental scanning tunneling microscopy (STM) measurements. [11] We found that the activation energy for the Ge ad-atom surface intermixing E i predicted by the SW potential (1.36 eV) is significantly lower than that obtained with Tersoff (1.97 eV). Using ab initio and NEBM, Ge/Si intermixing activation energies between 0.8 and 1.6 eV were reported for particular surface intermixing pathways. [9,19] Thus, SW potential predictions compare better with ab initio results than Tersoff potential does. We observed that most Ge ad-atoms intermixed into the topmost Si substrate (dimers) layer, though they were occasionally adsorbed in the second layer. SW potential predicts that around 90% of the intermixed Ge ad-atoms end up in the dimers layer at high temperatures, and close to 100% at low temperatures. With Tersoff potential, 85% of intermixed Ge ad-atoms end up in the topmost surface layer, being the remaining 15% in the second substrate layer, independently of temperature. The amount of Ge atoms adsorbed in deeper layers was insignificant for any of the two potentials.
Frequencies for diffusion are higher than for surface intermixing in both potentials, although frequencies for diffusion across dimers and surface intermixing approach at high temperatures. The most frequent mechanism is parallel diffusion at all temperatures. We obtained similar values for E d⟂ and E d∕∕ in the case of Tersoff, but E d∕∕ is significantly lower than E d⟂ for SW. Our results with the SW potential are in good agreement with the previously cited theoretical calculations and experiments (see Table 4). However, Tersoff potential overestimates E d∕∕ and does not predict the high anisotropy between the parallel and perpendicular components of diffusion.
From the inspection of the pathways followed by Ge ad-atoms in our CMD simulations, we identified the adsorption sites, that is, the local minimum energy positions that the Ge atom visits as it moves on the Si(001) 2 × 1 reconstructed surface. Some of these positions (called A, B, B 0 , C, D, H, M, P, and T) were already reported for the adsorption of Si [59][60][61] and Ge ad-atoms [10,17,18,22] on a Si(001) 2 × 1 surface. Other positions (that we call D Top , N, and Z) have not been reported previously, to our knowledge. All adsorption sites are represented in Figure 5. The six positions A, Table 4. Activation energies (in eV) for surface intermixing, and diffusion along and across dimer rows.

Present work
Other studies SW Tersoff Ab initio Lattice-gas STM Refs. [17,18] Ref. [8] Ref. [19] Ref. [9] Ref. [20] Ref. [22] Ref. [  At the topmost layer in the X-direction, Si atoms were paired to form the dimers of the 2 × 1 reconstruction of the (001) surface. Si atoms at the four bottom layers in the X-direction were held fixed at their bulk positions. To obtain the minimum energy configuration for each adsorption site at 0 K, one Ge ad-atom was placed on the adsorption position, and the system was relaxed using the conjugate-gradient minimization procedure.
The calculated energies of Ge ad-atom adsorption sites are represented in Table 5, along with results from other authors. Energies were measured with respect to the most stable position in every case. For SW the global minimum corresponds to site B in the trough and, among the positions on top of the dimer rows, the most stable is M, being H, P, and D sites only slightly higher in energy. Our results are consistent with the previous study of Roland and Gilmer with SW, [10] although energy values are not exactly the same, probably because these authors employed alternative SW parametrizations for Si and Ge. Calculations with the Tersoff potential also predict that site B is the most stable. [22] Ab initio [8,[16][17][18] and hybrid quantum and molecular mechanics (QM/MM) [62] calculations identify the lowest energy state of the Ge ad-atom on the dimer rows at site M. Although for the SW potential the global minimum energy position is in the trough (at site B), it predicts that sites M and H are the most stable configurations on top of the dimer rows. Taking into account all the results presented so far, we only consider the SW potential for the rest of our simulation study, since SW potential describes the thermodynamical properties of the material, and the anisotropy of surface diffusion in better agreement with ab initio, and Tersoff potential does not improve the description of adsorption sites with respect to SW.
The overall description of Ge surface diffusion and intermixing given by SW is in agreement with ab initio, but the predicted adsorption sites are not. To elucidate the contributions of different atomistic pathways to the overall behavior described by MD calculations, we analyze the mechanisms underlying Ge surface Ref. [22] Ref. [16] Refs. [17,18] Ref. [8] Ref.
[62]  diffusion and intermixing and determine the analogies and differences with ab initio observations. With respect to surface intermixing, we identified two main pathways (sketched in Figure 6): in the first one (I1), the Ge ad-atom, initially located on top of one dimer at site D, pushes one of the Si atoms in the dimer toward the trough, and occupies the initial position of the Si atom, which in turn ends up in site B 0 ; in the second pathway (I2), the Ge ad-atom, which initially occupies site B in the trough, moves downward to position B 0 , and pushes up a neighboring dimer Si atom that moves to site D. We observed that mechanism I1 is more frequent. Unless the deposited Ge atom lay directly over the dimer rows, this mechanism requires a previous diffusion from the most stable configuration (B) to position D. The energy diagrams determined by NEBM associated with both surface intermixing pathways (Figure 6c,d) reveal energy barriers of 1.37 and 1.50 eV for mechanisms I1 and I2, respectively, from the most stable configuration (B). These energies are very close to the effective intermixing energy obtained in dynamic simulations (Table 4). Surface intermixing pathways proposed in ab initio studies are similar to those identified in our simulations, but the starting position for Ge surface intermixing is site M instead of sites D or B. [9,19] We observed that Ge atoms in the through may occasionally become incorporated into the Si subsurface layer, by directly exchanging positions with Si atoms in such a layer, although this is an infrequent event specially at low temperatures. With respect to surface diffusion, Ge ad-atoms deposited on a Si(001) 2 × 1 substrate follow different pathways, that we disentangled into those contributing to parallel or perpendicular diffusion. We determined three main diffusion paths parallel to dimer rows and two paths perpendicular to dimer rows, which are depicted in Figure 7. Corresponding energy diagrams evaluated using NEBM are plotted in Figure 8.   T1 and site B for T2), are 0.77 eV to jump from the dimer row to the trough (T1), and 1.06 eV to jump from the trough onto the dimer row, breaking-up the Si dimer below (T2). Note that sites C and D behave as saddle points for parallel diffusion, but they are potential wells for perpendicular diffusion, allowing for switching from perpendicular to parallel mechanisms.
Although the deposited Ge atom was initialized at a random (Y, Z) position above the substrate, we noticed that it tends to bind to atoms located in the dimer rows (as they are the first atoms that it encounters in its way down to the substrate), becoming frequently adsorbed in M position. Ge ad-atom surface diffusion is a combination of all the parallel and perpendicular paths described above. Ge ad-atoms adsorbed in the dimer rows typically diffuse on top of them (parallel diffusion) for a few jumps, and then they hop into the trough (perpendicular diffusion). Ge ad-atoms diffuse along the trough (parallel diffusion) until they occasionally jump again on top of the dimer rows (perpendicular diffusion). The effective energy barriers for parallel and perpendicular diffusion obtained from the Arrhenius plots of Figure 4 are not coincident with the diffusion barriers of any particular path, but they encompass all these mechanisms weighted by their occurrence probability.
Static calculations using ab initio techniques [8,18,20] described surface diffusion pathways that agree with our MD description of parallel paths over the dimer row (P2 and P3), and also with the perpendicular path between contiguous dimers (T1), being the lowest energy path that along the top of dimer rows through site D (P2). However, ab initio studies did not offer a parallel path along the dimer rows (P1), nor a perpendicular path above dimers (T2), probably due to site B not being among the lowest energy minima sites.

Dynamics of Ge Ad-Atoms on Ge Monolayers Grown over Si(001)
Due to the lattice mismatch between Si and Ge, as the heteroepitaxy of Ge on Si proceeds, layers of stressed Ge form on the Si substrate. To study the behavior of Ge ad-atoms deposited over previously grown Ge layers on Si, we carried out an analogous study to the one presented so far, but now considering that atoms in the top layers are Ge instead of Si. The topmost Ge layer has also the 2 × 1 reconstruction. We refer to this substrate with m Ge monolayers on Si(001) as a Ge m /Si. We limited our study to one (Ge 1 /Si) and two (Ge 2 /Si) monolayers of Ge over Si(001) because, according to experiments and theoretical calculations, at least 2-3 wetting layers of Ge on Si can be grown before 3D islands nucleate. [63] Ge ad-atom adsorption sites in these systems are almost coincident with the ones found on the Si(001) 2 × 1 surface, except for Z and A positions. Site Z is not observed as a local energy minimum on the Ge 1,2 /Si surfaces, neither A in the Ge 2 /Si surface. Energies of the Ge ad-atom adsorption positions are shown in Table 6. Results obtained for the pure Si substrate are indi- Table 6. Ge ad-atom energies, in eV, at adsorption positions on pure Si and Ge 1,2 /Si reconstructed (001) substrates, relative to the most stable configuration in each case. cated again for the sake of comparison. Adsorption energies on the Ge 1 /Si and Ge 2 /Si surfaces are almost identical except for B 0 and N which are higher for the Ge 2 /Si surface. As it occurred for Si, site B is also the global energy minimum on the Ge 1,2 /Si surfaces. Relative energy values are, in general, slightly higher in the case of the Ge 1,2 /Si surfaces than in Si, except for C position, whose relative energy is significantly lower on the Ge 1,2 /Si surfaces. Being C a quasi-saddle point for diffusion along the trough, the corresponding energy barrier lowers in the Ge 1,2 /Si substrates case, as we shall see. The most stable position on top of the dimer rows is again site M, and it is almost degenerate in energy with P and H adsorption sites.
To characterize diffusion of Ge ad-atoms on Ge 1,2 /Si surfaces, we performed CMD simulations at temperatures between 0.45T Si m,SW and 0.65T Si m,SW for the Ge 1 /Si surface, and between 0.4T Si m,SW and 0.6T Si m,SW for the Ge 2 /Si surface. These temperatures are lower than the ones considered for the pure Si substrate because the Ge 1,2 /Si surfaces degrade more easily. This is due to the stress in Ge monolayers and the lower melting temperature of Ge respect to Si. In Ge 1,2 /Si systems, we cannot properly speak Figure 9. Arrhenius plot of the frequencies obtained in CMD simulations with the SW potential for parallel and perpendicular diffusion on pure Si (squares), Ge 1 /Si (circles), and Ge 2 /Si (triangles) substrates. Solid lines are best fits to Equation (6). Orange and blue color corresponds to diffusion along and across dimer rows, respectively. Cross marks represent the calculations on Ge 1 /Si substrate when stopping conditions were included. Table 7. Parameters obtained from the Arrhenius plot of Figure 9 for different substrates.
Si about surface intermixing events, because when the Ge ad-atom exchanges its position with a Ge atom in the top substrate monolayer, the diffusing species is again Ge, so the simulation system is equivalent to the one we had prior to the exchange. This system allows the comparison between the cumulative summation of MSD i from individual simulations halted after Ge exchanges, with the total MSD obtained from a long simulation that was not stopped when Ge exchange occurred. Diffusion frequencies of Ge on a Ge 1 /Si surface from both types of simulations are analogous at all simulated temperatures, as shown in Figure 9 (crosses and circles). This fact validates our methodology for the extraction of the MSD slope from the cumulative summation of individual MSD i . Results for Ge ad-atom frequencies of parallel and perpendicular diffusion hops on Ge 1,2 /Si surfaces are displayed in Figure 9. Table 7 shows the frequency prefactors as well as the activation energies. Previous results on a pure Si substrate are also included to facilitate the comparison. Diffusion frequencies are higher on Ge 1,2 /Si than on pure Si, and the activation energies are lower, in agreement with experimental observations. [11,64] This indicates that stressed Ge monolayers accelerate atom dynamics, as predicted in theoretical calculations. [64][65][66] Our extracted frequency prefactors 0 are close to 10 13 s −1 , which is the standard value used in most of the theoretical and experimental studies. [67][68][69] Frequencies for parallel diffusion hops in Ge 1 /Si and Ge 2 /Si surfaces are similar, while frequencies for perpendicular diffusion hops are almost one order of magnitude larger in the Ge 2 /Si case. Since activation energies are almost the same, this difference is due to the larger frequency prefactor for perpendicular diffusion in the Ge 2 /Si surface. This increase in the frequency prefactor is probably due to an increase in entropy [64] associated with a greater number of possible hopping paths when there are two Ge monolayers, as we will show later.
In Ge 1,2 /Si substrates, we observed that Ge ad-atom diffusion across dimer rows is mostly mediated by the exchange of a Ge adatom with a Ge atom of the reconstructed topmost surface layer. For the Ge 1 /Si surface this atomic exchange (Figure 10a) occurs in the same way as the Ge surface intermixing mechanism I2 in the Si substrate, but in this case it displaces a Ge atom of the dimer. This mechanism also occurs in the case of the Ge 2 /Si surface, but it involves a higher energy barrier (1.05 eV). For this surface, the main exchange mechanism becomes the one depicted in Figure 10b, in which the Ge ad-atom displaces diagonally through an intermediate local minimum (site BP) to position P. Figure 11 shows the energy diagrams for the main Ge adatom diffusion mechanisms along and across dimer rows on the  (Figure 11b) include T1 and T2 pathways, similar to those found on the pure Si substrate. Energy barriers calculated from the most stable position in each path (site M for T1 and site B for T2), are 0.60 eV to jump from the dimer row to the trough (T1), and 1.12 eV to jump from the trough onto the dimer row, breaking-up the Ge dimer below (T2). We also included pathways T3(Ge 1 /Si) and T3(Ge 2 /Si), which correspond to the exchange mediated mechanisms shown in Figure 10a,b, respectively.

Discussion and Conclusions
We found that the SW potential, using parametrizations of Balamane et al. [30] for Si and Posselt et al. [34] for Ge, with the mixing rules proposed by Gilmer and Grabow, [36] resulted in the best compromise for reproducing basic bulk material properties (l o , E coh , T m , and elastic parameters) and atomic surface properties. Tersoff potential has some drawbacks: the temperature scale is uncertain because it does not correctly predict the melting temperatures of Si and Ge, and the Ge ad-atom surface diffusion anisotropy is not captured.
Previous studies of Ge surface diffusion and intermixing on Si using ab-initio and CMD are limited to static calculations. To our knowledge, no full dynamics studies have been carried out, probably because these two mechanisms are interdependent and it was not possible to characterize them individually using the standard methods of homoepitaxial systems. The diffusing Ge ad-atom frequently exchanges its position with a Si surface atom, and the resulting system differs from the one prior to that event. Consequently the Ge ad-atom diffusion cannot be followed for simulation times long enough to extract meaningful MSD. To obtain surface intermixing time constants, we initially attempted to follow the evolution of hundreds Ge ad-atoms simultaneously deposited on the surface of a large simulation cell. Since the surface intermixing frequency is several orders of magnitude lower than that of surface diffusion, Ge atoms (separated even several nanometers apart) were likely to find each other before surface intermixing occurred. Thus, our initial attempt always led to the formation of Ge-Ge ad-dimers, which altered the frequency of surface intermixing events. We devised a novel methodology (CAIS) to characterize Ge ad-atom surface diffusion and surface intermixing that overcomes this interdependency in heterosystems: our simulation scheme separates surface diffusion and intermixing mechanisms by running one single ad-atom at a time on a small cell for a limited amount of time; the collective treatment of many individual simulations (that can be run in parallel as they are independent to each other) provide statistically significant data. Even in cells with a single Ge ad-atom, the free diffusion time was short because we stopped the simulation when we identified surface intermixing or degradation events, as they introduced artifacts on the evaluation of the MSD. Statistics, in this case, were improved by the cumulative summation of the MSD i extracted from all the short-time simulations. CAIS differs from "parallel replica" [58] and "temperature-accelerated dynamics" [70] methods, as these use parallelism to extend or extrapolate the time scale for simulations at low temperatures, while our scheme employs parallelism to improve statistics in the determination of parameters in short time-scale simulations.
Here CAIS is employed to characterize surface diffusion and intermixing of Ge ad-atoms on Si-Ge systems with SW potential, but this scheme can be transferred to any heterosystem and any interatomic potential. It can also be used to characterize the effective transition barriers between any two different states A and B in homo-or hetero-systems. To adapt our simulation scheme to other systems, the details for the identification of A and B must be specifically stated. A and B can be defined as specific local minima or in a more general way, as we have done in this work (A: Ge ad-atom diffusing on the Si surface, and B: Ge atom in a surface position-intermixed Ge atom). In addition, it is important www.advancedsciencenews.com www.advtheorysimul.com to identify when the system evolves to a state different from those to be characterized because from that instant the dynamics correspond to a different system. In our simulation, we identify when surface degradation occurred to halt the simulation. By extracting times at which the transition between A and B occurs, the associated time constant and the effective activation energy and prefactor for the transition can be evaluated assuming the validity of transition-state theory and first-order kinetics.
We obtained effective activation energies and frequency prefactors for the diffusion and surface intermixing mechanisms that came from the own dynamics of the system without resorting to static calculations, where particular diffusion and surface intermixing pathways are assumed a priori. The extracted parameters encompass all the possible atomic mechanisms, each one weighted by its own occurrence probability. In agreement with static ab initio calculations, we observed that surface diffusion is strongly anisotropic. We obtained that diffusion parallel to the dimer rows is two orders of magnitude more frequent than perpendicular diffusion, being surface intermixing the least frequent mechanism.
By analyzing the positions visited by the Ge ad-atom during its free movement on the Si surface, we extensively searched for ad-atom adsorption sites. This allowed us to disentangle surface diffusion and intermixing atomic mechanisms, to determine energy barriers for particular atomic pathways using NEBM, and to compare them with ab initio results. Ab initio calculations correctly predict the dimer buckling experimentally observed on the Si (001) reconstructed surface, [68,71,72] and identify the most stable Ge adsorption site at M on the dimer rows. [8,[16][17][18]62] However, SW potential does not reproduce dimer buckling and predicts that site B in the trough has the lowest energy. According to SW potential, the Ge ad-atom at site M is bonded to two dimers from the same row and to a second-layer Si atom, while the analysis of the charge density distribution in ab initio calculations does not show a directional bond to this second-layer Si atom. [17,18,20] This is a limitation of the SW potential (and also of the Tersoff potential), related to the way empirical potentials describe bonding geometry, [17,18] favoring the tetrahedral one and giving repulsive three-body interactions for this configuration. Nevertheless, SW potential determines that the M site is the most stable configuration on top of the dimer rows, resulting in diffusion pathways and energy barriers with respect to site M in qualitative agreement with ab initio calculations. The obtained effective activation energies for surface diffusion and intermixing are comparable to the energy barriers extracted from ab initio calculations, which indicates that SW potential gives a good overall description of the system, although atomic details may differ.
We explored the energetics of Ge ad-atom on Ge monolayers grown over Si. We found that all dynamic mechanisms are faster than on a pure Si substrate, with a higher anisotropy between parallel and perpendicular diffusion for the Ge 1 /Si surface than for the Ge 2 /Si surface. We noticed that perpendicular diffusion is mediated by the atomic exchange between the Ge ad-atom and a Ge atom on the topmost reconstructed surface layer. Stressed Ge layers decrease energy barriers and increase diffusion frequencies. The difference in the activation energies on Si and Ge 1,2 /Si substrates must be taken into account in higher-level simulations, where successive Ge layers are deposited as heteroepitaxial growth proceeds.
Our CMD study of Ge surface diffusion and intermixing on Si(001) 2 × 1 substrates is relevant because we validate the empirical potential and provide effective parameters derived from the own dynamics of the system. This fundamental study could have been done, at least partially, using ab initio molecular-dynamics techniques, as simulation cells are relatively small and simulation times relatively short. However, ab initio molecular dynamics would be restricted to a very narrow temperature range (at temperatures high enough to accelerate dynamics but low enough to ensure surface stability), much narrower than that accessible with CMD. The error would be large when extracting the slope in the Arrhenius plot from few points concentrated in a very narrow temperature window. To simulate the dynamic deposition of many large Ge layers, one must resort to CMD. To be able to scale deposition rate to the surface diffusion rate it is important to understand how the empirical potential describe the system and how surface diffusion and intermixing occur according to it. To achieve realistic deposition rates and sizes, higher level modeling techniques (KMC or continuous methods) must be used; they require effective diffusion and surface intermixing parameters (rather than very detailed atomic mechanisms) aimed to describe the global system behavior.
In summary, we showed that SW potential with the parametrization of Balamane et al. [30] for Si and Posselt et al. [34] for Ge using Gilmer-Grabow mixing rules [36] provides a good overall description of the Ge ad-atom diffusion and surface intermixing processes on a Si(001) reconstructed surface. Although the SW potential does not exactly reproduce all the atomistic details of ab initio calculations, the understanding of the atomic description provided by it is essential to interpret the mechanisms that lead to different structures (perfect wetting layers, islands, dislocations, etc.) in deposition simulations reachable by CMD approaches. We developed a new methodology to concurrently extract surface diffusion and intermixing effective parameters; we applied it to the study of Ge on Si, but it can be used in any heteroepitaxial system. With this methodology, we extracted effective barriers in dynamic simulations without assuming preestablished pathways. Our exploration of Ge ad-atom diffusion on Ge monolayers grown over Si showed that dynamics are accelerated with respect to the case of the pure Si substrate, and that diffusion across dimer rows involves the atomic exchange of the Ge ad-atom with a Ge atom on the surface.