Multiscale Kinetic Monte Carlo Simulation of Self

nanomaterials
Article
Multiscale Kinetic Monte Carlo Simulation of Self-Organized
Growth of GaN/AlN Quantum Dots
Jorge A. Budagosky 1 and Alberto García-Cristóbal 2, *
1 Nanotechnology on Surfaces and Plasma Group, Materials Science Institute of Seville (CSIC-US),
C/ Américo Vespucio 49, 41092 Seville, Spain
2 Instituto de Ciencia de los Materiales (ICMUV), Universidad de Valencia, C/ Catedràtic José Beltrán 2,
46980 Paterna, Spain
* Correspondence: [email protected]
Abstract: A three-dimensional kinetic Monte Carlo methodology is developed to study the strained
epitaxial growth of wurtzite GaN/AlN quantum dots. It describes the kinetics of effective GaN
adatoms on an hexagonal lattice. The elastic strain energy is evaluated by a purposely devised proce-
dure: first, we take advantage of the fact that the deformation in a lattice-mismatched heterostructure
is equivalent to that obtained by assuming that one of the regions of the system is subjected to
a properly chosen uniform stress (Eshelby inclusion concept), and then the strain is obtained by
applying the Green’s function method. The standard Monte Carlo method has been modified to
implement a multiscale algorithm that allows the isolated adatoms to perform long diffusion jumps.
With these state-of-the art modifications, it is possible to perform efficiently simulations over large
areas and long elapsed times. We have taylored the model to the conditions of molecular beam
epitaxy under N-rich conditions. The corresponding simulations reproduce the different stages of
the Stranski–Krastanov transition, showing quantitative agreement with the experimental findings

concerning the critical deposition, and island size and density. The influence of growth parameters,
Citation: Budagosky, J.A.; such as the relative fluxes of Ga and N and the substrate temperature, is also studied and found
García-Cristóbal, A. Multiscale to be consistent with the experimental observations. In addition, the growth of stacked layers of
Kinetic Monte Carlo Simulation of quantum dots is also simulated and the conditions for their vertical alignment and homogenization
Self-Organized Growth of GaN/AlN are illustrated. In summary, the developed methodology allows one to reproduce the main features
Quantum Dots. Nanomaterials 2022, of the self-organized quantum dot growth and to understand the microscopic mechanisms at play.
12, 3052. https://doi.org/10.3390/
nano12173052
Keywords: kinetic Monte Carlo; heteroepitaxy; Stranski-Krastanov growth mode; strain relaxation;
Academic Editors: Walter III-N semiconductors; gallium nitride; nucleation; self-organized quantum dots
Lacarbonara and Giovanni Formica
PACS: 64.75 Yz; 68.55 A-; 77.80 bn; 77.84 Bw; 81.16 Rf
Received: 10 August 2022
Accepted: 27 August 2022
Published: 2 September 2022
Publisher’s Note: MDPI stays neutral 1. Introduction

with regard to jurisdictional claims in
Semiconductor heteroepitaxy involves, in most cases, a lattice mismatch between the
published maps and institutional affil-
epitaxial layer and the substrate that results in the appearance of strain in the system. Such
iations.
strained heteroepitaxy is an area of materials research of enormous interest because of both
its exploitation in a variety of electronic and optoelectronic devices and its academic value
as a playground for the study of the interplay of strain-induced thermodynamical and
Copyright: © 2022 by the authors.
kinetic effects [1]. One key issue with such heteroepitaxy is that the strain introduced must
Licensee MDPI, Basel, Switzerland. always find a way to be relieved as the thickness of the epitaxial film increases. In the
This article is an open access article so-called Stranski–Krastanov (SK) growth mode, this strain relief occurs, after a critical
distributed under the terms and thickness is reached, by elastic relaxation of the stressed film via the formation of coherent
conditions of the Creative Commons three-dimensional (3D) nanometric islands [1]. It was recognized very early that, if the
Attribution (CC BY) license (https:// materials are appropriately selected, the overgrown structures can accommodate a few
creativecommons.org/licenses/by/ confined electronic states, and therefore can be used as quantum dots (QDs) [2]. Because of
4.0/). its origin, they are often referred to as self-organized or self-assembled quantum dots. One
Nanomaterials 2022, 12, 3052. https://doi.org/10.3390/nano12173052 https://www.mdpi.com/journal/nanomaterials

Nanomaterials 2022, 12, 3052 2 of 21
advantage of this approach is the ease of the fabrication of high density ensembles of QDs,
a necessary condition for their use in optoelectronic devices [2]. On the other hand, due to
the stochastic nature of the growth process, the resulting islands present a certain size and
position dispersion. This inhomogeneity makes difficult to control the desired properties
necessary for device implementation. The most studied examples of strained heteroepitaxy
are the InAs/GaAs and Ge/Si systems (with lattice mismatch of 6% and 4%, respectively),
but there are many more combinations of potential technological interest, such as InAs/InP,
CdSe/ZnSe, PbSe/PbTe, and GaN/AlN, to name a few.
The parameters and mechanisms that underlie the self-assembly of 3D islands have
long been the subject of great interest. From the experimental point of view, the goal
has been to master the specific conditions that allow to select the growth mode, and
provide control over the size and density of the islands. In this respect, experimental
data indicate that the growth of self-organized QDs is sensitive to parameters such as
substrate temperature, fluxes of the constituent materials, growth interruption time, etc.
Although these studies provide a valuable empirical knowledge, the microscopic origin
of the established general trends is not fully understood, owing to the complexity of
the physics of strained layer growth. A complementary approach is the development
of numerical simulations. At present, the most common and versatile growth modeling
strategy is based on the kinetic Monte Carlo (KMC) method. Its advantage lies in that,
due to its stochastic nature, it is able to cover the experimentally relevant growth times
and system sizes, while having access in principle to the full atomistic structure. Since
the growth of self-assembled QDs in the SK mode is driven by the elastic relaxation, the
incorporation of the strain in the KMC methodology is mandatory. It has to be noted
though that the computational cost in simulating epitaxial systems with strain is more
expensive than for systems without strain because the elastic displacement field has to
be dynamically updated following the changes in the growing film. The KMC method
has been implemented with a varying degree of sophistication for the simulation of self-
assembled epitaxial growth [3–10].
Most of the existing work on the growth of self-assembled QDs, both experimental
and simulation, is devoted to the study of InAs/GaAs and Si/Ge systems with cubic
crystalline structures. There has followed, however, a remarkable progress with the strained
heteroepitaxial growth of group-III nitrides with wurtzite structure, which have a wide
range of applications in the field of optoelectronics and high power electronics [11]. The
GaN/AlN heterostructure, normally grown along the polar axis of the crystal structure,
is a typical example; it is characterized by a lattice mismatch of the order of 2.4% and has
been demonstrated to present a SK growth transition leading to wurtzite GaN QDs. In the
literature, GaN/AlN quantum dots have been obtained by three major epitaxial growth
methods: plasma-assisted molecular beam epitaxy [12], ammonia-assisted molecular beam
epitaxy [13], and metalorganic vapor-phase epitaxy [14,15]. In this paper, we focus on the
GaN/AlN QDs obtained by plasma-assisted molecular beam epitaxy (PA-MBE) at CEA
Grenoble, since their growth has been thoroughly analyzed by a range of experimental
techniques [12,16–25]. Concerning the KMC simulations, there exist just a few studies
devoted to the growth of III-N wurtzite-type materials. They are, however, mainly focused
on GaN homoepitaxy [26–31]. To the best of our knowledge, nobody has reported KMC
simulations of self-organized growth of III-N QDs so far.
In this paper, we propose a 3D lattice kinetic Monte Carlo model to study the growth
of self-organized GaN/AlN QDs. Our aim is to perform simulations on length scales
close to one hundred nanometers and time scales of tens of seconds, by including the
evolving strain distribution in the system. This goal poses a strong computational challenge
that we have tackled by: (i) simulating effective GaN adatoms, instead of following in
detail the separate behavior of real Ga and N atoms, (ii) implementing a fast multiscale
algorithm that allows the adatoms to perform long diffusion jumps across multiple sites [32],
and (iii) developing an efficient procedure to obtain the elastic strain energy distribution
based on continuum elasticity theory. In order to validate the proposed model, we have
tailored it for its comparison against the experimental resuls on GaN/AlN QDs grown
by PA-MBE under N-rich conditions. The corresponding simulations of the evolution
of the surface morphology reproduce the different stages of the SK transition from an
initial two-dimensional (2D) layer to the full development of 3D QDs, and demonstrate
quantitative agreement with the experimental findings concerning the critical deposition,
and island size and density. The influence of relevant growth parameters, such as growth
temperature and the relative fluxes of Ga and N, is also studied and found to be consistent
with the experimental observations. Additionally, we have briefly simulated the growth of
multiple stacked QD layers, and demonstrated how the strain field can lead to both vertical
alignment and increase the uniformity of the QD ensembles of the upper layers. These
findings are in qualitative agreement with the reported experimental results [22].
The paper is organized as follows. In Section 2, we describe the theoretical background
and the procedure used in our simulations. In Section 3, we discuss the results obtained,
including a comparison with available experiments. Finally, Section 4 contains the overall
discussion and the conclusions extracted from the work.
2. Description of the Modeling Strategy

