Modeling of Al and Ga Droplet Nucleation during Droplet Epitaxy or Droplet Etching

The temperature dependent density of Al and Ga droplets deposited on AlGaAs with molecular beam epitaxy is studied theoretically. Such droplets are important for applications in quantum information technology and can be functionalized e.g., by droplet epitaxy or droplet etching for the self-assembled generation of quantum emitters. After an estimation based on a scaling analysis, the droplet densities are simulated using first a mean-field rate model and second a kinetic Monte Carlo (KMC) simulation basing on an atomistic representation of the mobile adatoms. The modeling of droplet nucleation with a very high surface activity of the adatoms and ultra-low droplet densities down to 5 × 106 cm−2 is highly demanding in particular for the KMC simulation. Both models consider two material related model parameters, the energy barrier ES for surface diffusion of free adatoms and the energy barrier EE for escape of atoms from droplets. The rate model quantitatively reproduces the droplet densities with ES = 0.19 eV, EE = 1.71 eV for Al droplets and ES = 0.115 eV for Ga droplets. For Ga, the values of EE are temperature dependent indicating the relevance of additional processes. Interestingly, the critical nucleus size depends on deposition time, which conflicts with the assumptions of the scaling model. Using a multiscale KMC algorithm to substantially shorten the computation times, Al droplets up to 460 °C on a 7500 × 7500 simulation field and Ga droplets up to 550 °C are simulated. The results show a very good agreement with the experiments using ES = 0.19 eV, EE = 1.44 eV for Al, and ES = 0.115 eV, EE = 1.24 eV (T≤ 300 °C) or EE = 1.24 + 0.06 (T[°C] − 300)/100 eV (T>300 °C) for Ga. The deviating EE is attributed to a re-nucleation effect that is not considered in the mean-field assumption of the rate model.


