Order of Epitaxial Self-Assembled Quantum Dots: Linear Analysis
Lawrence H. Friedman
November 4, 2018
Dept. of Engineering Science and Mechanics, Pennsylvania State University, 212 Earth and Engineering Science
Building, University Park, Pennsylvania 16802
lfriedman@psu.edu
keywords: quantum dots, strained films, epitaxial growth, semiconductors
Abstract
Epitaxial self-assembled quantum dots (SAQDs) are of interest for nanostructured optoelectronic and electronic
devices such as lasers, photodetectors and nanoscale logic. Spatial order and size order of SAQDs are important
to the development of usable devices. It is likely that these two types of order are strongly linked; thus, a study of
spatial order will also have strong implications for size order. Here a study of spatial order is undertaken using a linear
analysis of a commonly used model of SAQD formation based on surface diffusion. Analytic formulas for film-height
correlation functions are found that characterize quantum dot spatial order and corresponding correlation lengths that
quantify order. Initial atomic-scale random fluctuations result in relatively small correlation lengths (about two dots)
when the effect of a wetting potential is negligible; however, the correlation lengths diverge when SAQDs are allowed
to form at a near-critical film height. The present work reinforces previous findings about anisotropy and SAQD order
and presents as explicit and transparent mechanism for ordering with corresponding analytic equations. In addition,
SAQD formation is by its nature a stochastic process, and various mathematical aspects regarding statistical analysis
of SAQD formation and order are presented.
1 Introduction
Epitaxial self-assembled quantum dots (SAQDs) represent an important step in the advancement of semiconductor
fabrication at the nanoscale that will allow breakthroughs in optoelectronics and electronics. [1, 2, 3, 4, 5, 6, 7, 8, 9,
10, 11, 12] Most frequent optoelectronic applications are high efficiency lasers with exotic wavelengths or photode-
tectors. [1, 3, 4, 5, 6, 7, 8, 10, 11, 12] SAQDs are the result of a transition from 2D growth to 3D growth in strained
epitaxial films such as SixGe1−x/Si and InxGa1−xAs/GaAs. This process is known as Stranski-Krastanow growth or
Volmer-Webber growth. [13, 1, 14, 15].
In applications, order is a key factor. There are two types of order, spatial and size. Spatial order refers to
the regularity of SAQD dot placement, and it is necessary for nano-circuitry applications. Size order refers to the
uniformity of SAQD size which determines the voltage and/or energy level quantization of SAQDs. It is reasonable
to expect that these type of order are linked, and it is important to understand the factors that determine SAQD order.
Further understanding should help in the design and simulation of both spontaneous “bottom up” self-assembly and
directed or guided self-assembly to enhance SAQD order. [16, 17, 18, 19, 20, 21, 22, 23]Here, an elaboration of and
further application of a linear analysis of SAQD order [24] is presented. The work reported here forms the basis of a
non-linear theory and modeling of SAQD order that will be reported in future work.
In [24] it was reported that one could calculate a correlation function using a linearized model of SAQD formation.
This correlation function included two correlation lengths that could be used to describe SAQD order. It was also found
that one effect of a hypothesized wetting potential was to enhance SAQD order when growth occurs near the critical
film height for 3D growth. Here, these results are expanded to create a more rigorous linearized theory of SAQD order
that will inform non-linear theories. In particular, the model is generalized to any model that combines local energy
effects such as surface energy density and non-local elastic destabilization, and the procedure for predicting order
based on any linear theory with peak wavelengths is presented. The hypothesized effect of elastic anisotropy in [24]
is verified with calculations using linear anisotropic elasticity theory. [25, 26] Details such as statistical fluctuation
and convergence are also addressed along with a discussion of the possible forms of linear anisotropic terms in SAQD
growth kinetics, and the effect of an atomic-scale cutoff in the continuum theory is addressed. Finally, the order
enhancing effect of growing near the critical threshold is explored in more detail using calculations appropriate to
Ge/Si SAQDs.
In the literature, two modes of SAQD formation are generally discussed, the thermal nucleation mode and the nu-
cleationless mode. [27, 28, 29] In the thermal nucleation mode, a 2D film surface is metastable, and the formation of
individual quantum dots is thermally activated. [27]. This growth mode leads to the formation of individual quantum
dots as uncorrelated or loosely correlated discrete events at essentially random locations. In the nucleationless mode,
the 2D film surface transitions from stable (or metastable) to unstable. In this mode, dots form everywhere at once
appearing at first as a cross-hatched ripple-like disturbance on the 2D film surface and then maturing into recogniz-
able individual dots.[27, 30, 28, 31, 32] These two modes are probably connected via an encompassing conceptual
and mathematical model1, and perhaps some of what is observed experimentally is in fact a hybrid mechanism. In
agreement with intuition, it appears that the nucleationless mode leads to a more ordered dot pattern than the thermal
nucleation mode that is dominated by randomness. 2 Thus, the presented analysis applies to the nucleationless mode.
There are various implementations of nucleationless growth models [28, 37, 38, 39, 40, 18, 34], although, there is
also a great deal of commonality among these models. In particular, they all include a non-local elastic effect and local
surface energies and/or local wetting energies. Here, a linear analysis of quantum dot order resulting from this class of
model is presented. Particular note is taken of the effects of stochastic initial conditions crystal anisotropy in general,
elastic anisotropy in particular, and the effect of varying film height as a control parameter as first introduced in [33].
A simple model similar to [28, 37, 38, 40, 18] is presented to produce numerical examples and explore the effects
of the average film height. Concurrently, a more abstract and general model is presented and analyzed that includes
non-local elastic strain effects, and a local combined surface and wetting energy. The linear model with stochastic
initial conditions and deterministic film height evolution will pave the way for more sophisticated analysis involving a
non-linear model of stochastic film height evolution.
As previously stated, one of the goals in the present work is to further explore the role of the wetting poten-
tial during growth near the stability threshold in film height. A wetting potential has been included in the analysis
and simulations in [38, 33, 37, 28]. Although somewhat controversial, the wetting potential plays an important phe-
nomenological role. It ensures that growth takes place in the Stranski-Krastanow mode: that a 3D unstable growth
occurs only after a critical layer thickness is achieved, and that a residual wetting layer persists. The physical origins
and consequences of the wetting potential are discussed in [41, 28]. The analysis presented here is usable in models
that neglect the wetting potential by simply setting it to zero. Another possibility is simply that the wetting potential
is simply an approximation to the stabilizing effect of intermixing. [42] That said, if the wetting potential is real, the
present analysis shows that it is beneficial to SAQD order to grow near the critical layer thickness.
The presented analytic formulas and linear analysis are intended to complement existing numerical models of
SAQD order. [43, 37, 44, 45] and to form a basis for future non-linear analytic analysis of SAQD order. The current
findings agree with previous work on the beneficial effects of elastic anisotropy to enhance in-plane order.
The linear analysis, of course, represents a simplification of the film evolution, and it applies only to the initial
stages of SAQD formation when the nominally flat film surface becomes unstable and transitions to three-dimensional
growth. However, the small surface fluctuation stage of SAQD growth determines the initial seeds of order or disorder
in an SAQD array; thus, the small fluctuation stage should have an important influence on the final outcome. At later
stages when surface fluctuations are large, there is a natural tendency of SAQDS to either order or ripen [33, 37, 46,
39, 47] Ordering systems tend to evolve slowly due to critical slowing down [39], while ripening tends to diminish
order further. [37] Thus, it is possible that the linear model could, in fact, yield good predictions of SAQD order.
The simplification and linearizion facilitates the development of analytic solutions that are most transparent, easily
portable to multiple material systems and have no effective limit on system size. Finally, it is virtually impossible to
have a thorough understanding of the full non-linear model without first having a thorough understanding of the linear
behavior.
The remainder of the paper is organized as follows. Section 2 presents the physical assumptions and mathematical
1It is likely that there a transition from stable, to metastable and finally to unstable. The analysis presented in [33] would appear to support such
a view where the film height acts as the control parameter driving the transition. There is also some controversy regarding whether all dot growth is
nucleationless or not. [34, 32]
2Compare various figures in [29, 35, 14, 31, 36].
approximations used to model film growth. Section 3 discusses the stochastic initial conditions and the resulting
correlation functions and correlation lengths. Section 4, presents a procedure for estimating SAQD order with an
application to Ge dots on a Si substrate. Section 5 presents conclusions, while Appendices A-F present additional
calculational details.
2 Modeling
The formation of SAQDs is modeled as a deterministic surface diffusion process with stochastic initial conditions. The
resulting equations and ultimately the sought after correlation functions are different depending on whether the film
surface is treated as one-dimensional isotropic, two-dimensional isotropic or two-dimensional anisotropic. The 1D
and 2D isotropic cases are discussed first, and then the essential differences of the 2D anisotropic model are presented.
The stochastic initial conditions need to be expressed in terms of the correlation functions that are also use to analyze
order; consequently, the discussion of the initial conditions is deferred to Sec. 3.2.
It should be noted that the results presented here are fairly general. There has been a good deal of recent work
refining the modeling of nucleationless growth processes to incorporate various phenomenological aspects of SAQD
growth. For example, the inclusion of orientation-dependent surface energy [38], strain-dependent surface energy [34]
and explicit modeling of atomic species segregation and film-substrate inter diffusion. [48] Two models are presented
here. One is a simple concrete example. It is the simplest model one can use including elastic effects surface energy
and wetting energy. The second model is more abstract and describes the general case of a local potential energy
that depends on both the film height and film height gradient. One effect that is not examined here is that of mixed
4-fold and two-fold symmetry. Such a mixing can occur due to diffusional anisotropy or surface energy anisotropy.
(Sec. 2.2.1.2 and Appendix D). However, a similar analysis procedure should work for these cases as well. The general
procedure for possible application to other models is discussed in Sec. 3.5.
The following discussion will use abstract vector notation, e.g. k instead of ki, etc. Also, because it is sometimes
computational expedient to perform one-dimensional modeling [24, 39, 17, 42], the case of a one dimensional surface
with two dimensional volume is discussed along with the case of an isotropic 2D surface. To facilitate this combined
discussion, the dimensionality of the surface will be denoted as d. In Secs. 3.3 and 3.4, d = 1, 2 will be substituted as
appropriate. Finally, much of the calculation involves reciprocal space. The convention used for the Fourier transforms
f(x) =
ddk eik·xfk, and fk = (2π)
ddx e−ik·xf(x)
following the example of [28].
2.1 1D and 2D Isotropic model
This discussion pertains to both the 1D model and the 2D isotropic model. The formation of SAQDs is modeled as a
surface diffusion process where the film height is a function of the lateral position and time. The system is treated as
deterministic with stochastic initial conditions. First, the general non-linear governing equations are presented. Then,
the linearized form is presented. Finally, the key behavior is reviewed.
The mathematical model uses film height, H(x, t) as the dependent variable and the horizontal position x and
time t as the independent variables. The film height evolves over time due to surface diffusion driven by a diffusion
potential µ(x, t) and a flux of new material Q. The surface velocity is thus
vn = nz∂tH = −∇S · D∇Sµ(x, t) +Q (1)
where nz is the vertical component of the surface normal nz = [1 + (∇H)2]−1/2, ∇S is the surface gradient, D is the
diffusivity, and ∇S · is the surface divergence.
2.1.1 Energetics
The diffusion potential µ(x, t) must produce Stranski-Krastanow growth. Thus, it must contain an elastic term that
destabilizes film growth, a surface energy term that stabilizes planar growth and a wetting energy that ensures a wetting
layer. The diffusion potential can be derived from a total free energy.
F = Felast + Fsurf. + Fwet
volume
dV ω +
surface
dAsurf. γ +
dAW (H)
where ω is the elastic energy density, γ is the surface energy density, W (H) is the wetting energy density. The
last integral corresponds to Fwet, and whether the integral should be taken over the film surface or the substrate is
ambiguous. The “simple” model (Sec. 2.1.1.1) assumes that the integral is over the substrate, while the “general”
model (Sec. 2.1.1.2) can accommodate both cases.
2.1.1.1 simple form The simplest possible model results if the integral corresponding to Fwet is taken over the
lateral positions x rather than over the actual free-surface. In concrete terms, one can use dV = d2xdz and dAsurf. =
1 + (∇H(x))2
to obtain the expression,
volume
d2xdz ω[H](x, z) +
x-plane
1 + (∇H(x))2
γ +W (H(x))
, (2)
where the “ω[H]” indicates that the elastic energy density is a non-local functional of the film height,H. The diffusion
potential µ can be found, similar to [15], by differentiating F with respect to the surface motion (Appendix A.1),
µ(x) = ΩδF/δH(x). Doing so for Eq. (2) (Appendix A.2),
µ(x) = Ω [ω(x)− γκ(x) +W ′ (H(x))] . (3)
where Ω is the atomic volume, ω(x) is the elastic energy density at the film surface (implicitly ω[H] (x,H(x))),
κ = ∇ ·
∇H(x)
1 + (∇H(x))2
]−1/2}
is the total surface curvature, and W ′(H) = ∂H(x)W (H(x)) is the
derivative of W (H(x)) evaluated at x.
2.1.1.2 general form It should be noted that Eq. (3) is not the same diffusion potential used in [38]. The wetting
potential used there can be derived by taking W (H) as an energy density of the free surface, not a density in the
x-plane. Expressions like Eq. (3) and Eq. (1) in [38] are part of a larger class of surface evolution models with more
or less the same linear behavior.
The surface and wetting energy can be combined and incorporated into a more general form, with a total free
energy Fsw and a free energy density Fsw(H,∇H) that depends on the film heightH(x) and the film height slope or
orientation ∇H(x). The total free energy is thus
F = Felast. + Fsw (4)
volume
d2xdz ω[H](x, z) +
x−plane
d2xFsw (H(x),∇H(x)) .
Fsw may not necessarily be the sum of separate surface energy and wetting energy contributions. It need only be a
local function ofH and ∇H. The corresponding diffusion potential is
µ(x) = Ω
ω(x) + F (10)sw (x)−∇ · F
sw (x)
, (5)
where F (mn)sw indicates the mth derivative with respect to H and the nth derivative with respect to ∇H. F
sw (x) =
∂H(x)Fsw (H(x),∇H(x)) and each vector component of F
sw (x) is
F(01)sw (x)
= ∂[∇H(x)]
Fsw (H(x),∇H(x)).
One can obtain the results of the simple model (Eqs. (2) and (3)) by setting
Fsw =
1 + (∇H(x))2
γ +W (H(x)) . (6)
A diffusion potential like Eq. (1) in [38] can be obtained by setting
Fsw =
1 + (∇H(x))2
[γ (∇H(x)) +W (H(x))] .
This is different from Eq. (6) in two ways. First, the surface energy density depends on the surface orientation. Second,
the Jacobian, J =
1 + (∇H(x))2
multiplies both the surface energy density and the wetting potential. Despite
these differences, the common form of the diffusion potential (Eq. (5)) among different models suggests that they
might all lead to similar linearized forms and behavior.
2.1.1.3 Linearization The diffusion potential is now linearized about the average film height H̄. In general, one
can control the amount of deposited material, and thus the average film height H̄. It is therefore useful to decompose
H(x) into the spatially averaged mean value and fluctuations about the average. Similar to [28],
H = H̄+ h(x, t). (7)
In the present calculation, H̄ is specified as constant in time. This assumption corresponds physically to a fast deposi-
tion and then an anneal. It is not too difficult to generalize to a time dependent H̄, but that is beyond the scope of this
manuscript. In [38, 49], deposition and evaporation is explicitly modeled.
All terms in µ(x, t) are now kept to only linear order in h(x, t). The elastic energy density ω is a non-local
functional of h(x, t) [40]; however, the equations generating ω(x) are translationally invariant. Thus, it is convenient
to use reciprocal space for the linearization. The curvature is trivially linearized as κ(x) → ∇2h(x) in real space or
κk → −k2hk in reciprocal space. The linearized elastic strain energy ω can be found in reciprocal space as in [15]
to be ωk = −2M(1 + ν)�2mhk, where M = E/(1 − ν) is the biaxial modulus, E is the Young modulus, ν is the
Poisson ratio, and �m is the film-substrate mismatch strain. This formula neglects possible differences in elastic moduli
between the film and substrate as in [28], but a similar method of analysis should apply to that case as well. Linearizing
Eqs. (3) and (5) in reciprocal space, µk is proportional to hk with a proportionality coefficient that depends on k and
µlin,k = f(k, H̄)hk (8)
where f(k, H̄) for three different isotropic cases, corresponding to Eqs. (3) and (5), and an abstracted general form, is
given by
f(k, H̄) =
−2M(1 + ν)�2mk + γk2 +W ′′(H̄)
; case a (Eq. (3))
−2M(1 + ν)�2mk + F 02k2 + F 20
; case b (Eq. (5))
−ak + bk2 + c ; case c (general)
. (9)
Due to isotropy, f(k, H̄) is independent of the direction of k, and only the wave number, k = ‖k‖, appears in the
right hand side. F (20)sw is the second derivative of Fsw with respect to H, and F
sw the second derivative of Fsw with
respect to ∇H. F (20)sw and F
sw depend on H̄ only; thus they are constants in the present analysis. See Appendix B.2
for more precise definitions and the derivation of f(k, H̄). Using Eq. (6), produces F (02)sw = γ and F
sw = W ′′(H̄)
which is identical to the simple case of Eq. (9), a. Case c, labeled as “general” where a, b, and c depend implicitly on
H̄ shows that f(k, H̄) for cases a and b have the same relatively simple form. It also emphasizes the dynamic effects
as opposed to the physical causes. There is a destabilizing term, −ak, a short wavelength cutoff term, bk2, and a term
that stabilizes the entire spectrum, c.
Despite the label “general,” there are of course limitations to the application of Eqs. (8) and (9). For example,
there has been recent work on the effects of strain-dependent surface energies. [34] The second form can not represent
such an effect because the derivation assumes that the surface energy only depends on local quantities, (H and ∇H)
whereas the strain effect is non-local. However, it is reasonable to conjecture that a more detailed analysis of the
effects of a strain dependent surface energy term would produce a coefficient function f(k, H̄) not very different from
the case c “general” form of Eq. (9). Thus, the following analysis may very well apply to this more exotic model, but
more study is needed to be certain.
2.1.2 Dynamics
As discussed in Sec. 2.1.1, the dynamics are derived assuming no flux of new material (Q = 0) and keeping only
terms to linear order in the height fluctuation, h(x, t). Under these assumptions, Eq. (1) can be decomposed into a
Table 1: Characteristic wave-numbers, characteristic times and associated dimensionless variables for the three cases
addressed in Eq. (9).
kc tc α β
case a 2M(1+ν)�
16DΩM4(1+ν)4�8m
γW ′′(H̄)
4M2(1+ν)2�4m
case b 2M(1+ν)�
(F (02)sw )
16DΩM4(1+ν)4�8m
F (02)sw F
4M2(1+ν)2�4m
case c a/b b3/(DΩa4) k/kc cb/a2
trivial equation for H̄ and an equation for the film height fluctuation by inserting Eq. (7).
dH̄/dt = 0 (10)
∂th(x) = −∇ · D∇µlin(x) (11)
where µlin(x) is the inverse Fourier transform of Eqs. (8) and (9), and it depends implicitly on the average film height
H̄. Note that the time dependence is implicitly while the coordinate dependence is explicit. The explicit coordinate
dependence serves to distinguish Assuming that the diffusivity D is constant, the Fourier transform of Eq. (11) gives
the linearized differential equation for the evolution of each Fourier component.
∂thk = −Dk2µk = −Dk2f(k, H̄)hk. (12)
Solving Eq. (12),
hk(t) = hk(0)e
σkt; (13)
σk = −Dk2f(k, H̄). (14)
The surface evolves in reciprocal space as an initial condition, hk(0) multiplied by an envelope function, eσkt. For
most values of H̄, this envelope function has a peak. As time passes, this peak narrows and can be approximated by a
gaussian. To analyze this behavior, appropriate dimensionless variables are defined. Then, the stability of the film is
discussed. Finally, σk is expanded about its peak to aid analytic calculations.
The time dependent behavior of the film height fluctuations is facilitated by using a characteristic wave number,
characteristic time and related dimensionless variables. For the “general” case c of Eq. (9), the characteristic wavenum-
ber is kc = a/b, and the characteristic time is tc = 1/(DΩbk4c ) = b3/(DΩa4). These characteristic dimensions can be
used to define a dimensionless wave vector, α = k/kc and a dimensionless wetting parameter β = c/(bk2c ) = cb/a
One can also define a dimensionless time, τ = t/tc. To obtain the corresponding characteristic scales for cases a and
b, one merely has to plug in the appropriate substitutes for a, b and c and follow the pattern. For example, for case
a, make the substitution a → Ω2M(2 + ν)�2m, etc. Table 1 summarizes these values for all three cases. For all three
cases, f(k, H̄) and the growth constant σk reduce to the following forms:
f(k, H̄) = f(kcα, H̄) = Ωbk2c
−α+ α2 + β
σk = σkcα = t
α− α2 − β
, (16)
where α = ‖α‖ = k/kc is the dimensionless wave number. These forms are plotted in Figs. 1a and 2. Fig. 1a shows
f(k, H̄)/Ωbk2c vs. α for an isotropic or one dimensional surface. Figs. 2 shows tcσk vs. α for a 2D anisotropic
surface (Sec. 2.2). However, the curves marked 0◦ are identical to the dispersion relation for a 1D or 2D isotropic
surface (compare Eqs. (9) and (23)).
2.1.3 Peaks
The peak growth rate and the corresponding wavenumber k can be found from Eq. (16). σk has a peak at k0 = kcα0
where
9− 32β
. (17)
Expanding σk about this peak to second order in k − k0,
σk ≈ σ0 −
σ2(k − k0)2
! ! 0.25!
increasing θk
Figure 1: Dimensionless diffusion potential prefactors vs. dimensionless wave number. (a) The one dimensional or
isotropic case with β = 0.3. (b) The elastically isotropic case with anisotropy �A = 0.1 (see Eq. (22)).
increasing θk
22.5◦
increasing θk
22.5◦
Figure 2: Dimensionless growth constant vs. dimensionless wave number. Curves are plotted for the elastically
anisotropic case, but the curves marked 0◦ are the same as for the isotropic cases. In (a), β = 0. In (b) β = 0.2.
Figure 3: Exponential Envelope eσkt as function of α for β = 0.208 and t/tc = 100. (a) 2D isotropic surface. (b) 2D
anisotropic surface with � = 0.1236.
The two constants are
t−1c α
0 (α0 − 2β) , (18)
σ2 = t
c (3α0 − 4β) . (19)
Inserting this approximation for σk into Eq. (13),
hk(t) = hk(0)e
σ0te−
2σ2t(k−k0)
. (20)
The individual initial surface fluctuation components grow with a gaussian shaped envelope. An example of this
envelope is plotted in Fig. 3(a). Notice that in two dimensions, the envelope forms a ring as the peak is about the
wave-number k0 but not about any particular point in the k-plane.
2.1.4 Stability and wetting potential
Stranski-Krastanow growth is marked by a transition from stable two-dimensional growth to unstable three-dimensional
growth once a critical height Hc is reached. [1] Eqs. (17), (18) and (20) are useful for analyzing the transition from
stable to unstable growth. In order for this transition to occur, there must be some stabilizing term in the diffusion
potential. In the present model, this means that there must be some surface energy-like term that varies strongly with
film height. This condition equates to stating that W ′′(H̄) or F 20sw or c (Eq. (9)) must be rather large if H̄ < Hc.
However, as H̄ increases, these terms are reduced. Finally, when H̄ > Hc, this term is no longer capable of stabilizing
the film against fluctuations of all possible wavelengths.
The critical value Hc can be found using the analysis from [33]. By inspection of Eqs. (8), (9) and (12), modes
with f > 0 increase the total free energy F as they grow; thus, they are stable and decay with time. Modes with
f < 0 decrease the total free energy F as they grow; thus, they are unstable and grow with time. This growth and
decay rule is easily verified by inspection of Eq. (14). Thus, stable growth occurs when f(k, H̄) > 0 for all values of
k, and unstable growth occurs when f(k, H̄) < 0 for some values of k. Thus, the transition from stable to unstable
growth occurs when the minimum value of f(k, H̄) just becomes negative. Using the same dimensional analysis as in
the previous section and following the discussion of [33], one finds that the minimum value, fmin = Ωbk2c (β − 1/4),
occurs at kmin/kc = αmin = 1/2. fmin first becomes negative, and the transition to unstable growth occurs when
the dimensionless wetting parameter (Table 1) drops to a critical value, β = 1/4 . β > 1/4 stable 2D growth, and
β < 1/4 unstable 3D growth. It is reasonable to suppose that W (H̄), W ′′(H̄), and thus β are positive monotonically
decreasing functions of H̄ so that the interface becomes less important for large values of H̄. For example, in [50] it is
assumed that W (H) = B/H, where B is constant. When β → 0, corresponding to large H̄, the case discussed in [28]
is obtained. A similar analysis can be done for cases b and c once one specifies how the terms F (20)sw and F
sw or a, b
and c depend on H̄.
Using a guessed form for a wetting potential, one can find the critical film heightHc by setting β = 1/4 . Applying
this condition to case a in Eq. (3)
W ′′(Hc) = γk2c/4.
Using the wetting potential of [50] as an example, W (H) = B/H,
Hc = 3
8B/(γk2c ) =
8Bγ/(2M(1 + ν)�2m)2. (21)
Conversely, one can fit a wetting potential to an observed or reasonable critical layer thickness from the same condition.
Using the example wetting potential from [50],
(2M(1 + ν)�2m)
as stated in [50].3
2.2 2D Anisotropic case
Crystal anisotropy leads to a dispersion relation σk that is both quantitatively and qualitatively different from the
isotropic case. Here the effect of elastic anisotropy is discussed in most detail. Other sources of anisotropy are
the surface and wetting energies. For example, in [38] the surface energy density is orientation dependent which
introduces a possible anisotropy in the dispersion relation. Possible sources of anisotropy are an anisotropic elastic
stiffness tensor, an orientation dependent surface energy or wetting potential or anisotropic diffusion. As discussed
below, the form of anisotropy to linear order in the height fluctuation, h, is somewhat restricted. Results are presented
for 4-fold symmetric surfaces, that is surfaces that have invariant dynamic evolution laws when rotated by 90◦. Possible
complications arising from 2-fold symmetric anisotropic terms (with 180◦ rotational symmetry) are also discussed. As
for the isotropic case, first the energetics are discussed, then the dynamics, and finally the expansion about the peaks
in the dispersion relation, σk.
2.2.1 Energetics
The discussion of energetics will first treat the effects of elastic anisotropy and then anisotropy resulting from surface
or wetting like terms.
3This result from [50] corresponds to the choice Fsw(H,∇H) =
1 + (∇H)2γ+W (H). However, the numerical model in [50] appears to
use Fsw(H,∇H) =
1 + (∇H)2 [γ(∇H) +W (H)]. This difference should lead to a slightly different critical film height in their numerical
model from the one that they predicted (Eq. (21)).
Figure 4: Plot of Eθk/(M�
m) for various materials. Symbols indicate values calculated using Appendix C. Solid lines
are the interpolation (Eq. (22)) using the values from Table 2.
Table 2: Elastic constants [51] and calculated values (see Appendix C) for various materials of interest at T = 300K.
c11 c12 c44 M
1011 ergcm3 10
11 erg
cm3 10
11 erg
cm3 10
11 erg
Ge 12.60 4.40 6.77 13.93 2.16 1.906 0.1176
Si 16.60 6.40 7.96 18.07 2.22 1.997 0.1005
InAs 8.34 4.54 3.95 7.94 2.70 2.09 0.226
GaAs 11.90 5.34 5.96 12.45 2.15 1.87 0.1302
2.2.1.1 Elastic anisotropy One would like to obtain a simple symbolic expression for the elastic energy density at
the free surface, ωk, to first order in hk for the elastically anisotropic case. Similar discussions can be found in [25, 26].
For the isotropic case, ωk = −2M(1 + ν)�2mhk. For the anisotropic case,
ωk = −Eθkkhk
where the prefactor Eθk is the decrease in elastic energy at the surface per unit wave number (k → 1) and unit amplitude
(hk → 1) . It is not constant, but instead depends on the θk, the angle that k makes with the x−direction. The case of a
cube-symmetry elastic stiffness tensor such as for Si is considered where one must specify three elastic constants c11,
c12 and c44. [51]. Growth on a (100) surface will produce an elastic energy prefactor Eθk that is four-fold symmetric
(symmetric upon rotations by 90◦). A procedure similar to [25, 26] based on a first order perturbation analysis is
followed (Appendix C). A relatively simple interpolation formula [24] is hypothesized and then verified numerically.
The interpolation procedure, suggested in [24] uses the lowest possible order expansion in sin(θk) and cos(θk)
that has the appropriate four-fold symmetry and then interpolates between θk = 0◦ and θk = 45◦. Thus,
Eθk = E0◦
1− �A sin2 (2θk)
where �A = (E0◦ − E45◦)/E0◦ is an anisotropy factor. This lowest order form turns out to be a very good fit to
numerical calculations (Fig. 4). Table 2 gives values of E0◦ and �A for some systems of interest. In the elastically
isotropic case, E0◦ = E45◦ = 2M(1 + ν) so that �A = 0.
There are two important differences from the elastically isotropic case. The first is obvious, that Eθk depends on
angular orientation, θk. The second is that the peak value of ωk is not the same as that for the elastically isotropic case
because in general, E0◦ 6= 2M(1 + ν). In [24], where the purpose was simply to investigate the mechanism by which
elastic anisotropy effects order, this second difference was neglected.
2.2.1.2 Surface and Wetting Energy Anisotropy The surface energy and wetting potential can be additional
sources of anisotropy if they depend on the surface orientation so that γ → γ(∇H) or W (H) → W (H,∇H) (for
example, [52, 38]). Then, to first order in h ,
µsurf.,k = Ω
γk2 + k · γ̃′′ · k
where γ̃′′ is the (2×2) matrix or Hessian matrix that results from taking the second derivatives of γ(∇H) with respect
to the two components of ∇H (Appendix B.1). Similarly
µwet,k = Ω
W (20) + k · W̃(02) · k
where W (20) and W̃(02) are the second derivatives of W (H,∇H) with respect to H and ∇H (Appendix B.1). For
both µsurf.,k and µwet,k, the first term is isotropic, and the second term contains any possible anisotropy.
The rank of the γ̃′′ and W̃(02) matrices greatly restricts the possible forms of the additional anisotropy. These
(2 × 2) matrices must be either two-fold symmetric or perfectly isotropic. Thus, if the surface energy and wetting
potential are four-fold symmetric as Eθk is, then γ̃
′′ → γ′′, a scalar, and W̃(02) → W (02), a scalar, and neither one
contributes any additional anisotropy. They do, however, help to stabilize or further destabilize the 2D surface as they
add terms proportional to k2. The effect of these additional terms is indistinguishable from the effect of varying the
value of the surface energy density, γ. [52, 31]
It should be noted that the (100) surface of a diamond or zinc-blend structures allows for anisotropy that is only 2-
fold symmetric (rotations by 180◦). Thus, they could “break” the four-fold symmetry that occurs when one considers
the elastic anisotropy alone. However, this “broken” symmetry is somewhat dubious because even the diamond and
zinc-blend structures have a screw symmetry (rotations by 90◦ and translation in the [100] direction by half a lattice
vector). Thus, if for example, W (H,∇H) is anisotropic with two-fold symmetry to linear order, there must be a
fast oscillation with changes in the film height H. In Appendix D, a similar term related to anisotropic diffusion is
discussed. There does not appear to be any evidence for this two-fold symmetry in the case of (100) surfaces of IV/IV
systems such as Ge/Si, but in III-V/III-V systems the four-fold symmetry of the (100) surface may indeed be “broken”
in this way corresponding to either a surface energy anisotropy or a diffusional anisotropy. [53, 54]. Further analysis
of such terms in any more detail would greatly complicate the present discussion, so it is left for future work. Most
of the modeling literature avoids this complication by not including the symmetry-breaking of the zinc-blend surface,
for example [25, 26, 38].
One can perform a similar analysis of the combined surface and wetting potential, Fsw(H,∇H) (case b). To linear
order the resulting anisotropic diffusion potential is (Appendix B.2)
µsw,k = Ω
F (20)sw + k · F̃
sw · k
Again, F̃(02)sw is a rank 2 tensor, and all of the same symmetry considerations apply here as well.
Because the two-fold symmetry anisotropic terms are excluded from the current discussion, and isotropic terms
simply “renormalize” the effective of surface energy, there will be no further consideration of anisotropy resulting
from the surface energy or wetting potential in this discussion. Further calculations will proceed assuming that the
surface energy density, γ, nor the wetting potential,W (H), depend on ∇H or similarly that Fsw(H,∇H) has a purely
isotropic dependence on ∇H. This assumption can be made without affecting any of the qualitative results.
2.2.1.3 total diffusion potential Having dispensed with the discussion of the various sources of anisotropy, the
total diffusion potential is stated for the case of 4-fold symmetric elastic anisotropy and a completely isotropic surface
energy and wetting potential. µk = f(k, H̄) with
f(k, H̄) =
1− �A sin2(2θk)
k + γk2 +W ′′(H̄)
; case a (Eq. (3))
1− �A sin2(2θk)
k + F (02)sw k2 + F
; case b (Eq. (5))
1− �A sin2(2θk)
+ bk2 + c ; case c (general)
. (23)
Table 3: Characteristic wave-numbers, characteristic times and associated dimensionless variables for the three cases
addressed in Eq. (9)
kc tc α β
case a E0◦/γ γ3/(DΩE40◦) k/kc γW ′′(H̄)/E20◦
case b E0◦/F
/(DΩE40◦) k/kc F
sw /E20◦
case c a/b b3/(DΩa4) k/kc cb/a2
2.2.2 Dynamics
The dynamics is governed by surface diffusion, just as for the fully isotropic case. It is assumed that the diffusivity
is isotropic as was done for the surface energy and the wetting energies; thus, all anisotropy in the film evolution
dynamics comes from elastic effects alone. The possibility and effects of an anisotropic diffusion potential is discussed
in Appendix D (also see [54]). The time dependence of the surface perturbations simply follows Eqs. (13) and (14),
but with Eq. (23) used for f(k, H̄). As for the isotropic case, appropriate characteristic wave numbers (kc) and
time scales (tc) can be found for each of the three cases along with the associated dimensionless wave vector α and
dimensionless wetting parameter β. These are listed in Table 3. The dispersion relation, σk can be expressed in terms
of these dimensionless variables (α and β), giving
σk = σkcα = t
1− �A sin2(2θk)
− α2 − β
. (24)
The stability behavior is essentially the same as for the isotropic case with a transition occurring at β = 1/4 corre-
sponding to H̄ = Hc.
2.2.3 Expansion about peaks
σk has 4 peaks at (k, θk) = (k0, π[n− 1]/2) with k0 = kcα0 (Eq. (17)) and n = 1 . . . 4. In vector form, there are four
peaks at
kn = k0 (cos(π(n− 1)/2)i + sin(π(n− 1)/2)j) .
Similar to the isotropic case, σk can be expanded about individual peaks so that in the vicinity of peak n, σk ≈ σn
σn = σ0 −
σ‖(k − k0)2 −
0(θk − nπ/2)
where σ0 is given by Eq. (18), σ‖ = σ2 given by Eq. (19), and
σ⊥ = 8�Aα0t
In terms of the vector components parallel and perpendicular to kn, k‖ and k⊥ respectively,
σn = σ0 −
σ‖(k‖ − k0)2 −
k‖ = cos[π(n − 1)/2]kx + sin[π(n − 1)/2]ky , and k⊥ = − sin[π(n − 1)/2]kx + cos[π(n − 1)/2]ky . The time
evolution of hkin the vicinity of one of the kn is
hk(t) ≈ hk(0)et(σ0−
2σ2(k‖−k0)
2− 12σ⊥k
3 Correlation Functions
Correlation functions and associated constants such as correlation lengths can be very useful for characterizing order.
In particular, the autocorrelation function (Eq. (25)) and its Fourier transform (Eq. (26)) also known as the spectrum
function can give a very good characterization of dot order (Figs. 6a and c and 5b, e and h). The autocorrelation
function is denoted CA(∆x) where ∆x is the difference vector between two points in the x−plane. The spectrum
function is a function of k, and it is denoted CAk . The goal here is to be able to predict these two functions and to
h(nm)
h(nm)
h(nm)
describe them quantitatively in a manner that can be used to characterize SAQD order with just a few numbers. The au-
tocorrelation function is the result of a spatial average over one experiment or one simulation (numerical experiment).
It is regular and repeatable because it is closely tied to the correlation function and spectrum function that results from
an ensemble average (Eqs. X and X). These are denoted as C(∆x) and the spectrum Ck respectively. Note that the
ensemble averaged functions do not have a superscript “A.” These ensemble average correlation functions are useful
in the analysis of stochastic ordinary and partial differential equations. [55, 56]. From a strictly technical viewpoint,
the spatial average and the ensemble average are not exactly the same; however, they are closely enough connected
that it is reasonable to use one as a substitute for the other (Sec. 3.1 and Appendix 3).
In the following, the analysis of SAQD order via autocorrelation and correlation functions is discussed (Sec. 3.1).
Then, the stochastic initial conditions are discussed (Sec. 3.2). Then, the prediction of the Fourier transforms of the
correlation functions is discussed (Sec. 3.3). The real-space correlation functions are presented (Sec. 3.4). Finally,
there are some notes regarding generalizing the analysis method to any dispersion relation that has peaks (Sec. 3.5),
for example, peaks related to broken four-fold symmetry or growth on a miscut substrate.
3.1 Correlation Functions and SAQD order
Auto-correlation functions are well-suited for investigating SAQD order. The autocorrelation function is defined as
CA(∆x) =
d2x′ h(∆x + x′)h(x′)∗. (25)
Its Fourier transform sometimes called the spectrum [56], spectrum function or power spectrum is
CAk =
(2π)d
d2∆x e−ik·∆xC(∆x) =
(2π)d
, (26)
where A is the projected area of the film in the x − y−plane. A periodic array of SAQDs leads to a periodic auto-
correlation function. A nearly periodic array leads to a range-limited periodic auto-correlation function. The ensemble-
mean of these autocorrelation functions can be calculated, and it is a good predictor of a SAQD order.
3.1.1 Periodic array
Consider a perfectly periodic height fluctuation corresponding to a perfect lattice of SAQDs,
h(x) =
exp [ikn · (x− xO)] (27)
plus higher order harmonic, where the dots have a height proportional to h0, N is the degree of symmetry, probably,
4-fold or 6-fold, xO is a random origin offset.
kn = k0
2π(n− 1)
i + sin
2π(n− 1)
In a linear analysis, the higher order harmonics do not come into play, so they are neglected here. In reciprocal space,
e−ikn·xOδd(k− kn)
plus higher order harmonic. The autocorrelation function is found by plugging Eq. (27) into Eq. (25) and simplifying,
CA(∆x) =
)2 N∑
exp [ikn ·∆x] (28)
plus higher order harmonic. In finding Eq. (28), the relation∫
d2x′ ei(km−kn)·x
= Aδkmkn = (2π)
dδd(km − kn) (29)
has been used. δkk′ is the Kronecker Delta, and δd(k− k′) is the Dirac Delta. Eq. (29) will be helpful whenever it is
necessary to take an areal average or sum over Dirac Delta functions. In reciprocal space,
CAk =
(2π)2
m,n=1
δ2(k− km)δ2(k− kn)
δ2(k− ki) (30)
plus higher order harmonics, where δd(k − kn) = (A/(2π)d)δkkn .4. Thus, the order of the SAQD lattice manifests
itself as periodic functions in real-space (Eq. (28)) and sharp peaks in reciprocal space (Eq. (30)).
3.1.2 Nearly Periodic array
A nearly periodic arrays shows deviation from perfect order. This deviation is shows itself by broadening of the peaks
of the spectrum function, CAk , and by range limited periodicity of the real-space autocorrelation function, C
A(∆x).
These two measure of disorder are naturally related.
The disorder in lateral dot size ∆size and spacing, ∆spacing are related to each other and to the broadening of the
peaks in CAk (Fig. 6.a and c). Prior to ripening, the size and spatial order should be related, as the volume of a dot
should be proportional to the amount of nearby material. If the SAQDs have nearly uniform size and spacing (peak-
to-peak distance) L0, the reciprocal space autocorrelation function will be tightly clustered around the wavenumber
characterizing the dot spacing k0 = 2π/L0. There are a number of such peaks depending on the system symmetry
(Fig. 6.a and c), but consider just one. Since the order is not perfect, the peak will have a finite width. Consequently,
there will be a scatter in the dot size. Since L0 = 2π/k0, the scatter in dot spacing (∆spacing) is related to the scatter in
Fourier components (∆k). Taking the derivative of the spacing-wavenumber relation and rearranging,
∆spacing
It is reasonable to expect that the fractional disorder in size (∆size /Lsize) is given by a similar (if not exactly the same)
number.
Another way to view spatial order (periodicity) is not by dot-dot distances, but the distance over which the dot array
can be considered periodic. This limited periodicity is evident in the film height autocorrelation function (Eq. (25) and
Figs. 5.b, e and h). Consider two distant dots. Their position will be completely uncorrelated, so it will be completely
random as to whether one position corresponds to a peak or a valley. Thus, for a large differences in position the
autocorrelation function tends to zero.
CA(∆xlarge) = 0
Similarly, the mean-square fluctuation of the film height can be large so that
CA(∆x = 0)� 0.
The distance over which the autocorrelation function, CA(∆x) decays to 0 is the correlation length, Lcor. Thus, Lcor
is a reasonable measure of spatial order.
The two measures of order ∆spacing and Lcor are intrinsically linked. The well known rule of Fourier transforms
states that the product of the real-space and reciprocal space widths must be greater than or equal to unity. Thus,
∆kLcor ≥ 1, or ∆spacing ≥ 2πL20/Lcor. Similarly, one can expect that ∆size ∼ L2size/Lcor. Thus, assuming that dot size
is governed by the amount of nearby material, small dispersions in dot size are only possible if there is long correlation
length.
3.1.3 Ensemble Correlation Functions / ergodicity
SAQDs are seeded by random fluctuations. Consequently, each experiment or simulation must be treated as just
one possible realization, and the autocorrelation function will be different for each realization. Thus, for analytic
4Eq. (29) has been used to help with summation
predictions, one must rely on ensemble averages. In [24], it was assumed that the ensemble average correlation
function was a good description of a SAQD order, an assumption that was born out by numerical calculations. Now,
this relation is put on a more solid ground. In particular, it is found that the ensemble correlation functions provide
good estimates of the auto correlation function and spectrum function produced by any particular realization. First,
it is shown that the mean value of the film-height fluctuation is zero. Then the method to calculate the ensemble-
averaged autocorrelation function and spectrum function is presented. Additional mathematical details are presented
in Appendix 3.
3.1.3.1 Mean fluctuation It is fairly straightforward to show that the ensemble mean film-height fluctuation is
zero. The governing dynamics (Eq. (12)) is invariant upon the substitution h(x, t) → −h(x, t). Thus, assuming that
one does not bias the initial conditions the mean fluctuations must be zero for all time,
〈h(x, t)〉 = 〈−h(x, t)〉 = 0, and 〈hk(t)〉 = 0.
This is a common situation, and it is most appropriate to characterize the film height fluctuations using the two-point
correlation function (or simply “the correlation function”). [55]
3.1.3.2 Correlation Function The autocorrelation function can be estimated by its ensemble average. Further-
more, this ensemble average is equivalent to the correlation function that can be easily calculated analytically. These
relations are first discussed for the real-space correlation functions and then their Fourier transforms. First, the statis-
tical properties of the autocorrelation function are discussed. Then the statistical properties of the spectrum function.
Finally, the method to The main results are reported here, and details of derivations are reported in Appendix E.
First it is noted that the autocorrelation function averaged over all realizations is equal to the ensemble correlation
function. 〈
CA(∆x)
= C(∆x), where C(∆x) = 〈h(∆x)h(0)〉 , (31)
where 〈. . . 〉 indicate an ensemble average. Eq. (31) assumes that the model of film-growth is translationally invariant.5
This relationship is fortunate, in that it allows one to predict the “typical” autocorrelation function using analytic tools
that apply only to ensemble averages.
Second, it is noted that as the area that is used to calculate the autocorrelation function becomes large, the autocor-
relation function tends towards it mean value,
CA(∆x) ≈ C(∆x) +O[A−1/2], (32)
where O[A−1/2] indicates statistical fluctuations about the mean value that become smaller and smaller as the area in
an experiment or the simulation area in a numerical experiment becomes larger. These fluctuations or noise die out as
A1/2. For example, the autocorrelation functions in Figs. 5.e and h are very close to the ensemble average autocorre-
lation functions Figs. 5.f and h, but have random fluctuations that are most visible far from the origin. This property,
that averaging over a parameter such as position is equivalent to averaging over all realizations, is known as ergodicity.
Individual realizations are tightly distributed about a “typical” behavior. This tight distribution lends credibility to the
notion that one can have representative experiments or simulations. Unfortunately, the “demonstration” of Eq. (32) in
Appendix E is not as general as one might like. Rigorously, it applies when the Fourier components of film height (hk)
are independent and normally distributed; however, it is reasonable to conjecture that a relationship like Eq. (32) holds
whenever the statistical distribution of film heights is suitably bounded as the boundedness of CAk plays an important
role in the derivations.
In reciprocal space, one finds that the ensemble-mean spectrum function is〈
= Ck, (33)
where Ck is defined as the prefactor appearing in the reciprocal-space two-point correlation function.
Ckk = 〈hkh∗k〉 = Ckδ
d(k− k′) = Ck
(2π)d
δkk′ , (34)
5A quick survey of literature will find that, virtually all published continuum models of SAQD formation are translationally invariant.
where Eq. (29) has been used. This form for the two-point correlation function in reciprocal space occurs if and only
if the system is translationally invariant. Eq. (33) is valuable because one can solve for Ck analytically in the linear
case or using various analytic approximations in the non-linear case. Unlike the autocorrelation function, the spectrum
function fluctuates greatly about its mean. In fact, the fluctuations are about 100% (Appendix E.2). These large
fluctuations result in the commonly observed speckle pattern for the spectrum function CAk (Figs. 6.a and c). Contrast
this pattern with ensemble-mean spectrum function Ck shown in Figs. 6.b and d. These speckles can be removed by a
smoothing operation, and a relation similar to Eq. (32) results (Appendix E.2.2). Finally, it should be noted that just
as CAk is the Fourier transform of C
A(∆x), Ck is the Fourier transform of C(∆x) (Appendix E.1).
3.2 Stochastic Initial Conditions
To model or simulate the formation of SAQDs, it is absolutely essential to include some sort of stochastic effect. An
initially flat film h(x, 0) = 0 is in unstable equilibrium. Thus, to seed the formation of quantum dots, it is necessary
to perturb the flat surface. The simplest method to do this is to use stochastic initial conditions with deterministic
evolution. One can tenuously suppose that white noise initial conditions do not “bias” the ultimate evolution of the
film. [57] Thus, the initial conditions are taken from an ensemble with zero mean,
〈h(x, 0)〉 = 0. (35)
and a spatial correlation function,
C(x,x′, 0) = 〈h(x, 0)h(x′, 0)∗〉 = ∆2δd (x− x′) , (36)
where the brackets 〈. . . 〉 indicate an ensemble average, ∆ is the noise amplitude, and δd(x) is the d−dimensional
Dirac Delta function. White noise conditions have an infinite amplitude which is not physical. Thus, a minimum
modification can be made to “cut off” the infinite fluctuations.
C(x,x′, 0) =
(2πb20)
(x− x′)2
In the limit b0 → 0, this correlation function reverts to the white noise correlation functions.
In reciprocal space,
Ckk′(0) = 〈hk(0)h∗k′(0)〉
= (2π)−2d
ddx′ e(−ik·x+ik
′·x′)C(x,x′, 0)
(2π)d
δd(k− k′)
Letting b0 → 0, the white noise reciprocal space correlation function is obtained. Thus, the initial spectrum function
Ck(0) =
(2π)d
The atomic-scale has a small and short-lived influence on the final film morphology (Appendix F), but the cutoff
procedure is useful for choosing a reasonable value of ∆2. It seems reasonable to choose ∆2 so that the initial r.m.s.
fluctuation
C(0, 0) = 〈h(0, 0)h(0, 0)∗〉1/2is one monolayer (1 ML). Also, choosing b0 = 1 ML as the atomic scale
cutoff is
∆2 = (2π)d/2(1 ML)2+d, (38)
where the natural unit 1 ML is, of course, material dependent.
Using stochastic initial conditions, one can integrate individual initial conditions to obtain representative samples
and then average over many realizations, the Monte Carlo approach, or one can calculate analytically, the statistical
measures of the ensemble. The ensemble statistical measures are strongly related to the statistical measures of order
for an individual realization, so the second approach is opted for here. Thus, the predicted SAQD order is ultimately
stated in terms of ensemble correlation functions.
3.3 Reciprocal Space Correlation Functions
The reciprocal space correlation function, Ckk′ , and spectrum function, Ck, are calculated for the 1D and 2D isotropic
case and then for the 2D anisotropic case. Generally Ck includes the length scales introduced in Sec. 2.1.3 as well as
the atomic scale cutoff b0.
Ckk′ = 〈hk(t)h∗k′(t)〉 = e
(σk+σk′ )t 〈hk(0)hk′(0)∗〉
(2π)2
e(σk+σk′ )t−
δ2(k− k′). (39)
Without much error, b0 can be neglected in the exponential (Appendix F). Using Eq. (34), the spectrum function is
then identified as
(2π)d
e2σkt. (40)
Ck is now calculated for each model: 1D isotropic, 2D isotropic and 2D anisotropic.
3.3.1 one-dimensional
The one dimensional surface is the simplest, so it is treated first. The spectrum function is simply
e2σ0t−
2 (2σ2t)(k−k0)
Ck has a peak at k = ±k0i. One can easily read off the correlation length as
Lcor =
2σ2t = k
2(3α0 − 4β)(t/tc). (41)
so that
e2σ0t−
cor(k−k0)
This approximation is valid when k0Lcor � 1. In terms of kx,
e2σ0t
cor(kx−k0)
cor(kx+k0)
3.3.2 2D isotropic
The 2D isotropic case is very similar;
(2π)2
e2σ0t−
cor(k−k0)
, (42)
where Lcor is the same as in Eq. (41). It has maximum that forms a ring in the k−plane as graphed in Fig. 6.b.
3.3.3 anisotropic
The anisotropic spectrum function is
(2π)2
e2σ0t
‖(k‖−k0)
2− 12L
⊥ , (43)
where
2σ‖t = k
(6α0 − 8β)(t/tc), (44)
2σ⊥t = k
16�α0(t/tc), (45)
k‖ = cos[π(n− 1)/2]kx + sin[π(n− 1)/2]ky , and k⊥ = − sin[π(n− 1)/2]kx + cos[π(n− 1)/2]ky and it is graphed
in Fig. 6.d. This approximation is valid when k0L‖ � 1 and k0L⊥ � 1.
Figure 6: CAk and Ck for Ge/Si as discussed in Sec. 4. (a,b) 2D isotropic surface. Eq. (42) is used for Ck. (c,d) 2D
anisotropic surface. Eq. (43) is used for Ck.
Figure 7: 1D isotropic surface in real space for Ge/Si as discussed in Sec. 4. (a) Example of h(x) plotted over a length
of 8Lcor. (b) corresponding reals space correlation functions plotted for range ±4Lcor. Filled plot is an example of
CA(∆x). Solid line isC(∆x) (Eq. (46)).
Loosely speaking, one can argue that the isotropic case is similar to letting �A → 0 in Eq. (45) so that the
perpendicular correlation length is always 0 regardless of time. A more conservative approach would be to argue that
L⊥ ≈ 2π/k0 for the isotropic model via inspection of Figs. 6(a) and (b). Even still, the more conservative result
guarantees that the perpendicular correlation length will always be the same as the dot spacing; thus, it will always
limit SAQD order to the first nearest neighbor at best.
3.4 Real Space Correlation Functions
The real space correlation functionsC(∆x) are now calculated for the 1D and 2D isotropic cases and the 2D elastically
anisotropic case.
3.4.1 one-dimensional
In one dimension,
C(∆x) =
dkx e
ikx∆xCk
e2σ0t−
2/L2cor2 cos (k0x) . (46)
Thus, C(∆x) has a damped periodicity indicating that it is imperfectly periodic (Fig. 7).
3.4.2 2D isotropic
In two dimensions with elastic isotropy,
C(∆x) =
d2k eik·∆xCk
(2π)2
e2σ0t
dk kei(k∆x cos(θk−θ∆x)e−
cor(k−k0)
Performing the angular integration first,
C(∆x) =
e2σ0t
dk kJ0(k∆x)e
− 12L
cor(k−k0)
where J0 is the zeroth Bessel function. In general, this integral is best performed numerically; however, it can be
solved in two important cases: ∆x→ 0 and Lcor →∞ (corresponding to long times). In the first case,
C(∆x) =
e2σ0t
dk ke−
cor(k−k0)
Under the same conditions that Eq. (42) is valid (k0Lcor � 1), the lower limit of the integral can be approximated as
−∞ so that
C(∆x = 0) =
∆2k0√
2πLcor
e2σ0t. (47)
This function gives the mean square surface height fluctuation. In the second case where Lcor →∞, e−
cor(k−k0)
(2π)1/2L−1cor δ(k − k0), so that
C(∆x) =
∆2k0√
2πLcor
e2σ0tJ0 (k0∆x) . (48)
This correlation function is the most ordered case for a 2D isotropic surface. It is graphed in Fig. 5c.
3.4.3 anisotropic
To find the real-space correlation function for the elastically anisotropic case, it is best to find the contribution from
each peak and then sum so that
C(∆x) =
(2π)2
e2σ0t
Cn(∆x) (49)
where
Cn(∆x) =
d2k eik·xe−
‖(k‖−k0)
2− 12L
∆x can be decomposed into the directions parallel and perpendicular to kn, so that ∆x‖ = cos(π(n − 1)/2)∆x +
sin(π(n− 1)/2)∆y and ∆x⊥ = − sin(π(n− 1)/2)∆x+ cos(π(n− 1)/2)∆y. Thus,
Cn(∆x) =
dk‖ e
ik‖∆x‖− 12L
‖(k‖−k0)
dk⊥ e
ik⊥∆x⊥− 12L
2 (x2‖/L2‖+x2⊥/L2⊥)eik0x‖ .
Plugging into Eq. (49),
C(∆x) =
πL‖L⊥
e2σ0t
2 (x2/L2‖+y2/L2⊥) cos(k0x) + e
− 12 (x2/L2⊥+y2/L2‖) cos(k0y)
. (50)
3.5 Generalizability
The dynamics and analysis used here were for a specific model, but the general procedure for analyzing the order
resulting from a linearized model should hold for any model with well-separated peaks in the dispersion relation, σk.
The procedure to follow is:
1. Generate the dispersion relation, σk as some function of k.
2. Find the peaks in the dispersion relation, kn, (n = 1 . . . N )
3. Expand about the peaks to generate the peak values, σn, and local Hessian matrix,(
∂ki∂kj
The spectrum function is then approximately
Ck(t) ≈
(2π)2
e2σnt exp
t (k− kn) · H̃n · (k− kn)
. (51)
4. Find the Eigenvalues of the local Hessian matrix, (Hn)I and (Hn)II . They should be negative, if there is a
peak at kn
5. Use the eigenvalues to determine the correlation lengths, (Ln)I =
2 |(Hn)I | t and (Ln)II =
2 |(Hn)II | t.
The real-space correlation function is
C(∆x, t) ≈
(Hn)I (Hn)II
e2σnt exp
x · H̃−1n · x
eikn·x. (52)
The “goodness” of these approximate forms requires that (Ln)
I and (Ln)
II be much less than the spacing
between peaks in the correlation function so that the gaussians do not overlap greatly. A reasonable test for
this no-overlap condition is ‖kn‖ (Ln)I � 1 and ‖kn‖ (Ln)II � 1, assuming that the peaks are not large in
number or very closely spaced.
4 Order Predictions
The real-space correlation function formulas (Eqs. (46), (47), and (50)) and correlation length formulas (Eqs. (41), (44)
and (45)) can now be used to estimate the order of SAQDs. Ge on Si is chosen for this example because this system has
received the most attention from theoretical work [58, 38, 31, 18, 39, 41, 27, 25, 26, and others], and it is the simplest
since it involves the diffusion of a single species. The procedure described below tries to predict the amount of order
when an initial atomic-scale fluctuation becomes “large”. “Large” is taken to be greater than atomic-scale. Beyond
this point, one would expect non-linear terms to become important. An example is presented for Ge on Si at 600K to
compare and contrast the 2D anisotropic results with the 1D isotropic and 2D isotropic results. The predictions are
also compared with a linear numerical calculation on a discrete reciprocal-space grid to test the approximations made
and to illustrate the relation between the surface profile (h(x)), the example autocorrelation functions (CA(x) and
CAk ) and the ensemble correlation functions (C(∆x) and Ck). Figs. 6, 7 and 5 show these results. Finally, the relation
between average film height and order is investigated.
4.1 Ge at 600K
The formulations for the three discussed cases are implemented for Ge/Si at 600K. The correlation lengths are esti-
mated for the end of the linear regime where fluctuations become large (greater than atomic scale). First, appropriate
physical constants are used to give the corresponding correlation length and correlation functions vs. time. These in-
clude an initial average film height H̄ and a white noise amplitude ∆ (Eq. (38)). These initial conditions approximate
a film at the beginning of an anneal that immediately follows a rapid deposition. The time tlarge is found by solving for
the time where the mean-square fluctuations are atomic scale,
h(x, t)2
= C(∆x = 0) = 1 ML2. At this point, the
correlation lengths are calculated.
Physical constants for the 2D anisotropic calculation are taken as follows. The elastic constants for Ge at 600 K
are c11 = 1.199 × 1012, c12 = 4.01 × 1011(from cS = 3.991), c44 = 6.73. [51] Using aGe = 0.5658nm and
aSi = 0.5431nm, it is found that �m = 0.0418. Using the procedure from (Appendix C), M = 1.332× 1012dyn/cm2.
E0◦ = 4.96 × 109erg/cm3, and E45◦ = 4.35 × 109erg/cm3 , giving �A = 0.1236. The atomic volume is Ω =
2.27 × 10−23 cm3. The estimated surface energy density is γ = 1927 erg/cm2. The wetting potential is estimated
by picking a plausible critical surface height, Hc ≈ 4 ML = 1.132 nm and setting W (H) = E20◦H3c/(8γH) =
2.315 × 10−6/H erg/cm2. The resulting characteristic wave number is kc = 0.257 nm−1. The initial film height is
taken to be H̄ = Hc + 0.25 ML = 1.203 nm and then allowed to evolve naturally. Thus, β = 0.208, α0 = 0.5658,
k0 = 0.1456 nm−1, σ0 = 0.1192/tc, σ‖ = 0.864/(k2c tc), σ⊥ = 0.559/(k
c tc), L‖ = 0.744k
0 (t/tc)
1/2, and
L⊥ = 0.599k
0 (t/tc)
1/2. The unspecified diffusivity has been absorbed into the characteristic time tc. From Eq. (38),
∆2 = 0.0403 nm4, and Eq. (50) gives
C(0) =
1.223× 10−3tc/t
e0.02385t/tc nm2.
The initial infinitely rough surface undergoes a smoothing described by the tc/t factor. Then the surface roughens due
to the exponential. The initial divergent roughness is an artifact of the non-physical white noise with the atomic scale
cutoff b0 neglected (Appendix F). The time for the fluctuations to become “large” again are found by setting
C(0) = h2large (53)
where hlarge = 1 ML = 0.283 nm. The solutions are t1 = 0.01527tc or t2 = 430tc. The first solution is discarded
since it is due to the non-physical white noise. At tlarge = t2, L‖ = 105.8 nm, and L⊥ = 85.2 nm. Taking L⊥ as
more limiting, the correlation spans about n = k0L⊥/π = 3.95 islands across. The corresponding reciprocal space
(Eq. (43)) and real-space correlation function (Eq. (50)) are shown in Figs. 6.d and 5.f respectively.
A corresponding numerical experiment is performed. A periodic surface of size l = 96(2π/k0) is used. Random
initial conditions consistent with Eq. (38) are used for k−space points on a square grid bounded by kx, ky = ±2k0.
The relation between discrete and continuous Fourier components is used, (hk)discrete = [(2π)d/A]hk. Eqs. (13)
and (14) are used without any additional approximation to find hk at time t = tlarge. The resulting CAk , a portion of
the height profile h(x) and CA(∆x) are plotted in Figs. 6(c), 5(d) and 5(f) respectively.
Similar calculations can be performed for the one-dimensional and two-dimensional elastically isotropic cases.
Isotropic values used previously [58, 24] are about E = 1.361 × 1012 dyn/cm2 and ν = 0.198 giving M = E/(1 −
ν) = 1.697 × 1012 dyn/cm2 and E = 2M(1 + ν) = 7.10 × 109 erg/cm3. Using the same critical surface height,
Hc = 4 ML, W (H) = 4.74×10−6/H erg/cm2. The resulting characteristic wave number is kc = 0.368 nm−1. If the
film is grown to H̄ = Hc+0.25 ML = 1.203 nm and then allowed to evolve naturally, β = 0.208; thus, α0 = 0.5658,
k0 = 0.208 nm−1, σ0 = 0.1192/tc, σ2 = 0.864/(k2c tc), and Lcor = 0.744k
0 (t/tc)
1/2. In one dimension, Eq. (46)
is used to find the mean square height fluctuation. Using Eq. (38) with d = 1, ∆2 = 0.0568 nm3, and
C(0, t) = 0.01271(t/tc)
−1/2e0.0238t/tc .
Setting C(0, t) = (1 ML)2 = 0.0801 nm2, t1 = 0.0252tc, and t2 = 186.9tc. At t2, Lcor = 48.8 nm, and n =
k0Lcor/π = 3.24, so about 3 dots in a row should be well correlated. The corresponding numerical calculation of
size l = 96(2π/k0) is performed. A portion of h(x), CA(∆x) and C(∆x) are shown in Fig. 7. In two dimensions,
Eq. (47) is used to find
h(x, t)2
C(0, t) = 9.40× 10−4(t/tc)−1/2e0.0238t/tc .
Setting C(0, t) = 0.0801 nm2, t1 = 1.376 × 10−4tc, and t2 = 306tc. At t2, Lcor = 62.4 nm, and n = k0Lcor/π =
4.14, and correlation is expected to extend about 4 dots. However, it should be noted that this correlation is not
lattice-like. Corresponding numerical results and ensemble correlation functions are shown in Figs. 6 and 5.a-c.
4.2 General case of β
In [24] it was suggested that allowing the film to evolve with β close to the stability threshold could enhance the SAQD
correlation. It is interesting to note what happens for different values of β. Similar analytic and numerical calculations
are performed for the large film-height limit, β = 0, for the 2D anisotropic Ge/Si surface. For β = 0, tlarge = 40.3tc,
L⊥ = 30.0 nm, and n = k0L⊥/π = 1.84, so one to two dots in a row are expected to be well correlated. h(x)
and real-space correlation functions are shown in Figs. 5g-i. The range of order is significantly less than for the case
β = 0.208 (Sec. 4.1). For Si/Ge at 600K, the 2D anisotropic predictions for tlarge and L⊥ are shown in Fig. 8. In
general, the closer β is to the critical value 0.25, the longer the correlation length. One can manipulate equation (53)
to find that tlarge/tc varies approximately but not exactly as (β − 1/4)−1 × ln[h2large
�A/(∆2k2c )]. Consequently,
L⊥ ∼ (β − 1/4)−1/2. Furthermore, the appearance of hlarge and ∆2 inside the logarithm shows that the final order
estimates are not overly sensitive to the guesses for ∆2 and h2large. The divergence of L⊥ with β − 1/4 is initially
encouraging, but it is clear that for the parameters used for Ge/Si, subatomic control of the film height is needed to
yield significantly enhanced long range correlations. Also as one approaches this threshold, one can probably expect
thermal activation to nucleate subcritical SAQDs whose effect on supercritically formed SAQDs is uncertain. There
should be some interesting phenomena at the theH → Hc.
l arge c
l a r g e
k 0L �
� 0 �
l a r g e � �
Figure 8: tlarge and L⊥ vs. β for Si/Ge using the 2D anisotropic model as described in Sec. 4. Units are normalize to
characteristic time tc and predicted number of correlated dots (n = k0L⊥/π).
5 Discussion/Conclusions
The order of epitaxial self-assembled quantum dots during initial stages of growth has been studied using a common
model of surface diffusion with stochastic initial conditions. It has been shown that correlation functions of small
surface height fluctuations can be predicted analytically using corresponding ensemble average correlation functions.
These correlation functions are characterized by correlation lengths that can be predicted by analytic formulas given
certain reasonable assumptions about the diffusion potential and the height and lateral scale of initial atomic scale
random fluctuations. Thus, the linear model of film surface height evolution via surface diffusion has enabled analytic
predictions of epitaxial SAQD order that are valid for small film height fluctuations. To what extent the initial degree
of order persists into later stages of growth remains to be studied, but the order of initial stages should certainly have
a strong influence on final outcomes. Furthermore, the linear analysis should provide insight into the less tractable
non-linear behavior. These predictions of SAQD order have been used to investigate the role of crystal anisotropy and
initial film height.
Crystal anisotropy has been shown to play an important role in enhancing SAQD order as observed in previous
numerical simulations continuum and atomistic numerical simulations. [43, 37, 44, 45] If a four-fold symmetry is
assumed for the governing dynamics, the effect of crystal anisotropy to linear order is felt through elastic anisotropy
alone. It is shown that elastic anisotropy is required to produce a lattice-like structure of SAQDs. The enhanced spatial
order should in turn lead to enhanced size order, a consequence that must be confirmed with non-linear studies, but
appears to be true based on the present available literature.
The role of initial film height has been shown to greatly influence order. Growth near the critical film height for
dot formation can enhance order. This order enhancement comes from increasing the duration of the linear small-
fluctuation stage of growth. In fact, the predicted correlation lengths diverge when the initial film height approaches
the critical film height from above. Achieving large correlation lengths in this manor is of course practically limited
by ability to control film heights to subatomic accuracy. Additionally, one should be careful when interpreting the
continuum model in such a context, as the effect of atomic discreteness might be greater at the transition film height.
Finally, it is likely that additional randomizing effects of thermal activation will effectively cut off this divergence
when the critical film height is approached from below during deposition.
Finally, the presented method may be useful as a first step in the analysis of methods to enhance SAQD order. It is
reasonable to suppose that under some circumstances initial growth stages will be very important while for others they
will not. For example, prior work on vertical stacking appears to confirm the presented ordering mechanism. [44].
Vertical stacking not only achieves vertical correlation of dots, but each layer is more ordered horizontally than the
one below. Additionally, a “growth window” was found, whereby to achieve enhanced order, the evolution of each
layer be terminated before ripening begins. The reported simulation [44] supports the following scenario for SAQD
order development. Order is enhanced during the small fluctuation stage as described here. Once the fluctuations are
sufficiently large, the seeded dots evolve towards their equilibrium shapes. Finally, the dots begin to ripen and order
diminishes. Order is transfered via strain to the next layer so that the next layer gets a head start on its initial ordering.
Thus, the multiple layers of dots effectively draws out the linear growth stage. It may be possible to modify the present
model to predict the correlation length of each SAQD layer.
A Diffusion Potential
The diffusion potential is calculated in terms of the film height H that is a function of the in plane coordinates x =
xi + yj. The elastic and surface energy portions of the diffusion potential can be found in [15]
µelast(x) = Ωω(x), and µsurf = −Ωγκ(x),
where Ω is the atomic volume, ω(x) is the elastic energy density at the film surface, γ is the surface energy density,
and κ is the total surface curvature. However, other calculations need to be included:
1. µwet for the two wetting potential cases, Eq. (3) and (5),
2. and µsurf and µwet when the surface energy density γ and wetting energy density W also depend on surface
orientation.
Before these case are addressed, a general form for the diffusion potential is justified.
A.1 General Form µ = ΩδF/δH(x)
The diffusion potential, µ(x), is the change in free energy, F , when a particle is added at a position, x. Note that µ(x)
and F are relative energies. They can be used to compare the binding energy of one site on the surface in comparison
with another site, but should not be interpreted as an absolute binding energy or total formation energy of the surface.
If a particle has a volume Ω, then the diffusion potential at x is related to the variation of free energy with volume,
δF = Ω−1
ddxµ(x)δV (x), (54)
where δV (x) is the volume variation at x. Calculating δV (x), V =
ddxH(x).Therefore, δV (x) = δH(x).
Substituting into δF (Eq. (54)), δF = Ω−1
ddxµ(x)δH(x) or µ(x) = ΩδF/δH(x).
A.2 Simple Model
Starting from Eq. (2), µ(x) is found by taking the variational derivative,
µelast.(x) = Ω
δH(x)
volume
ddxdz ω[H](x, z) = Ωω (x)
where the “[H]” indicates that the elastic energy, ω, is a nonlocal functional of the film height H, and ω(x) =
ω[H] (x,H(x)), the elastic energy density evaluated above lateral position x at the free surface (z = H(x)). See [15]
for details of the derivation. The surface energy diffusion potential is
µsurf.(x,t) = Ω
δH(x)
1 + (∇H(x))2
= −Ω∇ ·
1 + (∇H(x))2
γ = −Ωγκ(x).
The wetting energy diffusion potential is
µwet(x) = Ω
δH(x)
ddxW (H(x))
= ΩW ′(H(x))
Putting these three terms together, one obtains Eq. (3)
A.3 General Model
Consider the general form for the combined surface energy and wetting potential,
Fsw =
ddxFsw(H(x),∇H(x))
as in Eq. (4) so that the free energy is an integral over the x−plane of an energy density that depends on H(x) and
∇H(x) locally. The corresponding diffusion potential is
µ(x) = Ω
δH(x)
F (10)sw (H(x),∇H(x))−∇ · F
sw (H(x),∇H(x))
B Linearized Diffusion Potential and Anisotropy
The linearized diffusion potential µlin, k is found by finding µ(x) to first order in height fluctuations (h), to get µlin(x)
and then taking the Fourier transform to get µlin,k. The linearization of the simple isotropic diffusion potential corre-
sponding to Eqs. (2) and (3) was discussed in Sec. 2.1.1.1. Here, the more general diffusion potential corresponding
to Eqs (4) and (5) is linearized and then applied to the anisotropic simple model and the anisotropic general model.
Only the surface and wetting parts of the diffusion potential are discussed in this appendix. See ref. [15], Sec. 2.2.1.1
and Appendix C for discussion of µelast..
B.1 Linearizing the simple model
Consider a wetting potential and diffusion potential that both depend on the film height gradient ∇H, γ → γ(∇H)
and W (H) → W (H,∇H). Starting from Eq. (6) and expanding to second order in the film height fluctuation using
H(x) = H̄+ h(x) (Eq. (7)),
1 + (∇H)2
]−1/2
γ(∇H) =
(∇h)2 + . . .
γ + γ′ ·∇h+ γ̃′′ : ∇h∇h+ . . .
= γ + γ′ ·∇h−
γ (∇h)2 + γ̃′′ : ∇h∇h+O[h3]
where γ is γ(0), and the primes indicate the derivatives with respect to the surface height gradient.
γ′ = ∂∇Hγ(∇H)|∇H=0 , and γ̃
′′ = ∂∇H∂∇Hγ(∇H)|∇H=0 .
Taking the derivative with respect to ∇h results in a tensor of rank equal to the order of the derivative because ∇h is
a vector (tank 1 tensor). Taking the variational derivative, µsurf.(x) = ΩδFsurf./δh(x),
µsurf., lin(x) = Ω
γ∇2h(x)− γ̃′′ : ∇∇h(x)
The term with γ′ vanishes because it is the divergence of a constant (∇ · γ′). Taking the inverse Fourier transform,
µsurf., lin,k = Ω
−γk2 + k · γ̃′′ · k
hk. (55)
The first term is isotropic. The second term is parameterized by a rank 2 symmetric tensor.
Going through the same process, one finds essentially the same result for an orientation dependent wetting energy.
The step details are so close to the details for linearizing the more general form, Fsw(H,∇H), they are deferred to
(Appendix B.2). One finds that
µwet,lin,k = Ω
W (20) + k · W̃(02) · k
. (56)
where W (mn) = ∂mH∂
∇HW (H,∇H)|H=H̄,∇H=0 is the m
th and nth derivative of the wetting energy density with
respect toH and ∇H evaluated for a perfectly flat film of height H̄. W (mn) is a tensor of rank n.
B.1.1 isotropic case
In the isotropic case, γ̃′′ → γ′′Ĩ, where Ĩ is the identity operator, and γ′′ is a scalar. Similarly, W̃(02) →W (02)Ĩ. One
thus gets for the combined surface and wetting parts of the diffusion potential,
µsw,lin,k = Ω
−γ + γ′′ +W (02)
k2 +W (20)
Thus, in the isotropic case, the linear order effect of introducing a surface orientation to either the surface energy or
the wetting potential is simply to change the apparent surface energy density by γ → γ − γ′′ −W (02).
B.1.2 anisotropic case
The surface and wetting parts of the diffusion potential (Eqs. (55) and (56)) can admit only a limited anisotropy.
They both contain rank 2 symmetric tensors, γ̃′′ and W̃(02) in the x−plane. For a two-dimensional surface, this
means that they can either have two-fold-symmetric (rotations by 180◦) anisotropy or none at all. Thus, for the case
considered in Sec. 2.2.1.2, four-fold-symmetric anisotropy , the surface and wetting parts of the diffusion potential
must be completely isotropic. As discussed in Sec. 2.2.1.2, the (100) surface of zinc-blend structures, such as the
mentioned Ge, Si, InAs and GaAs present a rather complicated situation. For simplicity, it is assumed here that the
surface and wetting energies are at least four-fold symmetric. Consequently, they are completely isotropic.
Finally, it should be noted that if Fsw depends on higher order derivatives, then the discussion is greatly compli-
cated and a larger class of anisotropic terms is admissible. For example, whenFsw → Fsw(H,∇H,∇∇H,∇∇∇H, . . . )
is expanded aboutH(x) = H̄ to quadratic order in h, it would contains tensors of rank 6 and maybe even higher.
B.2 Linearizing the general model
The elastic part of the linearized diffusion potential was discussed in Sec. 2.2.1.1 and Appendix C . Eq. (56) can be
found by using all of the following steps with the substitution Fsw →W . The surface-wetting part of the diffusion po-
tential µ(x) is found by expanding Fsw to second order in the film-height fluctuation, h, and then taking the variational
derivative. Expanding Fsw about h = 0 and ∇h = 0,
Fsw(H̄+ h,∇h) = F (00)sw + F
sw h+ F
sw ·∇h+ hF
sw ·∇h . . .
· · ·+
F (20)sw h
F̃(02)sw : ∇h∇h+O[h
Note that in this expansion, all the F (mn)sw terms are constant with respect to h and depend implicitly on the average
film height, H̄. The first index indicates the mth derivative with respect to h. The second index indicates the nth
derivative with respect to ∇h. The derivatives are evaluated for a perfectly flat surface of height H̄. Thus,
F (mn)sw = ∂
∇HFsw (H,∇H)|H=H̄,∇H=0 .
Since ∇h is a vector in the x−plane, F (mn)sw is a tensor of rank n. Taking the variational derivative of Fsw =∫
ddxFsw(H,∇H) and keeping terms to order h1,
δh(x)
= F (10)sw −∇ · F
sw + F
sw h−∇ ·
F̃(02)sw ·∇h
Note that the F (00)sw term vanishes because it is constant, and the F
sw term vanishes upon simplification. Additionally,
the F (10)sw can be neglected if one enforces the condition that the film-height fluctuations do not add or subtract material
from the surface, namely that
ddx δh(x, t) = 0. Alternatively, one can discard it in anticipation of taking the gradient
of the diffusion potential, since it is a constant. The term ∇ · F(01)sw = 0 for the same reasons, or because F
sw is a
constant. Multiplying through by the atomic volume,
µlin(x) = Ω
F (20)sw h− F̃
sw : ∇∇h
. (57)
B.2.1 isotropic case
In the isotropic case, F̃(02)sw must be proportional to the identity so that F̃
sw = F
sw Ĩ; thus,
µsw,lin(x) = Ω
F (20)sw h(x)− F
2h(x)
Taking the inverse Fourier transform of this equation,
µsw,lin,k = Ω
F (20)sw + F
This gives case b in Eq. (9).
B.2.2 anisotropic case
If the surface is anisotropic, then F̃(02)sw in Eq. (57) is a rank 2 symmetric tensor in the x−plane. Thus, it can have two
distinct eigenvalues, and automatically has 2-fold rotational symmetry (rotations by 180◦). If any other symmetry is
assumed such as 4-fold symmetry (rotations by 90◦), then F̃(02)sw must be fully isotropic. Taking the inverse Fourier
transform,
µsw,lin,k = Ω
F (20)sw + k · F̃
sw · k
In Eq. (23), case b, it is assumed that there is four-fold symmetry, resulting in a surface-wetting part of the diffusion
potential that is completely isotropic.
C Elastic Anisotropy
In principal, the anisotropic elastic energy ωk is found in the same fashion as the isotropic elastic energy. [15] The
flat film, initially in a state of biaxial stress, is perturbed by a small periodic surface fluctuation of amplitude h0. An
appropriate elastic field is added to satisfy the perturbed traction-free boundary condition at the free surface. Finally,
the elastic energy is evaluated at the free surface to first order in h0. The coefficient h0 is the sought after ωk. The
equations themselves are cumbersome and best solved using a numeric implementation, so an abstract procedure for
calculating ωk is outlined here. ωk is found for k = 1 but arbitrary θk.
Let the surface have a height variation
h(x) = h0e
To first order in h0, the surface normal is
n(x) = −ikh0eikxi + k.
The elastic energy needs to be calculated to first order in h0. To find the elastic energy, it is necessary to find the
perturbing elastic field to first order in h0.
The initial unperturbed stress state is
σ̃m =
 σm 0 00 σm 0
0 0 0
where σm =
c11 + c12 − 2c212/c11
�m. Note that this stress state is isotropic in the x−y-plane and thus independent
of rotations about the vertical axis. Under this stress state, a flat surface is traction-free. With the height perturbation,
the traction is
tj = (n · σ̃m)j = −ikh0M�mδj1e
ikx. (58)
Next to find the perturbing elastic fields. These are not isotropic in the x − y−plane, and it is necessary to take
into account the angle. First, the 3× 3× 3× 3 elastic stiffness tensor cijkl is constructed for the cube orientation from
the compact 9× 9 matrix cij . The tensor representation aids in rotation. The stiffness tensor is then passively rotated
in the x− y−plane by an angel θk,
cijkl(θk) =
m,n,p,q=1
R(θk)imR(θk)jnR(θk)kpR(θk)lqcmnpq
where
R(θk) =
 cos(θk) sin(θk) 0− sin(θk) cos(θk) 0
0 0 1
This passive rotation of cijkl is equivalent to actively rotating the wave vector k = ki by θk.
The appropriate form for the perturbing displacement field is found. Assume a displacement of the form
ui(x, y, z) = Uie
k(ix+κz),
where κ can have a complex value. The elastic equilibrium equations are
i,k,l=1
cijkl(θk)
ul = 0; j = 1 . . . 3.
Cjl(θk, κ)Ul
k2ek(ix+κz) = 0 (59)
where
Cjl(θk, κ) =
i,k=1
cijkl(θk)(iδi1 + δi3κ)(iδk1 + δk3κ).
Factoring out k2ek(ix+κz), the part in parenthesis must be identically zero.
To obtain a non-trivial solution, the determinant of Cjl(θk, κ) to zero. Six complex values of κ are found. The
values of κ with Re[κ] < 0 are discarded since the corresponding displacements blow up as z → −∞. Each of
the remaining values κ = κp with p = 1 . . . 3 is substituted back into Cjl(θk, κ), and Eq. (59) is solved to find the
corresponding eigenvectors, Upl . The total displacement is thus
ul(x, y, z) = i�mh0
k(ix+κpz),
where it is assumed that the perturbing elastic displacement field is proportional to h0and σm, and the factor of i is
put in for convenience. The coefficients Ap can be found from the traction-free boundary condition at the free surface.
The traction formula is
i,k,l=1
nicijkl(θk)
ul(x, y, z) = ik�mh0
i,k,l,p=1
nicijkl(θk)ApU
l (iδk1 + κ
pδk3)e
k(ix+κpz) (60)
The traction is already proportional to h0. Thus, all terms in the sum must be kept to zeroth order in h0 so that
h(x) = 0, and n(x) = k.
Thus, plugging z = 0 to Eq. (60),
tj = ik�mh0
(ic3j1l(θk) + κ
pc3j3l(θk))ApU
ikx. (61)
Since the total traction (Eqs. (58) and (61)) must be zero, the coefficients Ap are found from
KjpAp = Rj ,
where
Kjp =
(ic3j1l(θk) + κ
pc3j3l(θk))U
Rj = Mδj1
for j = 1 . . . 3. It is worth noting that only for the symmetry directions, θk = 0◦ and θk = 45◦ is the strain purely
plane-strain as it is for the elastically isotropic case.
The elastic energy at the film surface is found to order O(h0). If the stress and strain are expanded to first order in
h0, σ̃ = σ̃0 + σ̃1, and �̃ = �̃0 + �̃1, then
�̃ : c̃ : �̃ =
σ̃0 : �0 + σ̃0 : �̃1 +O(h
Thus,
U = U0 +M�m ((�1)11 + (�1)22)
(�1)11 =
= −�mkh0
(�1)22 = ∂u2/∂y = 0. Thus,
U = U0 − Eθkkh0e
where
Eθk = M�
where Apand U
1 are implicitly functions of θk. This procedure has been used to find the values of E0◦ and E45◦ for
Table. 2 and Sec. 4.
D Diffusional Anisotropy
In general, the surface diffusivity can depend on the film height H(x) and the surface orientation ∇H(x) so that the
surface current is
JS(x) = D̃(H(x),∇H(x)) ·∇sµ(x)
where ∇s is the surface gradient, and D̃ is a rank 2 tensor in the two-dimensional space tangent to the film surface at
x. Linearizing the surface current about a flat surface,
JS(x) = D̃(H̄) ·∇µlin(x)
where the diffusivity must be evaluated for h = 0 and ∇h = 0, since µlin(x) is already proportional to h(x). The lin-
earized diffusivity is a symmetric rank 2 tensor in the x−plane. Thus, it is similar to F̃sw discussed in Appendix B.2.2.
It is automatically either two-fold symmetry (rotations by 180◦) or it is completely isotropic. In Eq. (23), four-fold
symmetry of the surface is assumed. Thus, the diffusivity must be completely isotropic; D̃ → D, a scalar. Sec-
tion 2.2.1.2 and Appendix B.2.2 contain discussions of the symmetry properties of the various rank 2 tensors that
appear in the linear evolution equations. A limited case of diffusional anisotropy has been modeled via kinetic Monte
Carlo technique. [54]
E Correlation Functions
E.1 Mean Values
Equations (31) and (33) are central to the presented analysis. Here, they are derived. The two-point correlation func-
tions for a stochastic system are introduced. Then, the average of the autocorrelation function is taken and expressed
in terms of the two-point correlation functions. Finally, this average is simplified using the translational invariance of
the system (governing equations and ensemble of initial conditions).
The two-point real-space space correlation function is
C(x,x′) = 〈h(x)h(x′)∗〉 ,
and the reciprocal space correlation function is
Ckk′ = 〈hkh∗k′〉 .
These are related by the double Fourier transform,
Ckk′ = 1(2π)2d
ddxddx′ e−ik·x+ik
′·x′C(x,x′); (62)
C(x,x′) =
ddkddk′ eik·x−ik
′·x′Ckk′ . (63)
These ensemble correlation functions can be used to give the ensemble-mean autocorrelation function and spec-
trum function. In real space,
CA(∆x)
d2x′ 〈h(∆x + x′)h(x′)〉
d2x′ C(∆x + x′,x′). (64)
(2π)d
〈hkh∗k〉 =
(2π)d
Ckk. (65)
Fortunately, the translational invariance of the system simplifies these relations. Inspecting the governing equations
and invoking the translational invariance of the stochastic initial conditions, the resulting ensemble and its statistical
measures must also be translationally invariant. Thus under the translation by x′,
C(∆x + x′,x′) = C(∆x,0) = C(∆x), (66)
so that the independent variable is reduced to just the difference vector ∆x = x − x′. This relation can be used to
simplify both the real and reciprocal space relations.
The real space relation simplifies as follows.Inserting Eq. (66) into Eq. (64),
CA(∆x)
d2x′ C(∆x,0) = C(∆x). (67)
The reciprocal space relation (Eq. (62)) simplifies to
Ckk′ = Ckδ
2(k− k′) = Ck
(2π)d
δkk′ , (68)
where
(2π)d
d2∆x e−ik·∆xC(∆x).
One can see immediately from Eq. (67) that Ck is the Fourier transform of
CA(∆x)
= C(∆x), or one can plug
Eq. (68) into Eq. (65), to get
= Ck.
E.2 Variance and Convergence
The ergodic hypothesis is that an average with respect to a parameter such as position or time tends towards an
ensemble average. In this case,
CAk ≈
= Ck, (69)
and CA(∆x) ≈
CA(∆x)
= C(∆x).
when the surface area is very large. The ensemble average is a good substitute if the variance about the average
vanishes as the substrate area A becomes large. It is found that in reciprocal space,
Var(CAk ) =
= C2k. (70)
Thus, the ergodic hypothesis does not hold for CAk . In practice, C
k is a speckled version of Ck (Fig. 6) However, if
one smooths CAk by averaging over a small patch in reciprocal space of size ksmooth = 1/∆s, so that
CAk (∆s) =
)d/2 ∫
ddk′ e−
′−k)2CAk′ , (71)
then Var
CAk (∆s)
diminishes as 1/A. For sufficiently large ∆s,〈
CAk (∆s)
≈ Ck, (72)
CAk (∆s)
πd/2∆ds
C2k. (73)
Thus, the ergodic hypothesis (Eq. (69)) only holds for a smoothed version of CAk .
In real space,
CA(∆x)
CA(∆x)
CA(∆x)
(2π)d
e2ik·∆xC2k + C
, (74)
where the integral is bounded (finite) provided that either t > 0 or the atomic scale cutoff b0 > 0. Thus, the ergodic
hypothesis holds for the real space autocorrelation function.
E.2.1 Eq. (70)
First,
CAk C
is calculated.
CAk C
(2π)d
〈hkh∗khk′h
k′〉 .
Assume that he distribution of hk is gaussian. Also, assume that h(x) is real so that hkh−k = |hk|
2. Then,〈
= Ck1Ck2δ
d(k1 − k4)δd(k2 − k3) . . .
. . . +Ck1Ck2δ
d(k1 + k3)δ
d(k2 + k4) . . .
. . . +Ck1Ck3δ
d(k1 − k2)δd(k3 − k4).
Thus,
CAk C
(2π)d
δd(k− k′)
. . .
. . . +C2k
δd(k + k′)
+ CkCk′
δd(0)
. (75)
= C2k
δkk′ + δk(−k′)
+ CkCk′ , (76)
where Eq. (29) has been used liberally. Setting k = k′, results in Eq. (70).
E.2.2 Eq. (73)
Now consider CAk smoothed over a length ∆s (Eq. (71)). The mean value is
CAk (∆s)
)d/2 ∫
ddk′ e−
′−k)2 〈CAk′〉 .
)d/2 ∫
ddk′ e−
′−k)2Ck′ .
For sufficiently small ksmooth, (sufficiently large ∆s), Eq. (72) results.
The variance of CAk (∆s) is now calculated. First, it is necessary to calculate
CAk (∆s)
CAk (∆s)
ddk′ e−
′−k)2 . . .
. . . ×
ddk′′ e−
′′−k)2 〈CAk′CAk′′〉 .
Using Eq. (75) and Eq. (29) as needed,
CAk (∆s)
ddk′ddk′′ e−
′−k)2e−
′′−k)2
(2π)d
. . .
· · · ×
δd(k′ − k′′)
+ C2k
δd(k′ + k′′)
+ CkCk′
δd(0)
′−k)2C2k′ + e
− 12 ∆
s[(k′−k)2+(k′+k)2]C2k′
. . .
· · ·+
)d/2 ∫
ddk′ e−
′−k)2Ck′
The first integral is bounded (finite) because Ck is bounded. Let its finite value be denoted I . The second integral is
simply
CAk (∆s)
.Thus,
Var(CAk (∆s)) =
∆2ds I
a finite value that decreases as A−1 as required for the ergodic hypothesis to hold. For sufficiently small ksmooth (large
∆s), I ≈ (π/∆2s)d/2C2k, and Eq. (73) results. It should also be noted that the large ∆s required for this approximation
also creates a more stringent requirement that A be large.
E.2.3 Eq. (74)
Now, consider the real space auto-correlation function. First,
CA(∆x)CA(∆x)
is needed.
CA(∆x)CA(∆x)
ddkddk′ ei(k+k
′)·∆x 〈CAk CAk′〉
Proceeding in a fashion similar to the previous section (making use of Eqs. (75) and (29) as needed) ,
CA(∆x)CA(∆x)
(2π)2d
ddkddk′ ei(k+k
′)·∆x
δd(k− k′)
. . .
· · ·+ C2k
δd(k + k′)
+ CkCk′
δd(0)
(2π)d
e2ik·∆xC2k + C
. . .
· · ·+
ddk eik·∆xCk
ddk′ eik
′·∆xCk
(2π)d
e2ik·∆xC2k + C
CA(∆x)
Thus, Eq. (74) results. For the variance to be vanishing, the integral in Eq. (74) must be bounded (finite). If time,
t > 0, the exponential in Eq. (77) guarantees that the integral is bounded. For time t = 0, the integral is only bounded
if the atomic scale cutoff b0 > 0.
F Atomic Scale Cutoff
Starting from Eq. (39),
(2π)d
e2σkt−
. (77)
The effect of the small scale cutoff is both small and short-lived, as it only works to suppress fluctuations with large
wavenumbers. The most important fluctuations have wavenumbers between 0 and 2kc. Thus, the typical size of the
cutoff term is about b20k
c . If a typical dot size or spacing size 10 nm, and a typical atomic scale is 10
−1 nm, a typical
value for this term is about 10−3 − 10−2. To calculate the effect of the cutoff, it can absorbed into the time-dependent
part with the substitution
so that its effect lasts only as long as a perturbation with atomic scale curvature (κ = b0). Thus, Eq. (40) is a good
approximation.
Acknowledgement
Thanks to L. Fang and C. Kumar for useful comments during the writing of this article.
References
[1] D. Bimberg, M. Grnudmann, and N. N. Ledentsov. Quantum Dot Heterostructures. John Wiley & Sons, 1999.
[2] O. P. Pchelyakov, Yu. B. Bolkhovityanov, A. V. Dvurechenski, L. V. Sokolov, A. I. Nikiforov, A. I. Yakimov,
and B. Voigtländer. SiliconGermanium nanostructures with quantum dots: Formation mechanisms and electrical
properties. Semiconductors, 34(11):122947, 2000. [doi:10.1134/1.1325416].
[3] M. Grundmann. The present status of quantum dot lasers. Physica E, 5:167, 2000. [doi:10.1016/S1386-
9477(99)00041-7].
[4] Pierre M. Petroff, Axel Lorke, and Atac Imamoglu. Epitaxially self-assembled quantum dots. Physics Today,
pages 46–52, May 2001.
[5] Hui-Yun Liu, Bo Xu, Yong-Qiang Wei, Ding Ding, Jia-Jun Qian, Qin Han, Ji-Ben Liang, and Zhan-Guo Wang.
High-power and long-lifetime InAs/GaAs quantum-dot laser at 1080 nm. Applied Physics Letters, 79(18):2868–
70, 2001. [doi:10.1063/1.1415416].
[6] F. Heinrichsdorff, M.H. Mao, N. Kirstaedter, A. Krost, D. Bimberg, A. O. Kosogov, and P. Werner. Room-
temperature continuous-wave lasing from stacked InAs/GaAs quantum dots grown by metalorganic chemical
vapor deposition. Applied Physics Letters, 71(1):22–4, 1997. [doi:doi:10.1063/1.120556].
[7] D. Bimberg, N.N. Ledentsov, and J.A. Lott. Quantum-dot vertical-cavity surface-emitting laser. MRS Bulletin,
27(7):531–7, 2002.
[8] N. N. Ledentsov. Long-wavelength quantum-dot lasers on GaAs substrates: From media to device con-
cepts. IEEE Journal of Selected Topics in Quantum Electronics, 8(5):1015–23, September/October 2002.
[doi:10.1109/JSTQE.2002.804236].
[9] M Friesen, P Rugheimer, D. E. Savage, M. G. Lagally, D. W. van der Weide, R Joynt, and M. A. Eriksson.
Practical design and simulation of silicon-based quantum-dot qubits. Physical Review B, 67(12):121301 (R),
2003. [doi:10.1103/PhysRevB.67.121301].
[10] Yi-Chang Cheng, San-Te (Cing-Ming) Yang, Jyh-Neng Yang, Liann-Be Chang, and Li-Zen Hsieh. Fabrication of
a far-infrared photodetector based on InAs/GaAs quantum-dot superlattices. Optical Engineering, 42(1):11923,
2003. [doi:doi:10.1117/1.1525277].
[11] R. Krebs, S. Deubert, J.P. Reithmaier, and A. Forchel. Improved performance of MBE grown quantum-dot
lasers with asymmetricdots in a well design emitting near 1.3 µm. Journal of Crystal Growth, 251:7427, 2003.
[doi:10.1016/S0022-0248(02)02385-0].
[12] Hiroyuki Sakaki. Progress and prospects of advanced quantum nanostructures and roles of molecular beam
epitaxy. Journal of Crystal Growth, 251:9–16, 2003. [doi:10.1016/S0022-0248(03)00831-5].
[13] B. J. Spencer, P. W. Voorhees, and S. H. Davis. Morphological instability in epitaxially strained dislocation-free
films. Physical Review Letters, 67(26):3696–3699, 1991. [doi:10.1103/PhysRevLett.67.3696].
[14] Karl Brunner. Si/ge nanostructures. Reports on Progress in Physics, 65(1):27–72, 2002. [doi:10.1088/0034-
4885/65/1/202].
[15] L. B. Freund and S. Suresh. Thin Film Materials: Stress, Defect Formation and Surface Evolution, chapter 8.
Cambridge University Press, 2003.
[16] S. Yu Shiryaev, E. Verstlund Pedersen, F. Jensen, J. Wulff Petersen, J. Lundsgaard Hansen, and A. Nylandsted
Larson. Dislocation patterning - a new tool for spatial manipulation of Ge islands. Thin solid films, 294(1-
2):311–314, 1997. [doi: 10.1016/S0040-6090(96)09240-1].
[17] C. Kumar and L. H. Friedman. Simulation of thermal field directed self assembly of epitaxial self-assembled Ge
quantum dots. Journal of Applied Physics, in press.
[18] Lawrence H. Friedman and Jian Xu. Feasibility study for thermal-field directed self-assembly of heteroepitaxial
quantum dots. Applied Physics Letters, 88:093105, 2006. [doi:10.1063/1.2179109].
[19] S. Krishna, D. Zhu, J. Xu, and P. Bhattacharya. Structural and luminescence characteristics of cycled sub-
monolayer InAs/GaAs quantum dots with room-temperature emission at 1.3 µm. Journal of Applied Physics,
86:6135–8, 1999. [doi:10.1063/1.371664].
[20] R. Hull, J.L. Gray, M. Kammler, T. Vandervelde, T. Kobayashi, P. Kumar, T. Pernell, J.C. Bean, J.A. Floro,
and F.M. Ross. Precision placement of heteroepitaxial semiconductor quantum dots. Materials Science and
Engineering B, 101:1–8, 2003. [doi:10.1016/S0921-5107(02)00680-3].
[21] O. Guise, Jr. J. T. Yates, J. Levy, J. Ahner, V. Vaithyanathan, and D. G. Schlom. Patterning of sub-10nm ge
islands on si(100) by direct self-assembly. Applied Physics Letters, 87:171902, 2005. [doi:10.1063/1.2112198].
[22] X. Niu, R. Vardavas, R. E. Caflisch, and C. Ratsch. Level set simulation of directed self-assembly during epitaxial
growth. Physical Review B, 74(19):193403, Nov 2006. [doi:10.1103/PhysRevB.74.193403.
[23] Z. M. Zhao, T. S. Yoon, W. Feng, B. Y. Li, J. H. Kim, J. Liu, O. Hulko, Y. H. Xie, H. M. Kim, K. B.
Kim, H. J. Kim, K. L. Wang, C. Ratsch, R. Caflisch, D. Y. Ryu, and T. P. Russell. The challenges in
guided self-assembly of ge and inas quantum dots on si. THIN SOLID FILMS, 508(1-2):195–199, Jun 2006.
[doi:10.1016/j.tsf.2005.08.407].
[24] Lawrence H. Friedman. Anisotropy and order of epitaxial self-assembled quantum dots. Physical Review B, in
press.
[25] Y. Obayashi and K. Shintani. Directional dependence of surface morphological stability of heteroepitaxial layers.
Journal of Applied Physics, 84(6):3141, 1998. [doi:10.1063/1.368468].
[26] C. S. Ozkan, W. D. Nix, and H. J. Gao. Stress-driven surface evolution in heteroepitaxial thin films: Anisotropy
of the two-dimensional roughening mode. JOURNAL OF MATERIALS RESEARCH, 14(8):3247–3256, Aug
1999. [doi:10.1557/JMR.1999.043].
[27] J. Tersoff and F. K. LeGoues. Competing relaxation mechanisms in strained layers. Physical Review Letters,
72(22):3570–3573, May 1994. [doi:10.1103/PhysRevLett.72.3570].
[28] B. J. Spencer, P. W. Voorhees, and S. H. Davis. Morphological instability in epitaxially strained dislocation-
free solid films: Linear stability theory. Journal of Applied Physics, 73(10):4955–4970, 1993. [doi:
10.1063/1.353815].
[29] J. M. Baribeau, X. Wu, N. L. Rowell, and D. J. Lockwood. Ge dots and nanostructures grown epitaxially
on si. JOURNAL OF PHYSICS-CONDENSED MATTER, 18(8):R139–R174, Mar 2006. [doi:10.1088/0953-
8984/18/8/R01].
[30] D. J. Srolovitz. On the stability of surfaces of stressed solids. Acta Metallurgica, 37(2):621–625, 1989.
[doi:10.1016/0001-6160(89)90246-0].
[31] H. J. Gao and W. D. Nix. Surface roughening of heteroepitaxial thin films. ANNUAL REVIEW OF MATERIALS
SCIENCE, 29:173–209, 1999. [doi:0.1146/annurev.matsci.29.1.173].
[32] P. Sutter and M. G. Lagally. Nucleationless three-dimensional island formation in low-misfit heteroepitaxy.
Physical Review Letters, 84(20):4637, 2000. [doi:10.1103/PhysRevLett.84.4637.
[33] A. A. Golovin, S. H. Davis, and P. W. Voorhees. Self-organization of quantum dots in epitaxially strained solid
films. Physical Review E, 68:056203, 2003. [doi:10.1103/PhysRevE.68.056203].
[34] A. Ramasubramaniam and V. B. Shenoy. Growth and ordering of si-ge quantum dots on strain patterned sub-
strates. JOURNAL OF ENGINEERING MATERIALS AND TECHNOLOGY-TRANSACTIONS OF THE ASME,
127(4):434–443, Oct 2005. [doi:10.1115/1.1924559].
[35] I. Berbezier, A. Ronda, F. Volpi, and A. Portavoce. Morphological evolution of SiGe layers. Surface Science,
531:231–243, 2003. [doi:10.1016/S0039-6028(03)00488-6].
[36] J. R. R. Bortoleto, H. R. Gutierrez, M. A. Cotta, J. Bettini, L. P. Cardoso, and M. M. G. de Carvalho. Spatial order-
ing in InP/InGaP nanostructures. Applied Physics Letters, 82(20):3523–3525, 2003. [doi:10.1063/1.1572553].
[37] P. Liu, Y. W. Zhang, and C. Lu. Formation of self-assembled heteroepitaxial islands in elastically anisotropic
films. Physical Review B, 67:165414, 2003. [doi: 10.1103/PhysRevB.67.165414].
[38] Y.W. Zhang, A.F. Bower, and P. Liu. Morphological evolution driven by strain induced surface diffusion. Thin
solid films, 424:9–14, 2003. [doi:10.1016/S0040-6090(02)00897-0].
[39] Yu U. Wang, Yongmei M. Jin, and Armen G. Khachaturyan. Phase field microelasticity modeling of surface
instability of heteroepitaxial thin films. Acta Materialia, 52:81–92, 2004. [doi:10.1016/j.actamat.2003.08.027].
[40] W. T. Tekalign and B. J. Spencer. Evolution equation for a thin epitaxial film on a deformable substrate. Journal
of Applied Physics, 96(10):5505–5512, 2004. [doi:10.1063/1.1766084].
[41] M. J. Beck, A. van de Walle, and M. Asta. Surface energetics and structure of the ge wetting layer on si(100).
Physical Review B, 70(20):205337, Nov 2004. [doi:10.1103/PhysRevB.70.205337].
[42] Y. H. Tu and J. Tersoff. Origin of apparent critical thickness for island formation in heteroepitaxy. Physical
Review Letters, 93(21):216101, Nov 2004. [doi:10.1103/PhysRevLett.93.216101.
[43] V. Holy, G. Springholz, M. Pinczolits, and G. Bauer. Strain induced vertical and lateral correlations in quantum
dot superlattices. Physical Review Letters, 83(2):356–359, 1999. [doi:10.1103/PhysRevLett.83.356].
[44] P. Liu, Y. W. Zhang, and C. Lu. Three-dimensional finite-element simulations of the self-organized growth of
quantum dot superlattices. Physical Review B, 68:195314, 2003. [doi:10.1103/PhysRevB.68.195314].
[45] G. Springholz, M. Pinczolits, V. Holy, S. Zerlauth, I. Vavra, and G. Bauer. Vertical and lateral ordering in
self-organized quantum dot superlattices. Physica E, 9:149–163, 2001. [doi:10.1016/S1386-9477(00)00189-2.
[46] P. Liu, Y. W. Zhang, and C. Lu. Coarsening kinetics of heteroepitaxial islands in nucleationless stranski-krastanov
growth. Physical Review B, 68:035402, 2003. [doi:10.1103/PhysRevB.68.035402].
[47] F. M. Ross, J. Tersoff, and R. M. Tromp. Coarsening of self-assembled ge quantum dots on Si(001). Physical
Review Letters, 80(5):984–7, 1998. [doi:10.1103/PhysRevLett.80.984].
[48] J. Tersoff. Kinetic surface segregation and the evolution of nanostructures. Applied Physics Letters, 83(2):353–
355, 2003. [doi:doi:10.1063/1.1592304].
[49] A. Ramasubramaniam and V. B. Shenoy. A spectral method for the nonconserved surface evolution of nanocrys-
talline gratings below the roughening transition. Journal of Applied Physics, 97(11):114312, 2005. [doi:
10.1063/1.1897837].
[50] Y. W. Zhang and A. F. Bower. Three-dimensional analysis of shape transitions in strained-heteroepitaxial islands.
Applied Physics Letters, 78(18):2706–2708, 2001. [doi:10.1063/1.1354155].
[51] L. E. Vorbyev. Handbook Series On Semiconductor Parameters, volume 1. World Scientific, 1996.
[52] A. A. Golovin, M. S. Levine, T. V. Savina, and S. H. Davis. Faceting instability in the presence of wet-
ting interactions: A mechanism for the formation of quantum dots. Physical Review B, 70:235342, 2004.
[doi:10.1103/PhysRevB.70.235342].
[53] B L Liang, Zh M Wang, Yu I Mazur, V V Strelchuck, K Holmes, J H Lee, and G J Salamo. Ingaas quantum
dots grown on b-type high index gaas substrates: surface morphologies and optical propertiesmorphologies and
optical properties. Nanotechnology, 17(11):2736–2740, 2006. [doi:10.1088/0957-4484/17/11/004].
[54] M. Meixner, R. Kunert, and E. Scholl. Control of strain-mediated growth kinetics of self-assembled semicon-
ductor quantum dots. Physical Review B, 67:195301, 2003. [doi: 10.1103/PhysRevB.67.195301].
[55] Robert Zwanzig. Nonequilbrium Statistical Mechanics. Oxford University Press, New York, 2001.
[56] C. W. Gardiner. Handbook of Stochastic Methods for Physics Chemistry and the Natural Sciences. Springer,
New York, 3rd edition, 2004.
[57] M. C. Cross and P. C. Hohenberg. Pattern formation outside equilibrium. Reviews of Modern Physics, 65(3):851–
1112, 1993. [doi:10.1103/RevModPhys.65.851].
[58] B. J. Spencer, S. H. Davis, and P. W. Voorhees. Morphological instability in epitaxially strained dislocation-free
solid films: Nonlinear evolution. Physical Review B, 47(15):9760, 1993. [doi: 10.1103/PhysRevB.47.9760].
	Introduction
	Modeling
	1D and 2D Isotropic model
	Energetics
	simple form
	general form
	Linearization
	Dynamics
	Peaks
	Stability and wetting potential
	2D Anisotropic case
	Energetics
	Elastic anisotropy
	Surface and Wetting Energy Anisotropy
	total diffusion potential
	Dynamics
	Expansion about peaks
	Correlation Functions
	Correlation Functions and SAQD order
	Periodic array
	Nearly Periodic array
	Ensemble Correlation Functions / ergodicity
	Mean fluctuation
	Correlation Function
	Stochastic Initial Conditions
	Reciprocal Space Correlation Functions
	one-dimensional
	2D isotropic
	anisotropic
	Real Space Correlation Functions
	one-dimensional
	2D isotropic
	anisotropic
	Generalizability
	Order Predictions
	Ge at 600K
	General case of 
	Discussion/Conclusions 
	Diffusion Potential
	General Form =F/H(x)
	Simple Model
	General Model
	Linearized Diffusion Potential and Anisotropy
	Linearizing the simple model
	isotropic case
	anisotropic case
	Linearizing the general model
	isotropic case
	anisotropic case
	Elastic Anisotropy
	Diffusional Anisotropy
	Correlation Functions
	Mean Values
	Variance and Convergence
	Eq. (??)
	Eq. (??)
	Eq. (??)
	Atomic Scale Cutoff

  Epitaxial self-assembled quantum dots (SAQDs) are of interest for
nanostructured optoelectronic and electronic devices such as lasers,
photodetectors and nanoscale logic. Spatial order and size order of SAQDs are
important to the development of usable devices. It is likely that these two
types of order are strongly linked; thus, a study of spatial order will also
have strong implications for size order. Here a study of spatial order is
undertaken using a linear analysis of a commonly used model of SAQD formation
based on surface diffusion. Analytic formulas for film-height correlation
functions are found that characterize quantum dot spatial order and
corresponding correlation lengths that quantify order. Initial atomic-scale
random fluctuations result in relatively small correlation lengths (about two
dots) when the effect of a wetting potential is negligible; however, the
correlation lengths diverge when SAQDs are allowed to form at a near-critical
film height. The present work reinforces previous findings about anisotropy and
SAQD order and presents as explicit and transparent mechanism for ordering with
corresponding analytic equations. In addition, SAQD formation is by its nature
a stochastic process, and various mathematical aspects regarding statistical
analysis of SAQD formation and order are presented.

Introduction
Epitaxial self-assembled quantum dots (SAQDs) represent an important step in the advancement of semiconductor
fabrication at the nanoscale that will allow breakthroughs in optoelectronics and electronics. [1, 2, 3, 4, 5, 6, 7, 8, 9,
10, 11, 12] Most frequent optoelectronic applications are high efficiency lasers with exotic wavelengths or photode-
tectors. [1, 3, 4, 5, 6, 7, 8, 10, 11, 12] SAQDs are the result of a transition from 2D growth to 3D growth in strained
epitaxial films such as SixGe1−x/Si and InxGa1−xAs/GaAs. This process is known as Stranski-Krastanow growth or
Volmer-Webber growth. [13, 1, 14, 15].
In applications, order is a key factor. There are two types of order, spatial and size. Spatial order refers to
the regularity of SAQD dot placement, and it is necessary for nano-circuitry applications. Size order refers to the
uniformity of SAQD size which determines the voltage and/or energy level quantization of SAQDs. It is reasonable
to expect that these type of order are linked, and it is important to understand the factors that determine SAQD order.
Further understanding should help in the design and simulation of both spontaneous “bottom up” self-assembly and
directed or guided self-assembly to enhance SAQD order. [16, 17, 18, 19, 20, 21, 22, 23]Here, an elaboration of and
further application of a linear analysis of SAQD order [24] is presented. The work reported here forms the basis of a
non-linear theory and modeling of SAQD order that will be reported in future work.
In [24] it was reported that one could calculate a correlation function using a linearized model of SAQD formation.
This correlation function included two correlation lengths that could be used to describe SAQD order. It was also found
that one effect of a hypothesized wetting potential was to enhance SAQD order when growth occurs near the critical
film height for 3D growth. Here, these results are expanded to create a more rigorous linearized theory of SAQD order
that will inform non-linear theories. In particular, the model is generalized to any model that combines local energy
effects such as surface energy density and non-local elastic destabilization, and the procedure for predicting order
based on any linear theory with peak wavelengths is presented. The hypothesized effect of elastic anisotropy in [24]
is verified with calculations using linear anisotropic elasticity theory. [25, 26] Details such as statistical fluctuation
and convergence are also addressed along with a discussion of the possible forms of linear anisotropic terms in SAQD
growth kinetics, and the effect of an atomic-scale cutoff in the continuum theory is addressed. Finally, the order
enhancing effect of growing near the critical threshold is explored in more detail using calculations appropriate to
Ge/Si SAQDs.
In the literature, two modes of SAQD formation are generally discussed, the thermal nucleation mode and the nu-
cleationless mode. [27, 28, 29] In the thermal nucleation mode, a 2D film surface is metastable, and the formation of
individual quantum dots is thermally activated. [27]. This growth mode leads to the formation of individual quantum
dots as uncorrelated or loosely correlated discrete events at essentially random locations. In the nucleationless mode,
the 2D film surface transitions from stable (or metastable) to unstable. In this mode, dots form everywhere at once
appearing at first as a cross-hatched ripple-like disturbance on the 2D film surface and then maturing into recogniz-
able individual dots.[27, 30, 28, 31, 32] These two modes are probably connected via an encompassing conceptual
and mathematical model1, and perhaps some of what is observed experimentally is in fact a hybrid mechanism. In
agreement with intuition, it appears that the nucleationless mode leads to a more ordered dot pattern than the thermal
nucleation mode that is dominated by randomness. 2 Thus, the presented analysis applies to the nucleationless mode.
There are various implementations of nucleationless growth models [28, 37, 38, 39, 40, 18, 34], although, there is
also a great deal of commonality among these models. In particular, they all include a non-local elastic effect and local
surface energies and/or local wetting energies. Here, a linear analysis of quantum dot order resulting from this class of
model is presented. Particular note is taken of the effects of stochastic initial conditions crystal anisotropy in general,
elastic anisotropy in particular, and the effect of varying film height as a control parameter as first introduced in [33].
A simple model similar to [28, 37, 38, 40, 18] is presented to produce numerical examples and explore the effects
of the average film height. Concurrently, a more abstract and general model is presented and analyzed that includes
non-local elastic strain effects, and a local combined surface and wetting energy. The linear model with stochastic
initial conditions and deterministic film height evolution will pave the way for more sophisticated analysis involving a
non-linear model of stochastic film height evolution.
As previously stated, one of the goals in the present work is to further explore the role of the wetting poten-
tial during growth near the stability threshold in film height. A wetting potential has been included in the analysis
and simulations in [38, 33, 37, 28]. Although somewhat controversial, the wetting potential plays an important phe-
nomenological role. It ensures that growth takes place in the Stranski-Krastanow mode: that a 3D unstable growth
occurs only after a critical layer thickness is achieved, and that a residual wetting layer persists. The physical origins
and consequences of the wetting potential are discussed in [41, 28]. The analysis presented here is usable in models
that neglect the wetting potential by simply setting it to zero. Another possibility is simply that the wetting potential
is simply an approximation to the stabilizing effect of intermixing. [42] That said, if the wetting potential is real, the
present analysis shows that it is beneficial to SAQD order to grow near the critical layer thickness.
The presented analytic formulas and linear analysis are intended to complement existing numerical models of
SAQD order. [43, 37, 44, 45] and to form a basis for future non-linear analytic analysis of SAQD order. The current
findings agree with previous work on the beneficial effects of elastic anisotropy to enhance in-plane order.
The linear analysis, of course, represents a simplification of the film evolution, and it applies only to the initial
stages of SAQD formation when the nominally flat film surface becomes unstable and transitions to three-dimensional
growth. However, the small surface fluctuation stage of SAQD growth determines the initial seeds of order or disorder
in an SAQD array; thus, the small fluctuation stage should have an important influence on the final outcome. At later
stages when surface fluctuations are large, there is a natural tendency of SAQDS to either order or ripen [33, 37, 46,
39, 47] Ordering systems tend to evolve slowly due to critical slowing down [39], while ripening tends to diminish
order further. [37] Thus, it is possible that the linear model could, in fact, yield good predictions of SAQD order.
The simplification and linearizion facilitates the development of analytic solutions that are most transparent, easily
portable to multiple material systems and have no effective limit on system size. Finally, it is virtually impossible to
have a thorough understanding of the full non-linear model without first having a thorough understanding of the linear
behavior.
The remainder of the paper is organized as follows. Section 2 presents the physical assumptions and mathematical
1It is likely that there a transition from stable, to metastable and finally to unstable. The analysis presented in [33] would appear to support such
a view where the film height acts as the control parameter driving the transition. There is also some controversy regarding whether all dot growth is
nucleationless or not. [34, 32]
2Compare various figures in [29, 35, 14, 31, 36].
approximations used to model film growth. Section 3 discusses the stochastic initial conditions and the resulting
correlation functions and correlation lengths. Section 4, presents a procedure for estimating SAQD order with an
application to Ge dots on a Si substrate. Section 5 presents conclusions, while Appendices A-F present additional
calculational details.
2 Modeling
The formation of SAQDs is modeled as a deterministic surface diffusion process with stochastic initial conditions. The
resulting equations and ultimately the sought after correlation functions are different depending on whether the film
surface is treated as one-dimensional isotropic, two-dimensional isotropic or two-dimensional anisotropic. The 1D
and 2D isotropic cases are discussed first, and then the essential differences of the 2D anisotropic model are presented.
The stochastic initial conditions need to be expressed in terms of the correlation functions that are also use to analyze
order; consequently, the discussion of the initial conditions is deferred to Sec. 3.2.
It should be noted that the results presented here are fairly general. There has been a good deal of recent work
refining the modeling of nucleationless growth processes to incorporate various phenomenological aspects of SAQD
growth. For example, the inclusion of orientation-dependent surface energy [38], strain-dependent surface energy [34]
and explicit modeling of atomic species segregation and film-substrate inter diffusion. [48] Two models are presented
here. One is a simple concrete example. It is the simplest model one can use including elastic effects surface energy
and wetting energy. The second model is more abstract and describes the general case of a local potential energy
that depends on both the film height and film height gradient. One effect that is not examined here is that of mixed
4-fold and two-fold symmetry. Such a mixing can occur due to diffusional anisotropy or surface energy anisotropy.
(Sec. 2.2.1.2 and Appendix D). However, a similar analysis procedure should work for these cases as well. The general
procedure for possible application to other models is discussed in Sec. 3.5.
The following discussion will use abstract vector notation, e.g. k instead of ki, etc. Also, because it is sometimes
computational expedient to perform one-dimensional modeling [24, 39, 17, 42], the case of a one dimensional surface
with two dimensional volume is discussed along with the case of an isotropic 2D surface. To facilitate this combined
discussion, the dimensionality of the surface will be denoted as d. In Secs. 3.3 and 3.4, d = 1, 2 will be substituted as
appropriate. Finally, much of the calculation involves reciprocal space. The convention used for the Fourier transforms
f(x) =
ddk eik·xfk, and fk = (2π)
ddx e−ik·xf(x)
following the example of [28].
2.1 1D and 2D Isotropic model
This discussion pertains to both the 1D model and the 2D isotropic model. The formation of SAQDs is modeled as a
surface diffusion process where the film height is a function of the lateral position and time. The system is treated as
deterministic with stochastic initial conditions. First, the general non-linear governing equations are presented. Then,
the linearized form is presented. Finally, the key behavior is reviewed.
The mathematical model uses film height, H(x, t) as the dependent variable and the horizontal position x and
time t as the independent variables. The film height evolves over time due to surface diffusion driven by a diffusion
potential µ(x, t) and a flux of new material Q. The surface velocity is thus
vn = nz∂tH = −∇S · D∇Sµ(x, t) +Q (1)
where nz is the vertical component of the surface normal nz = [1 + (∇H)2]−1/2, ∇S is the surface gradient, D is the
diffusivity, and ∇S · is the surface divergence.
2.1.1 Energetics
The diffusion potential µ(x, t) must produce Stranski-Krastanow growth. Thus, it must contain an elastic term that
destabilizes film growth, a surface energy term that stabilizes planar growth and a wetting energy that ensures a wetting
layer. The diffusion potential can be derived from a total free energy.
F = Felast + Fsurf. + Fwet
volume
dV ω +
surface
dAsurf. γ +
dAW (H)
where ω is the elastic energy density, γ is the surface energy density, W (H) is the wetting energy density. The
last integral corresponds to Fwet, and whether the integral should be taken over the film surface or the substrate is
ambiguous. The “simple” model (Sec. 2.1.1.1) assumes that the integral is over the substrate, while the “general”
model (Sec. 2.1.1.2) can accommodate both cases.
2.1.1.1 simple form The simplest possible model results if the integral corresponding to Fwet is taken over the
lateral positions x rather than over the actual free-surface. In concrete terms, one can use dV = d2xdz and dAsurf. =
1 + (∇H(x))2
to obtain the expression,
volume
d2xdz ω[H](x, z) +
x-plane
1 + (∇H(x))2
γ +W (H(x))
, (2)
where the “ω[H]” indicates that the elastic energy density is a non-local functional of the film height,H. The diffusion
potential µ can be found, similar to [15], by differentiating F with respect to the surface motion (Appendix A.1),
µ(x) = ΩδF/δH(x). Doing so for Eq. (2) (Appendix A.2),
µ(x) = Ω [ω(x)− γκ(x) +W ′ (H(x))] . (3)
where Ω is the atomic volume, ω(x) is the elastic energy density at the film surface (implicitly ω[H] (x,H(x))),
κ = ∇ ·
∇H(x)
1 + (∇H(x))2
]−1/2}
is the total surface curvature, and W ′(H) = ∂H(x)W (H(x)) is the
derivative of W (H(x)) evaluated at x.
2.1.1.2 general form It should be noted that Eq. (3) is not the same diffusion potential used in [38]. The wetting
potential used there can be derived by taking W (H) as an energy density of the free surface, not a density in the
x-plane. Expressions like Eq. (3) and Eq. (1) in [38] are part of a larger class of surface evolution models with more
or less the same linear behavior.
The surface and wetting energy can be combined and incorporated into a more general form, with a total free
energy Fsw and a free energy density Fsw(H,∇H) that depends on the film heightH(x) and the film height slope or
orientation ∇H(x). The total free energy is thus
F = Felast. + Fsw (4)
volume
d2xdz ω[H](x, z) +
x−plane
d2xFsw (H(x),∇H(x)) .
Fsw may not necessarily be the sum of separate surface energy and wetting energy contributions. It need only be a
local function ofH and ∇H. The corresponding diffusion potential is
µ(x) = Ω
ω(x) + F (10)sw (x)−∇ · F
sw (x)
, (5)
where F (mn)sw indicates the mth derivative with respect to H and the nth derivative with respect to ∇H. F
sw (x) =
∂H(x)Fsw (H(x),∇H(x)) and each vector component of F
sw (x) is
F(01)sw (x)
= ∂[∇H(x)]
Fsw (H(x),∇H(x)).
One can obtain the results of the simple model (Eqs. (2) and (3)) by setting
Fsw =
1 + (∇H(x))2
γ +W (H(x)) . (6)
A diffusion potential like Eq. (1) in [38] can be obtained by setting
Fsw =
1 + (∇H(x))2
[γ (∇H(x)) +W (H(x))] .
This is different from Eq. (6) in two ways. First, the surface energy density depends on the surface orientation. Second,
the Jacobian, J =
1 + (∇H(x))2
multiplies both the surface energy density and the wetting potential. Despite
these differences, the common form of the diffusion potential (Eq. (5)) among different models suggests that they
might all lead to similar linearized forms and behavior.
2.1.1.3 Linearization The diffusion potential is now linearized about the average film height H̄. In general, one
can control the amount of deposited material, and thus the average film height H̄. It is therefore useful to decompose
H(x) into the spatially averaged mean value and fluctuations about the average. Similar to [28],
H = H̄+ h(x, t). (7)
In the present calculation, H̄ is specified as constant in time. This assumption corresponds physically to a fast deposi-
tion and then an anneal. It is not too difficult to generalize to a time dependent H̄, but that is beyond the scope of this
manuscript. In [38, 49], deposition and evaporation is explicitly modeled.
All terms in µ(x, t) are now kept to only linear order in h(x, t). The elastic energy density ω is a non-local
functional of h(x, t) [40]; however, the equations generating ω(x) are translationally invariant. Thus, it is convenient
to use reciprocal space for the linearization. The curvature is trivially linearized as κ(x) → ∇2h(x) in real space or
κk → −k2hk in reciprocal space. The linearized elastic strain energy ω can be found in reciprocal space as in [15]
to be ωk = −2M(1 + ν)�2mhk, where M = E/(1 − ν) is the biaxial modulus, E is the Young modulus, ν is the
Poisson ratio, and �m is the film-substrate mismatch strain. This formula neglects possible differences in elastic moduli
between the film and substrate as in [28], but a similar method of analysis should apply to that case as well. Linearizing
Eqs. (3) and (5) in reciprocal space, µk is proportional to hk with a proportionality coefficient that depends on k and
µlin,k = f(k, H̄)hk (8)
where f(k, H̄) for three different isotropic cases, corresponding to Eqs. (3) and (5), and an abstracted general form, is
given by
f(k, H̄) =
−2M(1 + ν)�2mk + γk2 +W ′′(H̄)
; case a (Eq. (3))
−2M(1 + ν)�2mk + F 02k2 + F 20
; case b (Eq. (5))
−ak + bk2 + c ; case c (general)
. (9)
Due to isotropy, f(k, H̄) is independent of the direction of k, and only the wave number, k = ‖k‖, appears in the
right hand side. F (20)sw is the second derivative of Fsw with respect to H, and F
sw the second derivative of Fsw with
respect to ∇H. F (20)sw and F
sw depend on H̄ only; thus they are constants in the present analysis. See Appendix B.2
for more precise definitions and the derivation of f(k, H̄). Using Eq. (6), produces F (02)sw = γ and F
sw = W ′′(H̄)
which is identical to the simple case of Eq. (9), a. Case c, labeled as “general” where a, b, and c depend implicitly on
H̄ shows that f(k, H̄) for cases a and b have the same relatively simple form. It also emphasizes the dynamic effects
as opposed to the physical causes. There is a destabilizing term, −ak, a short wavelength cutoff term, bk2, and a term
that stabilizes the entire spectrum, c.
Despite the label “general,” there are of course limitations to the application of Eqs. (8) and (9). For example,
there has been recent work on the effects of strain-dependent surface energies. [34] The second form can not represent
such an effect because the derivation assumes that the surface energy only depends on local quantities, (H and ∇H)
whereas the strain effect is non-local. However, it is reasonable to conjecture that a more detailed analysis of the
effects of a strain dependent surface energy term would produce a coefficient function f(k, H̄) not very different from
the case c “general” form of Eq. (9). Thus, the following analysis may very well apply to this more exotic model, but
more study is needed to be certain.
2.1.2 Dynamics
As discussed in Sec. 2.1.1, the dynamics are derived assuming no flux of new material (Q = 0) and keeping only
terms to linear order in the height fluctuation, h(x, t). Under these assumptions, Eq. (1) can be decomposed into a
Table 1: Characteristic wave-numbers, characteristic times and associated dimensionless variables for the three cases
addressed in Eq. (9).
kc tc α β
case a 2M(1+ν)�
16DΩM4(1+ν)4�8m
γW ′′(H̄)
4M2(1+ν)2�4m
case b 2M(1+ν)�
(F (02)sw )
16DΩM4(1+ν)4�8m
F (02)sw F
4M2(1+ν)2�4m
case c a/b b3/(DΩa4) k/kc cb/a2
trivial equation for H̄ and an equation for the film height fluctuation by inserting Eq. (7).
dH̄/dt = 0 (10)
∂th(x) = −∇ · D∇µlin(x) (11)
where µlin(x) is the inverse Fourier transform of Eqs. (8) and (9), and it depends implicitly on the average film height
H̄. Note that the time dependence is implicitly while the coordinate dependence is explicit. The explicit coordinate
dependence serves to distinguish Assuming that the diffusivity D is constant, the Fourier transform of Eq. (11) gives
the linearized differential equation for the evolution of each Fourier component.
∂thk = −Dk2µk = −Dk2f(k, H̄)hk. (12)
Solving Eq. (12),
hk(t) = hk(0)e
σkt; (13)
σk = −Dk2f(k, H̄). (14)
The surface evolves in reciprocal space as an initial condition, hk(0) multiplied by an envelope function, eσkt. For
most values of H̄, this envelope function has a peak. As time passes, this peak narrows and can be approximated by a
gaussian. To analyze this behavior, appropriate dimensionless variables are defined. Then, the stability of the film is
discussed. Finally, σk is expanded about its peak to aid analytic calculations.
The time dependent behavior of the film height fluctuations is facilitated by using a characteristic wave number,
characteristic time and related dimensionless variables. For the “general” case c of Eq. (9), the characteristic wavenum-
ber is kc = a/b, and the characteristic time is tc = 1/(DΩbk4c ) = b3/(DΩa4). These characteristic dimensions can be
used to define a dimensionless wave vector, α = k/kc and a dimensionless wetting parameter β = c/(bk2c ) = cb/a
One can also define a dimensionless time, τ = t/tc. To obtain the corresponding characteristic scales for cases a and
b, one merely has to plug in the appropriate substitutes for a, b and c and follow the pattern. For example, for case
a, make the substitution a → Ω2M(2 + ν)�2m, etc. Table 1 summarizes these values for all three cases. For all three
cases, f(k, H̄) and the growth constant σk reduce to the following forms:
f(k, H̄) = f(kcα, H̄) = Ωbk2c
−α+ α2 + β
σk = σkcα = t
α− α2 − β
, (16)
where α = ‖α‖ = k/kc is the dimensionless wave number. These forms are plotted in Figs. 1a and 2. Fig. 1a shows
f(k, H̄)/Ωbk2c vs. α for an isotropic or one dimensional surface. Figs. 2 shows tcσk vs. α for a 2D anisotropic
surface (Sec. 2.2). However, the curves marked 0◦ are identical to the dispersion relation for a 1D or 2D isotropic
surface (compare Eqs. (9) and (23)).
2.1.3 Peaks
The peak growth rate and the corresponding wavenumber k can be found from Eq. (16). σk has a peak at k0 = kcα0
where
9− 32β
. (17)
Expanding σk about this peak to second order in k − k0,
σk ≈ σ0 −
σ2(k − k0)2
! ! 0.25!
increasing θk
Figure 1: Dimensionless diffusion potential prefactors vs. dimensionless wave number. (a) The one dimensional or
isotropic case with β = 0.3. (b) The elastically isotropic case with anisotropy �A = 0.1 (see Eq. (22)).
increasing θk
22.5◦
increasing θk
22.5◦
Figure 2: Dimensionless growth constant vs. dimensionless wave number. Curves are plotted for the elastically
anisotropic case, but the curves marked 0◦ are the same as for the isotropic cases. In (a), β = 0. In (b) β = 0.2.
Figure 3: Exponential Envelope eσkt as function of α for β = 0.208 and t/tc = 100. (a) 2D isotropic surface. (b) 2D
anisotropic surface with � = 0.1236.
The two constants are
t−1c α
0 (α0 − 2β) , (18)
σ2 = t
c (3α0 − 4β) . (19)
Inserting this approximation for σk into Eq. (13),
hk(t) = hk(0)e
σ0te−
2σ2t(k−k0)
. (20)
The individual initial surface fluctuation components grow with a gaussian shaped envelope. An example of this
envelope is plotted in Fig. 3(a). Notice that in two dimensions, the envelope forms a ring as the peak is about the
wave-number k0 but not about any particular point in the k-plane.
2.1.4 Stability and wetting potential
Stranski-Krastanow growth is marked by a transition from stable two-dimensional growth to unstable three-dimensional
growth once a critical height Hc is reached. [1] Eqs. (17), (18) and (20) are useful for analyzing the transition from
stable to unstable growth. In order for this transition to occur, there must be some stabilizing term in the diffusion
potential. In the present model, this means that there must be some surface energy-like term that varies strongly with
film height. This condition equates to stating that W ′′(H̄) or F 20sw or c (Eq. (9)) must be rather large if H̄ < Hc.
However, as H̄ increases, these terms are reduced. Finally, when H̄ > Hc, this term is no longer capable of stabilizing
the film against fluctuations of all possible wavelengths.
The critical value Hc can be found using the analysis from [33]. By inspection of Eqs. (8), (9) and (12), modes
with f > 0 increase the total free energy F as they grow; thus, they are stable and decay with time. Modes with
f < 0 decrease the total free energy F as they grow; thus, they are unstable and grow with time. This growth and
decay rule is easily verified by inspection of Eq. (14). Thus, stable growth occurs when f(k, H̄) > 0 for all values of
k, and unstable growth occurs when f(k, H̄) < 0 for some values of k. Thus, the transition from stable to unstable
growth occurs when the minimum value of f(k, H̄) just becomes negative. Using the same dimensional analysis as in
the previous section and following the discussion of [33], one finds that the minimum value, fmin = Ωbk2c (β − 1/4),
occurs at kmin/kc = αmin = 1/2. fmin first becomes negative, and the transition to unstable growth occurs when
the dimensionless wetting parameter (Table 1) drops to a critical value, β = 1/4 . β > 1/4 stable 2D growth, and
β < 1/4 unstable 3D growth. It is reasonable to suppose that W (H̄), W ′′(H̄), and thus β are positive monotonically
decreasing functions of H̄ so that the interface becomes less important for large values of H̄. For example, in [50] it is
assumed that W (H) = B/H, where B is constant. When β → 0, corresponding to large H̄, the case discussed in [28]
is obtained. A similar analysis can be done for cases b and c once one specifies how the terms F (20)sw and F
sw or a, b
and c depend on H̄.
Using a guessed form for a wetting potential, one can find the critical film heightHc by setting β = 1/4 . Applying
this condition to case a in Eq. (3)
W ′′(Hc) = γk2c/4.
Using the wetting potential of [50] as an example, W (H) = B/H,
Hc = 3
8B/(γk2c ) =
8Bγ/(2M(1 + ν)�2m)2. (21)
Conversely, one can fit a wetting potential to an observed or reasonable critical layer thickness from the same condition.
Using the example wetting potential from [50],
(2M(1 + ν)�2m)
as stated in [50].3
2.2 2D Anisotropic case
Crystal anisotropy leads to a dispersion relation σk that is both quantitatively and qualitatively different from the
isotropic case. Here the effect of elastic anisotropy is discussed in most detail. Other sources of anisotropy are
the surface and wetting energies. For example, in [38] the surface energy density is orientation dependent which
introduces a possible anisotropy in the dispersion relation. Possible sources of anisotropy are an anisotropic elastic
stiffness tensor, an orientation dependent surface energy or wetting potential or anisotropic diffusion. As discussed
below, the form of anisotropy to linear order in the height fluctuation, h, is somewhat restricted. Results are presented
for 4-fold symmetric surfaces, that is surfaces that have invariant dynamic evolution laws when rotated by 90◦. Possible
complications arising from 2-fold symmetric anisotropic terms (with 180◦ rotational symmetry) are also discussed. As
for the isotropic case, first the energetics are discussed, then the dynamics, and finally the expansion about the peaks
in the dispersion relation, σk.
2.2.1 Energetics
The discussion of energetics will first treat the effects of elastic anisotropy and then anisotropy resulting from surface
or wetting like terms.
3This result from [50] corresponds to the choice Fsw(H,∇H) =
1 + (∇H)2γ+W (H). However, the numerical model in [50] appears to
use Fsw(H,∇H) =
1 + (∇H)2 [γ(∇H) +W (H)]. This difference should lead to a slightly different critical film height in their numerical
model from the one that they predicted (Eq. (21)).
Figure 4: Plot of Eθk/(M�
m) for various materials. Symbols indicate values calculated using Appendix C. Solid lines
are the interpolation (Eq. (22)) using the values from Table 2.
Table 2: Elastic constants [51] and calculated values (see Appendix C) for various materials of interest at T = 300K.
c11 c12 c44 M
1011 ergcm3 10
11 erg
cm3 10
11 erg
cm3 10
11 erg
Ge 12.60 4.40 6.77 13.93 2.16 1.906 0.1176
Si 16.60 6.40 7.96 18.07 2.22 1.997 0.1005
InAs 8.34 4.54 3.95 7.94 2.70 2.09 0.226
GaAs 11.90 5.34 5.96 12.45 2.15 1.87 0.1302
2.2.1.1 Elastic anisotropy One would like to obtain a simple symbolic expression for the elastic energy density at
the free surface, ωk, to first order in hk for the elastically anisotropic case. Similar discussions can be found in [25, 26].
For the isotropic case, ωk = −2M(1 + ν)�2mhk. For the anisotropic case,
ωk = −Eθkkhk
where the prefactor Eθk is the decrease in elastic energy at the surface per unit wave number (k → 1) and unit amplitude
(hk → 1) . It is not constant, but instead depends on the θk, the angle that k makes with the x−direction. The case of a
cube-symmetry elastic stiffness tensor such as for Si is considered where one must specify three elastic constants c11,
c12 and c44. [51]. Growth on a (100) surface will produce an elastic energy prefactor Eθk that is four-fold symmetric
(symmetric upon rotations by 90◦). A procedure similar to [25, 26] based on a first order perturbation analysis is
followed (Appendix C). A relatively simple interpolation formula [24] is hypothesized and then verified numerically.
The interpolation procedure, suggested in [24] uses the lowest possible order expansion in sin(θk) and cos(θk)
that has the appropriate four-fold symmetry and then interpolates between θk = 0◦ and θk = 45◦. Thus,
Eθk = E0◦
1− �A sin2 (2θk)
where �A = (E0◦ − E45◦)/E0◦ is an anisotropy factor. This lowest order form turns out to be a very good fit to
numerical calculations (Fig. 4). Table 2 gives values of E0◦ and �A for some systems of interest. In the elastically
isotropic case, E0◦ = E45◦ = 2M(1 + ν) so that �A = 0.
There are two important differences from the elastically isotropic case. The first is obvious, that Eθk depends on
angular orientation, θk. The second is that the peak value of ωk is not the same as that for the elastically isotropic case
because in general, E0◦ 6= 2M(1 + ν). In [24], where the purpose was simply to investigate the mechanism by which
elastic anisotropy effects order, this second difference was neglected.
2.2.1.2 Surface and Wetting Energy Anisotropy The surface energy and wetting potential can be additional
sources of anisotropy if they depend on the surface orientation so that γ → γ(∇H) or W (H) → W (H,∇H) (for
example, [52, 38]). Then, to first order in h ,
µsurf.,k = Ω
γk2 + k · γ̃′′ · k
where γ̃′′ is the (2×2) matrix or Hessian matrix that results from taking the second derivatives of γ(∇H) with respect
to the two components of ∇H (Appendix B.1). Similarly
µwet,k = Ω
W (20) + k · W̃(02) · k
where W (20) and W̃(02) are the second derivatives of W (H,∇H) with respect to H and ∇H (Appendix B.1). For
both µsurf.,k and µwet,k, the first term is isotropic, and the second term contains any possible anisotropy.
The rank of the γ̃′′ and W̃(02) matrices greatly restricts the possible forms of the additional anisotropy. These
(2 × 2) matrices must be either two-fold symmetric or perfectly isotropic. Thus, if the surface energy and wetting
potential are four-fold symmetric as Eθk is, then γ̃
′′ → γ′′, a scalar, and W̃(02) → W (02), a scalar, and neither one
contributes any additional anisotropy. They do, however, help to stabilize or further destabilize the 2D surface as they
add terms proportional to k2. The effect of these additional terms is indistinguishable from the effect of varying the
value of the surface energy density, γ. [52, 31]
It should be noted that the (100) surface of a diamond or zinc-blend structures allows for anisotropy that is only 2-
fold symmetric (rotations by 180◦). Thus, they could “break” the four-fold symmetry that occurs when one considers
the elastic anisotropy alone. However, this “broken” symmetry is somewhat dubious because even the diamond and
zinc-blend structures have a screw symmetry (rotations by 90◦ and translation in the [100] direction by half a lattice
vector). Thus, if for example, W (H,∇H) is anisotropic with two-fold symmetry to linear order, there must be a
fast oscillation with changes in the film height H. In Appendix D, a similar term related to anisotropic diffusion is
discussed. There does not appear to be any evidence for this two-fold symmetry in the case of (100) surfaces of IV/IV
systems such as Ge/Si, but in III-V/III-V systems the four-fold symmetry of the (100) surface may indeed be “broken”
in this way corresponding to either a surface energy anisotropy or a diffusional anisotropy. [53, 54]. Further analysis
of such terms in any more detail would greatly complicate the present discussion, so it is left for future work. Most
of the modeling literature avoids this complication by not including the symmetry-breaking of the zinc-blend surface,
for example [25, 26, 38].
One can perform a similar analysis of the combined surface and wetting potential, Fsw(H,∇H) (case b). To linear
order the resulting anisotropic diffusion potential is (Appendix B.2)
µsw,k = Ω
F (20)sw + k · F̃
sw · k
Again, F̃(02)sw is a rank 2 tensor, and all of the same symmetry considerations apply here as well.
Because the two-fold symmetry anisotropic terms are excluded from the current discussion, and isotropic terms
simply “renormalize” the effective of surface energy, there will be no further consideration of anisotropy resulting
from the surface energy or wetting potential in this discussion. Further calculations will proceed assuming that the
surface energy density, γ, nor the wetting potential,W (H), depend on ∇H or similarly that Fsw(H,∇H) has a purely
isotropic dependence on ∇H. This assumption can be made without affecting any of the qualitative results.
2.2.1.3 total diffusion potential Having dispensed with the discussion of the various sources of anisotropy, the
total diffusion potential is stated for the case of 4-fold symmetric elastic anisotropy and a completely isotropic surface
energy and wetting potential. µk = f(k, H̄) with
f(k, H̄) =
1− �A sin2(2θk)
k + γk2 +W ′′(H̄)
; case a (Eq. (3))
1− �A sin2(2θk)
k + F (02)sw k2 + F
; case b (Eq. (5))
1− �A sin2(2θk)
+ bk2 + c ; case c (general)
. (23)
Table 3: Characteristic wave-numbers, characteristic times and associated dimensionless variables for the three cases
addressed in Eq. (9)
kc tc α β
case a E0◦/γ γ3/(DΩE40◦) k/kc γW ′′(H̄)/E20◦
case b E0◦/F
/(DΩE40◦) k/kc F
sw /E20◦
case c a/b b3/(DΩa4) k/kc cb/a2
2.2.2 Dynamics
The dynamics is governed by surface diffusion, just as for the fully isotropic case. It is assumed that the diffusivity
is isotropic as was done for the surface energy and the wetting energies; thus, all anisotropy in the film evolution
dynamics comes from elastic effects alone. The possibility and effects of an anisotropic diffusion potential is discussed
in Appendix D (also see [54]). The time dependence of the surface perturbations simply follows Eqs. (13) and (14),
but with Eq. (23) used for f(k, H̄). As for the isotropic case, appropriate characteristic wave numbers (kc) and
time scales (tc) can be found for each of the three cases along with the associated dimensionless wave vector α and
dimensionless wetting parameter β. These are listed in Table 3. The dispersion relation, σk can be expressed in terms
of these dimensionless variables (α and β), giving
σk = σkcα = t
1− �A sin2(2θk)
− α2 − β
. (24)
The stability behavior is essentially the same as for the isotropic case with a transition occurring at β = 1/4 corre-
sponding to H̄ = Hc.
2.2.3 Expansion about peaks
σk has 4 peaks at (k, θk) = (k0, π[n− 1]/2) with k0 = kcα0 (Eq. (17)) and n = 1 . . . 4. In vector form, there are four
peaks at
kn = k0 (cos(π(n− 1)/2)i + sin(π(n− 1)/2)j) .
Similar to the isotropic case, σk can be expanded about individual peaks so that in the vicinity of peak n, σk ≈ σn
σn = σ0 −
σ‖(k − k0)2 −
0(θk − nπ/2)
where σ0 is given by Eq. (18), σ‖ = σ2 given by Eq. (19), and
σ⊥ = 8�Aα0t
In terms of the vector components parallel and perpendicular to kn, k‖ and k⊥ respectively,
σn = σ0 −
σ‖(k‖ − k0)2 −
k‖ = cos[π(n − 1)/2]kx + sin[π(n − 1)/2]ky , and k⊥ = − sin[π(n − 1)/2]kx + cos[π(n − 1)/2]ky . The time
evolution of hkin the vicinity of one of the kn is
hk(t) ≈ hk(0)et(σ0−
2σ2(k‖−k0)
2− 12σ⊥k
3 Correlation Functions
Correlation functions and associated constants such as correlation lengths can be very useful for characterizing order.
In particular, the autocorrelation function (Eq. (25)) and its Fourier transform (Eq. (26)) also known as the spectrum
function can give a very good characterization of dot order (Figs. 6a and c and 5b, e and h). The autocorrelation
function is denoted CA(∆x) where ∆x is the difference vector between two points in the x−plane. The spectrum
function is a function of k, and it is denoted CAk . The goal here is to be able to predict these two functions and to
h(nm)
h(nm)
h(nm)
describe them quantitatively in a manner that can be used to characterize SAQD order with just a few numbers. The au-
tocorrelation function is the result of a spatial average over one experiment or one simulation (numerical experiment).
It is regular and repeatable because it is closely tied to the correlation function and spectrum function that results from
an ensemble average (Eqs. X and X). These are denoted as C(∆x) and the spectrum Ck respectively. Note that the
ensemble averaged functions do not have a superscript “A.” These ensemble average correlation functions are useful
in the analysis of stochastic ordinary and partial differential equations. [55, 56]. From a strictly technical viewpoint,
the spatial average and the ensemble average are not exactly the same; however, they are closely enough connected
that it is reasonable to use one as a substitute for the other (Sec. 3.1 and Appendix 3).
In the following, the analysis of SAQD order via autocorrelation and correlation functions is discussed (Sec. 3.1).
Then, the stochastic initial conditions are discussed (Sec. 3.2). Then, the prediction of the Fourier transforms of the
correlation functions is discussed (Sec. 3.3). The real-space correlation functions are presented (Sec. 3.4). Finally,
there are some notes regarding generalizing the analysis method to any dispersion relation that has peaks (Sec. 3.5),
for example, peaks related to broken four-fold symmetry or growth on a miscut substrate.
3.1 Correlation Functions and SAQD order
Auto-correlation functions are well-suited for investigating SAQD order. The autocorrelation function is defined as
CA(∆x) =
d2x′ h(∆x + x′)h(x′)∗. (25)
Its Fourier transform sometimes called the spectrum [56], spectrum function or power spectrum is
CAk =
(2π)d
d2∆x e−ik·∆xC(∆x) =
(2π)d
, (26)
where A is the projected area of the film in the x − y−plane. A periodic array of SAQDs leads to a periodic auto-
correlation function. A nearly periodic array leads to a range-limited periodic auto-correlation function. The ensemble-
mean of these autocorrelation functions can be calculated, and it is a good predictor of a SAQD order.
3.1.1 Periodic array
Consider a perfectly periodic height fluctuation corresponding to a perfect lattice of SAQDs,
h(x) =
exp [ikn · (x− xO)] (27)
plus higher order harmonic, where the dots have a height proportional to h0, N is the degree of symmetry, probably,
4-fold or 6-fold, xO is a random origin offset.
kn = k0
2π(n− 1)
i + sin
2π(n− 1)
In a linear analysis, the higher order harmonics do not come into play, so they are neglected here. In reciprocal space,
e−ikn·xOδd(k− kn)
plus higher order harmonic. The autocorrelation function is found by plugging Eq. (27) into Eq. (25) and simplifying,
CA(∆x) =
)2 N∑
exp [ikn ·∆x] (28)
plus higher order harmonic. In finding Eq. (28), the relation∫
d2x′ ei(km−kn)·x
= Aδkmkn = (2π)
dδd(km − kn) (29)
has been used. δkk′ is the Kronecker Delta, and δd(k− k′) is the Dirac Delta. Eq. (29) will be helpful whenever it is
necessary to take an areal average or sum over Dirac Delta functions. In reciprocal space,
CAk =
(2π)2
m,n=1
δ2(k− km)δ2(k− kn)
δ2(k− ki) (30)
plus higher order harmonics, where δd(k − kn) = (A/(2π)d)δkkn .4. Thus, the order of the SAQD lattice manifests
itself as periodic functions in real-space (Eq. (28)) and sharp peaks in reciprocal space (Eq. (30)).
3.1.2 Nearly Periodic array
A nearly periodic arrays shows deviation from perfect order. This deviation is shows itself by broadening of the peaks
of the spectrum function, CAk , and by range limited periodicity of the real-space autocorrelation function, C
A(∆x).
These two measure of disorder are naturally related.
The disorder in lateral dot size ∆size and spacing, ∆spacing are related to each other and to the broadening of the
peaks in CAk (Fig. 6.a and c). Prior to ripening, the size and spatial order should be related, as the volume of a dot
should be proportional to the amount of nearby material. If the SAQDs have nearly uniform size and spacing (peak-
to-peak distance) L0, the reciprocal space autocorrelation function will be tightly clustered around the wavenumber
characterizing the dot spacing k0 = 2π/L0. There are a number of such peaks depending on the system symmetry
(Fig. 6.a and c), but consider just one. Since the order is not perfect, the peak will have a finite width. Consequently,
there will be a scatter in the dot size. Since L0 = 2π/k0, the scatter in dot spacing (∆spacing) is related to the scatter in
Fourier components (∆k). Taking the derivative of the spacing-wavenumber relation and rearranging,
∆spacing
It is reasonable to expect that the fractional disorder in size (∆size /Lsize) is given by a similar (if not exactly the same)
number.
Another way to view spatial order (periodicity) is not by dot-dot distances, but the distance over which the dot array
can be considered periodic. This limited periodicity is evident in the film height autocorrelation function (Eq. (25) and
Figs. 5.b, e and h). Consider two distant dots. Their position will be completely uncorrelated, so it will be completely
random as to whether one position corresponds to a peak or a valley. Thus, for a large differences in position the
autocorrelation function tends to zero.
CA(∆xlarge) = 0
Similarly, the mean-square fluctuation of the film height can be large so that
CA(∆x = 0)� 0.
The distance over which the autocorrelation function, CA(∆x) decays to 0 is the correlation length, Lcor. Thus, Lcor
is a reasonable measure of spatial order.
The two measures of order ∆spacing and Lcor are intrinsically linked. The well known rule of Fourier transforms
states that the product of the real-space and reciprocal space widths must be greater than or equal to unity. Thus,
∆kLcor ≥ 1, or ∆spacing ≥ 2πL20/Lcor. Similarly, one can expect that ∆size ∼ L2size/Lcor. Thus, assuming that dot size
is governed by the amount of nearby material, small dispersions in dot size are only possible if there is long correlation
length.
3.1.3 Ensemble Correlation Functions / ergodicity
SAQDs are seeded by random fluctuations. Consequently, each experiment or simulation must be treated as just
one possible realization, and the autocorrelation function will be different for each realization. Thus, for analytic
4Eq. (29) has been used to help with summation
predictions, one must rely on ensemble averages. In [24], it was assumed that the ensemble average correlation
function was a good description of a SAQD order, an assumption that was born out by numerical calculations. Now,
this relation is put on a more solid ground. In particular, it is found that the ensemble correlation functions provide
good estimates of the auto correlation function and spectrum function produced by any particular realization. First,
it is shown that the mean value of the film-height fluctuation is zero. Then the method to calculate the ensemble-
averaged autocorrelation function and spectrum function is presented. Additional mathematical details are presented
in Appendix 3.
3.1.3.1 Mean fluctuation It is fairly straightforward to show that the ensemble mean film-height fluctuation is
zero. The governing dynamics (Eq. (12)) is invariant upon the substitution h(x, t) → −h(x, t). Thus, assuming that
one does not bias the initial conditions the mean fluctuations must be zero for all time,
〈h(x, t)〉 = 〈−h(x, t)〉 = 0, and 〈hk(t)〉 = 0.
This is a common situation, and it is most appropriate to characterize the film height fluctuations using the two-point
correlation function (or simply “the correlation function”). [55]
3.1.3.2 Correlation Function The autocorrelation function can be estimated by its ensemble average. Further-
more, this ensemble average is equivalent to the correlation function that can be easily calculated analytically. These
relations are first discussed for the real-space correlation functions and then their Fourier transforms. First, the statis-
tical properties of the autocorrelation function are discussed. Then the statistical properties of the spectrum function.
Finally, the method to The main results are reported here, and details of derivations are reported in Appendix E.
First it is noted that the autocorrelation function averaged over all realizations is equal to the ensemble correlation
function. 〈
CA(∆x)
= C(∆x), where C(∆x) = 〈h(∆x)h(0)〉 , (31)
where 〈. . . 〉 indicate an ensemble average. Eq. (31) assumes that the model of film-growth is translationally invariant.5
This relationship is fortunate, in that it allows one to predict the “typical” autocorrelation function using analytic tools
that apply only to ensemble averages.
Second, it is noted that as the area that is used to calculate the autocorrelation function becomes large, the autocor-
relation function tends towards it mean value,
CA(∆x) ≈ C(∆x) +O[A−1/2], (32)
where O[A−1/2] indicates statistical fluctuations about the mean value that become smaller and smaller as the area in
an experiment or the simulation area in a numerical experiment becomes larger. These fluctuations or noise die out as
A1/2. For example, the autocorrelation functions in Figs. 5.e and h are very close to the ensemble average autocorre-
lation functions Figs. 5.f and h, but have random fluctuations that are most visible far from the origin. This property,
that averaging over a parameter such as position is equivalent to averaging over all realizations, is known as ergodicity.
Individual realizations are tightly distributed about a “typical” behavior. This tight distribution lends credibility to the
notion that one can have representative experiments or simulations. Unfortunately, the “demonstration” of Eq. (32) in
Appendix E is not as general as one might like. Rigorously, it applies when the Fourier components of film height (hk)
are independent and normally distributed; however, it is reasonable to conjecture that a relationship like Eq. (32) holds
whenever the statistical distribution of film heights is suitably bounded as the boundedness of CAk plays an important
role in the derivations.
In reciprocal space, one finds that the ensemble-mean spectrum function is〈
= Ck, (33)
where Ck is defined as the prefactor appearing in the reciprocal-space two-point correlation function.
Ckk = 〈hkh∗k〉 = Ckδ
d(k− k′) = Ck
(2π)d
δkk′ , (34)
5A quick survey of literature will find that, virtually all published continuum models of SAQD formation are translationally invariant.
where Eq. (29) has been used. This form for the two-point correlation function in reciprocal space occurs if and only
if the system is translationally invariant. Eq. (33) is valuable because one can solve for Ck analytically in the linear
case or using various analytic approximations in the non-linear case. Unlike the autocorrelation function, the spectrum
function fluctuates greatly about its mean. In fact, the fluctuations are about 100% (Appendix E.2). These large
fluctuations result in the commonly observed speckle pattern for the spectrum function CAk (Figs. 6.a and c). Contrast
this pattern with ensemble-mean spectrum function Ck shown in Figs. 6.b and d. These speckles can be removed by a
smoothing operation, and a relation similar to Eq. (32) results (Appendix E.2.2). Finally, it should be noted that just
as CAk is the Fourier transform of C
A(∆x), Ck is the Fourier transform of C(∆x) (Appendix E.1).
3.2 Stochastic Initial Conditions
To model or simulate the formation of SAQDs, it is absolutely essential to include some sort of stochastic effect. An
initially flat film h(x, 0) = 0 is in unstable equilibrium. Thus, to seed the formation of quantum dots, it is necessary
to perturb the flat surface. The simplest method to do this is to use stochastic initial conditions with deterministic
evolution. One can tenuously suppose that white noise initial conditions do not “bias” the ultimate evolution of the
film. [57] Thus, the initial conditions are taken from an ensemble with zero mean,
〈h(x, 0)〉 = 0. (35)
and a spatial correlation function,
C(x,x′, 0) = 〈h(x, 0)h(x′, 0)∗〉 = ∆2δd (x− x′) , (36)
where the brackets 〈. . . 〉 indicate an ensemble average, ∆ is the noise amplitude, and δd(x) is the d−dimensional
Dirac Delta function. White noise conditions have an infinite amplitude which is not physical. Thus, a minimum
modification can be made to “cut off” the infinite fluctuations.
C(x,x′, 0) =
(2πb20)
(x− x′)2
In the limit b0 → 0, this correlation function reverts to the white noise correlation functions.
In reciprocal space,
Ckk′(0) = 〈hk(0)h∗k′(0)〉
= (2π)−2d
ddx′ e(−ik·x+ik
′·x′)C(x,x′, 0)
(2π)d
δd(k− k′)
Letting b0 → 0, the white noise reciprocal space correlation function is obtained. Thus, the initial spectrum function
Ck(0) =
(2π)d
The atomic-scale has a small and short-lived influence on the final film morphology (Appendix F), but the cutoff
procedure is useful for choosing a reasonable value of ∆2. It seems reasonable to choose ∆2 so that the initial r.m.s.
fluctuation
C(0, 0) = 〈h(0, 0)h(0, 0)∗〉1/2is one monolayer (1 ML). Also, choosing b0 = 1 ML as the atomic scale
cutoff is
∆2 = (2π)d/2(1 ML)2+d, (38)
where the natural unit 1 ML is, of course, material dependent.
Using stochastic initial conditions, one can integrate individual initial conditions to obtain representative samples
and then average over many realizations, the Monte Carlo approach, or one can calculate analytically, the statistical
measures of the ensemble. The ensemble statistical measures are strongly related to the statistical measures of order
for an individual realization, so the second approach is opted for here. Thus, the predicted SAQD order is ultimately
stated in terms of ensemble correlation functions.
3.3 Reciprocal Space Correlation Functions
The reciprocal space correlation function, Ckk′ , and spectrum function, Ck, are calculated for the 1D and 2D isotropic
case and then for the 2D anisotropic case. Generally Ck includes the length scales introduced in Sec. 2.1.3 as well as
the atomic scale cutoff b0.
Ckk′ = 〈hk(t)h∗k′(t)〉 = e
(σk+σk′ )t 〈hk(0)hk′(0)∗〉
(2π)2
e(σk+σk′ )t−
δ2(k− k′). (39)
Without much error, b0 can be neglected in the exponential (Appendix F). Using Eq. (34), the spectrum function is
then identified as
(2π)d
e2σkt. (40)
Ck is now calculated for each model: 1D isotropic, 2D isotropic and 2D anisotropic.
3.3.1 one-dimensional
The one dimensional surface is the simplest, so it is treated first. The spectrum function is simply
e2σ0t−
2 (2σ2t)(k−k0)
Ck has a peak at k = ±k0i. One can easily read off the correlation length as
Lcor =
2σ2t = k
2(3α0 − 4β)(t/tc). (41)
so that
e2σ0t−
cor(k−k0)
This approximation is valid when k0Lcor � 1. In terms of kx,
e2σ0t
cor(kx−k0)
cor(kx+k0)
3.3.2 2D isotropic
The 2D isotropic case is very similar;
(2π)2
e2σ0t−
cor(k−k0)
, (42)
where Lcor is the same as in Eq. (41). It has maximum that forms a ring in the k−plane as graphed in Fig. 6.b.
3.3.3 anisotropic
The anisotropic spectrum function is
(2π)2
e2σ0t
‖(k‖−k0)
2− 12L
⊥ , (43)
where
2σ‖t = k
(6α0 − 8β)(t/tc), (44)
2σ⊥t = k
16�α0(t/tc), (45)
k‖ = cos[π(n− 1)/2]kx + sin[π(n− 1)/2]ky , and k⊥ = − sin[π(n− 1)/2]kx + cos[π(n− 1)/2]ky and it is graphed
in Fig. 6.d. This approximation is valid when k0L‖ � 1 and k0L⊥ � 1.
Figure 6: CAk and Ck for Ge/Si as discussed in Sec. 4. (a,b) 2D isotropic surface. Eq. (42) is used for Ck. (c,d) 2D
anisotropic surface. Eq. (43) is used for Ck.
Figure 7: 1D isotropic surface in real space for Ge/Si as discussed in Sec. 4. (a) Example of h(x) plotted over a length
of 8Lcor. (b) corresponding reals space correlation functions plotted for range ±4Lcor. Filled plot is an example of
CA(∆x). Solid line isC(∆x) (Eq. (46)).
Loosely speaking, one can argue that the isotropic case is similar to letting �A → 0 in Eq. (45) so that the
perpendicular correlation length is always 0 regardless of time. A more conservative approach would be to argue that
L⊥ ≈ 2π/k0 for the isotropic model via inspection of Figs. 6(a) and (b). Even still, the more conservative result
guarantees that the perpendicular correlation length will always be the same as the dot spacing; thus, it will always
limit SAQD order to the first nearest neighbor at best.
3.4 Real Space Correlation Functions
The real space correlation functionsC(∆x) are now calculated for the 1D and 2D isotropic cases and the 2D elastically
anisotropic case.
3.4.1 one-dimensional
In one dimension,
C(∆x) =
dkx e
ikx∆xCk
e2σ0t−
2/L2cor2 cos (k0x) . (46)
Thus, C(∆x) has a damped periodicity indicating that it is imperfectly periodic (Fig. 7).
3.4.2 2D isotropic
In two dimensions with elastic isotropy,
C(∆x) =
d2k eik·∆xCk
(2π)2
e2σ0t
dk kei(k∆x cos(θk−θ∆x)e−
cor(k−k0)
Performing the angular integration first,
C(∆x) =
e2σ0t
dk kJ0(k∆x)e
− 12L
cor(k−k0)
where J0 is the zeroth Bessel function. In general, this integral is best performed numerically; however, it can be
solved in two important cases: ∆x→ 0 and Lcor →∞ (corresponding to long times). In the first case,
C(∆x) =
e2σ0t
dk ke−
cor(k−k0)
Under the same conditions that Eq. (42) is valid (k0Lcor � 1), the lower limit of the integral can be approximated as
−∞ so that
C(∆x = 0) =
∆2k0√
2πLcor
e2σ0t. (47)
This function gives the mean square surface height fluctuation. In the second case where Lcor →∞, e−
cor(k−k0)
(2π)1/2L−1cor δ(k − k0), so that
C(∆x) =
∆2k0√
2πLcor
e2σ0tJ0 (k0∆x) . (48)
This correlation function is the most ordered case for a 2D isotropic surface. It is graphed in Fig. 5c.
3.4.3 anisotropic
To find the real-space correlation function for the elastically anisotropic case, it is best to find the contribution from
each peak and then sum so that
C(∆x) =
(2π)2
e2σ0t
Cn(∆x) (49)
where
Cn(∆x) =
d2k eik·xe−
‖(k‖−k0)
2− 12L
∆x can be decomposed into the directions parallel and perpendicular to kn, so that ∆x‖ = cos(π(n − 1)/2)∆x +
sin(π(n− 1)/2)∆y and ∆x⊥ = − sin(π(n− 1)/2)∆x+ cos(π(n− 1)/2)∆y. Thus,
Cn(∆x) =
dk‖ e
ik‖∆x‖− 12L
‖(k‖−k0)
dk⊥ e
ik⊥∆x⊥− 12L
2 (x2‖/L2‖+x2⊥/L2⊥)eik0x‖ .
Plugging into Eq. (49),
C(∆x) =
πL‖L⊥
e2σ0t
2 (x2/L2‖+y2/L2⊥) cos(k0x) + e
− 12 (x2/L2⊥+y2/L2‖) cos(k0y)
. (50)
3.5 Generalizability
The dynamics and analysis used here were for a specific model, but the general procedure for analyzing the order
resulting from a linearized model should hold for any model with well-separated peaks in the dispersion relation, σk.
The procedure to follow is:
1. Generate the dispersion relation, σk as some function of k.
2. Find the peaks in the dispersion relation, kn, (n = 1 . . . N )
3. Expand about the peaks to generate the peak values, σn, and local Hessian matrix,(
∂ki∂kj
The spectrum function is then approximately
Ck(t) ≈
(2π)2
e2σnt exp
t (k− kn) · H̃n · (k− kn)
. (51)
4. Find the Eigenvalues of the local Hessian matrix, (Hn)I and (Hn)II . They should be negative, if there is a
peak at kn
5. Use the eigenvalues to determine the correlation lengths, (Ln)I =
2 |(Hn)I | t and (Ln)II =
2 |(Hn)II | t.
The real-space correlation function is
C(∆x, t) ≈
(Hn)I (Hn)II
e2σnt exp
x · H̃−1n · x
eikn·x. (52)
The “goodness” of these approximate forms requires that (Ln)
I and (Ln)
II be much less than the spacing
between peaks in the correlation function so that the gaussians do not overlap greatly. A reasonable test for
this no-overlap condition is ‖kn‖ (Ln)I � 1 and ‖kn‖ (Ln)II � 1, assuming that the peaks are not large in
number or very closely spaced.
4 Order Predictions
The real-space correlation function formulas (Eqs. (46), (47), and (50)) and correlation length formulas (Eqs. (41), (44)
and (45)) can now be used to estimate the order of SAQDs. Ge on Si is chosen for this example because this system has
received the most attention from theoretical work [58, 38, 31, 18, 39, 41, 27, 25, 26, and others], and it is the simplest
since it involves the diffusion of a single species. The procedure described below tries to predict the amount of order
when an initial atomic-scale fluctuation becomes “large”. “Large” is taken to be greater than atomic-scale. Beyond
this point, one would expect non-linear terms to become important. An example is presented for Ge on Si at 600K to
compare and contrast the 2D anisotropic results with the 1D isotropic and 2D isotropic results. The predictions are
also compared with a linear numerical calculation on a discrete reciprocal-space grid to test the approximations made
and to illustrate the relation between the surface profile (h(x)), the example autocorrelation functions (CA(x) and
CAk ) and the ensemble correlation functions (C(∆x) and Ck). Figs. 6, 7 and 5 show these results. Finally, the relation
between average film height and order is investigated.
4.1 Ge at 600K
The formulations for the three discussed cases are implemented for Ge/Si at 600K. The correlation lengths are esti-
mated for the end of the linear regime where fluctuations become large (greater than atomic scale). First, appropriate
physical constants are used to give the corresponding correlation length and correlation functions vs. time. These in-
clude an initial average film height H̄ and a white noise amplitude ∆ (Eq. (38)). These initial conditions approximate
a film at the beginning of an anneal that immediately follows a rapid deposition. The time tlarge is found by solving for
the time where the mean-square fluctuations are atomic scale,
h(x, t)2
= C(∆x = 0) = 1 ML2. At this point, the
correlation lengths are calculated.
Physical constants for the 2D anisotropic calculation are taken as follows. The elastic constants for Ge at 600 K
are c11 = 1.199 × 1012, c12 = 4.01 × 1011(from cS = 3.991), c44 = 6.73. [51] Using aGe = 0.5658nm and
aSi = 0.5431nm, it is found that �m = 0.0418. Using the procedure from (Appendix C), M = 1.332× 1012dyn/cm2.
E0◦ = 4.96 × 109erg/cm3, and E45◦ = 4.35 × 109erg/cm3 , giving �A = 0.1236. The atomic volume is Ω =
2.27 × 10−23 cm3. The estimated surface energy density is γ = 1927 erg/cm2. The wetting potential is estimated
by picking a plausible critical surface height, Hc ≈ 4 ML = 1.132 nm and setting W (H) = E20◦H3c/(8γH) =
2.315 × 10−6/H erg/cm2. The resulting characteristic wave number is kc = 0.257 nm−1. The initial film height is
taken to be H̄ = Hc + 0.25 ML = 1.203 nm and then allowed to evolve naturally. Thus, β = 0.208, α0 = 0.5658,
k0 = 0.1456 nm−1, σ0 = 0.1192/tc, σ‖ = 0.864/(k2c tc), σ⊥ = 0.559/(k
c tc), L‖ = 0.744k
0 (t/tc)
1/2, and
L⊥ = 0.599k
0 (t/tc)
1/2. The unspecified diffusivity has been absorbed into the characteristic time tc. From Eq. (38),
∆2 = 0.0403 nm4, and Eq. (50) gives
C(0) =
1.223× 10−3tc/t
e0.02385t/tc nm2.
The initial infinitely rough surface undergoes a smoothing described by the tc/t factor. Then the surface roughens due
to the exponential. The initial divergent roughness is an artifact of the non-physical white noise with the atomic scale
cutoff b0 neglected (Appendix F). The time for the fluctuations to become “large” again are found by setting
C(0) = h2large (53)
where hlarge = 1 ML = 0.283 nm. The solutions are t1 = 0.01527tc or t2 = 430tc. The first solution is discarded
since it is due to the non-physical white noise. At tlarge = t2, L‖ = 105.8 nm, and L⊥ = 85.2 nm. Taking L⊥ as
more limiting, the correlation spans about n = k0L⊥/π = 3.95 islands across. The corresponding reciprocal space
(Eq. (43)) and real-space correlation function (Eq. (50)) are shown in Figs. 6.d and 5.f respectively.
A corresponding numerical experiment is performed. A periodic surface of size l = 96(2π/k0) is used. Random
initial conditions consistent with Eq. (38) are used for k−space points on a square grid bounded by kx, ky = ±2k0.
The relation between discrete and continuous Fourier components is used, (hk)discrete = [(2π)d/A]hk. Eqs. (13)
and (14) are used without any additional approximation to find hk at time t = tlarge. The resulting CAk , a portion of
the height profile h(x) and CA(∆x) are plotted in Figs. 6(c), 5(d) and 5(f) respectively.
Similar calculations can be performed for the one-dimensional and two-dimensional elastically isotropic cases.
Isotropic values used previously [58, 24] are about E = 1.361 × 1012 dyn/cm2 and ν = 0.198 giving M = E/(1 −
ν) = 1.697 × 1012 dyn/cm2 and E = 2M(1 + ν) = 7.10 × 109 erg/cm3. Using the same critical surface height,
Hc = 4 ML, W (H) = 4.74×10−6/H erg/cm2. The resulting characteristic wave number is kc = 0.368 nm−1. If the
film is grown to H̄ = Hc+0.25 ML = 1.203 nm and then allowed to evolve naturally, β = 0.208; thus, α0 = 0.5658,
k0 = 0.208 nm−1, σ0 = 0.1192/tc, σ2 = 0.864/(k2c tc), and Lcor = 0.744k
0 (t/tc)
1/2. In one dimension, Eq. (46)
is used to find the mean square height fluctuation. Using Eq. (38) with d = 1, ∆2 = 0.0568 nm3, and
C(0, t) = 0.01271(t/tc)
−1/2e0.0238t/tc .
Setting C(0, t) = (1 ML)2 = 0.0801 nm2, t1 = 0.0252tc, and t2 = 186.9tc. At t2, Lcor = 48.8 nm, and n =
k0Lcor/π = 3.24, so about 3 dots in a row should be well correlated. The corresponding numerical calculation of
size l = 96(2π/k0) is performed. A portion of h(x), CA(∆x) and C(∆x) are shown in Fig. 7. In two dimensions,
Eq. (47) is used to find
h(x, t)2
C(0, t) = 9.40× 10−4(t/tc)−1/2e0.0238t/tc .
Setting C(0, t) = 0.0801 nm2, t1 = 1.376 × 10−4tc, and t2 = 306tc. At t2, Lcor = 62.4 nm, and n = k0Lcor/π =
4.14, and correlation is expected to extend about 4 dots. However, it should be noted that this correlation is not
lattice-like. Corresponding numerical results and ensemble correlation functions are shown in Figs. 6 and 5.a-c.
4.2 General case of β
In [24] it was suggested that allowing the film to evolve with β close to the stability threshold could enhance the SAQD
correlation. It is interesting to note what happens for different values of β. Similar analytic and numerical calculations
are performed for the large film-height limit, β = 0, for the 2D anisotropic Ge/Si surface. For β = 0, tlarge = 40.3tc,
L⊥ = 30.0 nm, and n = k0L⊥/π = 1.84, so one to two dots in a row are expected to be well correlated. h(x)
and real-space correlation functions are shown in Figs. 5g-i. The range of order is significantly less than for the case
β = 0.208 (Sec. 4.1). For Si/Ge at 600K, the 2D anisotropic predictions for tlarge and L⊥ are shown in Fig. 8. In
general, the closer β is to the critical value 0.25, the longer the correlation length. One can manipulate equation (53)
to find that tlarge/tc varies approximately but not exactly as (β − 1/4)−1 × ln[h2large
�A/(∆2k2c )]. Consequently,
L⊥ ∼ (β − 1/4)−1/2. Furthermore, the appearance of hlarge and ∆2 inside the logarithm shows that the final order
estimates are not overly sensitive to the guesses for ∆2 and h2large. The divergence of L⊥ with β − 1/4 is initially
encouraging, but it is clear that for the parameters used for Ge/Si, subatomic control of the film height is needed to
yield significantly enhanced long range correlations. Also as one approaches this threshold, one can probably expect
thermal activation to nucleate subcritical SAQDs whose effect on supercritically formed SAQDs is uncertain. There
should be some interesting phenomena at the theH → Hc.
l arge c
l a r g e
k 0L �
� 0 �
l a r g e � �
Figure 8: tlarge and L⊥ vs. β for Si/Ge using the 2D anisotropic model as described in Sec. 4. Units are normalize to
characteristic time tc and predicted number of correlated dots (n = k0L⊥/π).
5 Discussion/Conclusions
The order of epitaxial self-assembled quantum dots during initial stages of growth has been studied using a common
model of surface diffusion with stochastic initial conditions. It has been shown that correlation functions of small
surface height fluctuations can be predicted analytically using corresponding ensemble average correlation functions.
These correlation functions are characterized by correlation lengths that can be predicted by analytic formulas given
certain reasonable assumptions about the diffusion potential and the height and lateral scale of initial atomic scale
random fluctuations. Thus, the linear model of film surface height evolution via surface diffusion has enabled analytic
predictions of epitaxial SAQD order that are valid for small film height fluctuations. To what extent the initial degree
of order persists into later stages of growth remains to be studied, but the order of initial stages should certainly have
a strong influence on final outcomes. Furthermore, the linear analysis should provide insight into the less tractable
non-linear behavior. These predictions of SAQD order have been used to investigate the role of crystal anisotropy and
initial film height.
Crystal anisotropy has been shown to play an important role in enhancing SAQD order as observed in previous
numerical simulations continuum and atomistic numerical simulations. [43, 37, 44, 45] If a four-fold symmetry is
assumed for the governing dynamics, the effect of crystal anisotropy to linear order is felt through elastic anisotropy
alone. It is shown that elastic anisotropy is required to produce a lattice-like structure of SAQDs. The enhanced spatial
order should in turn lead to enhanced size order, a consequence that must be confirmed with non-linear studies, but
appears to be true based on the present available literature.
The role of initial film height has been shown to greatly influence order. Growth near the critical film height for
dot formation can enhance order. This order enhancement comes from increasing the duration of the linear small-
fluctuation stage of growth. In fact, the predicted correlation lengths diverge when the initial film height approaches
the critical film height from above. Achieving large correlation lengths in this manor is of course practically limited
by ability to control film heights to subatomic accuracy. Additionally, one should be careful when interpreting the
continuum model in such a context, as the effect of atomic discreteness might be greater at the transition film height.
Finally, it is likely that additional randomizing effects of thermal activation will effectively cut off this divergence
when the critical film height is approached from below during deposition.
Finally, the presented method may be useful as a first step in the analysis of methods to enhance SAQD order. It is
reasonable to suppose that under some circumstances initial growth stages will be very important while for others they
will not. For example, prior work on vertical stacking appears to confirm the presented ordering mechanism. [44].
Vertical stacking not only achieves vertical correlation of dots, but each layer is more ordered horizontally than the
one below. Additionally, a “growth window” was found, whereby to achieve enhanced order, the evolution of each
layer be terminated before ripening begins. The reported simulation [44] supports the following scenario for SAQD
order development. Order is enhanced during the small fluctuation stage as described here. Once the fluctuations are
sufficiently large, the seeded dots evolve towards their equilibrium shapes. Finally, the dots begin to ripen and order
diminishes. Order is transfered via strain to the next layer so that the next layer gets a head start on its initial ordering.
Thus, the multiple layers of dots effectively draws out the linear growth stage. It may be possible to modify the present
model to predict the correlation length of each SAQD layer.
A Diffusion Potential
The diffusion potential is calculated in terms of the film height H that is a function of the in plane coordinates x =
xi + yj. The elastic and surface energy portions of the diffusion potential can be found in [15]
µelast(x) = Ωω(x), and µsurf = −Ωγκ(x),
where Ω is the atomic volume, ω(x) is the elastic energy density at the film surface, γ is the surface energy density,
and κ is the total surface curvature. However, other calculations need to be included:
1. µwet for the two wetting potential cases, Eq. (3) and (5),
2. and µsurf and µwet when the surface energy density γ and wetting energy density W also depend on surface
orientation.
Before these case are addressed, a general form for the diffusion potential is justified.
A.1 General Form µ = ΩδF/δH(x)
The diffusion potential, µ(x), is the change in free energy, F , when a particle is added at a position, x. Note that µ(x)
and F are relative energies. They can be used to compare the binding energy of one site on the surface in comparison
with another site, but should not be interpreted as an absolute binding energy or total formation energy of the surface.
If a particle has a volume Ω, then the diffusion potential at x is related to the variation of free energy with volume,
δF = Ω−1
ddxµ(x)δV (x), (54)
where δV (x) is the volume variation at x. Calculating δV (x), V =
ddxH(x).Therefore, δV (x) = δH(x).
Substituting into δF (Eq. (54)), δF = Ω−1
ddxµ(x)δH(x) or µ(x) = ΩδF/δH(x).
A.2 Simple Model
Starting from Eq. (2), µ(x) is found by taking the variational derivative,
µelast.(x) = Ω
δH(x)
volume
ddxdz ω[H](x, z) = Ωω (x)
where the “[H]” indicates that the elastic energy, ω, is a nonlocal functional of the film height H, and ω(x) =
ω[H] (x,H(x)), the elastic energy density evaluated above lateral position x at the free surface (z = H(x)). See [15]
for details of the derivation. The surface energy diffusion potential is
µsurf.(x,t) = Ω
δH(x)
1 + (∇H(x))2
= −Ω∇ ·
1 + (∇H(x))2
γ = −Ωγκ(x).
The wetting energy diffusion potential is
µwet(x) = Ω
δH(x)
ddxW (H(x))
= ΩW ′(H(x))
Putting these three terms together, one obtains Eq. (3)
A.3 General Model
Consider the general form for the combined surface energy and wetting potential,
Fsw =
ddxFsw(H(x),∇H(x))
as in Eq. (4) so that the free energy is an integral over the x−plane of an energy density that depends on H(x) and
∇H(x) locally. The corresponding diffusion potential is
µ(x) = Ω
δH(x)
F (10)sw (H(x),∇H(x))−∇ · F
sw (H(x),∇H(x))
B Linearized Diffusion Potential and Anisotropy
The linearized diffusion potential µlin, k is found by finding µ(x) to first order in height fluctuations (h), to get µlin(x)
and then taking the Fourier transform to get µlin,k. The linearization of the simple isotropic diffusion potential corre-
sponding to Eqs. (2) and (3) was discussed in Sec. 2.1.1.1. Here, the more general diffusion potential corresponding
to Eqs (4) and (5) is linearized and then applied to the anisotropic simple model and the anisotropic general model.
Only the surface and wetting parts of the diffusion potential are discussed in this appendix. See ref. [15], Sec. 2.2.1.1
and Appendix C for discussion of µelast..
B.1 Linearizing the simple model
Consider a wetting potential and diffusion potential that both depend on the film height gradient ∇H, γ → γ(∇H)
and W (H) → W (H,∇H). Starting from Eq. (6) and expanding to second order in the film height fluctuation using
H(x) = H̄+ h(x) (Eq. (7)),
1 + (∇H)2
]−1/2
γ(∇H) =
(∇h)2 + . . .
γ + γ′ ·∇h+ γ̃′′ : ∇h∇h+ . . .
= γ + γ′ ·∇h−
γ (∇h)2 + γ̃′′ : ∇h∇h+O[h3]
where γ is γ(0), and the primes indicate the derivatives with respect to the surface height gradient.
γ′ = ∂∇Hγ(∇H)|∇H=0 , and γ̃
′′ = ∂∇H∂∇Hγ(∇H)|∇H=0 .
Taking the derivative with respect to ∇h results in a tensor of rank equal to the order of the derivative because ∇h is
a vector (tank 1 tensor). Taking the variational derivative, µsurf.(x) = ΩδFsurf./δh(x),
µsurf., lin(x) = Ω
γ∇2h(x)− γ̃′′ : ∇∇h(x)
The term with γ′ vanishes because it is the divergence of a constant (∇ · γ′). Taking the inverse Fourier transform,
µsurf., lin,k = Ω
−γk2 + k · γ̃′′ · k
hk. (55)
The first term is isotropic. The second term is parameterized by a rank 2 symmetric tensor.
Going through the same process, one finds essentially the same result for an orientation dependent wetting energy.
The step details are so close to the details for linearizing the more general form, Fsw(H,∇H), they are deferred to
(Appendix B.2). One finds that
µwet,lin,k = Ω
W (20) + k · W̃(02) · k
. (56)
where W (mn) = ∂mH∂
∇HW (H,∇H)|H=H̄,∇H=0 is the m
th and nth derivative of the wetting energy density with
respect toH and ∇H evaluated for a perfectly flat film of height H̄. W (mn) is a tensor of rank n.
B.1.1 isotropic case
In the isotropic case, γ̃′′ → γ′′Ĩ, where Ĩ is the identity operator, and γ′′ is a scalar. Similarly, W̃(02) →W (02)Ĩ. One
thus gets for the combined surface and wetting parts of the diffusion potential,
µsw,lin,k = Ω
−γ + γ′′ +W (02)
k2 +W (20)
Thus, in the isotropic case, the linear order effect of introducing a surface orientation to either the surface energy or
the wetting potential is simply to change the apparent surface energy density by γ → γ − γ′′ −W (02).
B.1.2 anisotropic case
The surface and wetting parts of the diffusion potential (Eqs. (55) and (56)) can admit only a limited anisotropy.
They both contain rank 2 symmetric tensors, γ̃′′ and W̃(02) in the x−plane. For a two-dimensional surface, this
means that they can either have two-fold-symmetric (rotations by 180◦) anisotropy or none at all. Thus, for the case
considered in Sec. 2.2.1.2, four-fold-symmetric anisotropy , the surface and wetting parts of the diffusion potential
must be completely isotropic. As discussed in Sec. 2.2.1.2, the (100) surface of zinc-blend structures, such as the
mentioned Ge, Si, InAs and GaAs present a rather complicated situation. For simplicity, it is assumed here that the
surface and wetting energies are at least four-fold symmetric. Consequently, they are completely isotropic.
Finally, it should be noted that if Fsw depends on higher order derivatives, then the discussion is greatly compli-
cated and a larger class of anisotropic terms is admissible. For example, whenFsw → Fsw(H,∇H,∇∇H,∇∇∇H, . . . )
is expanded aboutH(x) = H̄ to quadratic order in h, it would contains tensors of rank 6 and maybe even higher.
B.2 Linearizing the general model
The elastic part of the linearized diffusion potential was discussed in Sec. 2.2.1.1 and Appendix C . Eq. (56) can be
found by using all of the following steps with the substitution Fsw →W . The surface-wetting part of the diffusion po-
tential µ(x) is found by expanding Fsw to second order in the film-height fluctuation, h, and then taking the variational
derivative. Expanding Fsw about h = 0 and ∇h = 0,
Fsw(H̄+ h,∇h) = F (00)sw + F
sw h+ F
sw ·∇h+ hF
sw ·∇h . . .
· · ·+
F (20)sw h
F̃(02)sw : ∇h∇h+O[h
Note that in this expansion, all the F (mn)sw terms are constant with respect to h and depend implicitly on the average
film height, H̄. The first index indicates the mth derivative with respect to h. The second index indicates the nth
derivative with respect to ∇h. The derivatives are evaluated for a perfectly flat surface of height H̄. Thus,
F (mn)sw = ∂
∇HFsw (H,∇H)|H=H̄,∇H=0 .
Since ∇h is a vector in the x−plane, F (mn)sw is a tensor of rank n. Taking the variational derivative of Fsw =∫
ddxFsw(H,∇H) and keeping terms to order h1,
δh(x)
= F (10)sw −∇ · F
sw + F
sw h−∇ ·
F̃(02)sw ·∇h
Note that the F (00)sw term vanishes because it is constant, and the F
sw term vanishes upon simplification. Additionally,
the F (10)sw can be neglected if one enforces the condition that the film-height fluctuations do not add or subtract material
from the surface, namely that
ddx δh(x, t) = 0. Alternatively, one can discard it in anticipation of taking the gradient
of the diffusion potential, since it is a constant. The term ∇ · F(01)sw = 0 for the same reasons, or because F
sw is a
constant. Multiplying through by the atomic volume,
µlin(x) = Ω
F (20)sw h− F̃
sw : ∇∇h
. (57)
B.2.1 isotropic case
In the isotropic case, F̃(02)sw must be proportional to the identity so that F̃
sw = F
sw Ĩ; thus,
µsw,lin(x) = Ω
F (20)sw h(x)− F
2h(x)
Taking the inverse Fourier transform of this equation,
µsw,lin,k = Ω
F (20)sw + F
This gives case b in Eq. (9).
B.2.2 anisotropic case
If the surface is anisotropic, then F̃(02)sw in Eq. (57) is a rank 2 symmetric tensor in the x−plane. Thus, it can have two
distinct eigenvalues, and automatically has 2-fold rotational symmetry (rotations by 180◦). If any other symmetry is
assumed such as 4-fold symmetry (rotations by 90◦), then F̃(02)sw must be fully isotropic. Taking the inverse Fourier
transform,
µsw,lin,k = Ω
F (20)sw + k · F̃
sw · k
In Eq. (23), case b, it is assumed that there is four-fold symmetry, resulting in a surface-wetting part of the diffusion
potential that is completely isotropic.
C Elastic Anisotropy
In principal, the anisotropic elastic energy ωk is found in the same fashion as the isotropic elastic energy. [15] The
flat film, initially in a state of biaxial stress, is perturbed by a small periodic surface fluctuation of amplitude h0. An
appropriate elastic field is added to satisfy the perturbed traction-free boundary condition at the free surface. Finally,
the elastic energy is evaluated at the free surface to first order in h0. The coefficient h0 is the sought after ωk. The
equations themselves are cumbersome and best solved using a numeric implementation, so an abstract procedure for
calculating ωk is outlined here. ωk is found for k = 1 but arbitrary θk.
Let the surface have a height variation
h(x) = h0e
To first order in h0, the surface normal is
n(x) = −ikh0eikxi + k.
The elastic energy needs to be calculated to first order in h0. To find the elastic energy, it is necessary to find the
perturbing elastic field to first order in h0.
The initial unperturbed stress state is
σ̃m =
 σm 0 00 σm 0
0 0 0
where σm =
c11 + c12 − 2c212/c11
�m. Note that this stress state is isotropic in the x−y-plane and thus independent
of rotations about the vertical axis. Under this stress state, a flat surface is traction-free. With the height perturbation,
the traction is
tj = (n · σ̃m)j = −ikh0M�mδj1e
ikx. (58)
Next to find the perturbing elastic fields. These are not isotropic in the x − y−plane, and it is necessary to take
into account the angle. First, the 3× 3× 3× 3 elastic stiffness tensor cijkl is constructed for the cube orientation from
the compact 9× 9 matrix cij . The tensor representation aids in rotation. The stiffness tensor is then passively rotated
in the x− y−plane by an angel θk,
cijkl(θk) =
m,n,p,q=1
R(θk)imR(θk)jnR(θk)kpR(θk)lqcmnpq
where
R(θk) =
 cos(θk) sin(θk) 0− sin(θk) cos(θk) 0
0 0 1
This passive rotation of cijkl is equivalent to actively rotating the wave vector k = ki by θk.
The appropriate form for the perturbing displacement field is found. Assume a displacement of the form
ui(x, y, z) = Uie
k(ix+κz),
where κ can have a complex value. The elastic equilibrium equations are
i,k,l=1
cijkl(θk)
ul = 0; j = 1 . . . 3.
Cjl(θk, κ)Ul
k2ek(ix+κz) = 0 (59)
where
Cjl(θk, κ) =
i,k=1
cijkl(θk)(iδi1 + δi3κ)(iδk1 + δk3κ).
Factoring out k2ek(ix+κz), the part in parenthesis must be identically zero.
To obtain a non-trivial solution, the determinant of Cjl(θk, κ) to zero. Six complex values of κ are found. The
values of κ with Re[κ] < 0 are discarded since the corresponding displacements blow up as z → −∞. Each of
the remaining values κ = κp with p = 1 . . . 3 is substituted back into Cjl(θk, κ), and Eq. (59) is solved to find the
corresponding eigenvectors, Upl . The total displacement is thus
ul(x, y, z) = i�mh0
k(ix+κpz),
where it is assumed that the perturbing elastic displacement field is proportional to h0and σm, and the factor of i is
put in for convenience. The coefficients Ap can be found from the traction-free boundary condition at the free surface.
The traction formula is
i,k,l=1
nicijkl(θk)
ul(x, y, z) = ik�mh0
i,k,l,p=1
nicijkl(θk)ApU
l (iδk1 + κ
pδk3)e
k(ix+κpz) (60)
The traction is already proportional to h0. Thus, all terms in the sum must be kept to zeroth order in h0 so that
h(x) = 0, and n(x) = k.
Thus, plugging z = 0 to Eq. (60),
tj = ik�mh0
(ic3j1l(θk) + κ
pc3j3l(θk))ApU
ikx. (61)
Since the total traction (Eqs. (58) and (61)) must be zero, the coefficients Ap are found from
KjpAp = Rj ,
where
Kjp =
(ic3j1l(θk) + κ
pc3j3l(θk))U
Rj = Mδj1
for j = 1 . . . 3. It is worth noting that only for the symmetry directions, θk = 0◦ and θk = 45◦ is the strain purely
plane-strain as it is for the elastically isotropic case.
The elastic energy at the film surface is found to order O(h0). If the stress and strain are expanded to first order in
h0, σ̃ = σ̃0 + σ̃1, and �̃ = �̃0 + �̃1, then
�̃ : c̃ : �̃ =
σ̃0 : �0 + σ̃0 : �̃1 +O(h
Thus,
U = U0 +M�m ((�1)11 + (�1)22)
(�1)11 =
= −�mkh0
(�1)22 = ∂u2/∂y = 0. Thus,
U = U0 − Eθkkh0e
where
Eθk = M�
where Apand U
1 are implicitly functions of θk. This procedure has been used to find the values of E0◦ and E45◦ for
Table. 2 and Sec. 4.
D Diffusional Anisotropy
In general, the surface diffusivity can depend on the film height H(x) and the surface orientation ∇H(x) so that the
surface current is
JS(x) = D̃(H(x),∇H(x)) ·∇sµ(x)
where ∇s is the surface gradient, and D̃ is a rank 2 tensor in the two-dimensional space tangent to the film surface at
x. Linearizing the surface current about a flat surface,
JS(x) = D̃(H̄) ·∇µlin(x)
where the diffusivity must be evaluated for h = 0 and ∇h = 0, since µlin(x) is already proportional to h(x). The lin-
earized diffusivity is a symmetric rank 2 tensor in the x−plane. Thus, it is similar to F̃sw discussed in Appendix B.2.2.
It is automatically either two-fold symmetry (rotations by 180◦) or it is completely isotropic. In Eq. (23), four-fold
symmetry of the surface is assumed. Thus, the diffusivity must be completely isotropic; D̃ → D, a scalar. Sec-
tion 2.2.1.2 and Appendix B.2.2 contain discussions of the symmetry properties of the various rank 2 tensors that
appear in the linear evolution equations. A limited case of diffusional anisotropy has been modeled via kinetic Monte
Carlo technique. [54]
E Correlation Functions
E.1 Mean Values
Equations (31) and (33) are central to the presented analysis. Here, they are derived. The two-point correlation func-
tions for a stochastic system are introduced. Then, the average of the autocorrelation function is taken and expressed
in terms of the two-point correlation functions. Finally, this average is simplified using the translational invariance of
the system (governing equations and ensemble of initial conditions).
The two-point real-space space correlation function is
C(x,x′) = 〈h(x)h(x′)∗〉 ,
and the reciprocal space correlation function is
Ckk′ = 〈hkh∗k′〉 .
These are related by the double Fourier transform,
Ckk′ = 1(2π)2d
ddxddx′ e−ik·x+ik
′·x′C(x,x′); (62)
C(x,x′) =
ddkddk′ eik·x−ik
′·x′Ckk′ . (63)
These ensemble correlation functions can be used to give the ensemble-mean autocorrelation function and spec-
trum function. In real space,
CA(∆x)
d2x′ 〈h(∆x + x′)h(x′)〉
d2x′ C(∆x + x′,x′). (64)
(2π)d
〈hkh∗k〉 =
(2π)d
Ckk. (65)
Fortunately, the translational invariance of the system simplifies these relations. Inspecting the governing equations
and invoking the translational invariance of the stochastic initial conditions, the resulting ensemble and its statistical
measures must also be translationally invariant. Thus under the translation by x′,
C(∆x + x′,x′) = C(∆x,0) = C(∆x), (66)
so that the independent variable is reduced to just the difference vector ∆x = x − x′. This relation can be used to
simplify both the real and reciprocal space relations.
The real space relation simplifies as follows.Inserting Eq. (66) into Eq. (64),
CA(∆x)
d2x′ C(∆x,0) = C(∆x). (67)
The reciprocal space relation (Eq. (62)) simplifies to
Ckk′ = Ckδ
2(k− k′) = Ck
(2π)d
δkk′ , (68)
where
(2π)d
d2∆x e−ik·∆xC(∆x).
One can see immediately from Eq. (67) that Ck is the Fourier transform of
CA(∆x)
= C(∆x), or one can plug
Eq. (68) into Eq. (65), to get
= Ck.
E.2 Variance and Convergence
The ergodic hypothesis is that an average with respect to a parameter such as position or time tends towards an
ensemble average. In this case,
CAk ≈
= Ck, (69)
and CA(∆x) ≈
CA(∆x)
= C(∆x).
when the surface area is very large. The ensemble average is a good substitute if the variance about the average
vanishes as the substrate area A becomes large. It is found that in reciprocal space,
Var(CAk ) =
= C2k. (70)
Thus, the ergodic hypothesis does not hold for CAk . In practice, C
k is a speckled version of Ck (Fig. 6) However, if
one smooths CAk by averaging over a small patch in reciprocal space of size ksmooth = 1/∆s, so that
CAk (∆s) =
)d/2 ∫
ddk′ e−
′−k)2CAk′ , (71)
then Var
CAk (∆s)
diminishes as 1/A. For sufficiently large ∆s,〈
CAk (∆s)
≈ Ck, (72)
CAk (∆s)
πd/2∆ds
C2k. (73)
Thus, the ergodic hypothesis (Eq. (69)) only holds for a smoothed version of CAk .
In real space,
CA(∆x)
CA(∆x)
CA(∆x)
(2π)d
e2ik·∆xC2k + C
, (74)
where the integral is bounded (finite) provided that either t > 0 or the atomic scale cutoff b0 > 0. Thus, the ergodic
hypothesis holds for the real space autocorrelation function.
E.2.1 Eq. (70)
First,
CAk C
is calculated.
CAk C
(2π)d
〈hkh∗khk′h
k′〉 .
Assume that he distribution of hk is gaussian. Also, assume that h(x) is real so that hkh−k = |hk|
2. Then,〈
= Ck1Ck2δ
d(k1 − k4)δd(k2 − k3) . . .
. . . +Ck1Ck2δ
d(k1 + k3)δ
d(k2 + k4) . . .
. . . +Ck1Ck3δ
d(k1 − k2)δd(k3 − k4).
Thus,
CAk C
(2π)d
δd(k− k′)
. . .
. . . +C2k
δd(k + k′)
+ CkCk′
δd(0)
. (75)
= C2k
δkk′ + δk(−k′)
+ CkCk′ , (76)
where Eq. (29) has been used liberally. Setting k = k′, results in Eq. (70).
E.2.2 Eq. (73)
Now consider CAk smoothed over a length ∆s (Eq. (71)). The mean value is
CAk (∆s)
)d/2 ∫
ddk′ e−
′−k)2 〈CAk′〉 .
)d/2 ∫
ddk′ e−
′−k)2Ck′ .
For sufficiently small ksmooth, (sufficiently large ∆s), Eq. (72) results.
The variance of CAk (∆s) is now calculated. First, it is necessary to calculate
CAk (∆s)
CAk (∆s)
ddk′ e−
′−k)2 . . .
. . . ×
ddk′′ e−
′′−k)2 〈CAk′CAk′′〉 .
Using Eq. (75) and Eq. (29) as needed,
CAk (∆s)
ddk′ddk′′ e−
′−k)2e−
′′−k)2
(2π)d
. . .
· · · ×
δd(k′ − k′′)
+ C2k
δd(k′ + k′′)
+ CkCk′
δd(0)
′−k)2C2k′ + e
− 12 ∆
s[(k′−k)2+(k′+k)2]C2k′
. . .
· · ·+
)d/2 ∫
ddk′ e−
′−k)2Ck′
The first integral is bounded (finite) because Ck is bounded. Let its finite value be denoted I . The second integral is
simply
CAk (∆s)
.Thus,
Var(CAk (∆s)) =
∆2ds I
a finite value that decreases as A−1 as required for the ergodic hypothesis to hold. For sufficiently small ksmooth (large
∆s), I ≈ (π/∆2s)d/2C2k, and Eq. (73) results. It should also be noted that the large ∆s required for this approximation
also creates a more stringent requirement that A be large.
E.2.3 Eq. (74)
Now, consider the real space auto-correlation function. First,
CA(∆x)CA(∆x)
is needed.
CA(∆x)CA(∆x)
ddkddk′ ei(k+k
′)·∆x 〈CAk CAk′〉
Proceeding in a fashion similar to the previous section (making use of Eqs. (75) and (29) as needed) ,
CA(∆x)CA(∆x)
(2π)2d
ddkddk′ ei(k+k
′)·∆x
δd(k− k′)
. . .
· · ·+ C2k
δd(k + k′)
+ CkCk′
δd(0)
(2π)d
e2ik·∆xC2k + C
. . .
· · ·+
ddk eik·∆xCk
ddk′ eik
′·∆xCk
(2π)d
e2ik·∆xC2k + C
CA(∆x)
Thus, Eq. (74) results. For the variance to be vanishing, the integral in Eq. (74) must be bounded (finite). If time,
t > 0, the exponential in Eq. (77) guarantees that the integral is bounded. For time t = 0, the integral is only bounded
if the atomic scale cutoff b0 > 0.
F Atomic Scale Cutoff
Starting from Eq. (39),
(2π)d
e2σkt−
. (77)
The effect of the small scale cutoff is both small and short-lived, as it only works to suppress fluctuations with large
wavenumbers. The most important fluctuations have wavenumbers between 0 and 2kc. Thus, the typical size of the
cutoff term is about b20k
c . If a typical dot size or spacing size 10 nm, and a typical atomic scale is 10
−1 nm, a typical
value for this term is about 10−3 − 10−2. To calculate the effect of the cutoff, it can absorbed into the time-dependent
part with the substitution
so that its effect lasts only as long as a perturbation with atomic scale curvature (κ = b0). Thus, Eq. (40) is a good
approximation.
Acknowledgement
Thanks to L. Fang and C. Kumar for useful comments during the writing of this article.
References
[1] D. Bimberg, M. Grnudmann, and N. N. Ledentsov. Quantum Dot Heterostructures. John Wiley & Sons, 1999.
[2] O. P. Pchelyakov, Yu. B. Bolkhovityanov, A. V. Dvurechenski, L. V. Sokolov, A. I. Nikiforov, A. I. Yakimov,
and B. Voigtländer. SiliconGermanium nanostructures with quantum dots: Formation mechanisms and electrical
properties. Semiconductors, 34(11):122947, 2000. [doi:10.1134/1.1325416].
[3] M. Grundmann. The present status of quantum dot lasers. Physica E, 5:167, 2000. [doi:10.1016/S1386-
9477(99)00041-7].
[4] Pierre M. Petroff, Axel Lorke, and Atac Imamoglu. Epitaxially self-assembled quantum dots. Physics Today,
pages 46–52, May 2001.
[5] Hui-Yun Liu, Bo Xu, Yong-Qiang Wei, Ding Ding, Jia-Jun Qian, Qin Han, Ji-Ben Liang, and Zhan-Guo Wang.
High-power and long-lifetime InAs/GaAs quantum-dot laser at 1080 nm. Applied Physics Letters, 79(18):2868–
70, 2001. [doi:10.1063/1.1415416].
[6] F. Heinrichsdorff, M.H. Mao, N. Kirstaedter, A. Krost, D. Bimberg, A. O. Kosogov, and P. Werner. Room-
temperature continuous-wave lasing from stacked InAs/GaAs quantum dots grown by metalorganic chemical
vapor deposition. Applied Physics Letters, 71(1):22–4, 1997. [doi:doi:10.1063/1.120556].
[7] D. Bimberg, N.N. Ledentsov, and J.A. Lott. Quantum-dot vertical-cavity surface-emitting laser. MRS Bulletin,
27(7):531–7, 2002.
[8] N. N. Ledentsov. Long-wavelength quantum-dot lasers on GaAs substrates: From media to device con-
cepts. IEEE Journal of Selected Topics in Quantum Electronics, 8(5):1015–23, September/October 2002.
[doi:10.1109/JSTQE.2002.804236].
[9] M Friesen, P Rugheimer, D. E. Savage, M. G. Lagally, D. W. van der Weide, R Joynt, and M. A. Eriksson.
Practical design and simulation of silicon-based quantum-dot qubits. Physical Review B, 67(12):121301 (R),
2003. [doi:10.1103/PhysRevB.67.121301].
[10] Yi-Chang Cheng, San-Te (Cing-Ming) Yang, Jyh-Neng Yang, Liann-Be Chang, and Li-Zen Hsieh. Fabrication of
a far-infrared photodetector based on InAs/GaAs quantum-dot superlattices. Optical Engineering, 42(1):11923,
2003. [doi:doi:10.1117/1.1525277].
[11] R. Krebs, S. Deubert, J.P. Reithmaier, and A. Forchel. Improved performance of MBE grown quantum-dot
lasers with asymmetricdots in a well design emitting near 1.3 µm. Journal of Crystal Growth, 251:7427, 2003.
[doi:10.1016/S0022-0248(02)02385-0].
[12] Hiroyuki Sakaki. Progress and prospects of advanced quantum nanostructures and roles of molecular beam
epitaxy. Journal of Crystal Growth, 251:9–16, 2003. [doi:10.1016/S0022-0248(03)00831-5].
[13] B. J. Spencer, P. W. Voorhees, and S. H. Davis. Morphological instability in epitaxially strained dislocation-free
films. Physical Review Letters, 67(26):3696–3699, 1991. [doi:10.1103/PhysRevLett.67.3696].
[14] Karl Brunner. Si/ge nanostructures. Reports on Progress in Physics, 65(1):27–72, 2002. [doi:10.1088/0034-
4885/65/1/202].
[15] L. B. Freund and S. Suresh. Thin Film Materials: Stress, Defect Formation and Surface Evolution, chapter 8.
Cambridge University Press, 2003.
[16] S. Yu Shiryaev, E. Verstlund Pedersen, F. Jensen, J. Wulff Petersen, J. Lundsgaard Hansen, and A. Nylandsted
Larson. Dislocation patterning - a new tool for spatial manipulation of Ge islands. Thin solid films, 294(1-
2):311–314, 1997. [doi: 10.1016/S0040-6090(96)09240-1].
[17] C. Kumar and L. H. Friedman. Simulation of thermal field directed self assembly of epitaxial self-assembled Ge
quantum dots. Journal of Applied Physics, in press.
[18] Lawrence H. Friedman and Jian Xu. Feasibility study for thermal-field directed self-assembly of heteroepitaxial
quantum dots. Applied Physics Letters, 88:093105, 2006. [doi:10.1063/1.2179109].
[19] S. Krishna, D. Zhu, J. Xu, and P. Bhattacharya. Structural and luminescence characteristics of cycled sub-
monolayer InAs/GaAs quantum dots with room-temperature emission at 1.3 µm. Journal of Applied Physics,
86:6135–8, 1999. [doi:10.1063/1.371664].
[20] R. Hull, J.L. Gray, M. Kammler, T. Vandervelde, T. Kobayashi, P. Kumar, T. Pernell, J.C. Bean, J.A. Floro,
and F.M. Ross. Precision placement of heteroepitaxial semiconductor quantum dots. Materials Science and
Engineering B, 101:1–8, 2003. [doi:10.1016/S0921-5107(02)00680-3].
[21] O. Guise, Jr. J. T. Yates, J. Levy, J. Ahner, V. Vaithyanathan, and D. G. Schlom. Patterning of sub-10nm ge
islands on si(100) by direct self-assembly. Applied Physics Letters, 87:171902, 2005. [doi:10.1063/1.2112198].
[22] X. Niu, R. Vardavas, R. E. Caflisch, and C. Ratsch. Level set simulation of directed self-assembly during epitaxial
growth. Physical Review B, 74(19):193403, Nov 2006. [doi:10.1103/PhysRevB.74.193403.
[23] Z. M. Zhao, T. S. Yoon, W. Feng, B. Y. Li, J. H. Kim, J. Liu, O. Hulko, Y. H. Xie, H. M. Kim, K. B.
Kim, H. J. Kim, K. L. Wang, C. Ratsch, R. Caflisch, D. Y. Ryu, and T. P. Russell. The challenges in
guided self-assembly of ge and inas quantum dots on si. THIN SOLID FILMS, 508(1-2):195–199, Jun 2006.
[doi:10.1016/j.tsf.2005.08.407].
[24] Lawrence H. Friedman. Anisotropy and order of epitaxial self-assembled quantum dots. Physical Review B, in
press.
[25] Y. Obayashi and K. Shintani. Directional dependence of surface morphological stability of heteroepitaxial layers.
Journal of Applied Physics, 84(6):3141, 1998. [doi:10.1063/1.368468].
[26] C. S. Ozkan, W. D. Nix, and H. J. Gao. Stress-driven surface evolution in heteroepitaxial thin films: Anisotropy
of the two-dimensional roughening mode. JOURNAL OF MATERIALS RESEARCH, 14(8):3247–3256, Aug
1999. [doi:10.1557/JMR.1999.043].
[27] J. Tersoff and F. K. LeGoues. Competing relaxation mechanisms in strained layers. Physical Review Letters,
72(22):3570–3573, May 1994. [doi:10.1103/PhysRevLett.72.3570].
[28] B. J. Spencer, P. W. Voorhees, and S. H. Davis. Morphological instability in epitaxially strained dislocation-
free solid films: Linear stability theory. Journal of Applied Physics, 73(10):4955–4970, 1993. [doi:
10.1063/1.353815].
[29] J. M. Baribeau, X. Wu, N. L. Rowell, and D. J. Lockwood. Ge dots and nanostructures grown epitaxially
on si. JOURNAL OF PHYSICS-CONDENSED MATTER, 18(8):R139–R174, Mar 2006. [doi:10.1088/0953-
8984/18/8/R01].
[30] D. J. Srolovitz. On the stability of surfaces of stressed solids. Acta Metallurgica, 37(2):621–625, 1989.
[doi:10.1016/0001-6160(89)90246-0].
[31] H. J. Gao and W. D. Nix. Surface roughening of heteroepitaxial thin films. ANNUAL REVIEW OF MATERIALS
SCIENCE, 29:173–209, 1999. [doi:0.1146/annurev.matsci.29.1.173].
[32] P. Sutter and M. G. Lagally. Nucleationless three-dimensional island formation in low-misfit heteroepitaxy.
Physical Review Letters, 84(20):4637, 2000. [doi:10.1103/PhysRevLett.84.4637.
[33] A. A. Golovin, S. H. Davis, and P. W. Voorhees. Self-organization of quantum dots in epitaxially strained solid
films. Physical Review E, 68:056203, 2003. [doi:10.1103/PhysRevE.68.056203].
[34] A. Ramasubramaniam and V. B. Shenoy. Growth and ordering of si-ge quantum dots on strain patterned sub-
strates. JOURNAL OF ENGINEERING MATERIALS AND TECHNOLOGY-TRANSACTIONS OF THE ASME,
127(4):434–443, Oct 2005. [doi:10.1115/1.1924559].
[35] I. Berbezier, A. Ronda, F. Volpi, and A. Portavoce. Morphological evolution of SiGe layers. Surface Science,
531:231–243, 2003. [doi:10.1016/S0039-6028(03)00488-6].
[36] J. R. R. Bortoleto, H. R. Gutierrez, M. A. Cotta, J. Bettini, L. P. Cardoso, and M. M. G. de Carvalho. Spatial order-
ing in InP/InGaP nanostructures. Applied Physics Letters, 82(20):3523–3525, 2003. [doi:10.1063/1.1572553].
[37] P. Liu, Y. W. Zhang, and C. Lu. Formation of self-assembled heteroepitaxial islands in elastically anisotropic
films. Physical Review B, 67:165414, 2003. [doi: 10.1103/PhysRevB.67.165414].
[38] Y.W. Zhang, A.F. Bower, and P. Liu. Morphological evolution driven by strain induced surface diffusion. Thin
solid films, 424:9–14, 2003. [doi:10.1016/S0040-6090(02)00897-0].
[39] Yu U. Wang, Yongmei M. Jin, and Armen G. Khachaturyan. Phase field microelasticity modeling of surface
instability of heteroepitaxial thin films. Acta Materialia, 52:81–92, 2004. [doi:10.1016/j.actamat.2003.08.027].
[40] W. T. Tekalign and B. J. Spencer. Evolution equation for a thin epitaxial film on a deformable substrate. Journal
of Applied Physics, 96(10):5505–5512, 2004. [doi:10.1063/1.1766084].
[41] M. J. Beck, A. van de Walle, and M. Asta. Surface energetics and structure of the ge wetting layer on si(100).
Physical Review B, 70(20):205337, Nov 2004. [doi:10.1103/PhysRevB.70.205337].
[42] Y. H. Tu and J. Tersoff. Origin of apparent critical thickness for island formation in heteroepitaxy. Physical
Review Letters, 93(21):216101, Nov 2004. [doi:10.1103/PhysRevLett.93.216101.
[43] V. Holy, G. Springholz, M. Pinczolits, and G. Bauer. Strain induced vertical and lateral correlations in quantum
dot superlattices. Physical Review Letters, 83(2):356–359, 1999. [doi:10.1103/PhysRevLett.83.356].
[44] P. Liu, Y. W. Zhang, and C. Lu. Three-dimensional finite-element simulations of the self-organized growth of
quantum dot superlattices. Physical Review B, 68:195314, 2003. [doi:10.1103/PhysRevB.68.195314].
[45] G. Springholz, M. Pinczolits, V. Holy, S. Zerlauth, I. Vavra, and G. Bauer. Vertical and lateral ordering in
self-organized quantum dot superlattices. Physica E, 9:149–163, 2001. [doi:10.1016/S1386-9477(00)00189-2.
[46] P. Liu, Y. W. Zhang, and C. Lu. Coarsening kinetics of heteroepitaxial islands in nucleationless stranski-krastanov
growth. Physical Review B, 68:035402, 2003. [doi:10.1103/PhysRevB.68.035402].
[47] F. M. Ross, J. Tersoff, and R. M. Tromp. Coarsening of self-assembled ge quantum dots on Si(001). Physical
Review Letters, 80(5):984–7, 1998. [doi:10.1103/PhysRevLett.80.984].
[48] J. Tersoff. Kinetic surface segregation and the evolution of nanostructures. Applied Physics Letters, 83(2):353–
355, 2003. [doi:doi:10.1063/1.1592304].
[49] A. Ramasubramaniam and V. B. Shenoy. A spectral method for the nonconserved surface evolution of nanocrys-
talline gratings below the roughening transition. Journal of Applied Physics, 97(11):114312, 2005. [doi:
10.1063/1.1897837].
[50] Y. W. Zhang and A. F. Bower. Three-dimensional analysis of shape transitions in strained-heteroepitaxial islands.
Applied Physics Letters, 78(18):2706–2708, 2001. [doi:10.1063/1.1354155].
[51] L. E. Vorbyev. Handbook Series On Semiconductor Parameters, volume 1. World Scientific, 1996.
[52] A. A. Golovin, M. S. Levine, T. V. Savina, and S. H. Davis. Faceting instability in the presence of wet-
ting interactions: A mechanism for the formation of quantum dots. Physical Review B, 70:235342, 2004.
[doi:10.1103/PhysRevB.70.235342].
[53] B L Liang, Zh M Wang, Yu I Mazur, V V Strelchuck, K Holmes, J H Lee, and G J Salamo. Ingaas quantum
dots grown on b-type high index gaas substrates: surface morphologies and optical propertiesmorphologies and
optical properties. Nanotechnology, 17(11):2736–2740, 2006. [doi:10.1088/0957-4484/17/11/004].
[54] M. Meixner, R. Kunert, and E. Scholl. Control of strain-mediated growth kinetics of self-assembled semicon-
ductor quantum dots. Physical Review B, 67:195301, 2003. [doi: 10.1103/PhysRevB.67.195301].
[55] Robert Zwanzig. Nonequilbrium Statistical Mechanics. Oxford University Press, New York, 2001.
[56] C. W. Gardiner. Handbook of Stochastic Methods for Physics Chemistry and the Natural Sciences. Springer,
New York, 3rd edition, 2004.
[57] M. C. Cross and P. C. Hohenberg. Pattern formation outside equilibrium. Reviews of Modern Physics, 65(3):851–
1112, 1993. [doi:10.1103/RevModPhys.65.851].
[58] B. J. Spencer, S. H. Davis, and P. W. Voorhees. Morphological instability in epitaxially strained dislocation-free
solid films: Nonlinear evolution. Physical Review B, 47(15):9760, 1993. [doi: 10.1103/PhysRevB.47.9760].
	Introduction
	Modeling
	1D and 2D Isotropic model
	Energetics
	simple form
	general form
	Linearization
	Dynamics
	Peaks
	Stability and wetting potential
	2D Anisotropic case
	Energetics
	Elastic anisotropy
	Surface and Wetting Energy Anisotropy
	total diffusion potential
	Dynamics
	Expansion about peaks
	Correlation Functions
	Correlation Functions and SAQD order
	Periodic array
	Nearly Periodic array
	Ensemble Correlation Functions / ergodicity
	Mean fluctuation
	Correlation Function
	Stochastic Initial Conditions
	Reciprocal Space Correlation Functions
	one-dimensional
	2D isotropic
	anisotropic
	Real Space Correlation Functions
	one-dimensional
	2D isotropic
	anisotropic
	Generalizability
	Order Predictions
	Ge at 600K
	General case of 
	Discussion/Conclusions 
	Diffusion Potential
	General Form =F/H(x)
	Simple Model
	General Model
	Linearized Diffusion Potential and Anisotropy
	Linearizing the simple model
	isotropic case
	anisotropic case
	Linearizing the general model
	isotropic case
	anisotropic case
	Elastic Anisotropy
	Diffusional Anisotropy
	Correlation Functions
	Mean Values
	Variance and Convergence
	Eq. (??)
	Eq. (??)
	Eq. (??)
	Atomic Scale Cutoff