This section describes the specifics of the lattice kinetic Monte Carlo (KMC) method
that we have implemented for the simulation of epitaxial growth of wurtzite GaN self-
assembled quantum dot layers on (0001)AlN surfaces. First, we set up the model geometry
and comment on the selected events and relevant parameters, as well as the technical
solutions we have adopted to accelerate the computations. Afterwards, we present the
procedure followed to calculate and incorporate the strain energy of the lattice-mismatched
system into the heteroepitaxial growth simulations.
2.1. KMC Model and Implementation

The system under study consists of a GaN thin film growing on a substrate. The
substrate plays a passive role here, not being affected by the phenomena taking place at
the surface. Therefore, the simulation domain reduces in practice to the deposited GaN
film. In principle, two separate atomic species, Ga and N, should be simulated explicitly
within the KMC model. Such a detailed description of all possible processes would be
challenging, but one aspect of GaN growth provides a simplification. Under sufficiently
N-rich conditions to insure microscopic stoichiometry, the Ga adatom kinetics, rather than
the Ga-N reaction kinetics, is the rate-limiting step. The GaN growth can thus be effectively
regarded as a single species process involving the attachment and migration of Ga adatoms.
Although the N kinetics is not considered explicitly in this approach, it will be taken care of
in an averaged way by using effective parameters for the Ga adatom diffusion rate that
depend on the Ga/N flux ratio. For this reason, we use the term “effective GaN adatoms”
to refer to the simulated species. The GaN film is thus modeled as a set of effective GaN
adatoms located at the sites of a hexagonal lattice of parameters a∗ and c∗ ; see Figure 1.
The in-plane projection of the domain, shown in Figure 1b, has l sites on each side, so
that the total number of two-dimensional (2D) lattice sites is Nl = l × l, √
and the simulated
surface is S = Nl ∆x1 ∆x2 , where ∆x1 = a∗ and ∆x2 = a∗ cos(π/6) = 3a∗ /2. Solid-on-
solid restrictions are imposed (i.e., lateral overhangs and vacancies are not allowed) and
periodic boundary conditions are applied in both lateral directions. The GaN adatoms of
our simulation distribute dynamically over this lattice as explained below. It should be
noted though that our model is not strictly atomistic. The consideration of effective GaN
adatoms hopping on a lattice and the use of the effective parameters introduced below have
to be understood as a mean-field approach [33,34]. In all simulations we used c∗ = cAlN /2,
where cAlN = 0.4982 nm is the axial lattice parameter of the substrate wurtzite structure (c∗
being therefore the thickness of one monolayer, ML). The effective in-plane parameter is
taken to be a∗ = 0.502 nm. Our model does not contemplate the occurrence of substrate
surface reconstruction, vacancies or intermixing, which nevertheless are thought to play a
minor role in the growth of GaN self-assembled quantum dots.
Ga adatom
(a)
(i)
deposition (iii)
desorption
GaN film hI
(ii)
diffusion
c*
I a*
Substrate
(b)
Figure 1. (a) Cross section of the system under study. The actual domain of simulation is the indicated
GaN film, modeled as a set of adatoms located at predefined lattice sites, and represented here as
spheres in hexagonal unit cells of parameters a∗ and c∗ . The in-plane sites are labeled by integers
I, and the corresponding local heights are denoted by h I (see text). Also sketched are the possible
processes that a given adatom can experience: (i) deposition, (ii) diffusion, and (iii) desorption.
(b) Top view of the GaN film, with indication of the 2D hexagonal lattice, and illustration of the
possible diffusion jumps that an adatom at x I can make.
The instantaneous configuration of the film can be specified by an occupation array

N
{σI } I =l 1 , where I is a suitably defined integer index running over all 2D lattice sites, and
σI represents the number of adatoms piled up at site I. The local thickness of the GaN
layer at each site I is then h I = σI c∗ , see Figure 1a, and is used below to represent the
surface morphology. For later reference, we denote by x I the in-plane position vector of 2D
lattice site I. Within the KMC method, the configuration of the GaN film evolves in time by
stochastic changes of the occupation of the discrete lattice sites: In its standard version, at
each time step, upon the random selection of a site, a properly selected process (among a
catalog of possible events, see below) is performed there, and the configuration array {σI }
is updated accordingly. The corresponding time interval is calculated and added to the
elapsed simulation time [35,36].
For the simulation of the epitaxial growth we have included three types of pro-
cesses: deposition (adsorption), surface diffusion, and evaporation (desorption), which
are schematically indicated in Figure 1a. Each elementary event is associated with a rate,
as follows:
(i) During the growth simulation, the adatoms are randomly deposited onto the
surface, with a deposition rate rdep = FNl , where F represent the deposition flux (in ML/s).
(ii) The adatoms sitting on the surface are allowed to randomly hop to one of its six
nearest neighbor sites, see Figure 1b, with a jump rate defined by an Arrhenius-type ex-
pression:
E
rdif = ν0 exp − dif , (1)
kB T
where Edif is the diffusion energy barrier, k B the Boltzmann constant, and T the substrate
temperature. The factor ν0 , called attempt frequency, defines the time scale of the surface
processes, and is typically related with the vibrational frequency of a surface adatom. Here,
we take ν0 = 1013 s−1 .
(iii) The adatom desorption process is analogously described by the rate:

E
rdes = ν0 exp − des , (2)
kB T
where Edes is the corresponding energy barrier.