Introduction
Semiconductor quantum dots (QDs) are established as quantum emitters and represent essential building blocks for quantum information technology [1][2][3]. For their fabrication, the use of self-assembly techniques is a powerful approach providing versatile semiconductor nanostructures by molecular beam epitaxy (MBE) [4]. Two major paths are pursuit relying on different mechanisms for the self-assembly. In the Stranski-Krastanov growth, epitaxial layers of different lattice constant are deposited and the driving force for self-assembled formation of e.g., InAs/GaAs QDs or Ge/Si QDs is the minimization of the strain energy [5][6][7]. As a drawback of this technique, the resulting QDs are substantially strained. More flexible regarding the choice of the materials are droplet-based techniques [8] such as the droplet epitaxy [9][10][11][12][13] or droplet etching [14][15][16][17][18]. Here, the driving force for self-assembly is the minimization of the solid and liquid surface and interface energies. For droplet epitaxy, often Ga droplets are deposited on AlGaAs surfaces and subsequently recrystallized to form GaAs QDs [19,20] or quantum rings [21]. For droplet etching, usually Al droplets are deposited on AlGaAs, transformed into low-density nanoholes during post-growth annealing, and subsequently filled with GaAs to form GaAs QDs [22,23] or quantum dot molecules [24].
For both methods, droplet epitaxy and droplet etching, the desired feature density is equal to the density of the initial droplets. Thus, the precise control of the droplet density is essential for the creation of quantum emitters with well-defined properties. For instance, low density QDs [25] are relevant as single-photon sources for quantum information technology and high density dots [26] for device applications such as lasers or solar cells. The central process parameter controlling the droplet density is the substrate temperature T. In an earlier study [12], we have discussed the T-dependent density of Ga droplets, but only for low temperatures up to 300°C and using only basic approaches for modeling. Here, we study also Al droplets and significantly expand the temperature range to consider technologically relevant processes. Furthermore, we compare three different approaches for the modeling of the experimental data. In this sense, the present manuscript follows two goals: first, the better understanding of the experimental droplet densities and the determination of relevant material parameters and, second, the evaluation of different concepts for droplet modeling.
Popular approaches for modeling of nucleation on surfaces during crystal growth are coupled mean-field rate equations [12,[27][28][29] and kinetic Monte Carlo (KMC) simulations [30][31][32][33]. Rate equations describe the time-dependent average density of objects on the surface such as mobile adatoms (monomers) and clusters of various size. In principle, every cluster size requires the calculation of an individual equation which can be very timeconsuming for system containing large clusters. To overcome this issue, a critical cluster size i is introduced which represents the smallest stable cluster size [27][28][29]. Clusters with size below i are unstable and the monomer escape rate is larger than the rate of attachment. Now, all stable clusters can be averaged and do not require individual equations.
As an advantage in comparison to a mean-file rate model, a kinetic Monte Carlo simulation is based on an atomistic representation of the surface adatoms and allows a modeling also of the droplet size distribution. On the other side, the computation time of a KMC simulation substantially increases for large simulation fields (required for the present low droplet densities) and for a high rate of monomer diffusion processes. This limits the maximal process temperature accessible by this method. Starting with a conventional KMC approach, we have simulated the nucleation of Al droplets up to a temperature T of 300°C. By using a multiscale KMC algorithm [32], we have reduced the computation time by a factor of 30,000 in comparison to a conventional KMC approach and achieved very good agreement with experimental Al droplet densities up to a temperature of 460°C on a 7500 × 7500 simulation field and with Ga droplets up to 550°C.

Experimental Droplet Density
The sample fabrication using solid-source MBE on (001) GaAs wafers is described previously [34,35]. In brief, after deposition of an atomically flat AlGaAs layer, the As flux is reduced and the droplet material is deposited. For Ga, the material flux to the surface is F = 0.8 mL/s, the deposition time = 4 s, and the resulting coverage of droplet material on the substrate surface is Ft = 3.2 mL. For Al, the parameters are F = 0.4 mL/s, deposition time = 2.5 s, and coverage Ft = 1.0 mL. The droplets are formed in Volmer-Weber growth mode [36] driven by the minimization of the liquid and solid surface energies and of the liquid-solid interface energy. For higher process temperatures [35], the deposited droplets drill nanoholes into the substrate, which is called droplet etching [14][15][16][17][18]. Relevant processes are here the diffusion of As from the crystalline substrate into the droplet material and the removal of the droplet material by spreading over the substrate surface [37]. Importantly, since every droplet forms a nanohole [17], the initial droplet density can be determined from nanohole data. After growth, the surfaces are characterized using atomic force microscopy (AFM) and the droplet densities N D are determined. Figure 1 shows values of N D after deposition of Ga and Al as function of the temperature T during droplet deposition. The data are taken from previous publications [34,35].

Computation
The computations are performed on personal computers. We have started the simulations with the rate equations model using the interpreted programming language Python3. To reduce the computation time we have tested Julia (about two time faster than Python) and finally switched to the compiled languages Delphi (about 10 time faster than Python) and C++ (about 1.5 times faster than Delphi). The final computations of both the rate equations model as well as the kinetic Monte Carlo simulation are done parallel using two codes one in Delphi and one in C++. Typical computation times of the rate model range from a few minutes up to a few hours and depend on the simulated process temperature. The computation time is determined by the accuracy of the time interval for the numerical integration, where process conditions with a smaller monomer density require a smaller time interval. On the other side, the computation times of the KMC simulations are much longer and depend crucially on the simulated process temperature, since at a higher T the rate of surface activity increases and the reduced droplet densities require larger simulation fields. A conventional KMC simulation for a process temperature of T = 300°C using a 1400 × 1400 simulation field takes about one week, which we consider as a limit for a reasonable parameterization. A multiscale KMC approach (see below) computes the same process substantially faster in about 20 s. Here, the limit for Al droplets are simulations for T = 460°C using a 7500 × 7500 simulation field and for Ga droplets T = 550°C, which require more than one week.

Scaling Analysis
We start the analysis of the experimental droplet density using classical nucleation theory [27][28][29], which predicts the density of stable three-dimensional clusters by a scaling law For complete condensation of three-dimensional clusters the scaling energy becomes E = p(E S + E i /i), with the energy barrier E S for surface diffusion, p = i/(i + 2.5), the critical cluster size i, and the energy E i of a critical cluster.
A scaling analysis of the experimental droplet densities in Figure 1 yields a fitted slope of E = 0.74 eV for Al (T = 280. . . 570°C) and E = 0.35 eV for Ga droplets. The reduction of N D at T > 570°C for Al droplets was earlier [12] related to droplet coarsening by Ostwald ripening [29,38].
Earlier flux dependent data [12] indicate p = 0.5 and thus i = 2.5 for a low T = 200°C and Ga droplets. We will show below that this value of i is not valid for the higher temperatures discussed here. For a comparison with the rate model described below, an energy barrier E E for escape of atoms from a cluster is introduced with E E = 2E i /i [12]. Now, the scaling energy becomes E = p(E S + E E /2).

Mean-Field Rate Equations
Mean-field rate equations represent a common approach to model nucleation processes during crystal growth [12,[27][28][29]. The present model considers the average densities of three types of objects on the growing surface: (1) Mobile adatoms (monomers) with density N 1 . (2) Small droplet-like immobile clusters of density N s that are composed of s atoms and potentially unstable, i.e., below the critical cluster size i. (3) Droplets with density N D and an average volume of s atoms that are above the critical cluster size and, thus, stable. Droplets and clusters are approximated by a hemispherical shape.
Furthermore, three classes of processes are considered: (1) Arrival of the impinging material with a flux F on the surface and formation of monomers. (2) Attachment of mobile monomers to either other monomers, to clusters with size s, or to droplets. The corresponding rates are R A,1 , R A,s , and R A,D , respectively. In this scheme, collisions between mobile adatoms represent nucleation events. (3) Escape of monomers from clusters or droplets with the respective rates R E,s and R E,D .
The driving force for monomer mobility and their attachment to other objects is the surface diffusion where the hopping frequency is given by the diffusion coefficient [39], k B Boltzmann's constant, and h Planck's constant. We have to consider also the capture numbers [27,28,40] σ s , σ D = 2π(r/λ)k 1 (r/λ)/k 0 (r/λ), which reflects the depletion of the monomer density around monomers, clusters, or droplets, where k 0 , k 1 are modified Bessel functions, r is the radius, and λ −2 = (F/D) + 2σ 1 N 1 + σ D N D is the surface diffusion length [40]. We assume that clusters and droplets are hemispherically shaped with a radius of r = 3 3s/(2π). This allows to calculate the monomer attachment rates to other monomers (s = 1) and clusters with size s as R A,s = N 1 σ s N s D, and to droplets as Escape of atoms from clusters with size s is considered with rate R E,s = 2πrζ N s D E and from droplets with rate , ζ = exp(r c /r) describes the enhancement of the vapor pressure for small droplets due to the Gibbs-Thomson effect, r c = 2γV mol /(N A k B T), γ is the surface tension (0.89 N/m for Al and 0.67 N/m for Ga), V mol the molar volume (10.0 × 10 −6 m 3 /mol for Al and 11.8 × 10 −6 m 3 /mol for Ga), and N A the Avogadro constant.
In the following a mean-field model basing on a set of coupled rate-equations is introduced which describes the droplet nucleation and allows to calculate the time dependence of the average droplet density. Since the critical cluster size i depends on the model parameters and is not know initially, a maximum possible unstable cluster size j > i is introduced. The first equation describes the monomer density, which is balanced by the flux of impinging adatoms, monomer attachment to other monomers, clusters, and droplets, as well as by escape from clusters and droplets: The density of clusters with size of s = 2. . . j − 1 atoms evolves as: and that of clusters with s = j as: Finally, the droplet density follows: For the model calculations, the j + 1 individual rate equations (Equations (2)- (5)) are solved numerically (see Section 3). Figure 2 show an example for the simulated time evolution of various quantities during Al droplet deposition at F = 0.4 mL/s, deposition time 2.5 s, T = 500°C, E S = 0.19 eV, and E E = 1.71 eV. The monomer density N 1 initially increases, followed by a decrease due to nucleation and attachment events. The droplet density N D saturates very fast. We note that this is not the case at low temperatures. Both the surface diffusion length λ = 300. . . 350 nm and the capture number σ D = 1.4. . . 3.1 increase slightly during deposition. As an interesting result, the critical cluster size i increases from 10 to 14. The critical cluster is the smallest cluster size s which fulfills the condition R A,s ≥ R E,s . In the numerical calculations, j must be always larger then i. In addition, finally the droplet volume which increases nearly linear with time according to V = s Ft/N D . The rate model has two free model parameters, the energy barrier E S for surface diffusion of monomers and the energy barrier E E for escape of monomers from clusters or droplets. We start with the parameterization of the model for the Al droplet density using the experimental flux F = 0.4 mL/s and a deposition time of 2.5 s, which yields a droplet material coverage Ft = 1.0 mL. In a first step, we have calculated N D for varied E S and E E at a constant T = 500°C. Results are shown in Figure 3a and demonstrate that N D depends on E S and even stronger on E E . Interestingly, the expected decrease of N D with decreasing E S is observed only for very large E E = 10 eV, i.e., irreversible aggregation. Here, the mobile adatoms perform a faster migration at a lower E S and attach preferred to existing islands instead of nucleating new ones. In contrast to that, a lower E N ≤ 2.0 eV yields two regimes separated by a maximum of N D . We assume that E E influences N D by two competing effects. First, the nucleation rate will be reduced due to the escape of atoms from clusters with size below i. Second, the monomer density is increased by escape which yields an enhanced nucleation rate. Therefore, the regime at lower E S is diffusion controlled and N D decreases with decreasing E S , whereas the regime at higher E S is escape controlled and N D decreases with increasing E S .   A comparison with droplet densities obtained from AFM measurements establish that a range of pairs E S , E E provides agreement between model calculations and experiments (Figure 3a). A summary of the possible pairs is plotted in Figure 3b. The same procedure for T = 280°C gives a second range of pairs E S and E E (Figure 3b). As a key point, there is only one pair E S = 0.19 eV, E E = 1.71 eV which yields agreement for both temperatures and represents an unambiguous determination of both energy values. Figure 1 demonstrates the very good reproduction of the experimental droplet densities by the rate model using these parameters.
As a further interesting point, the rate model reproduces also the experimental reduction of N D at T > 570°C. This reduction was earlier [12] attributed to droplet coarsening by Ostwald ripening [29,38]. Since the present model does not consider droplet coarsening, the relevance of Ostwald ripening for droplet formation is now questionable.
In addition, Figure 3c shows that the temperature dependence of the critical cluster size i = 3. . . 39 at the end of deposition follows approximately i ∝ exp(T). Since a time and T-dependent i is not considered in the scaling model (Equation (1)), a scaling analysis of the droplet densities establishes as an only rough approximation. Nevertheless, from E = p(E S + E E /2) or p = E/(E S + E E /2), an effective critical cluster size i = 6 can be estimated.
The parameterization of the model for the Ga droplet density is more complex. Here, the experimental conditions are F = 0.8 mL/s and a deposition time = 4 s. Again, pairs E S , E E with agreement between calculated and experiments N D are determined (Figure 4a). However, in contrast to the Al droplets above, for the Ga droplet density there is no pair which provides agreement for the whole temperature range. Only the low temperature regime from 140. . . 300°C is reproduced by a common pair E S = 0.115 eV, E E = 1.51 eV and calculations using these parameters agree well with the AFM data in this regime (Figure 4b).
The behavior at higher temperatures cannot be reproduced using these energies (Figure 4b). We assume that above 300°C either E S or E E depends on T. Since a variation of E S has an only small effect, probably E E is the relevant T-dependent parameter. Assuming a constant E S = 0.115 eV, Figure 4c shows values of E E which provide agreement with the experimental N D at varied T. Figure 4c clearly establishes that above 300°C the experimental N D cannot be reproduced assuming a constant E E . This is also confirmed by the example in Figure 4b which agrees only at T = 500°C using E S = 0.115 eV and a constant E E = 1.71 eV. As a consequence, for a reproduction of the whole temperature range, the escape energy must be T-dependent with E E = 1.51 eV for T ≤ 300°C and E E = 1.51 + 0.1 (T[°C] − 300)/100 eV for T > 300°C. This choice for E E allows a very good reproduction of the Ga droplet density (Figures 1 and 4b).   An analysis of the scaling behavior makes sense only for T ≤ 300°C where E E is constant. Here, p = E/(E S + E E /2) yields an effective critical cluster size of about i = 2.

Kinetic Monte Carlo Simulations
The kinetic Monte Carlo (KMC) simulation model considers in an atomistic picture the activity of mobile atoms (monomers) on the growing surface and their nucleation to droplets (size s ≥ 2 atoms). The model assumes that only monomers are mobile. The surface is represented as a square simulation field with m n = m x × m y elements and cyclic boundary conditions. The positions of the respective monomers and droplets as well a the individual droplet sizes define the status of the growth process.
We start with a conventional KMC approach (called MC1), where individual atoms can arrive on the surface with flux F, hop to a nearest-neighbor surface site with rate D = ν exp[−E S /(k B T)], or escape from a hemispherically shaped droplet composed of s atoms with rate R E,s = 2πrζν exp[−E E /(k B T)], r is the droplet radius, and ζ the Gibbs-Thomson enhancement as is already defined above for the rate model. In the MC1 model, the various rates sum up to a total activity rate: with the monomer number n 1 and the numbers n s of droplets with size s. This allows the random selection of the next process as well as the calculation of the time interval up to the next event dt = 1/R tot . The arrival of an atom on the surface can result either in the formation of a new monomer, a nucleation event by a direct hit to another monomer, or in droplet growth by a direct hit. Monomer diffusion can cause a site-change, a nucleation event, or the attachment to a droplet. An escape from a droplet yields a new monomer at a random angle and distance r + 2 from the droplet center. A scheme of a single loop of the Monte Carlo model MC1 is shown in Figure 5.

RndSel atom
Delete atom

RndSel site
New status New monomer

Nucleation Attachment
New dimer Drop growth For a reasonable statistics we consider a minimum of 20 droplets per simulation field. This requires for Al droplets, e.g., a 1400 × 1400 field for T = 300°C, 4000 × 4000 for T = 400°C, and 10,000 × 10,000 for T = 500°C. Simulation runs for Al droplets at T = 300°C are performed on a 1400 × 1400 grid using the energy values E S = 0.19 eV and E E = 1.71 eV as determined by the rate model above. The simulated droplet density of 2.3 × 10 10 cm −2 ± 7% is significantly larger compared to the results of the rate model. This disagreement demonstrates that the determined energy values are model-dependent. A discussion will be given below.
We note that a precise parameterization or simulations for higher temperatures are not accessible with the MC1 model due to the very long computation time. On a personal computer, a single simulation run for T = 300°C using the above parameters takes about one week. This is caused by the large simulation fields in combination with the very high number of diffusion processes due to the small E S . Simulations for a higher temperature with even larger simulation fields and a much higher D would require a substantially longer simulation time.
To speed-up the Monte Carlo simulation, in the modified simulation model MC2 we replace monomer diffusion via a large number of nearest-neighbour hops by fewer jumps over longer distances. This approximation is established by DeVita et al. as a multiscale kinetic Monte Carlo algorithm [32]. In detail, Equation (6) indicates that for D m n F + ∑ s n s R E,s an individual monomer performs much more diffusion events compared to arrival plus escape. To reduce the number of simulation steps, all diffusion events within the time interval τ I = 1/(m n F + ∑ s n s R E,s ) up to the occurrence of the next arrival or escape are summarized. The monomer surface diffusion length is λ = √ τD according to the Einstein relation, with the diffusion time τ. In the time interval τ I the diffusing monomer travels a distance d I = √ τ I D. Or, in other words, the rate R I = 1/τ I = m n F + ∑ s n s R E,s for travelling a distance d I by diffusion is given by the time interval up to the next arrival or escape. If there are other objects at a distance smaller than d I , the diffusion can result in either a displacement of the monomer, in nucleation by collision with another monomer, or in attachment to a droplet. The probability p for a collision with another object depends on the circular segment r/(πd) covered by the object, where d is the distance and r is the radius of the object. This gives the nucleation rate R N,ij = rD/(πd 3 ij ) at which monomer i collides with a second monomer j, with the distance d ij between both. Accordingly, the rate of attachment to a droplet k is R A,ik = rD/(πd 3 ik ). Replacing the rate n 1 D for monomer diffusion in Equation (6) by the above modifications yields for the total activity rate: As for the MC1 model, R tot is used for the random selection of the next process and the calculation of dt = 1/R tot .
With the modified MC2 model, simulation runs are performed for Al droplets on a 1400 × 1400 grid using the energy values E S = 0.19 eV and E E = 1.71 eV as for the rate model and MC1. In the temperature range of 150°C ≤ T ≤ 300°C the simulated Al droplet densities agree within ±4% to the results of the MC1 simulation. This indicates that the approximations applied for the multiscale KMC algorithm MC2 are compatible with the conventional MC1 model. As a huge advantage, a corresponding simulation run takes about 20 s, which is about 30,000 times faster than the MC1 simulation. This extreme improvement of the computation time suggests the modifications applied for MC2 also for processes at higher temperatures.
A parameterization of the modified MC2 simulation for T = 300°C using a 1400 × 1400 field yields agreement with the experimental droplet density for E S = 0.19 eV and E E = 1.44 eV (Figure 6c). Similar to the rate model (Figure 3a), the influence of E S is only weak. Further simulations using the same E S and E E demonstrate also very good agreement with the experimental droplet densities at higher T (Figure 6c). Due to the more than exponential increase of the time needed for simulation runs as function of T (Figure 6d), the model is limited for Al droplets to T ≤ 460°C. For illustration, Figure 6a shows an AFM image of an AlGaAs surface after Al droplet deposition at T = 460°C. At this temperature, the deposited droplets are transformed into nanoholes during postgrowth annealing [35]. Figure 6b shows a corresponding surface with droplets simulated with MC2.
As an interesting point, the value of E E = 1.44 eV determined with MC2 is much smaller than the E E = 1.71 eV given by the rate model. In order to examine the origin of the deviating energy parameters, we study the relevant processes in more detail. The droplet density is balanced by nucleation (collisions between mobile monomers) and dimer break break-off (monomer escape from dimers). To separate both processes, we have performed simulations for irreversible aggregation, where escape processes are frozen by a very high E E = 10 eV and droplet nucleation is controlled only by E S . Here, for Al droplets, T = 460°C, and E S = 0.19 eV, the rate model yields N D = 2.8 × 10 10 cm −2 and the MC2 simulation a nearly identical N D = 2.6 × 10 10 cm −2 . The KMC model simulates nucleation basing on a more realistic atomistic picture, where diffusing monomers located on individual positions can collide. Here, the very small deviation between KMC and rate model results justifies the approximations made by the rate model, where the nucleation rate is calculated via N 2 1 σ 1 D assuming an average monomer density and a capture number which describes the monomer depletion around other objects. Now we switch to reversible aggregation with E E = 1.71 eV according to the rate model. Here the rate model yields N D = 2.4 × 10 8 cm −2 in agreement with the AFM data, whereas the MC2 simulation computes a much higher N D = 9.2 × 10 9 cm −2 . This example demonstrates that the difference between rate model and MC2 simulation is related to dimer break-off. The rate model calculates dimer break-off via 2πrζ N 2 D E and, after escape, the monomers increases the average monomer density on the whole surface. For the atomistic Monte Carlo simulation, after escape, a monomer is still in the close proximity of the other monomer with a high probability for a re-nucleation and, thus, a recovery of the dimer. This effect causes an effective reduction of the rate of dimer break-off events for the KMC simulation and finally results in much higher droplet densities compared to the rate model. We note, that we consider the atomistic approach of the KMC simulation to be more realistic.
Finally, we have simulated the density of Ga droplets using the MC2 model. Since for Al droplets the value of E S is model-independent, we use here E S = 0.115 eV as given by the rate model. Furthermore, the results of the rate model establish that for Ga the value of E E depends on the temperature (Figure 4). A parameterization of the MC2 model yields for T ≤ 300°C a constant value of E E = 1.24 eV and for T > 300°C a T-dependent E E = 1.24 + 0.06 (T[°C] − 300)/100 eV. These values of E S and E E allow a very good reproduction of the Ga droplet density by MC2 as is demonstrated in Figure 7. Like for the MC2 simulations of the Al droplet density, the values of E E are smaller than those of the mean-field model (E E = 1.51 eV for T ≤ 300°C and E E = 1.51 + 0.1 (T[°C] − 300)/100 eV for T > 300°C). We assume that the above re-nucleation effect also reduces the effective rate of escape processes from Ga droplets and requires a lower E E in comparison to the mean-field model. The accessible temperatures up to 550°C are higher compared to MC2 simulations of Al droplets which is caused by the lower densities of Ga droplets and the, thus, reduced simulation fields.

Conclusions
The deposition temperature controls the density of Ga and Al droplets during droplet epitaxy and etching over several orders of magnitude. This allows the self-assembled creation of quantum emitters with a tailored feature density. In particular ultra-low densities around 10 7 cm −2 are interesting for applications in quantum information technology, where individual quantum emitters can be easily addressed. The droplet density is balanced by nucleation events caused by collisions between diffusing adatoms and the escape of atoms from droplets. Both processes are characterized by the material dependent activation energies E S for surface diffusion and E E for escape.
The present manuscript studies the applicability of three different theoretical approaches for the modeling of the temperature-dependent density of Ga and Al droplets. An analysis basing on a scaling law requires the knowledge of the critical nucleus size i. Since for the droplets studied here, the value of i depends on temperature and varies with the deposition time, a simple scaling analysis is questionable. A model using mean-field rate equations is able to quantitatively reproduce the experimental Al droplet densities over the whole temperature range using E S = 0.19 eV and E E = 1.71 eV. This is an important result and demonstrates that the temperature dependent nucleation is balanced by the kinetics of droplet nucleation and the escape of atoms from droplets. Other processes such as Ostwald ripening as assumed earlier [12] are not relevant. Furthermore, the small value of E S indicates a very high activity of the mobile Al adatoms on the surface. In contrast to that, for Ga droplets the experimental droplet densities can be reproduced by E S = 0.115 eV, E E = 1.51 eV only at low temperatures up to 300°C, higher temperatures require a T-dependent E E . This result establishes that for T > 300°C additional processes become relevant that modify the binding energy of adatoms to the droplets.
In a further approach for modeling, the Al droplet density is studied using a kinetic Monte Carlo (KMC) simulation basing on an atomistic representation of the mobile adatoms on the surface. A conventional KMC model requires very long computation times for modeling high process temperatures with high surface activity and low droplet density. Using a multiscale KMC algorithm [32] to shorten the computation times, the KMC simulation quantitatively reproduces the Al droplet densities up to T = 460°C. Interestingly, the value of E S = 0.19 eV is equal to the result obtained using the rate model, but the value of E E = 1.44 eV is much smaller. This is attributed to the mean-field assumption of the rate model, where adatoms after escape from a droplet increase the average droplet density over the whole surface. In contrast to that, in the KMC model the escaped atoms are still in the close proximity of the droplet with a high probability of a recapturing. This reduces the effective rate of escape processes and, thus, requires a lower E E in comparison to the mean-field model.
We conclude that a KMC model is more accurate for modeling of droplet nucleation compared to a mean-field rate model. On the other side, the computation times of the rate model are substantially shorter which allows the modeling also of high process temperatures. Therefore, a modified rate model which considers the spatial distribution of adatoms after escape from a droplet is highly desirable.
Author Contributions: C.H. coordinated the study, designed the models, programmed in Delphi, performed simulation runs, and prepared the manuscript. S.F. assisted in the model design, programmed in C++, performed simulation runs, and assisted in the manuscript preparation. All authors have read and agreed to the published version of the manuscript.