The energy barriers of both diffusion and desorption processes are actually composed
of various contributions:
(0)
Edif = Edif + nEb + mEstep − Estr , (3a)
(0)
Edes = Edes + nEb − Estr . (3b)
(0) (0)
The terms Edif (Edes ) are the contribution to the diffusion (desorption) barriers due
to the binding to the underlying surface. For both diffusion and desorption barriers we
assumed a simple bond-counting rule, where the energy Eb is a measure of the binding
energy between lateral neighbor adatoms. Accordingly, in case a given adatom is bonded to
n neighbors, an increase +nEb of the diffusion and desorption barriers is introduced. In the
case of the hexagonal lattice, n runs from 0 (isolated adatom) to 6 (completely surrounded
adatom). In some instances, the selected diffusion jump would take place across a single
step, i.e., between sites with a difference in height of c∗ , as would occur at an island edge.
Single jumps across multiple steps are explicitly forbidden in our model. In the single
step jumps one has to take into account the (Ehrlich-)Schwöbel effect [30,37–39]. This is
effectively accounted for here by adding a term +mEstep , where Estep is the magnitude
of the Schwöbel barrier that the adatom has to surmount to cross the step edge, and m
is the number of in-plane jump directions in which there is a step. For simplicity, in our
simulations the additional step barrier is only included for downward step jumps. In the
growth of a material lattice-mismatched to the substrate (as in the case of GaN on AlN), the
growing layer is (inhomogeneously) strained and there is a corresponding accumulation of
elastic energy. Consequently, the surface adatoms will diffuse (desorb) more often from
the strained sites than from the relaxed ones [40]. As time passes by, this will lead to a
net flux of adatoms from the strained to the relaxed areas, which is at the origin of the
self-assembly of islands in these heteroepitaxial systems. To incorporate this effect into
our simulation procedure, we introduce a strain-induced barrier term − Estr , reflecting
the increased likelihood of diffusion and evaporation events at sites where Estr is larger.
The calculation of Estr will be discussed in Section 2.2. We note that both Edif and Edes
actually depend on the local environment (through n and m) and absolute position (through
space-dependent Estr ) of the surface adatom considered.
During the KMC simulations, the system just described is evolved in time by using
an event-based Bortz–Kalos–Lebowitz (BKL) algorithm [32,41]. The corresponding time
interval ∆t is calculated and added to the elapsed simulation time. In this work, compelled
by the general requirement to increase the efficiency of growth simulations in large surfaces
and for sufficiently long times, and the specific need to include in a simple way the effect of
strain on the surface kinetics, we have incorporated the following modifications over the
standard BKL algorithm:
• At every time step, the systematic search across the entire surface for the sites attainable
by the process selected by the BKL algorithm has a cost in CPU time scaling with
the surface size. To avoid this problem, we have used an inverted-list algorithm [42].
Since the updating of the inverted lists is a local procedure, the adatom search process
becomes now independent of the surface size. This approach is specially efficient in the
case of having a not too extensive catalog of possible events. In our case, in the absence
of Schwöbel and strain-induced barriers, the system would have a very simple rate
structure: There are only seven transitions for the diffusion and desorption processes,
corresponding to the possible bonds to nearest neighbors (n = 0, . . . , 6). To include
the (Ehrlich-)Schwöbel effect (in cases where the jump occurs across a step) and the
inhomogeneous strain-induced barrier Estr , without having to modify the structure of
inverted lists mentioned above, we use an acceptance-rejection algorithm [42].
• In growth regimes where surface diffusion plays a dominant role, most of the com-
putation time is spent in calculating the adatom diffusion random walks. However,
the associated jump rates for isolated and bonded adatoms can differ by orders of
magnitude, due to the relative factor exp − knETb . We use the multiscale kinetic Monte
B
Carlo method proposed in Ref. [32] that takes advantage of this disparity in the
adatom dynamics by allowing the isolated adatoms to make longer jumps than the
bonded ones. This reduces drastically the number of Monte Carlo steps required to
complete the simulation. For example, for a surface with Nl = 500 × 500 lattice sites,
a tenfold reduction in computation time can be easily achieved. The long jumps in
the multiscale model are an effective way of simulating by a single event a chain of
multiple short jumps. We mentioned earlier that the presence of a nonuniform surface
elastic energy is expected to lead, for sufficiently long times, to a preferential diffusion
towards more relaxed regions. Therefore, it is natural to introduce, for long jumps
(with associated long elapsed times), a bias in the jump direction that gives preference
to the arrival sites with the lowest values of Estr . In this work, we have implemented
that idea as follows: Once an isolated adatom (n = 0) has been selected, we scan its
neighborhod in search of obstacles (steps to upper terraces, other isolated adatoms or
clusters). This search sets the number of sites L (>1) for the long jump, and the jump
direction α (α = 1, . . . , 6) is chosen according to the following probability distribution:
6
pα = wα / ∑ w β , (4)
β =1
with
(α)
!
ηα Estep + k( L − 1)∆Estr
wα = exp − , (5)
kB T
where the term ηα Estep accounts for the possible presence of Schwöbel barriers: ηα = 1
(ηα = 0) if there is (there is not) a step when making a jump of L sites along direction
(α)
α. The strain-induced bias term k( L − 1)∆Estr is proportional to the elastic energy
(α)
difference ∆Estr between the target site, located L sites away along direction α, and
the initial site. The factor ( L − 1) introduces an explicit dependence of the bias on
the jump length L. The coefficient k is a free input parameter of the model. A similar
parameter was used in Ref. [9] to simulate superlattices of PbSe/PbEuTe QDs. In our
case, the optimum value k = 12 was determined so that the simulations reproduce the
general features of growth experiments in GaN/AlN systems.
In principle, one could try to determine the parameters of the KMC model on the basis
of the energetics obtained from ab initio total energy calculations [29,43]. However, this is
outside the scope of this work, and we will instead estimate the strain-independent model
(0) (0)
parameters (Edif , Edes , Eb , Estep ) by trying to reproduce various well established features of
the homoepitaxial growth of GaN (see Section 3.1).
Since this work is ultimately interested in the growth of QDs in the Stranski–Krastanov
mode, all the simulations begin with a first monolayer of GaN already formed. The output
of the simulation at later times is a certain array {h I }, which defines the inhomogeneous
thickness of the GaN film. The simulations can be analyzed qualitatively by direct obser-
vation of the resulting surface morphology, and quantitatively by computing, at selected
stages of the simulation, quantities such as the dot density or the dot height and diam-
eter distribution. A quantity of special interest that we extract from our simulations is
the surface step density Sd (associated to the degree of surface roughness), which has
been demonstrated to provide a direct representation of the specular electron diffraction
(RHEED) intensities recorded in a variety of growth experiments [33]. In the analysis
of RHEED experiments that we present below, we have represented the quantity 1 − Sd
instead of Sd , since the maxima (minima) in the RHEED traces occur at situations when the
density of steps at the surface of the GaN layer is minimal (maximal).
2.2. Calculation of the Strain-Induced Energy Barrier Estr

As advanced above, to simulate the growth process in GaN/AlN systems, due to
the mismatch between the lattice parameters of GaN and AlN, it is necessary to take
into account the strain field accumulated in the GaN epitaxial layer. Since morphological
changes are restricted to the surface of the material, we focus here on the calculation of
the elastic strain energy density at the free surface of GaN, which will enter into the KMC
calculations through the strain-induced barrier Estr .
At this point, there are alternative possibilities to calculate Estr . One approach is
based on the so-called ball-and-spring models, where an explicit physical displacement
of the adatoms from their initial lattice site is allowed, and a system of restoring forces is
implemented (off-lattice models) [10,44,45]. However, given the effective character of the
adatoms in our model, there is not much sense in pursuing an atomistic description of the
elastic energy. We have opted instead for another approach, in which the strain field in the
growing layer+substrate system is treated in the framework of the continuum theory of
elastic media, and afterwards used to determine Estr . To couple the discrete nature of our
original model to a continuous description of the system, it is necessary to define first a
smoothed version of the surface that can be used in the calculation of the strain energy. To
that end, we have used a coarse-graining strategy that uses an appropriate filter to obtain
a function h( x⊥ ) describing the continuous surface morphology from the values {h I }
defining the GaN layer thickness in the original discrete description (see Figure 2a) [46].
Once h( x⊥ ) is determined, the strain field at the surface and the corresponding strain
energy density U ( x⊥ ) are calculated in this paper as follows. The mismatched GaN layer
(together with the buried QD layers, if present) is considered as an inclusion in the matrix
given by the AlN substrate (see Figure 2b). According to Eshelby’s theory, the inclusion
regions act as the source of the strain, which extends throughout the system [47]. For inclu-
sions buried under a flat surface, the exact strain and stress distribution can in principle
be obtained by using the Green’s function method [48]. However, in our case we face the
need of implementing the stress-free boundary condition on the non-flat surface given
by h( x⊥ ). We have circumvented this problem by using an small-slope approximation
of the GaN surface topography, in which it is exploited that |∇⊥ h( x⊥ )| 1 to apply a
perturbation procedure and simplify the problem [49]. Experimental observations gen-
erally show that the diameter of the GaN/AlN QDs exceeds by one order of magnitude
their height, which justifies the use of this approach. The simplification conveyed by the
perturbative treatment allows one to approximate the system with inhomogeneous thick-
ness and stress-free boundary conditions by a system with a flat surface and a prescribed
distribution of surface tractions, determined by h( x⊥ ) and the lattice mismatch [49]. It is

then possible to use the elastic Green’s function for a transversely isotropic half-space [50].
By following the sketched procedure, we can obtain the strain, stress, and strain energy
density U ( x⊥ ) at the surface of an arbitrary system of small-slope (GaN) islands on a
(0001)-oriented semi-infinite (AlN) matrix. The details of our strain calculations will be
discussed more extensively in a forthcoming publication [51]. The required parameters
are the lattice mismatch of the GaN/AlN system, ( aAlN − aGaN )/aGaN = −0.024 and
(cAlN − cGaN )/cGaN = −0.039, and the AlN elastic constants.
Figure 2. Schematic illustration of an epitaxial GaN layer on a semi-infinite AlN substrate. Both the
original discrete description {h I } (a) and the smoothed version h( x⊥ ) (b) of the surface morphology
are shown. The topography in part (b) is overlaid by a color map representing the strain energy
Estr ( x⊥ ), as calculated within continuum elasticity theory (see text).
Finally, we define the energy Estr ( x⊥ ) by:
Estr ( x⊥ ) = Ω∗ U ( x⊥ ) f ( x⊥ ) , (6)
√
where Ω∗ = 3/2 a∗ 2 c∗ is the volume of the unit cell occupied by each GaN adatom, and
f ( x⊥ ) is a function that heuristically takes into account the formation of the GaN wetting
layer [52]:
1
f ( x⊥ ) = 1 − h i . (7)
h( x )
1 + exp 4 c∗⊥ − γ
The parameter γ = 1.5 has been chosen in order to obtain a critical thickness of 2 ML
for the SK transition, as occurs in the growth experiments that we compare below with our
simulations. In Figure 2b, we illustrate the calculation of Estr ( x⊥ ) for a typical situation.
The discrete values of the strain-induced barrier energy appearing in the KMC simulations
are then given by Estr ( I ) = Estr ( x⊥ = x I ).
Despite the efficiency of the procedure devised, the cost in CPU time to update
the surface strain energy density of the film at each time step would still be prohibitive.
However, major changes in the surface morphology occur at a time scale much greater than
that of the individual surface events. Therefore, in all the simulations we have updated the
strain-induced energy only at every nstr = 500(l/64) time steps. Tests using lower values
of nstr show no qualitative changes.
3. Simulation Results
In this section, we present the results of our numerical simulations. We first present a
strategy to calibrate some of the parameters of the model by a comparison of the simulations
with GaN homoepitaxial growth experiments. We will afterwards investigate the effect of
some relevant parameters (GaN coverage Θ, flux ratio φGa /φN and substrate temperature
T) on the morphology of a single GaN QD layer grown on AlN. Finally, we will study
the growth of stacks of GaN/AlN QD layers. All our simulations have been performed
on a lattice of Nl = l × l = 256 × 256 sites (corresponding to an area of S ∼ 128.4 nm
× 111.2 nm).
3.1. Calibration of the Activation Energies by Homoepitaxial Growth Simulations

It is well established experimentally that the relative amount of Ga and N on the
surface significantly influences the kinetics of MBE growth [20,21]. Although the details
behind that effect cannot be explicitly taken into account in our effective KMC model, we
explore here the possibility of indirectly accounting for it by proposing a dependence of
the activation energies on the ratio between the Ga and N fluxes. However, due to the
evaporation processes and decomposition of GaN at the high substrate temperatures at
which MBE growth takes place, the effective fluxes (φGa eff , φeff ) are not necessarily equal to
N
the nominal fluxes (φGa , φN ) (the output fluxes of the Knudsen cells of the MBE system). In
the N-rich regime that we are interested in, we can consider the N desorption as negligible,
so that φN eff = φ , and the growth of GaN is limited by the net flux of Ga [20,21,53]. In
N
Appendix A we propose a simple model to obtain φGa eff in terms of ( φ , φ ), as well as the
Ga N
temperature, which is valid only in the N-rich regime. In the simulations in this paper we
eff to the deposition flux F introduced in Section 2.1.
assimilate the effective flux φGa
To obtain an estimate of the activation energies of our KMC model, we first assume
them to be dependent on the effective flux ratio η ≡ φGa eff /φ through appropiate model
N
functions. The free parameters of these functions have been fixed by fitting the RHEED
traces measured during GaN homoepitaxial growth [20] against the corresponding tempo-
ral evolution of 1 − Sd calculated from our simulations. Figure 3 shows the RHEED traces
for three nominal ratios φGa /φN together with the time evolution of 1 − Sd calculated with
the fitted activation energies. Of course, in the simulations of homoepitaxial growth (in
our case GaN on GaN), the term − Estr is absent from the diffusion and desorption energy
barrier expressions (3). The approximate expressions for the activation energies that are
obtained after this calibration procedure are:
(0)
Edif ' 1.78 − 4.87 × 10−3 e4η eV ,
Eb ' 0.41 − 0.12 η eV ,
(0)
Edes ' 2.53 − 0.17 η eV ,
η 14.4
Estep ' 0.12 − 0.06 eV . (8)
0.7514.4 + η 14.4
As we noted above, due to the simplicity of the growth model, the activation energies
of the different processes should be regarded as effective parameters, i.e., parameters whose
values have to be understood as averages of the microscopic processes that we are not
considering explicitly.
3.2. Onset of the Stranski–Krastanov Growth Mode and Formation of Quantum Dots
In this section, we analyze the heteroepitaxial growth of a GaN film on an AlN
substrate as a function of the coverage Θ (amount of deposited GaN). The substrate
temperature is fixed to T = 730 ◦ C. The Ga and N fluxes are chosen so that φGa /φN = 0.8,
and the growth rate is ∼0.15 ML/s, i.e., the growth is performed under a slight excess of N.
RHEED intensity
1  Sd
GaN  0.89
RHEED intensity (arb.unts)
1  Sd (arb.unts)
GaN  0.82
GaN  0.68
0 5 10 15 20 25 30 35 40
Time (s)
Figure 3. Comparison between the RHEED traces measured during GaN homoepitaxial growth at
T = 730 ◦ C under different ratios φGa /φN in Ref. [20] with the step density 1 − Sd curves obtained
from simulations with the fitted activation energies. In all cases, φN = 0.28 ML/s.
Figure 4 shows the evolution of the GaN film morphology as the deposition progresses,
from a completely formed first monolayer up to a final GaN coverage of Θ = 5 ML. It can
be observed that the growth of the second monolayer takes place by the 2D layer-by-layer
mode, i.e., by the formation of 1 ML high nuclei that grow progressively into laterally
extended 2D islands (platelets) by capturing the adatoms that had been deposited on the
surface (see the results in Figure 4 for a coverage of 1.5 ML). These platelets coalesce until
the second monolayer is fully completed. However, from that point on the morphology
changes drastically: although initially 2D platelets start to grow as before, at a deposition
of 2.2–2.4 ML new clusters are formed on top of some of those platelets, and simultane-
ously, due to the strain effects, the lateral growth of the platelets becomes energetically
inconvenient and therefore strongly inhibited before coalescence becomes important. This
qualitative change, with the formation of stable 3D islands (2 or more ML high) on top of
the flat surface (clearly observed in Figure 4 for a coverage of 2.6 ML), marks the onset of the
Stranski–Krastanov transition [21]. In order to quantitatively assess these observations, we
show in Figure 5a the evolution of the calculated step density 1 − Sd as a function of Θ. The
oscillation between the first and second ML is indicative of the layer-by-layer growth mode.
Shortly after the completion of the second ML (identified by the disappearance of the oscil-
lation), the critical thickness is reached at Θ ∼ 2.25 ML, and the growth is changed to the
SK mode, with the formation of the first 3D islands (quantum dots), as commented above.
The SK transition does not manifest itself very clearly in this representation, although
a very small elbow associated to it can be hinted at 2.25–2.4 ML. For further deposition,
from 2.2 to 3 ML, the quantity 1 − Sd decreases monotonically, in that part of the process
the existing islands continue to grow, while at the same time new 3D islands begin to
form. This is best observed in Figure 5b, where we depict the evolution of the QD density
with increasing GaN coverage. The simulation results demonstrate that the QD density
increases rapidly from the onset of the SK transition until approximately 3.0 ML. From
that moment on and up to 5 ML, the number of QDs remains constant, and the growth
continues exclusively by increasing the volume of the existing QDs, as can be observed
in Figure 4. For the sake of comparison, we also present in Figure 5b the corresponding
experimental results of the QD density reported in Ref. [21]. We conclude that the results
from the numerical simulations concerning the different stages of the growth process and
the overall features of the quantum dot ensemble are generally in good agreement with the
growth experiments. There remain, however, some discrepancies that could be addressed
in further work. Most noticeably, the experiments demonstrate the appearance, with in-
creasing coverage after density saturation, of a bimodal height distribution, with fixed
height islands coexisting with larger islands ever growing in size. Our present model does
not reproduce such behavior.
Figure 4. Simulation of the process of growth of 5 ML of GaN on AlN. The various panels show
the GaN surface morphology at different instants (corresponding to the specified values of the GaN
coverage Θ). The color scale indicating the discrete height is independently chosen in each panel, but
the corresponding maximum height is specified in every case.
1.0
1Sd (arb. units)
0.9 2D 3D
0.8 (a) ~ 2.25 ML
1012
QDs density (cm-2)
1011
Experimental data
1010 (b) KMC simulation
109
1 2 3 4 5
GaN coverage  (ML)
Figure 5. (a) Evolution of the step density 1 − Sd as a function of Θ. The vertical line at ≈2.25 ML
indicates the critical thickness at which the 2D-3D SK transition occurs. (b) QD density as a function
of Θ. The results obtained from the KMC simulations are compared with the experimental data of
Ref. [21].
3.3. QD Growth as a Function of φGa /φN Ratio

One of the factors easily controllable in GaN/AlN heteroepitaxial growth experiments,
and known to influence the appearance of the SK transition and the characteristics of the
resulting QDs, is the ratio between the nominal fluxes of Ga and N, φGa /φN [18–20]. In this
section, we use our simulations to study the morphology of the QD layer as a function of
φGa /φN . The total amount of GaN deposited is Θ = 3 ML and the substrate temperature is
fixed to T = 730 ◦ C.
Figure 6 shows the film morphologies simulated for values of φGa between 0.15
and 0.28 ML/s. The N flux φN has been set at 0.28 ML/s. As can be observed, for φGa
eff /φ
less than 0.23 ML/s (φGa N ≈ 0.88), the SK transition competes with the kinetically
induced formation of extended 2D platelets: in this regime the mobility of the adatoms on
the surface, directly dependent on the ratio φGa eff /φ , is not high enough to promote the
N
nucleation of 3D precursor islands well defined and separated from each other. Therefore,
the system partially relaxes its elastic energy by means of a rough surface instead of
through a genuine 2D-3D transition. This behavior generally agrees with that observed in
growth experiments [18,19]. Furthermore, for simulations with increasing φGa /φN beyond
φGa ≥ 0.23, we observe a significant increase in the average island height, along with a
slight decrease in the diameter, together with a progressive size homogenization. The origin
(0)
of this trend lies in the decrease in both the diffusion barrier Edif and the binding energy
Eb (that can be associated to an increase in the effective mobility of the adatoms) as φGa
increases, so that the diffusion of the GaN adatoms across the surface is predominantly
determined by the inhomogeneities in the elastic energy. Therefore, the precursor 2D
islands slow down their lateral growth earlier than for growth under lower φGa /φN , and
the SK transition to 3D growth sets is clear. Our calculations (not shown) also demonstrate
that for Ga fluxes beyond 0.23 ML/s, the QD density decays linearly as we get closer to
the stoichiometry (which corresponds to φGa = φN + δφGaN ∼ 0.3 ML/s). Simulations
carried out considering a higher N flux of value φN = 0.34 ML/s (not shown here) reveal a
tendency of the GaN layer to relax its elastic energy via the formation of irregular islands
of small height and large diameter (instead of fully developed 3D dots), with essentially
constant density.
Figure 6. Surface morphology of a GaN epitaxial layer grown on an AlN substrate for simulations
with different nominal Ga fluxes φGa (specified in each panel of the figure). The N flux φN has been
set at 0.28 ML/s. The color scale indicating the discrete height is the same for all panels. A reference
height is provided in the panel for φGa = 0.28 ML/s.
3.4. QD Growth as a Function of Substrate Temperature

We have observed that variations in the mobility of the adatoms, as a function of
φGa /φN , can affect the growth (for a fixed temperature). On the other hand, since the
surface diffusion is a thermally activated process, significant changes in adatom mobility

can in principle be achieved by varying the temperature. It is relevant, therefore, to study
in detail the growth of QDs as a function of the substrate temperature. To that end, we have
performed simulations at various temperatures T between 610 ◦ C and 760 ◦ C. The total
amount of GaN deposited is fixed to 3 ML. The fluxes have been set to φN = 0.34 ML/s
and φGa = 0.25 ML/s, respectively. However, the nominal Ga flux φGa has been modified
when necessary (for T > 730 ◦ C) so as to compensate the increment of desorption with
temperature and thus maintain an approximately constant value of the growth rate around
∼0.25 ML/s. This allows one to ensure that the characteristics of the QD formation and
growth will depend exclusively on the adatom mobility activated by the temperature.
Figure 7 shows the different GaN layer morphologies obtained at various temperatures.
The first thing to notice is the absence of the SK transition for temperatures below 700 ◦ C,
attributed to the insufficient mobility of the adatoms. This finding agrees qualitatively with
the experiments of Ref. [21], where it is reported that the formation of QDs is suppressed for
temperatures below 680 ◦ C. The discrepancy in the temperature values can be tentatively
attributed to differences in other growth conditions between the simulation and experiment,
such as the N flux used, which is not specified in the report of the experiments. The
simulations also demonstrate that, as we increase the temperature from 700 ◦ C, there
is a reduction in the diameter of the QDs. This is accompanied by a clear decrease in
the QD density. This trend is best illustrated in Figure 8, which shows the values of
the QD density extracted from the simulations, together with experimental results from
Ref. [21]. As we have observed above, the island density tends to saturate at sufficiently high
GaN coverage. The experiments demonstrate that this saturation density depends on the
substrate temperature T, but when T > 730 ◦ C saturation occurs for GaN coverages larger
than 3.0 ML. Therefore, two different sets of experimental density values are displayed
in the figure: one corresponding to a fixed GaN coverage of 3 ML (as in the simulations)
and the other corresponding to the saturation density. The simulation results exhibit a
decrease in the island density with increasing temperature that lies between the two sets of
experimental results. This can be considered a good agreement on account that the growth
conditions established in the simulator may not be exactly the same as those used in the
experiment, as commented above.
Figure 7. Surface morphology of a GaN epitaxial layer grown on an AlN substrate as obtained by
simulations performed at different temperatures (specified in each panel of the figure). The color
scale indicating the discrete height is the same for all panels. A reference height is provided in the
panel for T = 760 ◦ C.
5
4
QD density (1011 cm-2)

3
1 density after 3 ML
saturation density
KMC simulation
680 700 720 740 760 780

Temperature (ºC)
Figure 8. QD density obtained from the KMC simulations as a function of temperature (blue circles).
For comparison, we also display two sets of experimental densities (red circles and green squares)
measured in Ref. [21] (see text).
3.5. Growth of Stacks of GaN/AlN QDs: Correlation Effects

The previous results have confirmed the utility of the model developed to investigate
the self-assembling process of a single layer of GaN QDs under different growth conditions.
In this section, we will show simulations of the growth of stacks of GaN/AlN QD layers.
For the modeling of the stacks, we proceed as follows: once the simulation of the first
GaN layer is finished, it is assumed that it is covered with an AlN layer (spacer) with a
chosen thickness d. This way, we obtain a new flat AlN surface on which the growth of the
next GaN layer can proceed. However, due to the presence of QDs embedded in the AlN
substrate, the surface will exhibit a nonuniform strain. This strain, calculated by means of
Eshelby’s inclusion theory [50,51], allows one to obtain the local AlN lattice parameters on
the free surface, which are used as input parameter in the calculation of the elastic energy
during the growth of the next QD layer. The process is repeated up to the completion
of the total number of QD layers. When the simulation of the last layer is finished, the
system is covered with an AlN capping layer having the same thickness as the spacers.
The above procedure allows us to compute the strain energy density U ( x) throughout the
entire volume of the heterostructure. This detailed description of the strain, which has not
been incorporated before in KMC growth simulations of GaN QD multilayers, is, however,
very useful to analyze the elastic interaction between the QDs through the barrier and its
dependence on the spacer thickness.
In the simulations, a total amount of 5 ML is deposited in each GaN layer. The growth
conditions were fixed at T = 730 ◦ C and φGa = 0.25 ML/s and φN = 0.31 ML/s. Although
real stacks can contain tens and even hundreds of layers [23,24,54], in this work, due
to computer memory limitations, we consider stacks of five QD layers each. We present
simulations of four stacks with different spacer thicknesses, labeled A (d = 41c∗ = 10.2 nm),
B (d = 21c∗ = 5.2 nm), C (d = 17c∗ = 4.2 nm) and D (d = 13c∗ = 3.2 nm). Figure 9 depicts,
next to each other, random cross sections of the four stacks. In the stack A, the different
layers seem to be totally independent. There is a great dispersion in sizes and a total
absence of vertical correlation. As the spacer thickness is decreased, it is already grasped in
stack B a certain tendency for the dots in successive layers to grow above each other, but an
important dispersion of sizes is still present. However, for stack C (spacing d = 4.2 nm) a
clear trend for the QDs to nucleate just above the QDs of the preceding layers is detected.
This behavior is more pronounced for the stack D (spacing d = 3.2 nm) where the spatial
correlation in the vertical direction and the effect of homogenization of sizes and density
is evident. This behavior should be even more pronounced with further increasing the
number of layers, as is observed in real stacks [23]. The figure also displays color maps of
the elastic energy density for the whole of the structure. It is apparent that the accumulation
of (tensile) strain energy in the regions of the AlN spacer above the buried dots can act as a
driving factor for the aligned growth of the QDs in the following layers.
Figure 9. Vertical cross sections through the simulated stacks consisting of five QD layers separated by
different spacer thicknesses: d = 10.2 nm (stack A), d = 5.2 nm (B), d = 4.2 nm (C), and d = 3.2 nm
(D). For each stack, the distribution of the elastic energy density U ( x) is shown by means of a
pixelated color map. The color scale has been truncated at 0.1 eV/nm3 , so that most the GaN volume,
which is strained beyond that value, appears as black colored.
4. Discussion
We have developed a (3D) lattice kinetic Monte Carlo model for the simulation of
the heteroepitaxial growth of GaN/AlN systems. The species simulated are effective GaN
adatoms. The model is standard in that it includes the relevant processes of deposition,
diffusion and desorption, and implements the BKL algorithm. In addition, it incorporates
several special characteristics that we highlight in the following. On account of the usual
wurtzite crystalline structure of the III-N semiconductors, the simulation domain is built
on an hexagonal lattice. The influence of the elastic strain inherent to lattice-mismatched
growth is incorporated through an inhomogeneous correction Estr ( I ) to the diffusion and
desorption energy barriers. This quantity defined at the lattice sites is related to the strain
energy density at the surface of the growing GaN film U ( x⊥ ). The latter quantity is in
turn calculated within linear elastic continuum theory by an efficient procedure purposely
devised and based on the combined use of the Eshelby inclusion concept and the Green’s
function method. This procedure accounts for the elastic anisotropy of the system and
allows to dynamically account for the deformation at the free surface of the irregular profile
of the growing film. The above strain energy considerations supersede previous works on
similar topics, which make drastic simplifications in various respects [8,9,55]. On the other
hand, our goal to study the self-organization effect imposes the need to perform simulations
over a large area domain on long time scales. In order to improve the computational speed
so as to fulfill the above requirements, we have introduced various modifications to the
standard BKL algorithm. The most important one has been to allow the isolated adatoms to
perform longer discrete jumps than the laterally bound adatoms (multiscale kinetic Monte
Carlo algorithm). The direction of these jumps has been selected according to the strain
landscape in which they take place. With these improvements, efficient realistic simulations
on surfaces with side lengths close to one hundred nanometers, over time intervals of tens
of seconds (multiplied by up to 5 in the case of stacks) are possible. We note that our
model assumes a flat substrate surface. In practice, most often the substrates present some
degree of deviation from flatness. For example, in the experiments of Ref. [21] to which we
compare our simulations, the AlN surface is characterized by about 30 nm wide terraces
and spiral hillocks. This structure however does not seem to affect the nucleation of the
dots which are distributed rather homogeneously. In some other works, vicinal substrates
with an intentional miscut are used, resulting in a step organization, made of regularly
ordered terraces, at the AlN surface [56]. When the terrace width is larger than the diffusion
length of the adatoms, their kinetics is essentially equal to that in a flat surface. This leads
to a homogeneous nucleation, with dots covering the entire surface. On the contrary, when
the diffusion length is comparable to the width of the terraces, the GaN nucleation at steps
is energetically favorable, leading to a clear alignment of the dots along the AlN steps. The
modeling of this specific behavior is interesting and could motivate future extensions of
our KMC model.
Let us now comment on what we consider the main limitations of the described
KMC model. In the first place, we note that it is not strictly atomistic since the simulated
units are effective GaN adatoms. This means that the parameters are also effective and, in
principle, their values must be obtained by trial and error. However, in our case we have
developed a systematic strategy to ascertain the optimal values adapted to our purposes.
Initially, simulations of homoepitaxial GaN PA-MBE growth have been quantitatively
compared with corresponding available experiments. This comparison has allowed one
to estimate the intrinsic parameters governing the growth dynamics on the GaN surface
in the N-rich regime, and its dependence on the nominal fluxes. Afterwards, the strain-
related parameters, notably the parameter k commented below, have been fixed to achieve
agreement with some key experimental results, such as the SK critical thickness, and
the correlation spacing in QD stacks. The other noteworthy point of the model is the
introduction of the free parameter k in Equation (5). As mentioned above, in the case of
isolated adatoms we have replaced the short jumps with longer jumps, which allow us
to considerably speed up the simulations. However, in doing this we are eliminating the
cumulative effect of the elastic inhomogeneities after many short jumps. To compensate
for it, we have added a bias in the direction of the long jumps. This bias is proportional to
the difference between the elastic energies of the initial and arrival sites and to the number
of lattice sites between them. The free parameter k does not alter the surface distribution
of Estr ( x⊥ ) calculated by the elastic continuum theory but just magnifies the amplitude
of the strain-induced modulation between the initial and final sites to compensate for the
loss of the cumulative effect of many short random jumps. A similar parameter has been
introduced in the KMC simulations of Ref. [9]. It is remarkable that a single parameter k is
sufficient for obtaining a reasonably good correspondence with the experimental data, as
commented below.
The capabilities of the developed KMC model have been demonstrated by comparing
its simulations to the comprehensive growth experiments of Ref. [21]. In particular, the
influence of macroscopic growth parameters such as the surface coverage, the deposition
fluxes and the substrate temperature has been investigated. We summarize here the main
findings obtained in the simulations that are in qualitative (and even quantitative in many
instances) agreement with the experiments: (i) The KMC simulations demonstrate clearly
that the GaN growth passes four stages: (1) initially, the growth is layer by layer; (2) sub-
sequently, 2D precursor islands form, which (3) later transform into genuine 3D islands
through the SK transition; (4) during further GaN growth, the 3D islands increase their size,
whereas their density remains constant. (ii) Within the model the adatom mobility, which is
coded through the diffusion rate, depends directly on the temperature and indirectly on the
impinging material fluxes. The simulations have demonstrated that these parameters (and
therefore the adatom mobility) strongly influence the SK transition. This is in agreement
with the experiments where it is observed that the surface morphology of GaN layers grown
by PA-MBE depends crucially on the Ga/N flux ratio and the substrate temperature. Finally,
we have presented simulations of the growth of stacked GaN QD layers separated by AlN
spacers, fully taking into account the elastic strain field in the multilayer system. The study
as a function of the spacer layer thickness reveals the important role of the in-plane lattice
parameter modulation in each spacer layer surface (induced by the buried QDs) in directing
the vertical alignment of the dots, together with their size and density homogenization.
Again, the simulations yield results in good qualitative agreement with experimental data
concerning the observed vertical correlations [23]. Still, the model fails to reproduce some
features clearly observed in the experiments. The most notable is the bimodal distribution
of sizes observed in these systems under certain growth conditions [21]. On the other
side, although the simulations results for the QD height and diameter are compatible with
those reported in the experiments, they do not clearly reproduce the typical hexagonal base
pyramid shapes found in real systems of self-assembled GaN QDs [17]. This is probably
due to the lack of an explicit atomistic description of the GaN wurtzite structure.
In summary, the KMC model developed here for the simulation of strained het-
eroepitaxial growth can be considered as an excellent compromise between continuum
models [57] and fully atomistic simulations [58,59]. We note that, even though the model
has been tailored to the study of PA-MBE growth of the GaN/AlN system under N-rich
conditions, the general set up can be modified or extended to have wider applicability. In
this respect, it could be easily adapted, e.g., to deal with the growth of GaN quantum dots
on non-polar surfaces [25].
Author Contributions: J.A.B. and A.G.-C. contributed to the modeling, analysis of the results, and
writing the manuscript. All authors have read and agreed to the published version of the manuscript.
Funding: The contribution of J. A. Budagosky to this work was made under the auspices of the EU
H2020 program through the ERC Starting Grant 3DScavengers (grant agreement 851929). The work
of A. García-Cristóbal was performed in the framework of the project I+D+I PID2020-112507GB-I0,
funded by MCIN/AEI/10.13039/501100011033 and by “ERDF A way of making Europe”, and under
Grant PROMETEO/2021/082 (ENIGMA), funded by Generalitat Valenciana.
Data Availability Statement: Not applicable.
Conflicts of Interest: The authors declare no conflict of interest.
Abbreviations
The following abbreviations are used in this manuscript:
KMC Kinetic Monte Carlo

SK Stranski–Krastanov
QD Quantum dot
PA-MBE Plasma-assisted molecular beam epitaxy
ML Monolayer
BKL Bortz–Kalos–Lebowitz
eff
Appendix A. Model for the Effective Flux φGa
eff that incorporates the influence of
We propose a simple model for the effective flux φGa
the thermal decomposition of GaN. If we would only take into account the Ga evaporation,
the net Ga flux should be less than the nominal one leaving the Knudsen cell. This effect
can be represented by a sticking factor α (α ∈ [0, 1]), which is the fraction of Ga atoms
impinging on the surface that will not evaporate and remain on the surface.
However, within the growth temperature range relevant here, the decomposition
of GaN previously formed on the surface represent an extra source of Ga that must be
incorporated into the global calculation of the effective flux. The rate of GaN decomposition
is the thermal decomposition rate of GaN in the vacuum, δφGaN , which can be described as
an Arrhenius type function:
EGaN
− kB T
δφGaN = νGaN e , (A1)
whose parameters, νGaN = 6.3 × 1013 s−1 and EGaN = 3.1 eV, are taken from Ref. [53].
In our model, we start from the above observations and define the growth rate of
GaN as:
νGaN = αφGa − δφGaN + αδφGaN . (A2)
The term δφGaN , given by (A1), represents the correction to the growth rate due to
the decomposition of GaN, whereas αφGa is the net flux of deposited Ga that remains on
the surface without evaporating, and αδφGaN is the extra Ga flux resulting from the GaN
decomposition which, given the N-rich environment, will remain on the surface. The
effective Ga flux on the surface is therefore:
eff
φGa = α(φGa + δφGaN ) . (A3)
The coefficient α itself depends in general on the nominal fluxes φGa and φN . We
assume this dependence to follow a distribution type function with the form:
1−A
α(φGa , φN ) = A + . (A4)
1 + e B(φGa −φN −δφGaN )
In order to determine the parameters A and B, we focus on the stoichiometry condition
at the surface which, for a given nominal N flux φN , is reached by increasing the nominal Ga
(s)
flux up to a value φGa equal to φN plus the correction corresponding to the Ga desorption
rate that, in these conditions, coincides with the GaN decomposition rate:
(s)
φGa = φN + δφGaN . (A5)
At the stoichiometry situation, the growth rate reaches the value:
(s)
νGaN = φN − δφGaN . (A6)
Substituting (A6) and (A5) into (A2), we obtain an expression for the coefficient α in
stoichiometry conditions:
φN
α(s) = . (A7)
φN + 2δφGaN

(s)
Now, by requiring α φGa , φN = α(s) , we determine A:
φN − 2δφGaN
A(φN ) = . (A8)
φN + 2δφGaN
Then, by imposing
the condition that the growth rate is maximum in stoichiometry
∂νGaN
conditions, i.e., ∂φGa (s) = 0, the parameter B is found to be:
φGa
φN
B(φN ) = . (A9)
δφGaN (φN + 2δφGaN )
Finally, with the coefficient α determined by the above procedure, the effective Ga
eff can be obtained by means of (A3) for any growth condition within the N-rich
flux φGa
regime. In order to test the validity of our model, we have compared in Figure A1 our
results for νGaN as a function of φGa with the growth rates measured experimentally in
Refs. [20,53]. An excellent agreement with the experiments is found even at temperatures
as high as 770 ◦ C. This shows that, at least up to this temperature, the desorption of N can
be considered negligible.
0.5
0.4 (a) TS = 730 ºC

GaN growth rate GaN (ML/s)
0.3
0.2
N flux N = 0.28 ML/s
0.1
0.0
0.5
0.4 (b) TS = 770 ºC
0.3
0.2
0.1
0.0
0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8
Nominal Ga flux Ga (ML/s)

Figure A1. Growth rate of GaN as a function of the impinging Ga flux, for two substrate temperatures
and two impinging N fluxes. The dots of (a,b) are experimental data taken from Refs. [20,53],
respectively. The lines represent fits calculated within the framework of the model described in
the text.
References
1. Shchukin, V.; Ledentsov, N.N.; Bimberg, D. Epitaxy of Nanostructures; Springer Science & Business Media: Berlin/Heidelberg,
Germany, 2004.
2. Bimberg, D.; Grundmann, M.; Ledentsov, N.N. Quantum Dot Heterostructures; John Wiley & Sons: Hoboken, NJ, USA, 1999.
3. Barabasi, A.L. Self-assembled island formation in heteroepitaxial growth. Appl. Phys. Lett. 1997, 70, 2565–2567. [CrossRef]
4. Schöll, E.; Bose, S. Kinetic Monte Carlo simulation of the nucleation stage of the self-organized growth of quantum dots. Solid-State
Electron. 1998, 42, 1587–1591. [CrossRef]
5. Rottler, J.; Maass, P. Second Layer Nucleation in Thin Film Growth. Phys. Rev. Lett. 1999, 83, 3490–3493. [CrossRef]
6. Meixner, M.; Schöll, E. Kinetically enhanced correlation and anticorrelation effects in self-organized quantum dot stacks. Phys.
Rev. B 2003, 67, 121202(R). [CrossRef]
7. Meixner, M.; Kunert, R.; Schöll, E. Control of strain-mediated growth kinetics of self-assembled semiconductor quantum dots.
Phys. Rev. B 2003, 67, 195301. [CrossRef]
8. Zhu, R.; Pan, E.; Chung, P. Fast multiscale kinetic Monte Carlo simulations of three-dimensional self-assembled quantum dot
islands. Phys. Rev. B 2007, 75, 205339. [CrossRef]
9. Mixa, M.; Holỳ, V.; Springholz, G.; Bauer, G. Kinetic Monte Carlo simulation of self-organized growth of PbSe/PbEuTe quantum
dot multilayers. Phys. Rev. B 2009, 80, 045325. [CrossRef]
10. Schulze, T.; Smereka, P. Kinetic Monte Carlo simulation of heteroepitaxial growth: Wetting layers, quantum dots, capping, and
nanorings. Phys. Rev. B 2012, 86, 235313. [CrossRef]
11. Nakamura, S.; Pearton, S.; Fasol, G. The Blue Laser Diode: The Complete Story; Springer Science & Business Media: Berlin/Heidelberg,
Germany, 2000.
12. Daudin, B.; Widmann, F.; Feuillet, G.; Samson, Y.; Arlery, M.; Rouvière, J.L. Stranski-Krastanov growth mode during the molecular
beam epitaxy of highly strained GaN. Phys. Rev. B 1997, 56, R7069–R7072. [CrossRef]
13. Damilano, B.; Grandjean, N.; Semond, F.; Massies, J.; Leroux, M. From visible to white light emission by GaN quantum dots on Si
(111) substrate. Appl. Phys. Lett. 1999, 75, 962–964. [CrossRef]
14. Miyamura, M.; Tachibana, K.; Arakawa, Y. High-density and size-controlled GaN self-assembled quantum dots grown by
metalorganic chemical vapor deposition. Appl. Phys. Lett. 2002, 80, 3937–3939. [CrossRef]
15. Bellmann, K.; Tabataba-Vakili, F.; Wernicke, T.; Strittmatter, A.; Callsen, G.; Hoffmann, A.; Kneissl, M. Desorption induced GaN
quantum dots on (0001) AlN by MOVPE. Phys. Status Solidi RRL 2015, 9, 526–529. [CrossRef]
16. Widmann, F.; Simon, J.; Daudin, B.; Feuillet, G.; Rouvière, J.; Pelekanos, N.; Fishman, G. Blue-light emission from GaN
self-assembled quantum dots due to giant piezoelectric effect. Phys. Rev. B 1998, 58, R15989–R15992. [CrossRef]
17. Arlery, M.; Rouvière, J.; Widmann, F.; Daudin, B.; Feuillet, G.; Mariette, H. Quantitative characterization of GaN quantum-dot
structures in AlN by high-resolution transmission electron microscopy. Appl. Phys. Lett. 1999, 74, 3287–3289. [CrossRef]
18. Bourret, A.; Adelmann, C.; Daudin, B.; Rouvière, J.L.; Feuillet, G.; Mula, G. Strain relaxation in (0001) AlN/GaN heterostructures.
Phys. Rev. B 2001, 63, 245307. [CrossRef]
19. Mula, G.; Adelmann, C.; Moehl, S.; Oullier, J.; Daudin, B. Surfactant effect of gallium during molecular-beam epitaxy of GaN on
AlN (0001). Phys. Rev. B 2001, 64, 195406. [CrossRef]
20. Adelmann, C.; Brault, J.; Jalabert, D.; Gentile, P.; Mariette, H.; Mula, G.; Daudin, B. Dynamically stable gallium surface coverages
during plasma-assisted molecular-beam epitaxy of (0001) GaN. J. Appl. Phys. 2002, 91, 9638–9645. [CrossRef]
21. Adelmann, C.; Daudin, B.; Oliver, R.; Briggs, G.; Rudd, R. Nucleation and growth of GaN/ AlN quantum dots. Phys. Rev. B 2004,
70, 125427. [CrossRef]
22. Gogneau, N.; Fossard, F.; Monroy, E.; Monnoye, S.; Mank, H.; Daudin, B. Effects of stacking on the structural and optical
properties of self-organized GaN/AlN quantum dots. Appl. Phys. Lett. 2004, 84, 4224–4226. [CrossRef]
23. Chamard, V.; Schülli, T.; Sztucki, M.; Metzger, T.; Sarigiannidou, E.; Rouvière, J.L.; Tolan, M.; Adelmann, C.; Daudin, B. Strain
distribution in nitride quantum dot multilayers. Phys. Rev. B 2004, 69, 125327. [CrossRef]
24. Cros, A.; Garro, N.; Llorens, J.; García-Cristóbal, A.; Cantarero, A.; Gogneau, N.; Monroy, E.; Daudin, B. Raman study and
theoretical calculations of strain in GaN quantum dot multilayers. Phys. Rev. B 2006, 73, 115313. [CrossRef]
25. Founta, S.; Coraux, J.; Jalabert, D.; Bougerol, C.; Rol, F.; Mariette, H.; Renevier, H.; Daudin, B.; Oliver, R.; Humphreys, C.; et al.
Anisotropic strain relaxation in a-plane GaN quantum dots. J. Appl. Phys. 2007, 101, 063541. [CrossRef]
26. Wang, K.; Singh, J.; Pavlidis, D. Theoretical study of GaN growth: A Monte Carlo approach. J. Appl. Phys. 1994, 76, 3502–3510.
[CrossRef]
27. Mao, H.; Lee, S.; Park, S.J. The Monte Carlo simulation of epitaxial growth of hexagonal GaN. Surf. Sci. 1999, 432, L617–L620.
[CrossRef]
28. Xu, D.; Zapol, P.; Stephenson, G.B.; Thompson, C. Kinetic Monte Carlo simulations of GaN homoepitaxy on c-and m-plane
surfaces. J. Chem. Phys. 2017, 146, 144702. [CrossRef]
29. Chugh, M.; Ranganathan, M. Lattice kinetic Monte Carlo simulation study of the early stages of epitaxial GaN (0001) growth.
Appl. Surf. Sci. 2017, 422, 1120–1128. [CrossRef]
30. Turski, H.; Krzyżewski, F.; Feduniewicz-Żmuda, A.; Wolny, P.; Siekacz, M.; Muziol, G.; Cheze, C.; Nowakowski-Szukudlarek, K.;
Xing, H.G.; Jena, D.; et al. Unusual step meandering due to Ehrlich-Schwoebel barrier in GaN epitaxy on the N-polar surface.
Appl. Surf. Sci. 2019, 484, 771–780. [CrossRef]
31. Erwin, S.C.; Lyons, J.L. Atomic Layer Epitaxy of III-Nitrides: A Microscopic Model of Homoepitaxial Growth. ACS Appl. Mater.
Interfaces 2020, 12, 49245–49251. [CrossRef]
32. DeVita, J.P.; Sander, L.M.; Smereka, P. Multiscale kinetic Monte Carlo algorithm for simulating epitaxial growth. Phys. Rev. B
2005, 72, 205421. [CrossRef]
33. Shitara, T.; Vvedensky, D.; Wilby, M.; Zhang, J.; Neave, J.; Joyce, B. Step-density variations and reflection high-energy electron-
diffraction intensity oscillations during epitaxial growth on vicinal GaAs (001). Phys. Rev. B 1992, 46, 6815–6824. [CrossRef]
34. Vvedensky, D.; Haider, N.; Shitara, T.; Smilauer, P. Evolution of surface morphology during epitaxial growth. Philos. Trans. R. Soc.
Lond. Ser. A Phys. Eng. Sci. 1993, 344, 493–505. [CrossRef]
35. Voter, A.F. Introduction to the kinetic Monte Carlo method. In Radiation Effects in Solids; NATO Science Series II: Mathematics,
Physics and Chemistry; Springer: Dordrecht, The Netherlands, 2007; Volume 235, pp. 1–23. [CrossRef]
36. Andersen, M.; Panosetti, C.; Reuter, K. A practical guide to surface kinetic Monte Carlo simulations. Front. Chem. 2019, 7, 202.
[CrossRef]
37. Ehrlich, G.; Hudda, F. Atomic view of surface self-diffusion: Tungsten on tungsten. J. Chem. Phys. 1966, 44, 1039–1049. [CrossRef]
38. Schwoebel, R.L.; Shipsey, E.J. Step motion on crystal surfaces. J. Appl. Phys. 1966, 37, 3682–3686. [CrossRef]
39. Kaufmann, N.A.; Lahourcade, L.; Hourahine, B.; Martin, D.; Grandjean, N. Critical impact of Ehrlich–Schwöbel barrier on GaN
surface morphology during homoepitaxial growth. J. Cryst. Growth 2016, 433, 36–42. [CrossRef]
40. Grandusky, J.; Jindal, V.; Raynolds, J.; Guha, S.; Shahedipour-Sandvik, F. Density functional calculations of the strain effects on
binding energies and adatom diffusion on (0001) GaN surfaces. Mater. Sci. Eng. B 2009, 158, 13–18. [CrossRef]
41. Bortz, A.B.; Kalos, M.H.; Lebowitz, J.L. A new algorithm for Monte Carlo simulation of Ising spin systems. J. Comput. Phys. 1975,
17, 10–18. [CrossRef]
42. Schulze, T.P. Efficient kinetic Monte Carlo simulation. J. Comput. Phys. 2008, 227, 2455–2462. [CrossRef]
43. Kratzer, P.; Penev, E.; Scheffler, M. First-principles studies of kinetics in epitaxial growth of III–V semiconductors. Appl. Phys. A
2002, 75, 79–88. [CrossRef]
44. Lam, C.H.; Lung, M.; Sander, L.M. Fast kinetic Monte Carlo simulation of strained heteroepitaxy in three dimensions. J. Sci.
Comput. 2008, 37, 73–88. [CrossRef]
45. Russo, G.; Smereka, P. A multigrid-Fourier method for the computation of elastic fields with application to heteroepitaxy.
Multiscale Model. Simul. 2006, 5, 130–148. [CrossRef]
46. Tandon, S.; Beleggia, M.; Zhu, Y.; De Graef, M. On the computation of the demagnetization tensor for uniformly magnetized
particles of arbitrary shape. Part II: Numerical approach. J. Magn. Magn. Mater. 2004, 271, 27–38. [CrossRef]
47. Eshelby, J.D. The determination of the elastic field of an ellipsoidal inclusion, and related problems. Proc. R. Soc. Lond. Ser. A
Math. Phys. Sci. 1957, 241, 376–396. [CrossRef]
48. Mura, T. Micromechanics of Defects in Solids; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1987.
49. Xiang, Y.; E, W. Nonlinear evolution equation for the stress-driven morphological instability. J. Appl. Phys. 2002, 91, 9414–9422.
[CrossRef]
50. Pan, E.; Yuan, F. Three-dimensional Green’s functions in anisotropic bimaterials. Int. J. Solids Struct. 2000, 37, 5329–5351.
[CrossRef]
51. Budagosky, J.A.; García-Cristóbal, A. Strain field of a wurtzite type inclusion in a semi-infinite matrix. 2022, to be submitted.
52. Eisenberg, H.R.; Kandel, D. Origin and properties of the wetting layer and early evolution of epitaxially strained thin films. Phys.
Rev. B 2002, 66, 155429. [CrossRef]
53. Fernández-Garrido, S.; Koblmüller, G.; Calleja, E.; Speck, J.S. In situ GaN decomposition analysis by quadrupole mass spectrome-
try and reflection high-energy electron diffraction. J. Appl. Phys. 2008, 104, 033541. [CrossRef]
54. Kalliakos, S.; Bretagnon, T.; Lefebvre, P.; Taliercio, T.; Gil, B.; Grandjean, N.; Damilano, B.; Dussaigne, A.; Massies, J. Photolumi-
nescence energy and linewidth in GaN/AlN stackings of quantum dot planes. J. Appl. Phys. 2004, 96, 180–185. [CrossRef]
55. Song, X.; Feng, H.; Liu, Y.M.; Yu, Z.Y.; Yin, H.Z. Kinetic Monte Carlo simulations of three-dimensional self-assembled quantum
dot islands. Chin. Phys. B 2013, 23, 016802. [CrossRef]
56. Brault, J.; Tanaka, S.; Sarigiannidou, E.; Rouvière, J.L.; Daudin, B.; Feuillet, G.; Nakagawa, H. Linear alignment of GaN quantum
dots on AlN grown on vicinal SiC substrates. J. Appl. Phys. 2003, 93, 3108–3110. [CrossRef]
57. Villain, J. Continuum models of crystal growth from atomic beams with and without desorption. J. Phys. I 1991, 1, 19–42.
[CrossRef]
58. Liang, K.; Sun, X.; Wu, G.; Zhang, L.; Liu, S.; Gan, Z. The investigation of molecular beam epitaxy growth of GaN by molecular
dynamics simulation. Comput. Mater. Sci. 2020, 173, 109426. [CrossRef]
59. Li, R.; Wu, G.; Liang, K.; Wang, S.; Sun, X.; Han, X.; Xue, L.; Li, H.; Liu, S. Effects of AlN substrate orientation on crystalline
quality of wurtzite GaN films investigated via molecular dynamics. Comput. Mater. Sci. 2022, 202, 110991. [CrossRef]

Multiscale Kinetic Monte Carlo Simulation of Self

Uploaded by

Copyright:

Available Formats

Multiscale Kinetic Monte Carlo Simulation of Self

Uploaded by

Document Information

Original Title

Copyright

Available Formats

Share this document

Share or Embed Document

Sharing Options

Did you find this document useful?

Is this content inappropriate?

Copyright:

Available Formats

Multiscale Kinetic Monte Carlo Simulation of Self

Uploaded by

Copyright:

Available Formats

nanomaterials

Publisher’s Note: MDPI stays neutral 1. Introduction

Nanomaterials 2022, 12, 3052. https://doi.org/10.3390/nano12173052 https://www.mdpi.com/journal/nanomaterials

2. Description of the Modeling Strategy

2.1. KMC Model and Implementation

The instantaneous configuration of the film can be specified by an occupation array

where Edes is the corresponding energy barrier.

2.2. Calculation of the Strain-Induced Energy Barrier Estr

distribution of surface tractions, determined by h( x⊥ ) and the lattice mismatch [49]. It is

Finally, we define the energy Estr ( x⊥ ) by:

3.1. Calibration of the Activation Energies by Homoepitaxial Growth Simulations

RHEED intensity (arb.unts)

0.8 (a) ~ 2.25 ML

3.3. QD Growth as a Function of φGa /φN Ratio

3.4. QD Growth as a Function of Substrate Temperature

surface diffusion is a thermally activated process, significant changes in adatom mobility

QD density (1011 cm-2)

680 700 720 740 760 780

3.5. Growth of Stacks of GaN/AlN QDs: Correlation Effects

KMC Kinetic Monte Carlo

At the stoichiometry situation, the growth rate reaches the value:

0.4 (a) TS = 730 ºC

0.4 (b) TS = 770 ºC

Nominal Ga flux Ga (ML/s)

You might also like